python - 全局分布拟合共享一些参数而不指定 python 中的 bin 大小

Question

我有几个数据集分别非常适合 vonMises 分布。我正在寻找一种方法，让他们共享 mu 但具有不同的 kappas，而不用关心 bins 的选择。

当一个人只想适应一个模型时，这很简单:scipy here与原始数据相符。但我一直在使用 symfit 寻找全局拟合或 lmfit ，或在某些帖子中(here 和 here)，在所有情况下我们都必须指定 x 坐标和 y 坐标，这意味着之前必须为分布选择一些 bin 大小。

这是一些仅适用于两个数据集的人工数据，可以用作我需要的示例，尽管使用 scipy 分别拟合每个数据集。 (请注意，我不需要关心垃圾箱的选举)。

import numpy as np
import scipy.stats as st
import matplotlib.pyplot as plt

# creating the data
mu1, mu2 = .05, -.05
sigma1, sigma2 = 3.1, 2.9
n1, n2 = 8000, 9000
y1 = np.random.vonmises(mu1, sigma1, n1)
y2 = np.random.vonmises(mu2, sigma2, n2)

# fitting
dist = st.vonmises
*args1, loc1, scale1 = dist.fit(y1, fscale=1)
*args2, loc2, scale2 = dist.fit(y2, fscale=1)

x1 = np.linspace(np.min(y1), np.max(y1), 200)
x2 = np.linspace(np.min(y2), np.max(y2), 200)

pdf_fitted1 = dist.pdf(x1, *args1, loc=loc1, scale=scale1)
pdf_fitted2 = dist.pdf(x2, *args2, loc=loc2, scale=scale2)

# plotting
plt.hist(y1, bins=40, density=True, histtype='step', color='#1f77b4')
plt.hist(y2, bins=40, density=True, histtype='step', color='#ff7f0e')
plt.plot(x1, pdf_fitted1, color='#1f77b4')
plt.plot(x2, pdf_fitted2, color='#ff7f0e')
plt.show()

如果有人可以提供帮助，我会很高兴，在此先感谢。如有任何答复或评论，我们将不胜感激。

Answer 1

谢谢你提出这个很好的问题。原则上你应该能够使用开箱即用的 symfit 来解决这个问题，但是你的问题揭示了一些小问题。 symfit 是一个不断变化的项目，因此不幸的是，这种情况时有发生。但我创建了一个应该适用于当前 master 分支的解决方法，我希望尽快发布一个新版本来解决这些问题。

原则上，这是您已经在 LikeLihood 中找到的全局拟合示例的组合。目标函数。使用 LogLikelihood，您不需要分箱，而是直接使用测量值。但是，LikeLihood 似乎还不能正确处理多分量模型，因此我包含了一个固定版本的 LogLikelihood。

import matplotlib.pyplot as plt
from symfit import (
    Fit, parameters, variables, exp, cos, CallableModel, pi, besseli
)
from symfit.core.objectives import LogLikelihood
from symfit.core.minimizers import *
from symfit.core.printing import SymfitNumPyPrinter

# symbolic bessel is not converted into numerical yet, this monkey-patches it.
def _print_besseli(self, expr):
    return 'scipy.special.iv({}, {})'.format(*expr.args)

SymfitNumPyPrinter._print_besseli = _print_besseli

# creating the data
mu1, mu2 = .05, -.05  # Are these supposed to be opposite sign?
sigma1, sigma2 = 3.5, 2.5
n1, n2 = 8000, 9000
np.random.seed(42)
x1 = np.random.vonmises(mu1, sigma1, n1)
x2 = np.random.vonmises(mu2, sigma2, n2)

# Create a model for `n` different datasets.
n = 2
x, *xs = variables('x,' + ','.join('x_{}'.format(i) for i in range(1, n + 1)))
ys = variables(','.join('y_{}'.format(i) for i in range(1, n + 1)))
mu, kappa = parameters('mu, kappa')
kappas = parameters(','.join('k_{}'.format(i) for i in range(1, n + 1)),
                    min=0, max=10)
mu.min, mu.max = - np.pi, np.pi  # Bound to 2 pi

# Create a model template, who's symbols we will replace for each component.
template = exp(kappa * cos(x - mu)) / (2 * pi * besseli(0, kappa))
model = CallableModel(
    {y_i: template.subs({kappa: k_i, x: x_i}) for y_i, x_i, k_i in zip(ys, xs, kappas)}
)
print(model)

class AlfredosLogLikelihood(LogLikelihood):
    def __call__(self, *args, **kwargs):
        evaluated_func = super(LogLikelihood, self).__call__(
            *args, **kwargs
        )

        ans = - sum([np.nansum(np.log(component))
                     for component in evaluated_func])
        return ans

fit = Fit(model, x_1=x1, x_2=x2, objective=AlfredosLogLikelihood)

x_axis = np.linspace(- np.pi, np.pi, 101)
fit_result = fit.execute()
print(fit_result)
x1_result, x2_result = model(x_1=x_axis, x_2=x_axis, **fit_result.params)

# plotting
plt.hist(x1, bins=40, density=True, histtype='step', color='#1f77b4')
plt.hist(x2, bins=40, density=True, histtype='step', color='#ff7f0e')
plt.plot(x_axis, x1_result, color='#1f77b4')
plt.plot(x_axis, x2_result, color='#ff7f0e')
plt.show()

输出如下:

[y_1(x_1; k_1, mu) = exp(k_1*cos(mu - x_1))/(2*pi*besseli(0, k_1)),
 y_2(x_2; k_2, mu) = exp(k_2*cos(mu - x_2))/(2*pi*besseli(0, k_2))]

Parameter Value        Standard Deviation
k_1       3.431673e+00 None
k_2       2.475649e+00 None
mu        1.030791e-02 None
Status message         b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
Number of iterations   13

我希望这能让您朝着正确的方向前进，感谢您揭露这个错误 ;)。