Я пытаюсь использовать PyMC 2.3 для получения оценки параметра составной модели. Под «составным» я подразумеваю случайную величину, которая следует за распределением, параметром которого является другая случайная величина. («вложенные» или «иерархические» иногда используются для обозначения этой ситуации, но я думаю, что они менее конкретны и в этом контексте вызывают больше путаницы).
Приведем конкретный пример. «Реальные» данные генерируются из составного распределения, которое представляет собой распределение Пуассона с параметром, распределенным как экспоненциальное. В обычном scipy данные генерируются следующим образом:
import numpy as np
from scipy.stats import distributions
np.random.seed(3) # for repeatability
nsamples = 1000
tau_true = 50
orig_lambda_sample = distributions.expon(scale=tau_true).rvs(nsamples)
data = distributions.poisson(orig_lambda_sample).rvs(nsamples)
Я хочу получить оценку параметра модели tau_true. Мой подход к моделированию этой проблемы в PyMC заключается в следующем:
tau = pm.Uniform('tau', 0, 100)
count_rates = pm.Exponential('count_rates', beta=1/tau, size=nsamples)
counts = pm.Poisson('counts', mu=count_rates, value=data, observed=True)
Обратите внимание, что я использую size=nsamples, чтобы иметь новую стохастическую переменную для каждой выборки.
Наконец, я запускаю MCMC как:
model = pm.Model([count_rates, counts, tau])
mcmc = pm.MCMC(model)
mcmc.sample(40000, 10000)
Модель сходится (хотя и медленно, > 10 ^ 5 итераций) к распределению с центром около 50 (tau_true). Однако определение 1000 стохастических переменных для соответствия одному распределению с одним параметром кажется излишним.
Есть ли лучший способ описать составную модель в PyMC?
PS Я также пробовал с более информативным априорным анализом (tau = pm.Normal('tau', mu=51, tau=1/2**2)), но результаты аналогичны, а сходимость по-прежнему медленная.
