参数估计

根据一系列独立样本，以及模型，找到参数

最大似然估计(Maximum likelihood estimation)

计算过程

1. 写出似然函数；
2. 如果无法直接求导的话，对似然函数取对数；
3. 求导数；
4. 求解模型中参数的最优值。

Example1：抽球

假设一个袋子装有白球与红球，比例未知，现在抽取10次（每次抽完都放回，保证事件独立性），假设抽到了7次白球和3次红球，在此数据样本条件下，可以采用最大似然估计法求解袋子中白球的比例（最大似然估计是一种“模型已定，参数未知”的方法）。当然，这种数据情况下很明显，白球的比例是70%，但如何通过理论的方法得到这个答案呢？一些复杂的条件下，是很难通过直观的方式获得答案的，这时候理论分析就尤为重要了，这也是学者们为何要提出最大似然估计的原因。我们可以定义从袋子中抽取白球和红球的概率如下：

为第一次采样，为第二次，为模型，模型参数

我们定义似然为：

两边取，取是为了将右边的乘号变为加号，方便求导，左边的通常称之为对数似然：

平均似然对数：

最大似然估计的过程就是找一个合适的，使得平均对数似然的值为最大。因此，可以得到以下最大估计的公式：

这里讨论的是2次采样的情况，拓展到多次采样的情况，次采样最大似然估计公式：

我们定义为模型（也就是之前公式中的），表示抽到白球的概率为，而抽到红球的概率为，因此10次抽取抽到白球7次的概率可以表示为：

将其描述为平均似然可得：

最大似然就是找到一个合适的，获得最大的平均似然。我们可以对平均似然的公式对求导，并令导数为0。

由此可得，当抽取白球的概率为0.7时，最可能产生10次抽取抽到白球7次的事件。

Example2：正态分布

假如有一组采样值，我们知道其服从正态分布，且标准差已知。当这个正态分布的期望为多少时，产生这个采样数据的概率为最大？

这个例子中正态分布就是模型，而期望就是前文提到的，似然函数如下：

正态分布的公式，当第一参数（期望）为0，第二参数（方差）为1时，分布为标准正态分布：

所以似然值为：

对上式求导可得：

最大后验概率估计(Maximum a posteriori)

假设有五个袋子，各袋中都有无限量的饼干(樱桃口味或柠檬口味)，已知五个袋子中两种口味的比例分别是

1.樱桃 100% 2.樱桃 75% + 柠檬 25% 3.樱桃 50% + 柠檬 50% 4.樱桃 25% + 柠檬 75% 5.柠檬 100%

如果只有上所述条件，从同一个袋子中连续拿到2个柠檬饼干，那这个袋子最有可能是上述五个的哪一个？

我们首先采用最大似然估计来解这个问题，写出似然函数。假设从袋子中能拿出柠檬饼干的概率为(我们通过这个概率来确定是从哪个袋子中拿出来的)，则似然函数可以写作

由于的取值是一个离散值，即上面描述中的0,25%，50%，75%，1。我们只需要评估一下这五个值哪个值使得似然函数最大即可，得到为袋子5。这里便是最大似然估计的结果。上述最大似然估计有一个问题，就是没有考虑到模型本身的概率分布，下面我们扩展这个饼干的问题。

假设拿到袋子1或5的概率都是0.1，拿到2或4的概率都是0.2，拿到3的概率是0.4，那同样上述问题的答案呢？这个时候就变MAP了（若本例中拿到五个袋子概率都是0.2，即为均匀分布时，MAL与MLE问题是等价；当不同时，即本例子，两问题不等价）。我们根据公式

写出我们的

根据题意的描述可知，的取值分别为0,25%, 50%, 75%, 1，的取值分别为0.1, 0.2, 0.4, 0.2, 0.1.分别计算出MAP函数的结果为：0, 0.0125, 0.125, 0.28125, 0.1.由上可知，通过MAP估计可得结果是从第四个袋子中取得的最高。

EM算法(Expectation–maximization Algorithm)

EM算法是一种迭代算法，用于含有隐变量的概率模型参数的极大似然估计，或极大后验概率估计。EM算法的每次迭代由两步组成：E步：求期望(Expectation)；M步求极大(Maximization)。

算法步骤

输入：观测变量数据，隐变量数据，联合分布，条件概率分布

输出：模型参数

1. 选择参数的初值，开始迭代

2. E步：记为第次迭代参数的估计值，在第次迭代的E步，计算

这里，是在给定观测数据和当前的参数估计下隐变量数据的条件概率分布。

1. M步：求使极大化的，确定第次迭代的参数的估计值

