python - 在python中为一对分布生成MLE

好的，所以我当前的曲线拟合代码有一个步骤，该步骤使用scipy.stats根据数据确定正确的分布，

distributions = [st.laplace, st.norm, st.expon, st.dweibull, st.invweibull, st.lognorm, st.uniform]
mles = []

for distribution in distributions:
    pars = distribution.fit(data)
    mle = distribution.nnlf(pars, data)
    mles.append(mle)

results = [(distribution.name, mle) for distribution, mle in zip(distributions, mles)]

for dist in sorted(zip(distributions, mles), key=lambda d: d[1]):
    print dist
best_fit = sorted(zip(distributions, mles), key=lambda d: d[1])[0]
print 'Best fit reached using {}, MLE value: {}'.format(best_fit[0].name, best_fit[1])


print [mod[0].name for mod in sorted(zip(distributions, mles), key=lambda d: d[1])]

其中data是数值列表。到目前为止，这对于拟合单峰分布非常有效，已在脚本中确认，该脚本可从随机分布中随机生成值，并使用curve_fit重新确定参数。

现在，我想使代码能够处理双峰分布，如以下示例所示:

是否可以从scipy.stats中获取一对模型的MLE，以便确定特定的分布对是否适合数据？

distributions = [st.laplace, st.norm, st.expon, st.dweibull, st.invweibull, st.lognorm, st.uniform]
distributionPairs = [[modelA.name, modelB.name] for modelA in distributions for modelB in distributions]

并使用这些对获取适合数据的那对分布的MLE值？

最佳答案

这不是一个完整的答案，但它可以帮助您解决问题。假设您知道您的问题是由两种密度产生的。
一种解决方案是使用 k均值或EM算法。

初始化。
您可以通过将每个观测值影响到一个或另一个密度来初始化算法。然后初始化两个密度(初始化密度的参数，本例中的参数之一是“高斯”，“拉普拉斯”，依此类推...
迭代。
然后，反复执行以下两个步骤:

步骤1.
假设每个点的影响都是正确的，则优化参数。现在，您可以使用任何优化求解器。此步骤为您提供了适合您数据的最佳两个密度(具有给定参数)的估计。

步骤2。
您可以根据最大的可能性将每个观测分类为一个密度或另一个密度。

您重复直到收敛。

此网页对此进行了很好的解释
https://people.duke.edu/~ccc14/sta-663/EMAlgorithm.html

如果您不知道生成数据的密度是多少，该问题将更加棘手。您必须处理受罚的分类问题，这有点困难。

这是一个简单的编码示例:您知道您的数据来自2个不同的高斯分布(您不知道每个密度生成多少个变量)。在您的情况下，您可以调整此代码以在每对可能的密度对上循环(通常更长，但是根据经验，我认为可以)

import scipy.stats as st
import numpy as np

#hard coded data generation
data = np.random.normal(-3, 1, size = 1000)
data[600:] = np.random.normal(loc = 3, scale = 2, size=400)

#initialization

mu1 = -1
sigma1 = 1

mu2 = 1
sigma2 = 1

#criterion to stop iteration
epsilon = 0.1
stop = False

while  not stop :
    #step1
    classification = np.zeros(len(data))
    classification[st.norm.pdf(data, mu1, sigma1) > st.norm.pdf(data, mu2, sigma2)] = 1

    mu1_old, mu2_old, sigma1_old, sigma2_old = mu1, mu2, sigma1, sigma2

    #step2
    pars1 = st.norm.fit(data[classification == 1])
    mu1, sigma1 = pars1
    pars2 = st.norm.fit(data[classification == 0])
    mu2, sigma2 = pars2

    #stopping criterion
    stop = ((mu1_old - mu1)**2 + (mu2_old - mu2)**2 +(sigma1_old - sigma1)**2 +(sigma2_old - sigma2)**2) < epsilon

#result
print("The first density is gaussian :", mu1, sigma1)
print("The first density is gaussian :", mu2, sigma2)
print("A rate of ", np.mean(classification), "is classified in the first density")

希望能帮助到你。

关于python - 在python中为一对分布生成MLE，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/42033893/