移入时pymc中离散形式的意外行为

Question

我试图模拟一个骰子掷100次，其中我的数据是所有掷骰子的总和（类似于Jaynes'Brandeis骰子的最大熵原理）。这是我后来第一次尝试接近一个装了骰子的人。你知道吗

我正在使用pymc2.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()

Answer 1

这是PyMC版本2.3.6中采样步骤的问题。它在版本2.3.2中的工作与预期一样。我在github中与Chris Fonnesbeck讨论了这个问题，他在PyMC的开发版本中纠正了这个问题。以后，请检查您的版本以及它的行为。你知道吗

Answer 2

不清楚为什么你一定会期望一个统一的分布。离散制服只是你的首选。该模型所具有的所有信息是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

因此，模型的行为似乎和我预期的一样。你知道吗