1. 重复第（2）步和第（3）步，直至收敛。

步骤说明

1. 参数的初值可以任意选择，但需注意EM算法对初值是敏感的

2. E步求。函数式中是未观测数据，是观测数据。注意的第1个变元表示要极大化的参数，第2个变元表示参数的当前估计值。每次迭代实际在求函数及其极大。

3. M步求极大化，得到，完成一次迭代，每次迭代使似然函数增大或达到局部极值

4. 给出停止迭代的条件，一般是对较小的正数，满足或则停止迭代。

Q函数

E步中的是EM算法的核心，称为函数。完全数据的对数似然函数关于在给定观测数据和当前参数下对未观测数据的条件概率分布的期望称为函数

Example

假设有3枚硬币，分别记作。这些硬币正面出现的概率分别是。进行如下掷硬币试验：先掷硬币，根据其结果选出硬币或硬币，正面选硬币，反面选；然后掷选出的硬币，掷硬币的结果，出现正面记作1，出现反面记作0；独立地重复n次试验（本例n = 10），结果如下

只能观测到掷硬币的结果，不能观测掷硬币过程。问如何估计三枚硬币正面出现的概率，即三硬币模型的参数。

三硬币模型可以写作：

这里，随机变量是观测变量，表示一次试验观测的结果是1或0；随机变量是隐变量，表示未观测到的掷硬币的结果；是模型参数。这一模型是以上数据的生成模型。注意，随机变量的数据可以观测，随机变量的数据不可观测。将观测数据表示为，未观测数据表示为，则观测数据的似然函数为

即

考虑求模型参数的极大似然估计，即

这个问题没有解析解，只有通过迭代的方法求解。EM算法就是可以用于求解这个问题的一种迭代算法。

EM算法首先选取参数的初值，记作，然后通过下面的迭代计算参数的估计值，直到收敛为止。第次迭代参数的估计值为。EM算法的第次迭代如下

E步：计算在模型参数下观测数据来自掷硬币的概率

M步：计算模型参数的新估计值

进行数字计算。假设模型参数的初值取为

由E步得到对于与均有

由M步迭代公式得到

再回到E步，得到

继续M步，得

收敛，于是得到模型的参数的极大似然估计

如果取初值，那么得到的模型参数的极大似然估计是。这就是说，EM算法与初值的选取有关，选择不同的初值可能得到不同的参数估计值。

Code实现

E-step：

import numpy as np
import math
pro_A, pro_B, por_C = 0.5, 0.5, 0.5
def pmf(i, pro_A, pro_B, por_C):
    pro_1 = pro_A * math.pow(pro_B, data[i]) * math.pow((1-pro_B), 1-data[i])
    pro_2 = pro_A * math.pow(pro_C, data[i]) * math.pow((1-pro_C), 1-data[i])
    return pro_1 / (pro_1 + pro_2)

M-step：

class EM:
    def __init__(self, prob):
        self.pro_A, self.pro_B, self.pro_C = prob
    # e_step
    def pmf(self, i):
        pro_1 = self.pro_A * math.pow(self.pro_B, data[i]) * math.pow((1-self.pro_B), 1-data[i])
        pro_2 = (1 - self.pro_A) * math.pow(self.pro_C, data[i]) * math.pow((1-self.pro_C), 1-data[i])
        return pro_1 / (pro_1 + pro_2)
    # m_step
    def fit(self, data):
        count = len(data)
        print('init prob:{}, {}, {}'.format(self.pro_A, self.pro_B, self.pro_C))
        for d in range(count):
            _ = yield
            _pmf = [self.pmf(k) for k in range(count)]
            pro_A = 1/ count * sum(_pmf)
            pro_B = sum([_pmf[k]*data[k] for k in range(count)]) / sum([_pmf[k] for k in range(count)])
            pro_C = sum([(1-_pmf[k])*data[k] for k in range(count)]) / sum([(1-_pmf[k]) for k in range(count)])
            print('{}/{}  pro_a:{:.3f}, pro_b:{:.3f}, pro_c:{:.3f}'.format(d+1, count, pro_A, pro_B, pro_C))
            self.pro_A = pro_A
            self.pro_B = pro_B
            self.pro_C = pro_C

测试

data=[1,1,0,1,0,0,1,0,1,1]
em = EM(prob=[0.5, 0.5, 0.5])
f = em.fit(data)
next(f)
# 第一次迭代
f.send(1)
# 第二次
f.send(2)
em = EM(prob=[0.4, 0.6, 0.7])
f2 = em.fit(data)
next(f2)
f2.send(1)
f2.send(2)