PyMC 的拓展

PyMC 变量

模型的上下文
在PyMC3中，我们通常会在模型的上下文内处理模型中所需的所有变量。如下所示：

import pymc3 as pm

model = pm.Model()
# 模型上下文
with model:
    parameter = pm.Exponential("poisson_param", 1.0)
    data_generator = pm.Poisson("data_generator", parameter)
model

与 PyMC 相比，PyMC3 中的变量必须定义在模型内。在给定的模型上下文中创建的任何变量都将自动分配给该模型。当然，你也可以尝试在模型上下文之外定义变量，观察程序是否报错。
通过 with 模型对象名称，我们可以在同一个模型的上下文中继续定义新的变量。
但是请注意，如果一个代码段运行了多次，程序就会报错。因为，多次运行，相当于你在一个模型中定义了多个名字相同的变量。解决方法是，从模型初始化的那段代码块开始，从头运行。

with model:
    data_plus_one = data_generator + 1
data_plus_one

一旦我们在上下文中定义了变量，我们就可以在上下文外使用这个变量了。如下，我们展示了 model 中的 parameter 对象中的某个属性。

# 具体含义后面会介绍，这里只需要知道它是 model 中的变量即可
parameter.tag.test_value

所有 PyMC3 变量都有一个初始值（即 test_value 属性），让我们输出它们看看：

print("parameter.tag.test_value =", parameter.tag.test_value)
print("data_generator.tag.test_value =", data_generator.tag.test_value)
print("data_plus_one.tag.test_value =", data_plus_one.tag.test_value)

test_value 仅仅用于模型之中，表示的是该变量所属分布函数的采样起点。我们可以通过修改参数 testval 值，来改变变量的初始状态。

with pm.Model() as model:
    parameter = pm.Exponential("poisson_param", 1.0, testval=0.5)

print("\nparameter.tag.test_value =", parameter.tag.test_value)

如果你使用的先验条件极不稳定，那么你就可以使用上面参数，来为这个分布函数的采样指定一个比较好的起点。

PyMC 中存在着两种类型的编程变量：随机型和确定型
随机型变量是一种不确定的变量。即使你知道这个变量的参数和组件的所有值，它仍然是随机的。此类别中包括类：Poisson、DiscreteUniform 和 Exponential。
确定型变量是已知变量参数和组成部分的非随机变量。 一种快速检测它为确定型变量的方法是：如果有一个变量 foo，我们知道 foo 的所有组件的值，我们就可以确定 foo 的值。那么 foo 就是一个确定型变量。

接下来，我们将对这两种变量进行具体的分析。
初始化随机型变量
下面，我们使用 pm.DiscreteUniform 函数来初始一个随机型变量，格式如下：
somevariable=pm.DiscreteUniform("discreteunivar",0,4)
在该函数中，传入的第一个参数是一个字符串，它表示的是变量的名字。 在接下来的分析中会讲到，变量名字将会被用来检索后验分布，所以最好使用描述性命名。
其他的参数随着变量的类别不同而不同。如，在上面的函数中，0,4 是离散均匀分布中随机变量的上界和下界。在 PyMC 文档 中，包含了其他随机变量的具体参数。
对于多变量的问题，与其构建一个随机变量的 Python 数组，不如在调用随机变量构建函数时，通过指定大小关键字，直接构建多个变量。让我们举个例子来对这句话进行阐述。
假设现在我们需要定义多个变量βi(i=1,...,N) 来进行建模，首先想到的方法就是逐个定义，代码如下：

beta_1 = pm.Uniform("beta_1", 0, 1)
beta_2 = pm.Uniform("beta_2", 0, 1)

为了避免这样的重复定义，我们可以通过一个参数，将这些变量一次性全部定义出来，并打包成一个变量，代码如下：

betas = pm.Uniform("betas", 0, 1, shape=N)

其中 shape 表示的就是变量的个数。
初始化确定性变量
由于我们模拟的大多数变量都是随机型变量。因此为了区分，我们一般使用 pm.Deterministic 函数对其进行封装。

deterministic_variable = pm.Deterministic("deterministic variable", some_function_of_variables)

我们只需要在 PyMC3 中调用 Deterministic 类，并传入我们想要传入的变量即可。在 PyMC 中变量即函数。因为每个变量，都会对应一个分布函数。因此上面的 some_function_of_variables 即为我们想要处理的变量。
如下所示，我们定义一个函数，然后将参数传入该函数，再把该函数传入 Deterministic 函数中，就能够生成最后的确定型变量。

