pymc3:具有多个观察变量的层次模型
pymc3: hierarchical model with multiple obsesrved variables
我有一个简单的分层模型,其中有很多个体,我有来自正态分布的小样本。这些分布的均值也服从正态分布。
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()
我期待 model_means
的后验结果看起来像我原来的均值分布。但它似乎收敛到 30
均值。如何从 pymc3 模型中恢复均值的原始标准差(在我的示例中为 12)?
这个问题让我纠结于 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()
我有一个简单的分层模型,其中有很多个体,我有来自正态分布的小样本。这些分布的均值也服从正态分布。
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()
我期待 model_means
的后验结果看起来像我原来的均值分布。但它似乎收敛到 30
均值。如何从 pymc3 模型中恢复均值的原始标准差(在我的示例中为 12)?
这个问题让我纠结于 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()