7

私は、正規分布からの小さなサンプルを持つ多数の個人を含む単純な階層モデルを持っています。これらの分布の平均も正規分布に従います。

import numpy as np

n_individuals = 200
points_per_individual = 10
means = np.random.normal(30, 12, n_individuals)
y = np.random.normal(means, 1, (points_per_individual, n_individuals))

PyMC3 を使用して、サンプルからモデル パラメーターを計算したいと考えています。

import pymc3 as pm
import matplotlib.pyplot as plt

model = pm.Model()
with model:
    model_means = pm.Normal('model_means', mu=35, sd=15)

    y_obs = pm.Normal('y_obs', mu=model_means, sd=1, shape=n_individuals, observed=y)

    trace = pm.sample(1000)

pm.traceplot(trace[100:], vars=['model_means'])
plt.show()

mcmc サンプル

の事後model_means分布が、元の平均分布のようになることを期待していました。しかし、それ30は手段の平均に収束するようです。pymc3 モデルから平均の元の標準偏差 (私の例では 12) を回復するにはどうすればよいですか?

4

1 に答える 1

7

この質問は、私が PyMC3 の概念に苦労していたものです。

n_individualsをモデル化するには観測確率変数が必要で、 をモデル化するyにはn_individual確率変数が必要meansです。hyper_meanこれらには、事前確率とhyper_sigmaそのパラメーターも必要です。sigmasの標準偏差の事前分布ですy。

import matplotlib.pyplot as plt

model = pm.Model()
with model:
    hyper_mean = pm.Normal('hyper_mean', mu=0, sd=100)
    hyper_sigma = pm.HalfNormal('hyper_sigma', sd=3)

    means = pm.Normal('means', mu=hyper_mean, sd=hyper_sigma, shape=n_individuals)
    sigmas = pm.HalfNormal('sigmas', sd=100)

    y = pm.Normal('y', mu=means, sd=sigmas, observed=y)

    trace = pm.sample(10000)

pm.traceplot(trace[100:], vars=['hyper_mean', 'hyper_sigma', 'means', 'sigmas'])
plt.show()

後部

于 2015-11-13T09:09:13.210 に答える