def subtract(x, y):
    return x - y


model = pm.Model()
with model:
    stochastic_1 = pm.Uniform("U_1", 0, 1)
    stochastic_2 = pm.Uniform("U_2", 0, 1)

    det_1 = pm.Deterministic("Delta", subtract(stochastic_1, stochastic_2))
det_1

上述所说的封装前置是一种简单的构建确定型变量的方法，但它并非唯一的方法。元素操作、加法、指数和类似的方式都可以产生确定型变量。例如，下面的代码也产生了一个确定型的变量。

# 该代码为上一个实验中的收发短信例子中的定义变量 λ1、λ2 和τ的代码片段
with pm.Model() as model:
    lambda_1 = pm.Exponential("lambda_1", 1.0)
    lambda_2 = pm.Exponential("lambda_2", 1.0)
    tau = pm.DiscreteUniform("tau", lower=0, upper=10)

new_deterministic_variable = lambda_1 + lambda_2
new_deterministic_variable

如上所示，如果变量 lambda_1 和 lambda_2 已知，则这两个参数相加的和 new_deterministic_variable 则已知。因此可以说 new_deterministic_variable 为确定型变量。

image.png

并且 PyMC3 的代码如下：

import numpy as np
# 在上个实验中，我们有74个数据点
n_data_points = 5
idx = np.arange(n_data_points)
with model:
    lambda_ = pm.math.switch(tau >= idx, lambda_1, lambda_2)

其实我们可以很清楚的明白，如果τ,λ1和λ2是已知的，那么λ 则也是完全已知的。因此λ 是一个确定性变量。
这里我们使用了 switch 函数使λ 在合适的时间从λ1变为λ2。该函数主要来自于 theano 的包。该工具包我们会在后面实验进行讨论。

实验二中遗留的问题

让我们利用我们今天学的 PyMC3 ，更加细致的讲解一下上一个试验中的一些代码。
既然说到了上一个收发短信的实验，那么接下来，我们就来看看λ1的先验分布函数到底长什么样？

from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt
import scipy.stats as stats
%matplotlib inline
figsize(12.5, 4)
# lambda_1 是在上上个代码块定义的
# 对 lambda_1 的分布函数进行采样，共采样20000次
samples = lambda_1.random(size=20000)

# 统计每次出现的频率
plt.hist(samples, bins=70, normed=True, histtype="stepfilled")
plt.title("Prior distribution for $\lambda_1$")
plt.xlim(0, 8)

上图中，我们可以清楚的看到λ1的先验概率。还记得在上一个实验中，我们定义了模型的变量后，接下来做了什么吗？接下来，我们需要将定义的模型与实际的数据结合起来，然后在对模型进行训练学习，最后得到每个参数的具体分布函数。
为了避免又去加载数据，这里我们直接用模拟一组真实数据为 data:

data = np.array([10, 25, 15, 20, 35])
data

接下来，我们使用的是定义随机变量的函数 Poisson 。并将真实数据加入模型中。而这个函数有个关键字参数：observed 。其实这个关键字的作用非常简单：将变量的当前值固定为真实数据 data，即修改随机型变量中的 test_value 的值。代码如下所示：

with model:
        # 第一个参数：定义的随机型变量的名字
        obs = pm.Poisson("obs", lambda_, observed=data)
print("value: ", obs.tag.test_value)

如上所示，其实我们将数据包含到模型中的方式就是：将随机型变量初始化为一个固定的值 data。
我们可以输出一下 model ，观察它到底存了哪些变量：

model

接下来，我们只需要进行训练即可。当然，由于下面的代码和实验 2 相同，因此，这里就不做赘述了。且关于模型学习的代码讲解，我们会在后面的实验中进行。

建模方法的总结

