如何实施最大似然估计类型2？

Question

如何实施最大似然估计类型2？

pythonscipysympybayesianpymc3

3

我正在尝试实现一个经验贝叶斯 ML-II（最大似然估计类型 II）方法，用于从历史数据中估计先验分布参数。

其中：

π(θ) 是先验分布的表达式
p(x|θ) 是数据分布的表达式
m(x) 是边缘分布的表达式

根据步骤，我需要首先进行积分，找到边缘分布的表达式，然后找到该表达式的极值以估计先验分布的参数。可以使用诸如 scipy.optimize 等方法来实现极值。那么问题是如何进行积分？

enter image description here

- abraxas

@merv 我们亲爱的Merv，谢谢你。我知道像pymc3这样的贝叶斯模型包。我问过他们，他们说经验贝叶斯不是一种通用方法，他们没有实现ML-II。特别是，这个ML-II不是一个MAP。看起来pymc3的MAP方法不执行ML-II。 - abraxas

也许 symfit 能够帮助你，看看这个文档中的例子。然后，你就能够使用 sympy 风格的解析表达式，并且使用 scipy 将它们拟合到你的数据中，但是不必与 scipy 交互。注意：我是 symfit 的作者。 - tBuLi

@tBuLi，我已经阅读了你的链接，看起来我仍然需要自己编写模型公式，这恰恰是我遇到的问题，我需要先找到边缘分布的表达式以进行整合。如果我要使用脑力积分来解决这个问题，那就有麻烦了。我希望能够编写一些通用代码来整合这个表达式。 - abraxas

我想确认一下，您的模型是否有一个分析表达式，但是希望有一些自动化的方法将其集成到边缘分布中，并将其拟合到数据集中？如果是这种情况，请告诉我，我将为您编写一个使用symfit的示例。 - tBuLi

@tBuLi，你的理解令我钦佩！事实上，估计先验分布参数的方法是对先验分布乘以似然函数进行积分，然后找出能够最大化边缘分布值的先验分布参数类型。在得到边缘分布的表达式后，将使用这些参数作为先验分布的参数进行最大后验估计。 - abraxas

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- tBuLi · Accepted Answer

这里使用symfit试着解决这个问题。我选择从一个没有协方差的二元正态分布中进行抽样作为示例。"最初的回答"。

import numpy as np
import matplotlib.pyplot as plt
from symfit import Model, Fit, Parameter, Variable, integrate, oo
from symfit.distributions import Gaussian
from symfit.core.objectives import LogLikelihood

# Make variables and parameters
x = Variable('x')
y = Variable('y')
m = Variable('m')
x0 = Parameter('x0', value=0.6, min=0.5, max=0.7)
sig_x = Parameter('sig_x', value=0.1)
y0 = Parameter('y0', value=0.7, min=0.6, max=0.9)
sig_y = Parameter('sig_y', value=0.05)

pdf = Gaussian(x=x, mu=x0, sig=sig_x) * Gaussian(x=y, mu=y0, sig=sig_y)
marginal = integrate(pdf, (y, -oo, oo), conds='none')
print(pdf)
print(marginal)

model = Model({m: marginal})

# Draw 10000 samples from a bivariate distribution
mean = [0.59, 0.8]
cov = [[0.11**2, 0], [0, 0.23**2]]
xdata, ydata = np.random.multivariate_normal(mean, cov, 10000).T

# We provide only xdata to the model
fit = Fit(model, xdata, objective=LogLikelihood)
fit_result = fit.execute()
print(fit_result)

xaxis = np.linspace(0, 1.0)
plt.hist(xdata, bins=100, density=True)
plt.plot(xaxis, model(x=xaxis, **fit_result.params).m)
plt.show()

这将为PDF和边际分布打印以下内容：

这将在PDF和边际分布中打印以下内容：

>>> exp(-(-x0 + x)**2/(2*sig_x**2))*exp(-(-y0 + y)**2/(2*sig_y**2))/(2*pi*Abs(sig_x)*Abs(sig_y))
>>> sqrt(2)*sig_y*exp(-(-x0 + x)**2/(2*sig_x**2))/(2*sqrt(pi)*Abs(sig_x)*Abs(sig_y))

And for the fit results:

Parameter Value        Standard Deviation
sig_x     1.089585e-01 7.704533e-04
sig_y     5.000000e-02 nan
x0        5.905688e-01 -0.000000e+00
Fitting status message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
Number of iterations:   9
Regression Coefficient: nan

你可以看到x0和sig_x已经被正确获得，但是关于y参数的信息无法获得。我认为在这个例子中这是有道理的，因为没有相关性，但是我会让你自己去处理那些细节问题;)。

"最初的回答"