好的,所以我当前的曲线拟合代码有一个步骤使用 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])]
其中数据是数值列表。到目前为止,这对于拟合单峰分布非常有效,在从随机分布中随机生成值并使用 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-mean 或 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/