最大似然估计

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

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

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把投掷硬币模拟出来：

?

1 2 3 4 5 6 7 8

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库把似然函数写出来：

?

1 2 3 4 5 6 7 8 9

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值：

?

1 2 3 4

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函数计算）：

?

1 2 3 4 5 6

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进行了极大似然估计，同时该方法可以应用于连续分布的情形，只要通过其概率密度函数得出其似然函数即可。