最大似然估計

felwell 2019-10-17

展開全文

　　極大似然估計（Maximum likelihood estimation, 簡稱MLE）是很常用的參數(shù)估計方法，極大似然原理的直觀想法是，一個隨機試驗如有若干個可能的結(jié)果A，B，C，...，若在一次試驗中，結(jié)果A出現(xiàn)了，那么可以認(rèn)為實驗條件對A的出現(xiàn)有利，也即出現(xiàn)的概率P(A)較大。也就是說，如果已知某個隨機樣本滿足某種概率分布，但是其中具體的參數(shù)不清楚，參數(shù)估計就是通過若干次試驗，觀察其結(jié)果，利用結(jié)果推出參數(shù)的大概值

EM算法：Expectation-Maximization：

最大似然：

已知：（1）樣本服從分布的模型，（2）觀測到的樣本

求解：模型的參數(shù)

總的來說：極大似然估計就是用來估計模型參數(shù)的統(tǒng)計學(xué)方法

最大似然數(shù)學(xué)問題（100名學(xué)生的身高問題）

樣本集X={x1,x2,…,xN} N=100

概率密度：p(xi|θ)抽到男生i（的身高）的概率

獨立同分布：同時抽到這100個男生的概率就是他們各自概率的乘積

θ是服從分布的參數(shù)

最大似然函數(shù)：l(a) = sum(logp(xi;a)) （對數(shù)是為了乘法轉(zhuǎn)加法）

什么樣的參數(shù) 能夠使得出現(xiàn)當(dāng)前這批樣本的概率最大

已知某個隨機樣本滿足某種概率分布，但是其中具體的參數(shù)不清楚，

參數(shù)估計就是通過若干次試驗，觀察其結(jié)果，利用結(jié)果推出參數(shù)的大概值。

　本文以一個簡單的離散型分布的例子，模擬投擲硬幣估計頭像(head)向上的概率。投擲硬幣落到地面后，不是head向上就是tail朝上，這是一個典型的伯努利實驗，形成一個伯努利分布，有著如下的離散概率分布函數(shù)：

其中，x等于1或者0，即結(jié)果，這里用1表示head、0表示tail。

對于n次獨立的投擲，很容易寫出其似然函數(shù)：

現(xiàn)在想用極大似然估計的方法把p估計出來。就是使得上面這個似然函數(shù)取極大值的情況下的p的取值，就是要估計的參數(shù)。

首先用Python把投擲硬幣模擬出來：

from scipy.stats import bernoulli

# 生成樣本

p_1 = 1.0 / 2 # 假設(shè)樣本服從p為1/2的bernouli分布

fp = bernoulli(p_1) # 產(chǎn)生伯努利隨機變量

xs = fp.rvs(100) # 產(chǎn)生100個樣本

print(xs[:30]) # 看看前面30個

# [0 1 1 1 1 0 0 0 0 1 0 1 0 0 0 0 1 1 1 1 0 1 1 0 0 1 0 0 0 1]

通過此模擬，使用sympy庫把似然函數(shù)寫出來：

import sympy

import numpy as np

# 估計似然函數(shù)

x, p, z = sympy.symbols('x p z', positive=True)

phi = p**x*(1-p)**(1-x) #分布函數(shù)

L = np.prod([phi.subs(x, i) for i in xs]) # 似然函數(shù)

print(L)

# p**52*(-p + 1)**48

從上面的結(jié)論可以看出，作100次伯努利實驗，出現(xiàn)positive、1及head的數(shù)目是52個，相應(yīng)的0也就是tail的數(shù)目是48個，比較接近我們設(shè)的初始值0.5即1.0/2（注意：現(xiàn)在我們假設(shè)p是未知的，要去估計它，看它經(jīng)過Python的極大似然估計是不是0.5！）。

下面，我們使用Python求解這個似然函數(shù)取極大值時的p值：

logL = sympy.expand_log(sympy.log(L))

sol = sympy.solve(sympy.diff(logL, p), p)

print(sol)

[13/25]

結(jié)果沒有什么懸念，13/25的值很接近0.5！

取對數(shù)后，上面Python的算法最后實際上是求解下式為0的p值：