移动间隔时 pymc 中 DiscreteUniform 的意外行为
Unexpected behaviour with DiscreteUniform in pymc when shifting interval
我正在尝试模拟 100 次掷骰子,其中我的 数据 是所有掷骰的总和(类似于 Jaynes 的 Brandeis 骰子的最大熵原理)。这是我后来第一次尝试接近已加载的骰子。
我正在使用 pymc 2.3
如果我用 DiscreteUniform('dice', 1, 6, size=N)
设置从 1 到 6 的骰子值,并设置一个总和值等于平均总和值 100*3.5=350,那么我得到均匀的后验分布,正如预期的那样。
但是如果我设置骰子值从0到5,总和等于100*2.5=250,分布不均匀。 0 值的采样率要低得多!因为我只是将值移动 1 个单位,所以我期望得到相同的结果。知道为什么它们不同吗?我做错了什么?
这是完整代码:
import numpy as np
import pymc as pm
import matplotlib.pyplot as plt
N = 100
shifts = (0, -1)
for shift in shifts:
obs_mean = 3.5+shift
obs_total = int(N*obs_mean)
sigma = 0.01*N
dice = pm.DiscreteUniform('dice', 1+shift, 6+shift, size=N)
@pm.deterministic
def calc_total(d=dice):
return np.sum(d)
total = pm.Normal('total', mu=calc_total, tau=1./sigma, observed=True, value=obs_total)
# package the full model in a dictionary
model1 = dict(dice=dice, calc_total=calc_total, total=total)
# run the basic MCMC:
S = pm.MCMC(model1)
S.sample(iter=100000, burn=10000)
dice_trace = S.trace('dice')[:]-shift
plt.hist(dice_trace.flat, bins=(0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5), normed=True, alpha=0.5)
plt.show()
编辑:根据评论,我做了一个更简单的模型:两个均匀分布,一个从 1 到 6,另一个从 0 到 5,然后是确定性函数 dice2
将其加 1,因此两个模型中的先验 dice2
相同,并且可能性仅取决于 dice2
,但是它们的后验分布不同。
另一个有趣的情况是当 shift 设置为 -7 时,结果只是反转骰子的符号,但结果不同。
import numpy as np
import pymc as pm
import matplotlib.pyplot as plt
N = 100
shifts = (0, -1)
for shift in shifts:
obs_mean = 3.5
obs_total = int(N*obs_mean)
sigma = 0.01*N
dice = pm.DiscreteUniform('dice', 1+shift, 6+shift, size=N)
@pm.deterministic
def dice2(d=dice):
return d-shift
@pm.deterministic
def calc_total(d=dice2):
return np.sum(d)
total = pm.Normal('total', mu=calc_total, tau=1./sigma, observed=True, value=obs_total)
# package the full model in a dictionary
model1 = dict(dice=dice, dice2=dice2, calc_total=calc_total, total=total)
# run the basic MCMC:
S = pm.MCMC(model1)
S.sample(iter=100000, burn=10000)
dice_trace = S.trace('dice2')[:]
plt.hist(dice_trace.flat, bins=(0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5), normed=True, alpha=0.5)
plt.show()
不清楚为什么您一定会期望均匀分布。离散制服只是你的先验。该模型的所有信息是 shift=0
的 350 和 shift=-1
的 250 的总和,并将生成具有此期望的参数估计。当我 运行 每个偏移值下的模型并查看轨迹时,我得到 shift=0
的以下分布(仅按唯一值查看摘要):
>>> pd.Series(shift_0.flat).value_counts()
2 1526136
4 1526011
3 1511494
5 1503698
6 1471922
1 1460739
对应于以下期望:
>>> pd.Series(shift_0.flat).mean() * 100
350.02311111111112
而 shift=-1
>>> pd.Series(shift_1.flat).value_counts()
1 1894489
2 1724072
3 1577420
4 1457896
5 1320425
0 1025698
dtype: int64
>>> pd.Series(shift_1.flat).mean() * 100
250.08703333333332
因此,该模型的行为似乎符合我的预期。
这是 PyMC 2.3.6 版采样步骤的问题。它在版本 2.3.2 中按预期工作。我在 github 中与 Chris Fonnesbeck 讨论了这个问题,他在 PyMC 的开发版本中更正了它。将来,请检查您拥有的版本及其行为方式。
我正在尝试模拟 100 次掷骰子,其中我的 数据 是所有掷骰的总和(类似于 Jaynes 的 Brandeis 骰子的最大熵原理)。这是我后来第一次尝试接近已加载的骰子。
我正在使用 pymc 2.3
如果我用 DiscreteUniform('dice', 1, 6, size=N)
设置从 1 到 6 的骰子值,并设置一个总和值等于平均总和值 100*3.5=350,那么我得到均匀的后验分布,正如预期的那样。
但是如果我设置骰子值从0到5,总和等于100*2.5=250,分布不均匀。 0 值的采样率要低得多!因为我只是将值移动 1 个单位,所以我期望得到相同的结果。知道为什么它们不同吗?我做错了什么?
这是完整代码:
import numpy as np
import pymc as pm
import matplotlib.pyplot as plt
N = 100
shifts = (0, -1)
for shift in shifts:
obs_mean = 3.5+shift
obs_total = int(N*obs_mean)
sigma = 0.01*N
dice = pm.DiscreteUniform('dice', 1+shift, 6+shift, size=N)
@pm.deterministic
def calc_total(d=dice):
return np.sum(d)
total = pm.Normal('total', mu=calc_total, tau=1./sigma, observed=True, value=obs_total)
# package the full model in a dictionary
model1 = dict(dice=dice, calc_total=calc_total, total=total)
# run the basic MCMC:
S = pm.MCMC(model1)
S.sample(iter=100000, burn=10000)
dice_trace = S.trace('dice')[:]-shift
plt.hist(dice_trace.flat, bins=(0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5), normed=True, alpha=0.5)
plt.show()
编辑:根据评论,我做了一个更简单的模型:两个均匀分布,一个从 1 到 6,另一个从 0 到 5,然后是确定性函数 dice2
将其加 1,因此两个模型中的先验 dice2
相同,并且可能性仅取决于 dice2
,但是它们的后验分布不同。
另一个有趣的情况是当 shift 设置为 -7 时,结果只是反转骰子的符号,但结果不同。
import numpy as np
import pymc as pm
import matplotlib.pyplot as plt
N = 100
shifts = (0, -1)
for shift in shifts:
obs_mean = 3.5
obs_total = int(N*obs_mean)
sigma = 0.01*N
dice = pm.DiscreteUniform('dice', 1+shift, 6+shift, size=N)
@pm.deterministic
def dice2(d=dice):
return d-shift
@pm.deterministic
def calc_total(d=dice2):
return np.sum(d)
total = pm.Normal('total', mu=calc_total, tau=1./sigma, observed=True, value=obs_total)
# package the full model in a dictionary
model1 = dict(dice=dice, dice2=dice2, calc_total=calc_total, total=total)
# run the basic MCMC:
S = pm.MCMC(model1)
S.sample(iter=100000, burn=10000)
dice_trace = S.trace('dice2')[:]
plt.hist(dice_trace.flat, bins=(0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5), normed=True, alpha=0.5)
plt.show()
不清楚为什么您一定会期望均匀分布。离散制服只是你的先验。该模型的所有信息是 shift=0
的 350 和 shift=-1
的 250 的总和,并将生成具有此期望的参数估计。当我 运行 每个偏移值下的模型并查看轨迹时,我得到 shift=0
的以下分布(仅按唯一值查看摘要):
>>> pd.Series(shift_0.flat).value_counts()
2 1526136
4 1526011
3 1511494
5 1503698
6 1471922
1 1460739
对应于以下期望:
>>> pd.Series(shift_0.flat).mean() * 100
350.02311111111112
而 shift=-1
>>> pd.Series(shift_1.flat).value_counts()
1 1894489
2 1724072
3 1577420
4 1457896
5 1320425
0 1025698
dtype: int64
>>> pd.Series(shift_1.flat).mean() * 100
250.08703333333332
因此,该模型的行为似乎符合我的预期。
这是 PyMC 2.3.6 版采样步骤的问题。它在版本 2.3.2 中按预期工作。我在 github 中与 Chris Fonnesbeck 讨论了这个问题,他在 PyMC 的开发版本中更正了它。将来,请检查您拥有的版本及其行为方式。