PyMC3煤矿灾害示例问题：添加第二个矿井

Question

我正在使用PyMC3，试图在PyMC3文档中安装mining disaster switchpoint model的修改版本。假设您有两个煤矿（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)

收益率与文档示例非常相似。以下是跟踪图：

在拟合保持矿井分离的模型时，我尝试了两种方法，由于不同的原因，这两种方法都是次优的。第一种方法是分别为每个矿井拟合两个数据可能性

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

第二个想法是使用与Hierarchical Linear Regression example类似的方法，并结合使用串联、索引和先验的形状输出，以拟合每个参数的向量和单个数据似然

模型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。然而，从后验角度看，它不适合切换点

总之，有没有一种方法可以使模型2A以允许计算WAIC的方式进行拟合，或者可以对模型2B进行任何更改以使其更符合后验条件

非常感谢你的帮助

Answer 1

我没有一个明确的答案，但这里有一些建议，应该可以帮助你把事情做好

首先，将ArviZ更新为其最新版本，从错误消息中可以看出，您的版本比第一个版本旧，并且支持多种可能性。即使看起来您正在使用PyMC3函数，PyMC3也会将其绘图和统计数据委托给ArviZ

然后，我建议大家看看ArviZ的教育资源。目前有一个open PR来增加关于这类问题的指导。这是笔记本的link。我认为它处于一个足够先进的状态，是有用的。如果不是，还有其他关于here on SO或PyMC3语篇1，2的问题。这些应该包括一些额外的例子

最后，这里是这些详细答案的关键思想。第一个关键点是没有一个正确的答案，根据你想问的问题，waic/loo可以用不同的方法计算。第二个关键想法是ArviZ让我们选择如何计算waic/loo以适应所有可能的问题，因此在多个可能性情况下，需要对log_likelihoods组中的数据进行后处理