Я пытаюсь приспособить обратное гамма-распределение к гистограмме.Я могу легко выбрать набор параметров хорошего соответствия, поиграв с двумя его параметрами (см. Оранжевую линию на Рис. 1).
Однако, когда я пытаюсь оценить потомков, используя конферансье, я просто выбираю начальные начальные точки цепи.
import numpy as np
import corner
import emcee
from scipy.special import gamma as gamma_function
import matplotlib.pyplot as plt
def _inverse_gamma_distribution(x0, alpha, beta):
return beta**alpha/gamma_function(alpha)* (1/x0)**(alpha + 1)* np.exp(-beta/x0)
flux = np.array([0., 0., 0.08753462, 0.48757902, 0.59385076, 0.32146012, 0.20280991, 0.06542532, 0.01888836, 0.00369042, 0.00133481, 0., 0., 0.10504154, 0.44777665])
bin_center = np.array([0.06898463, 0.12137053, 0.21353749, 0.37569469, 0.66099163, 1.16293883, 2.04605725, 3.59980263, 6.33343911, 11.1429583, 19.60475495, 0.06898463, 0.12137053, 0.21353749, 0.37569469])
error = np.array([0., 0., 0.03914667, 0.06965415, 0.0579539, 0.03214601, 0.01924986, 0.00824282, 0.00333902, 0.0011127, 0.0005045, 0., 0., 0.04288303, 0.0667506 ])
def lnprior(theta):
alpha, beta = theta
if 1.0 < alpha < 2.0 and 1.0 < beta < 2.0:
return 0.0
return -np.inf
def lnprob(theta, x, y, yerr):
lp = lnprior(theta)
if not np.isfinite(lp):
return -np.inf
return lp + lnlike(theta, x, y, yerr)
def lnlike(theta, x, y, yerr):
alpha, beta = theta
model = np.array(_inverse_gamma_distribution(x, alpha, beta))
return -0.5*np.sum((y - model)**2/yerr**2)
p0 = [1.5, 1.5]
ndim, nwalkers = 2, 100
pos = [np.array(p0) + 5e-1*np.random.randn(ndim) for i in range(nwalkers)]
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=(bin_center, flux, error))
sampler.run_mcmc(pos, 1000)
samples = sampler.chain[:, 500:, :].reshape((-1, ndim))
print("mean posterior: ", samples.T[0].mean(), samples.T[1].mean())
print("std posterior: ", samples.T[0].std(), samples.T[1].std())
fig = corner.corner(samples)
fig, ax = plt.subplots()
ax.set_xscale("log", nonposx='clip')
ax.set_yscale("log", nonposy='clip')
ax.set_ylim([0.0001, 10])
x = np.linspace(0.1, 20, 1000)
ax.errorbar(bin_center, flux, yerr=error, fmt=".")
ax.plot(x, _inverse_gamma_distribution(x, *p0))
for alpha, beta in samples[np.random.randint(len(samples), size=100)]:
ax.plot(x, _inverse_gamma_distribution(x, alpha, beta), "-", color="grey", alpha=0.1)
plt.show()
Среднее среднее значение и стандартное отклонение соответствуют значениям, которые я установил вpos = [np.array(p0) + 5e-1*np.random.randn(ndim) for i in range(nwalkers)] (то есть [р0 + - 0,5]).Изменение количества шагов и / или количества сожженных образцов не помогает.
Если я увеличу количество ходунков, он просто лучше сэмплирует входную позицию - они никогда не кажутсяиди куда угодно.
Это сообщение ссылка , но я не думаю, что оно применимо здесь - обратное гамма-распределение, возможно, являетсянемного раздражает, но проблема подбора кривой должна быть великолепно определена ..?
------------------- РЕДАКТИРОВАТЬ --------------------
Я знаю, что понимаю проблему: удаление точек данных, где flux равен нулю, заставляет пешеходов начать движение, и я получаю результат,Я нахожу это поведение немного странным, поэтому, чтобы обновить вопрос: , почему два отдельных выброса препятствуют движению конферансье?