小样本分析（三）

          我怀着日思夜想的问题，希望你能给我答案，你慷慨的给与我答案，又慷慨的给了我一个新的问题。

导读：

这篇是继前两篇小概率估算的后续，很多地方没有说，建议先粗略看一下前面的，了解一下估算的思路。本文是说当有一个经验知识之后（比如所有球里面红球占到0.01），利用熵最大假设和贝叶斯公式优化m次抽样发生n次时的概率估计。由于公式较多直接从word截图偷个闲。

熵最大时概率分布推导：

【注：拉格朗日乘数法我也忘记得差不多了，在网上找到的例子也主要是只有一个约束条件，这里到两个约束的拓展有点不确定，但是先做了，后面也得到了想要的结果，这里先不管了】

不难发现在 alpha 比较小的时候（比如-100），exp(alpha)无限接近0，方程有两个解alpha=0和alpha=-1/p。在我的计算机模拟中也是如此。

源代码：

# coding=utf-8__author__ = '01053185'import numpy as npimport mathdef error_f(alpha, p):    # 计算误差的绝对值 e^α+(1/α+p)∙(1-e^α)=0  方程    return abs(math.exp(alpha) + (1.0 / alpha + p) * (1 - math.exp(alpha)))class P_test():    def __init__(self):        self.alpha = 0        self.beta = 0        self.p_pre = 0        self.x_array = np.array([0.0] * 10000)        self.pro_x_array = np.array([0.0] * 10000)        self.step_x_array = np.array([0.0] * 10000)        self.p_estimate = 0  # 最后估计的概率值    def solve_function(self, p):        self.p_pre = p        alpha = -1000 + 10 ** -9        sign = 1.0        step = 1.0        error_now = error_f(alpha, self.p_pre)        T = 0        while (error_now > 10 ** -9):            temp_alpha = alpha + step * sign            temp_error = error_f(temp_alpha, self.p_pre)            if temp_error < error_now:                alpha = temp_alpha                error_now = temp_error                T += 1                if T == 2:                    step *= 2                continue            sign *= -1  # 反向尝试            temp_alpha = alpha + step * sign            temp_error = error_f(temp_alpha, self.p_pre)            if temp_error < error_now:                T = 0                alpha = temp_alpha                error_now = temp_error                continue            step *= 0.3        self.alpha = alpha        self.beta = 1 - math.log(1.0 / self.alpha * (math.exp(self.alpha) - 1))    def set_array(self):        start = 0.0        step = 10.0 ** -9        x_array = np.array([0.0] * 100000)  # 记录某个区间的中间点        pro_x_array = np.array([0.0] * 100000)  # 记录某个区间的概率        step_x_array = np.array([0.0] * 100000)  # 记录某个区间的宽度        step_i = 0        while start < 1:            # print step            temp_x = 0.5 * (start + min(start + step, 1))            step_end = min(start + step, 1)            temp_step_x = step_end - start            # 概率密度为：math.exp(self.alpha * temp_x - 1 + self.beta) 取其对数            temp_pro = self.alpha * temp_x - 1 + self.beta            if math.exp(temp_pro)*temp_step_x < 0.00001 and step<0.0005:                step *= 2                continue            else:                x_array[step_i] = temp_x                step_x_array[step_i] = temp_step_x                pro_x_array[step_i] = temp_pro            start =  step_end            step_i += 1        self.x_array = x_array[0:step_i]  # 每种情况的发生概率 以该区间的中值为准        self.pro_x_array = pro_x_array[0:step_i]  # 每种情况的对数概率值        self.step_x_array = step_x_array[0:step_i]  # 每种情况的取值区间大小    def get_pro(self, n, m):        # 计算不同的段x概率在m次试验中,发生n次的概率  使用对数避免精度损失        # self.pro_x_array * self.step_x_array 每种情况发生的原始分布        now_jiashe = np.log(self.x_array) * n + np.log(1 - self.x_array) * (m - n) + self.pro_x_array + np.log(self.step_x_array)  # 调整原始的经验分布        my_max = now_jiashe.max()  # 以最大值为1 进行计算        now_jiashe = now_jiashe - my_max        now_jiashe = np.exp(now_jiashe)/sum(np.exp(now_jiashe))  # 返回为不同情况的 真实的概率分布        return sum(now_jiashe * self.x_array)if __name__ == "__main__":    a = P_test()    a.solve_function(0.001)  #初始经验概率    # print sum(a.array_pro_pro)    a.set_array()    print a.get_pro(0,0) #0次抽样，0次发生    print a.get_pro(10000,10000)  #10000次抽样，10000次发生

阅码指南：

solve_function这个函数是利用随机行走的办法找方程的解，原方程有一个0伪解（alpha趋于0时，方程的误差也趋于零），这里通过初始设置避开了，找到那个非零解。

error_f 则是计算某个尝试解与方程的符合程度。这个里面真正的包含了需要求解的方程

set_array 这个函数则是知道概率密度分布函数后，将其划分为0-1区间各种不同概率的离散情况便于数值计算。这个划分（离散化）方法直接会影响到后面拟合概率的精度。精度的控制由 math.exp(temp_pro)*temp_step_x < 0.00001 and step<0.0005 这条判断语句限定，前者是对概率区间的初始概率进行限制，后者是对概率区间的大小进行限制。总之最后得到概率在0-1之间各种不同情况发生的概率.（精度原因都取了对数），也就是先验概率分布。

get_pro 则是结合m次抽样，发生n次之后这一事件，计算后验概率密度分布，最后返回该概率分布下的概率期望。

至于实现中的一些trick这里就不一一赘述了，希望看到此文的盆友能够在实际中应用一下，比较一下这一方法与你原有的优化方法优劣。

结语：

这个代码在开发的时候我犯了个小错误，导致0次抽样0次发生时返回概率与原假设不符，导致我把公式推导屡了两三遍，最后几近放弃，幸好鬼使神差的发现是个实现中的逻辑漏洞。至于此方法到底有多好，只有实际应用才能证明了。而这一方法凡是涉及概率估算的地方都是可以应用的。

各种错误（漏洞）在所难免，欢迎留言指正。至于小概率这一问题从问题的提出，到解决方法，从方法的理论到最后的实现在这三篇小概率的文章中算是基本都有了，算是讲了一个完整故事。好开心……

不过在前面曾今说过分布的熵其实代表了我们对估出概率的信心，由于这个值具体的实际应用并不多见，因此也没有给与太多重视，不过或许只是我们没有发现使用熵这一信息的方法而已。不过既然最后我们可以得到后验概率分布，就可以轻松地回答落在某概率区间概率（比如，我们可以说红球的概率取值在9%~10%的之间概率为90%，这时就对9.5%的这一估计结果算是比较确信了）。

另外在这里熵最大的分布服从指数衰减的规律（重尾分布），而另一个中随处可见的长尾分布也是经常令人想起指数函数，在物理认知中有一句话：封闭系统中熵最大的状态即为稳定平衡状态（可能熵增原理更容易被搜索到）。熵最大、长尾分布、重尾分布、稳定平衡态这些看似毫不相关的事件背后是否隐藏着不为人知关系。

国庆七天假，本来还想运用条件独立化来做一个贝叶斯分类器，来感受体会一下，结果玩提莫去了……