我正在玩 PyMC3，试图在 PyMC3 文档中适应采矿灾难切换点模型的修改版本。假设您有两个煤矿（mine1 和 mine2），每个煤矿在相同的年份范围内具有相似的灾难计数。
然而，mine1 在实施降低灾难次数的安全程序变更方面晚了 5 年：
import numpy as np
import matplotlib.pyplot as plt
mine1=np.array([0,4,5,4,0,1,4,3,4,0,6,3,3,4,0,2,6,3,3,5,4,5,3,1,4,4,1,5,5,3,4,2,5,2,2,3,
4,2,1,3,0,2,1,1,1,1,3,0,0,1,0,1,1,0,0,3,1,0,3,2,2,0,1,1,1,0,1,0,1,0,0,0,
2,1,0,0,0,1,1,0,2,3,3,1,0,2,1,1,1,1,2,4,2,0,0,1,4,0,0,0,1]);
mine2=np.array([3,3,4,0,2,6,2,3,4,3,7,4,1,5,4,1,5,5,3,4,1,6,2,2,2,4,4,0,4,0,3,3,1,0,3,2,
0,0,1,0,1,1,0,0,3,0,0,3,1,1,0,1,1,1,0,0,0,0,1,1,1,3,1,0,1,0,0,2,0,1,2,2,
0,0,3,3,0,2,3,2,4,2,0,0,1,3,0,0,1,2,0,1,1,0,0,2,0,2,0,0,0]);
both_mines = mine1+mine2;
years = np.arange(1849,1950);
fig, axs = plt.subplots(2);
axs[0].plot(years, both_mines,'ko');
axs[0].legend(['mines_summed'],loc='upper right');
axs[0].set_ylabel('disaster count')
axs[1].plot(years, mine1,'ro');
axs[1].plot(years, mine2,'bo');
axs[1].legend(['mine1','mine2'],loc='upper right');
axs[1].set_ylabel('disaster count')
我有兴趣测试更好的模型拟合是否来自对年度计数求和并将单个切换点拟合到这个总计数时间序列，或者将单独的模型拟合到两个矿山。
模型 1 - 跨矿总和的单一模型
import pymc3 as pm
with pm.Model() as model1:
switchpoint = pm.DiscreteUniform('switchpoint', lower=years.min(), upper=years.max());
early_rate = pm.Exponential('early_rate', 1)
late_rate = pm.Exponential('late_rate', 1)
rate = pm.math.switch(switchpoint >= years, early_rate, late_rate)
disasters_both_mines = pm.Poisson('disasters_both_mines', rate, observed=both_mines)
trace1 = pm.sample(10000,tune=2000);
pm.traceplot(trace1)
Yields 与文档示例非常相似。这是跟踪图：
在拟合使矿井分开的模型时，我尝试了两种由于不同原因都不是最佳的方法。第一个是拟合两个数据的可能性，分别为每个地雷。
模型 2a - 单独的地雷，两种可能性
with pm.Model() as model2a:
switchpoint_mine1 = pm.DiscreteUniform('switchpoint_mine1', lower=years.min(), upper=years.max());
switchpoint_mine2 = pm.DiscreteUniform('switchpoint_mine2', lower=years.min(), upper=years.max());
early_rate_sep = pm.Exponential('early_rate2', 1,shape=2)
late_rate_sep = pm.Exponential('late_rate2', 1,shape=2)
rate_mine1 = pm.math.switch(switchpoint_mine1>=years, early_rate_sep[0], late_rate_sep[0]);
rate_mine2 = pm.math.switch(switchpoint_mine2>=years, early_rate_sep[1], late_rate_sep[1]);
disasters_mine1 = pm.Poisson('disasters_mine1', rate_mine1, observed=mine1);
disasters_mine2 = pm.Poisson('disasters_mine2', rate_mine2, observed=mine2);
trace2a = pm.sample(10000,tune=2000);
pm.traceplot(trace2a);
合身看起来不错，似乎对开关点的差异很敏感。但是，我无法计算 WAIC 或 LOO 值，这意味着我无法比较模型 1 的拟合度。我猜是因为有两组观察结果？
例如
pm.waic(trace2a)
Traceback (most recent call last):
File "<ipython-input-270-122a6fb53049>", line 1, in <module>
pm.waic(trace2a)
File "<home dir>/opt/anaconda3/lib/python3.7/site-packages/pymc3/stats/__init__.py", line 24, in wrapped
return func(*args, **kwargs)
File "<home dir>/opt/anaconda3/lib/python3.7/site-packages/arviz/stats/stats.py", line 1164, in waic
raise TypeError("Data must include log_likelihood in sample_stats")
TypeError: Data must include log_likelihood in sample_stats
第二个想法是使用与分层线性回归示例类似的方法，并在先验上使用串联、索引和形状输出的组合，以拟合每个参数的向量和单个数据似然度。
模型 2b - 单独索引的地雷，单一似然函数
mine1_ind = np.ones(101,dtype=int)-1
mine2_ind = np.ones(101,dtype=int)*1
mine_ix = np.concatenate((mine1_ind,mine2_ind), axis=0);
concat_mines = np.concatenate((mine1,mine2), axis=0);
concat_years = np.transpose(np.concatenate((years,years), axis=0));
with pm.Model() as model2b:
switchpoint_mine1and2 = pm.DiscreteUniform('switchpoint_mine1and2', lower=years.min(), upper=years.max(),shape=2);
early_rate_mine1and2 = pm.Exponential('early_rate_mine1and2', 1,shape=2);
late_rate_mine1and2 = pm.Exponential('late_rate_mine1and2', 1,shape=2);
rate_mine1and2 = pm.math.switch(switchpoint_mine1and2[mine_ix]>=concat_years[mine_ix], early_rate_mine1and2[mine_ix], late_rate_mine1and2[mine_ix]);
disasters_mine1and2 = pm.Poisson('disasters_mine1and2', rate_mine1and2, observed=concat_mines);
trace2b = pm.sample(10000,tune=2000);
该模型适合并允许计算 WAIC。但是从后面看，它不适合切换点。
总而言之，有没有办法以允许计算 WAIC 的方式拟合 Model2a，或者是否可以对 Model2b 进行任何更改以使其更适合后验？
非常感谢您的帮助。