在贝叶斯建模的时候，思考数据是如何产生的能够加深我们对这种方法的理解。其实建模的本质，就是找到一个模型来描述这些数据的产生规律。现在让我们以短信收发为例，总结一下模型构建的步骤与思路。
在拿到数据后，我们首先开始思考，“什么随机变量最能描述这些统计数据呢？”。我们发现 Poisson 随机变量是一个很好的选择。因为，它能够很好地代表统计数据。于是，我们决定就用 Poisson 分布来模拟真实收数据的分布。
接下来，我们想，“如果短信接收数据服从 Poisson 分布的话，那么 Poisson 分布需要什么东西呢？”。OK，Poisson 分布需要一个参数λ。
那么我们知道λ 的具体值吗？不，实际上，不知道。我们猜测λ 有两个取值，一个对应于早期行为，另一个对应于后期行为。我们并不知道这个参数是在什么时候发生的变化。因此，我们还需要假设一个行为变化的时间点τ。
那么，对于这两个λ，我们应该选取一个怎样的分布呢？然后，我们想起了可以给正实数赋概率的指数分布函数。但是，指数分布函数也存在它自己的参数：α
那么，参数α 是什么样子呢？我们也不知道。此时我们可以继续为参数α 选择一个分布。但是，我们已经积累了很多未知变量了，现在还是停下来为妙。因此，我们觉得给 α 赋予一个固定的值，那么α 选择多少为宜呢？
因为我们知道λ 的范围在 10 到 30 之间，所以如果α 的取值太小是不适合的。类似的，取值太大，也不好。α 的最佳取值应该可以反映我们对λ 的判断，即此时，真实数据的均值应该和λ 的均值相等。这些，我们已经在实验2中，解释过了。
对于τ 什么时候发生，我们并没有什么概念。因此，我们假设了τ 是来自整个时段的一个离散平均分布。

在下图中，我们将上述过程进行了可视化，详细的表示了模型建立的过程。

image.png

不管是 PyMC 还是其他的概率编程语言，其设计理念都是为了讲述这一数据生成过程。

利用模型产生模拟数据

现在，我们将前面讨论的步骤反过来，就可以用来模拟真实数据集的产生了。
首先，让我们先确定行为发生变化的日期，这里我们随机假设一个日期：

# 设总共80天，利用 Numpy 随机确定行为变化的日期
tau = np.random.randint(0, 80)
print(tau)

从Exp(α) 分布中抽取λ1和λ2的值，随机指定两个参数的值：

# 这里的 alpha 为想模拟产生的数据的平均值的逆
# 可以根据实际情况，自己定义
alpha = 1./20.
lambda_1, lambda_2 = np.random.exponential(scale=1/alpha, size=2)
print(lambda_1, lambda_2)

接下来，我们只需要将λ 的值传入 Poisson 函数中即可。如下代码所示，前τ条数据，由参数λ=λ1的 Poisson 分布产生。第 τ 条到 第 80 条数据由λ=λ2的 Poisson 分布产生.

data = np.r_[stats.poisson.rvs(mu=lambda_1, size=tau),
             stats.poisson.rvs(mu=lambda_2, size=80 - tau)]
data

让我们对我们产生的数据继续可视化。

plt.bar(np.arange(80), data, color="#348ABD")
plt.bar(tau-1, data[tau - 1], color="r", label="user behaviour changed")
plt.xlabel("Time (days)")
plt.ylabel("count of text-msgs received")
plt.title("Artificial dataset")
plt.xlim(0, 80)
plt.legend()

从上图可以看出，我们模拟出来的数据并不像真实的观察数据，其实这是合理的。因为，模拟的数据和真实数据相同的可能性本来就很低。而我们使用 PyMC 来建模的目的，就是为了找到最优的参数，来最大化这种可能性。
虽然产生模拟数据的能力只是我们模型的一个小作用。但是，后面我们会发现，这个作用在贝叶斯推断方法中是很重要的。例如，在下面的代码中，我们可以产生更多的人工模拟数据。

def plot_artificial_sms_dataset():
    tau = stats.randint.rvs(0, 80)
    alpha = 1./20.
    lambda_1, lambda_2 = stats.expon.rvs(scale=1/alpha, size=2)
    data = np.r_[stats.poisson.rvs(mu=lambda_1, size=tau), stats.poisson.rvs(
        mu=lambda_2, size=80 - tau)]
    plt.bar(np.arange(80), data, color="#348ABD")
    plt.bar(tau - 1, data[tau-1], color="r", label="user behaviour changed")
    plt.xlim(0, 80)


figsize(12.5, 5)
plt.title("More example of artificial datasets")
for i in range(4):
    plt.subplot(4, 1, i+1)
    plot_artificial_sms_dataset()