如何使用python查找已知数据的多项分布参数？

Posted 2023-03-12

【中文标题】如何使用python查找已知数据的多项分布参数？

【英文标题】：How to find multinomial distribution parameter for a known data using python?

【发布时间】：2018-12-17 19:12:10

【问题描述】：

我刚开始学习 Python。这是一个数据框：

a=pd.DataFrame('A1':[0,1,2,3,2,1,6,0,1,1,7,10])

现在我认为这些数据服从多项分布。因此，12 个数字表示 12 个类别（类别 0、1、2...）的频率。例如，类别 0 的出现是 0。所以，我希望在给定这个数据的情况下找到多项式的所有参数。最后，我们得到了多项式的最佳参数（或者我们可以说每个数字的最佳概率）。例如，

category:    0,      1,     2,     3,      4...
weights:    0.001,  0.1,   0.2,   0.12,   0.2...

所以，我不需要测试数据来预测。这不是一个分类。作为一个新手，我什至不确定是否应该使用 scipy.stats.multinomial 或 sklearn 模型，或其他一些技术。那么，谁能给我一些帮助？

【参考方案1】：

最大似然估计（MLE）是获取分布参数的点估计的最重要的过程之一。这是你需要开始的。

解析解：

跨国分布是二项分布的扩展，MLE可以解析得到。请参阅此数学堆栈交换帖子 (MLE for Multinomial Distribution) 以获得完整的分析解决方案。该过程从定义一个似然函数开始，L(p) 以观察到的数据 x(i) 为条件，其中 p 和 x 是 k 类/类别 和 i= 0,1,...k 的概率和观察到的出现。在给定参数集 (p) 的情况下，它衡量观察一组 观测值 (x) 的可能性：

L(p) 等于：

主要思想是在参数（p）的范围内最大化似然函数值。给定总观察值n（即所有类别的出现总和），点估计等于：

a.values/a.values.sum()                        # point estimates for p = x/n

# array([[0.        ], [0.02941176], [0.05882353], [0.08823529], 
#        [0.05882353], [0.02941176], [0.17647059], [0.        ], 
#        [0.02941176], [0.02941176], [0.20588235], [0.29411765]])

数值解：

上面的结果也可以用scipy.optimize.minimize数值得到。请注意，L(p) 是阶乘项和指数项的乘积。阶乘项是一个常数，不依赖于参数值 (p)，因此不考虑进行优化。对于指数项，最好进行log变换来简化目标函数； MLE 的常见做法，因为 log 是 单调递增 函数。此外，由于scipy.optimize.minimize 用于最小化，我们将使用对数转换似然函数的负值。 注意，最大化函数值等于最小化其负值。

import pandas as pd
import numpy as np
import scipy.optimize as sciopt

# bounds for parameters to lie between (0,1), 
# absolute zero (0) for lower bound avoided as log takes an infinite value 
bnds = [(0.001e-12,1) for i in range(12)]

# Initializing parameters value for optimization
init_parameters = np.asarray([0.1 for i in range(12)])

# Negative Log Likelihood Function
neg_log_lik = lambda p: -np.sum([a.values[i]*np.log(p[i]) for i in range(12)])

# Constraint sum(p) = 1
cons = 'type': 'eq', 'fun': lambda p:  (sum([p[i] for i in range(12)]) - 1) 

# Minimizing neg_log_lik
results = sciopt.minimize(neg_log_lik, x0 = init_parameters, 
                          method='SLSQP', bounds= bnds, constraints= cons)

results.x                                    # point estimates for p

#   array([1.00000000e-15, 2.94179308e-02, 5.88243586e-02, 8.82394605e-02,
#          5.88243586e-02, 2.94059735e-02, 1.76454713e-01, 1.00000000e-15,
#          2.94134577e-02, 2.94135714e-02, 2.05849197e-01, 2.94156978e-01])

参考scipy.optimize.minimize有关上述实现的详细信息的文档。