Как реализовать тип оценки максимального правдоподобия 2? - PullRequest
Новая
       50

Как реализовать тип оценки максимального правдоподобия 2?

0 голосов
/ 03 марта 2019

Я пытаюсь реализовать эмпирический байесовский метод ML-II (метод оценки максимального правдоподобия, тип II) для оценки предыдущих параметров распределения по историческим данным

Где ：

  1. π (θ) является выражением для предыдущего распределения
  2. p (x | θ) является выражением для распределения данных
  3. m (x) является выражением для предельного распределения

В соответствии с шагами, мне нужно сначала интегрировать, чтобы найти выражение предельного распределения, а затем найти экстремальное значение этого выражения, чтобы оценить параметры предыдущего распределения.Экстремальные значения могут быть достигнуты с использованием таких методов, как scipy.optimize.Итак, вопрос в том, как мы интегрируем это?

enter image description here

1 Ответ

0 голосов
/ 06 марта 2019

Вот пример использования 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 и предельного распределения: 

>>> 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))

И для подгонкирезультаты: 

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

enter image description here

Вы можете видеть, что x0 и sig_x получены правильно, но никакой информации опараметр, чтобы сделать с y.Я думаю, что в этом примере это имеет смысл, поскольку здесь нет корреляции, но я оставлю вас в борьбе с такими деталями;).

...