PyMC3 -煤矿灾难示例-问题re:添加第二个矿井

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);

​

​

fit看起来不错，而且似乎对切换点的差异很敏感。然而，我不能计算出世界人工智能大会或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);

这个模型很适合，并且允许计算一个世界人工智能大会。然而，从后部看，它不能适应开关。

​

​

所以总而言之，有没有一种方法可以用一种方式来拟合Model2a，允许计算出一个世界人工智能大会，或者是否可以对Model2b进行任何更改，使其能够拟合更好的后验结果？

非常感谢您的帮助。

Answer 1

我没有一个明确的答案，但这里有一些建议，应该会帮助你让事情变得正常。

首先将ArviZ更新到其最新版本，从错误消息看，您的版本似乎比支持多重可能性的第一个版本旧。尽管看起来您正在使用PyMC3函数，但PyMC3将其绘图和统计信息委托给ArviZ。

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

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