Pymc 大小/索引问题
Pymc size / indexing issue
我正在尝试使用 pymc 2.3.5 为 Kruschke 的 "filtration-condensation experiment" 建模。 (numpy 1.10.1)
基本上有：
- 4组
- 每组40人
- 每个人有 64 次伯努利试验 (correct/incorrect)
我正在建模的内容：
每个人的结果都是二项分布（例如 64 人中有 45 人是正确的）。
我相信每个人的表现都是 Beta 分布。
这个 Beta 分布受个体所属群体的影响（通过参数 A=mu*kappa 和 B=(1-mu)*kappa）
我认为每个群体的影响有多强是 Gamma 分布（kappa 变量）
我相信每个组的平均值是 Beta 分布（mu 变量）
问题：
- 当我使用 "size=" 参数建模时，pymc get 丢失了
- 当我手动分离每个发行版时（无大小=），pymc 做得很好
我包含以下代码：
数据
import numpy as np
import seaborn as sns
import pymc as pm
from pymc.Matplot import plot as mcplot
%matplotlib inline
# Data
ncond = 4
nSubj = 40
trials = 64
N = np.repeat([trials], (ncond * nSubj))
z = np.array([45, 63, 58, 64, 58, 63, 51, 60, 59, 47, 63, 61, 60, 51, 59, 45,
61, 59, 60, 58, 63, 56, 63, 64, 64, 60, 64, 62, 49, 64, 64, 58, 64, 52, 64, 64,
64, 62, 64, 61, 59, 59, 55, 62, 51, 58, 55, 54, 59, 57, 58, 60, 54, 42, 59, 57,
59, 53, 53, 42, 59, 57, 29, 36, 51, 64, 60, 54, 54, 38, 61, 60, 61, 60, 62, 55,
38, 43, 58, 60, 44, 44, 32, 56, 43, 36, 38, 48, 32, 40, 40, 34, 45, 42, 41, 32,
48, 36, 29, 37, 53, 55, 50, 47, 46, 44, 50, 56, 58, 42, 58, 54, 57, 54, 51, 49,
52, 51, 49, 51, 46, 46, 42, 49, 46, 56, 42, 53, 55, 51, 55, 49, 53, 55, 40, 46,
56, 47, 54, 54, 42, 34, 35, 41, 48, 46, 39, 55, 30, 49, 27, 51, 41, 36, 45, 41,
53, 32, 43, 33])
condition = np.repeat([0,1,2,3], nSubj)
不起作用
# modeling
mu = pm.Beta('mu', 1, 1, size=ncond)
kappa = pm.Gamma('gamma', 1, 0.1, size=ncond)
# Prior
theta = pm.Beta('theta', mu[condition] * kappa[condition], (1 - mu[condition]) * kappa[condition], size=len(z))
# likelihood
y = pm.Binomial('y', p=theta, n=N, value=z, observed=True)
# model
model = pm.Model([mu, kappa, theta, y])
mcmc = pm.MCMC(model)
#mcmc.use_step_method(pm.Metropolis, mu)
#mcmc.use_step_method(pm.Metropolis, theta)
#mcmc.assign_step_methods()
mcmc.sample(100000, burn=20000, thin=3)
# outputs never converge and does vary in new simulations
mcplot(mcmc.trace('mu'), common_scale=False)
有效
z1 = z[:40]
z2 = z[40:80]
z3 = z[80:120]
z4 = z[120:]
Nv = N[:40]
mu1 = pm.Beta('mu1', 1, 1)
mu2 = pm.Beta('mu2', 1, 1)
mu3 = pm.Beta('mu3', 1, 1)
mu4 = pm.Beta('mu4', 1, 1)
kappa1 = pm.Gamma('gamma1', 1, 0.1)
kappa2 = pm.Gamma('gamma2', 1, 0.1)
kappa3 = pm.Gamma('gamma3', 1, 0.1)
kappa4 = pm.Gamma('gamma4', 1, 0.1)
# Prior
theta1 = pm.Beta('theta1', mu1 * kappa1, (1 - mu1) * kappa1, size=len(Nv))
theta2 = pm.Beta('theta2', mu2 * kappa2, (1 - mu2) * kappa2, size=len(Nv))
theta3 = pm.Beta('theta3', mu3 * kappa3, (1 - mu3) * kappa3, size=len(Nv))
theta4 = pm.Beta('theta4', mu4 * kappa4, (1 - mu4) * kappa4, size=len(Nv))
# likelihood
y1 = pm.Binomial('y1', p=theta1, n=Nv, value=z1, observed=True)
y2 = pm.Binomial('y2', p=theta2, n=Nv, value=z2, observed=True)
y3 = pm.Binomial('y3', p=theta3, n=Nv, value=z3, observed=True)
y4 = pm.Binomial('y4', p=theta4, n=Nv, value=z4, observed=True)
# model
model = pm.Model([mu1, kappa1, theta1, y1, mu2, kappa2, theta2, y2,
mu3, kappa3, theta3, y3, mu4, kappa4, theta4, y4])
mcmc = pm.MCMC(model)
#mcmc.use_step_method(pm.Metropolis, mu)
#mcmc.use_step_method(pm.Metropolis, theta)
#mcmc.assign_step_methods()
mcmc.sample(100000, burn=20000, thin=3)
# outputs converge and are not too much different in every simulation
mcplot(mcmc.trace('mu1'), common_scale=False)
mcplot(mcmc.trace('mu2'), common_scale=False)
mcplot(mcmc.trace('mu3'), common_scale=False)
mcplot(mcmc.trace('mu4'), common_scale=False)
mcmc.summary()
谁能给我解释一下为什么 mu[condition] 和 gamma[condition] 不起作用？ :)
我想问题是不将 theta 拆分为不同的变量，但无法理解为什么，也许有一种方法可以将形状参数传递给 size= on theta？
首先，我可以确认第一个版本不会导致稳定的结果。我无法确认的是第二个要好得多；我在第二个代码中也看到了非常不同的结果，对于不同的运行，第一个 mu 参数的值在 0.17 和 0.9 之间变化。
可以通过为马尔可夫链使用良好的起始值来解决收敛问题。这可以通过首先进行最大后验 (MAP) 估计，然后从那里启动马尔可夫链来完成。 MAP 步骤的计算成本低，并且会导致收敛的马尔可夫链，对于您的代码的两个变体都具有可重现的结果。作为参考和比较：我看到的四个 mu 参数的值在 0.94 / 0.86 / 0.72 和 0.71 左右。
您可以通过在使用 "model=pm.Model(..." 定义模型的行之后插入以下两行代码来进行 MAP 估计：
map_ = pm.MAP(model)
map_.fit()
Cameron Davidson-Pilon 的 Bayesian Methods for Hackers 以及有关 PyMC 的其他有用主题对这项技术进行了更详细的介绍。
我正在尝试使用 pymc 2.3.5 为 Kruschke 的 "filtration-condensation experiment" 建模。 (numpy 1.10.1) 基本上有：
- 4组
- 每组40人
- 每个人有 64 次伯努利试验 (correct/incorrect)
我正在建模的内容：
每个人的结果都是二项分布（例如 64 人中有 45 人是正确的）。
我相信每个人的表现都是 Beta 分布。
这个 Beta 分布受个体所属群体的影响（通过参数 A=mu*kappa 和 B=(1-mu)*kappa）
我认为每个群体的影响有多强是 Gamma 分布（kappa 变量）
我相信每个组的平均值是 Beta 分布（mu 变量）
问题：
- 当我使用 "size=" 参数建模时，pymc get 丢失了
- 当我手动分离每个发行版时（无大小=），pymc 做得很好
我包含以下代码：
数据
import numpy as np
import seaborn as sns
import pymc as pm
from pymc.Matplot import plot as mcplot
%matplotlib inline
# Data
ncond = 4
nSubj = 40
trials = 64
N = np.repeat([trials], (ncond * nSubj))
z = np.array([45, 63, 58, 64, 58, 63, 51, 60, 59, 47, 63, 61, 60, 51, 59, 45,
61, 59, 60, 58, 63, 56, 63, 64, 64, 60, 64, 62, 49, 64, 64, 58, 64, 52, 64, 64,
64, 62, 64, 61, 59, 59, 55, 62, 51, 58, 55, 54, 59, 57, 58, 60, 54, 42, 59, 57,
59, 53, 53, 42, 59, 57, 29, 36, 51, 64, 60, 54, 54, 38, 61, 60, 61, 60, 62, 55,
38, 43, 58, 60, 44, 44, 32, 56, 43, 36, 38, 48, 32, 40, 40, 34, 45, 42, 41, 32,
48, 36, 29, 37, 53, 55, 50, 47, 46, 44, 50, 56, 58, 42, 58, 54, 57, 54, 51, 49,
52, 51, 49, 51, 46, 46, 42, 49, 46, 56, 42, 53, 55, 51, 55, 49, 53, 55, 40, 46,
56, 47, 54, 54, 42, 34, 35, 41, 48, 46, 39, 55, 30, 49, 27, 51, 41, 36, 45, 41,
53, 32, 43, 33])
condition = np.repeat([0,1,2,3], nSubj)
不起作用
# modeling
mu = pm.Beta('mu', 1, 1, size=ncond)
kappa = pm.Gamma('gamma', 1, 0.1, size=ncond)
# Prior
theta = pm.Beta('theta', mu[condition] * kappa[condition], (1 - mu[condition]) * kappa[condition], size=len(z))
# likelihood
y = pm.Binomial('y', p=theta, n=N, value=z, observed=True)
# model
model = pm.Model([mu, kappa, theta, y])
mcmc = pm.MCMC(model)
#mcmc.use_step_method(pm.Metropolis, mu)
#mcmc.use_step_method(pm.Metropolis, theta)
#mcmc.assign_step_methods()
mcmc.sample(100000, burn=20000, thin=3)
# outputs never converge and does vary in new simulations
mcplot(mcmc.trace('mu'), common_scale=False)
有效
z1 = z[:40]
z2 = z[40:80]
z3 = z[80:120]
z4 = z[120:]
Nv = N[:40]
mu1 = pm.Beta('mu1', 1, 1)
mu2 = pm.Beta('mu2', 1, 1)
mu3 = pm.Beta('mu3', 1, 1)
mu4 = pm.Beta('mu4', 1, 1)
kappa1 = pm.Gamma('gamma1', 1, 0.1)
kappa2 = pm.Gamma('gamma2', 1, 0.1)
kappa3 = pm.Gamma('gamma3', 1, 0.1)
kappa4 = pm.Gamma('gamma4', 1, 0.1)
# Prior
theta1 = pm.Beta('theta1', mu1 * kappa1, (1 - mu1) * kappa1, size=len(Nv))
theta2 = pm.Beta('theta2', mu2 * kappa2, (1 - mu2) * kappa2, size=len(Nv))
theta3 = pm.Beta('theta3', mu3 * kappa3, (1 - mu3) * kappa3, size=len(Nv))
theta4 = pm.Beta('theta4', mu4 * kappa4, (1 - mu4) * kappa4, size=len(Nv))
# likelihood
y1 = pm.Binomial('y1', p=theta1, n=Nv, value=z1, observed=True)
y2 = pm.Binomial('y2', p=theta2, n=Nv, value=z2, observed=True)
y3 = pm.Binomial('y3', p=theta3, n=Nv, value=z3, observed=True)
y4 = pm.Binomial('y4', p=theta4, n=Nv, value=z4, observed=True)
# model
model = pm.Model([mu1, kappa1, theta1, y1, mu2, kappa2, theta2, y2,
mu3, kappa3, theta3, y3, mu4, kappa4, theta4, y4])
mcmc = pm.MCMC(model)
#mcmc.use_step_method(pm.Metropolis, mu)
#mcmc.use_step_method(pm.Metropolis, theta)
#mcmc.assign_step_methods()
mcmc.sample(100000, burn=20000, thin=3)
# outputs converge and are not too much different in every simulation
mcplot(mcmc.trace('mu1'), common_scale=False)
mcplot(mcmc.trace('mu2'), common_scale=False)
mcplot(mcmc.trace('mu3'), common_scale=False)
mcplot(mcmc.trace('mu4'), common_scale=False)
mcmc.summary()
谁能给我解释一下为什么 mu[condition] 和 gamma[condition] 不起作用？ :)
我想问题是不将 theta 拆分为不同的变量，但无法理解为什么，也许有一种方法可以将形状参数传递给 size= on theta？
首先，我可以确认第一个版本不会导致稳定的结果。我无法确认的是第二个要好得多；我在第二个代码中也看到了非常不同的结果，对于不同的运行，第一个 mu 参数的值在 0.17 和 0.9 之间变化。
可以通过为马尔可夫链使用良好的起始值来解决收敛问题。这可以通过首先进行最大后验 (MAP) 估计，然后从那里启动马尔可夫链来完成。 MAP 步骤的计算成本低，并且会导致收敛的马尔可夫链，对于您的代码的两个变体都具有可重现的结果。作为参考和比较：我看到的四个 mu 参数的值在 0.94 / 0.86 / 0.72 和 0.71 左右。
您可以通过在使用 "model=pm.Model(..." 定义模型的行之后插入以下两行代码来进行 MAP 估计：
map_ = pm.MAP(model)
map_.fit()
Cameron Davidson-Pilon 的 Bayesian Methods for Hackers 以及有关 PyMC 的其他有用主题对这项技术进行了更详细的介绍。