最大似然估计

时间:2023-02-13 16:00:30

  极大似然估计(Maximum likelihood estimation, 简称MLE)是很常用的参数估计方法,极大似然原理的直观想法是,一个随机试验如有若干个可能的结果A,B,C,...,若在一次试验中,结果A出现了,那么可以认为实验条件对A的出现有利,也即出现的概率P(A)较大。也就是说,如果已知某个随机样本满足某种概率分布,但是其中具体的参数不清楚,参数估计就是通过若干次试验,观察其结果,利用结果推出参数的大概值

EM算法:Expectation-Maximization:
最大似然:
    已知:(1)样本服从分布的模型, (2)观测到的样本
    求解:模型的参数
    总的来说:极大似然估计就是用来估计模型参数的统计学方法
    
    最大似然数学问题(100名学生的身高问题)
        样本集X={x1,x2,…,xN} N=100
        概率密度:p(xi|θ)抽到男生i(的身高)的概率
        独立同分布:同时抽到这100个男生的概率就是他们各自概率的乘积
        θ是服从分布的参数
        
        最大似然函数:l(a) = sum(logp(xi;a)) (对数是为了乘法转加法)
        什么样的参数 能够使得出现当前这批样本的概率最大
        已知某个随机样本满足某种概率分布,但是其中具体的参数不清楚,
        参数估计就是通过若干次试验,观察其结果,利用结果推出参数的大概值。

 本文以一个简单的离散型分布的例子,模拟投掷硬币估计头像(head)向上的概率。投掷硬币落到地面后,不是head向上就是tail朝上,这是一个典型的伯努利实验,形成一个伯努利分布,有着如下的离散概率分布函数:
最大似然估计

其中,x等于1或者0,即结果,这里用1表示head、0表示tail。

对于n次独立的投掷,很容易写出其似然函数:

最大似然估计

现在想用极大似然估计的方法把p估计出来。就是使得上面这个似然函数取极大值的情况下的p的取值,就是要估计的参数。

首先用Python把投掷硬币模拟出来:

from scipy.stats import bernoulli

# 生成样本
p_1 = 1.0 / 2  # 假设样本服从p为1/2的bernouli分布
fp = bernoulli(p_1)  # 产生伯努利随机变量
xs = fp.rvs(100)   # 产生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库把似然函数写出来:

import sympy
import numpy as np

# 估计似然函数
x, p, z = sympy.symbols('x p z', positive=True)
phi = p**x*(1-p)**(1-x)   #分布函数
L = np.prod([phi.subs(x, i) for i in xs])   # 似然函数
print(L)
# p**52*(-p + 1)**48

从上面的结论可以看出,作100次伯努利实验,出现positive、1及head的数目是52个,相应的0也就是tail的数目是48个,比较接近我们设的初始值0.5即1.0/2(注意:现在我们假设p是未知的,要去估计它,看它经过Python的极大似然估计是不是0.5!)。

下面,我们使用Python求解这个似然函数取极大值时的p值:

logL = sympy.expand_log(sympy.log(L))
sol = sympy.solve(sympy.diff(logL, p), p)
print(sol)
[13/25]

 

结果没有什么悬念,13/25的值很接近0.5!

取对数后,上面Python的算法最后实际上是求解下式为0的p值:

最大似然估计

上式留给网友自行推导,很多资料都可找到该式。这个式子,是著名的Logistic回归参数估计的极大似然估计算法的基础。

进一步,为了更加直观的理解投掷硬币的伯努利实验,我们给出以均值(均值为100*0.5=50)为中心对称的加总离散概率(概率质量函数(probability mass function),Python里面使用pmf函数计算):

from scipy.stats import binom
b = binom(100, .5)
# 以均值为中心对称的加总概率
g = lambda x: b.pmf(np.arange(-x, x) + 50).sum()
print(g(10))
0.9539559330706295

  对于上面的Python代码,可以通过下图更好地去理解:

最大似然估计

把这20个离散的概率全部显示出来,也可以看到在0.08左右取到它们的最大值

本文针对简单的离散概率质量函数的分布使用Python进行了极大似然估计,同时该方法可以应用于连续分布的情形,只要通过其概率密度函数得出其似然函数即可。