贝叶斯概率推断：短信数据推断行为

数据集是一个包含用户每天接收到短信条数的数据，利用用户短信数据来推断用户行为的变化．

%matplotlib inline
import numpy as np
from matplotlib import pyplot as plt
from IPython.core.pylabtools import figsize

直观感受下数据

# 加载数据
figsize(15,6)
sms_data = np.loadtxt('data/txtdata.csv')
# 样本数量
count_data = len(sms_data)
plt.bar(np.arange(count_data), sms_data)
plt.xlabel("Time (days)")
plt.ylabel("Text message received")
plt.title("Did the user's texting habits change over time?")
plt.grid(True)
plt.xlim(0, count_data)

(0, 74)

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-f0XRYezA-1601974878023)(output_2_1.png)]
　　仅从图中并不能很容易的观察到，这段时间内用户的行为发生了什么变化，对于这中离散的随机变量Poisson分布可以很好的模拟这种数据．假设iii天收到的短信条数CiC_iCi是服从参数λ\lambdaλ的泊松分布，记为：Ci～Poi(λ)C_i ～Poi(\lambda)Ci～Poi(λ)
　　确定了CiC_iCi的分布类型，但是参数λ\lambdaλ是不能确定的，因为我们并不能很直观的观察到它的取值．但是poisson分布的特点是：当λ\lambdaλ增大的时候，取大值的概率会增加，也就是说一天内收到信息较多的概率会增加．观察上图发现，40天之后，短信的数量是有明显的增加，这说明这期间λ\lambdaλ增加了．
　　为了模拟这个情况，我们要假设一个转折点τ\tauτ,参数λ\lambdaλ在τ\tauτ天之后取值开始变大，所以λ\lambdaλ的取值有两个，在τ\tauτ之前一个，之后一个．
　　　　　　　　　　　　　　λ1:t<τ\lambda_1: t < \tauλ1:t<τ
　　　　　　　　　　　　　　λ2:t>τ\lambda_2 : t > \tauλ2:t>τ
　　在贝叶斯推断下，需要对λ1\lambda_1λ1，λ2\lambda_2λ2分配一个相应的先验概率，但是如何分配一个好的先验概率呢？并不知道．怎么办？
　　在poisson分布中参数λ\lambdaλ是可以取任意正数的，刚好指数分布对与任意的正数都存在一个连续的概率密度函数，或许可以用指数分布来模拟参数λ\lambdaλ,但是，指数分布也需要一个参数α\alphaα,又多了一个未知变量．参数λ\lambdaλ服从参数为α\alphaα的指数分布．
　　　　　　　　　　　　　　λ1～Exp(α)\lambda_1～Exp(\alpha)λ1～Exp(α)
　　　　　　　　　　　　　　λ2～Exp(α)\lambda_2～Exp(\alpha)λ2～Exp(α)
　　α\alphaα是一个父变量，因为它会影响到参数λ\lambdaλ，那么问题来了，α\alphaα的值如何确定？在整个模型中，不希望这个参数被赋予太多的主观色彩，所以将这个参数设置为样本平均值的倒数，为什么呢？因为：模型中我们假设参数λ\lambdaλ是服从参数为α\alphaα的指数分布的，指数分布的期望就是参数的逆，即：1/α1/\alpha1/α.
　　
　　　　　　　　　　　　　　1N∑Ci≈E[λ∣α]=1α\dfrac{1}{N}\sum C_i \approx E[\lambda|\alpha] = \dfrac{1}{\alpha}N1∑Ci≈E[λ∣α]=α1

补充：对于参数α\alphaα如果你有更好的选择，可以使用两个不同的α\alphaα来模拟不同时期的参数λ\lambdaλ
　　对于参数τ\tauτ,很难选择合适的先验概率，我们假设每天的先验是相等的即：1/70：
　　　　　　　　　　　　　　τ～DiscreteUniform(1,70)\tau ～ DiscreteUniform(1,70)τ～DiscreteUniform(1,70)
　　　　　　　　　　　　　　$ P(\tau = k)=1/70, k = 1,2,…70$
　　假设了这么多，未知变量的整体先验分布是什么样的？下面使用pymc来模拟，未知变量的先验分布情况：pymc是一个贝叶斯分析库，现在已经有了新版本pymc3,新旧版本可同时安装．

import pymc as pm

# 超参数alpha,
# 设定参数alpha为样本平均值的逆
alpha = 1.0/sms_data.mean()
# 参数lambda_1,2服从参数为alpha的指数分布
lambda_1 = pm.Exponential('lambda_1', alpha)
lambda_2 = pm.Exponential('lambda_2', alpha)
# 参数tau的取值范围０～count_data
tau = pm.DiscreteUniform('tau', lower=0, upper=count_data)

tau.random()  # 0~74之间的随机正整数

array(13)

# lambda_函数返回的是一个跟sms_data 等长的lambda参数数组
# 告诉pyMc这是一个定性函数
@pm.deterministic
def lambda_(tau= tau, lambda_1 = lambda_1, lambda_2 = lambda_2):out = np.zeros(count_data)# 设每天收到的信息数量是服从泊松分布的# 在tau天之前泊松分布对应的参数lambda_1out[:tau] = lambda_1# 在tau天之后泊松分布对应的参数lambda_2out[tau:] = lambda_2return out

observation = pm.Poisson("obs", lambda_, value=sms_data, observed=True)
# 创建模型实例
model = pm.Model([observation, lambda_1, lambda_2, tau])

# 马尔科夫链蒙特卡洛
mcmc = pm.MCMC(model)
mcmc.sample(40000,10000)

 [-----------------100%-----------------] 40000 of 40000 complete in 12.4 sec

未知变量的后验分布

# lambda_1的后验分布,一个长度为30000的数组
lambda_1_samples = mcmc.trace('lambda_1')[:]
# lambda_2后验分布，一个长度为30000的数组
lambda_2_samples = mcmc.trace('lambda_2')[:]
# tau后验分布,一个长度为30000的数组
tau_samples = mcmc.trace('tau')[:]

参数可视化

# lambda_1,lambda_2,tau的后验分布直方图
figsize(15, 10)# lambda_1
ax = plt.subplot(311)
# 取消自动缩放
ax.set_autoscaley_on(False)
plt.hist(lambda_1_samples, histtype='stepfilled', bins=30, label="$\lambda_1$", normed=True)
plt.legend(loc="upper left")
plt.grid(True)
plt.title(r"""Posterior distributions of the variables$\lambda_1,\;\lambda_2,\;\tau$""")
# x轴坐标范围
plt.xlim([15, 30])
plt.xlabel("$\lambda_1$ value")# lambda_2
ax = plt.subplot(312)
ax.set_autoscaley_on(False)
plt.hist(lambda_2_samples, histtype='stepfilled', bins=30,  label=" $\lambda_2$",color='#7A68A6',normed=True)
plt.legend(loc="upper left")
plt.grid(True)
plt.xlim([15, 30])
plt.xlabel("$\lambda_2$ value")# tau
plt.subplot(313)w = 1.0 / tau_samples.shape[0] * np.ones_like(tau_samples)
plt.hist(tau_samples, bins=count_data, alpha=1,label=r"$\tau$",color="#467821",weights=w ,rwidth=2.)
plt.xticks(np.arange(count_data))
plt.grid(True)
plt.legend(loc="upper left")
plt.ylim([0, .75])
plt.xlim([35, len(sms_data) - 20])
plt.xlabel(r"$\tau$ (in days)")
plt.ylabel("probability");

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-TYOA1G3l-1601974878028)(output_10_0.png)]
　　从上图中观察，参数的合理值：λ1\lambda_1λ1大概为18,λ2\lambda_2λ2大概为23,两个参数的差别很明显，这说明在不同时期参数λ\lambdaλ确实发生了变化，这也说明用户的接收短信的行为也发生的变化．变量τ\tauτ返回的是一个离散变量，从图中看到在45天有超过６成的把握可以确定用户行为发生了改变．43,44也是潜在的转折点．

后验样本

在０～７０天中，期望每天收到的信息数量等价于参数λ\lambdaλ的期望，为什么呢?因为，模型是假设Ci～Poi(λ)C_i ～Poi(\lambda)Ci～Poi(λ)，poisson分布的期望值等于它的参数λ\lambdaλ.下面计算每天短信条数的期望值．

# 这个期望值我们假设是服从泊松分布的，分布的期望值＝参数lambda
# 如果天数在转折点tau之前，取值lambda_1,否则取lambda_2,然后在取平均值
expected_texts_per_day = np.zeros(count_data)
for day in range(0, count_data):ix = day < tau_samplesexpected_texts_per_day[day] = (lambda_1_samples[ix].sum()+lambda_2_samples[~ix].sum())/len(tau_samples)figsize(15,6)
plt.plot(range(count_data), expected_texts_per_day, lw=4, color="#E24A33",label="expected number ")
plt.xlim(0, count_data)
plt.xlabel("Day")
plt.ylabel("Expected # text-messages")
plt.title("Expected number of text-messages received")
plt.bar(np.arange(len(sms_data)),sms_data, color="#348ABD", label="observed texts per day")
plt.grid(True)
plt.legend(loc="best");

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-prAaBiYH-1601974878033)(output_11_0.png)]
　　观察上图中的结果发现，分析的结果很符合之前的估计，用户的行为确实发生了改变，而且变化是很突然的，所以可以推测情况产生的原因可能是：短信资费降低，或者逢年过节期间，或者天气提醒短信订阅等等．

如何确定两个λ\lambdaλ不同的

首先通过观察数据的图像，直观的根据先验信息判定λ1,λ2\lambda_1,\lambda_2λ1,λ2是不同的,因为后期收到短信的数量是有明显增加的．但是这样的先验估计可能存在严重的偏差．如何证实呢？通过参数λ\lambdaλ的后验分布来验证．
　　方法是计算出P(λ1<λ2∣data)P(\lambda_1 < \lambda_2|data)P(λ1<λ2∣data),即在获得参数后验分布情况的条件下，计算出λ1<λ2\lambda_1 < \lambda_2λ1<λ2的概率．如果这个概率接近50%，这仍不能确定我们的先验估计是正确的．如果概率值接近100%,那么可以确定λ1≠λ2\lambda_1 \neq \lambda_2λ1=λ2.先验估计正确．

#　通过lambda_1和lambda_2的后验分布确定，它们的值不同的概率
print("the probability :%.3f"%(lambda_1_samples < lambda_2_samples).mean())

the probability :1.000

很明显百分之百的把握λ1≠λ2\lambda_1 \neq \lambda_2λ1=λ2．
比较λ1,λ2\lambda_1,\lambda_2λ1,λ2差值为1,2,5,10的概率

# 两个值之间相差１，２，５，１０的概率
for d in [1, 2, 5, 10]:v = (abs(lambda_1_samples - lambda_2_samples)>=d).mean()print("the probability the difference is larger than %d : %f"%(d, v))

the probability the difference is larger than 1 : 1.000000
the probability the difference is larger than 2 : 1.000000
the probability the difference is larger than 5 : 0.519733
the probability the difference is larger than 10 : 0.000000

扩充到两个转折点

假设现在我们对一个转折点表示很怀疑，我们现在认为用户的行为发生了两次改变．扩充之后，用户的行为分为三个阶段，三个泊松分布对应三个λ1,λ2,λ3\lambda_1, \lambda_2, \lambda_3λ1,λ2,λ3，两个转折点τ1,τ2\tau_1,\tau_2τ1,τ2
　　　　　　　　　　　　　　　　　
　　　　　　　　　　　　　　　　　λ1:t<τ1\lambda_1 : t < \tau_1λ1:t<τ1
　　　　　　　　　　　　　　　　　λ2:τ1≤t≤τ2\lambda_2 : \tau_1 \leq t \leq \tau_2λ2:τ1≤t≤τ2
　　　　　　　　　　　　　　　　　λ3:t>τ2\lambda_3 : t > \tau_2λ3:t>τ2
　　　　　　　　　　　　　　　　　λ1～Exp(α)\lambda_1 ～ Exp(\alpha)λ1～Exp(α)
　　　　　　　　　　　　　　　　　λ2～Exp(α)\lambda_2 ～ Exp(\alpha)λ2～Exp(α)
　　　　　　　　　　　　　　　　　λ3～Exp(α)\lambda_3 ～ Exp(\alpha)λ3～Exp(α)
　　　　　　　　　　　　　　　　　τ1～DiscreteUniform(1,69)\tau_1 ～ DiscreteUniform(1,69)τ1～DiscreteUniform(1,69)
　　　　　　　　　　　　　　　　　τ2～DiscreteUniform(τ1,70)\tau_2 ～ DiscreteUniform(\tau_1,70)τ2～DiscreteUniform(τ1,70)

代码实现

# 超参数alpha,
# 设定参数alpha为样本平均值的逆(其实这个alpha,也可以设置为三个，每个对应不同的泊松分布)
# 为了方便起见，有不想参杂较多的主观色彩，仍采用样本均值的倒数
alpha = 1.0/sms_data.mean()
lambda_1 = pm.Exponential("lambda_1",alpha)
lambda_2 = pm.Exponential("lambda_2",alpha)
lambda_3 = pm.Exponential("lambda_3",alpha)
tau_1 = pm.DiscreteUniform("tau_1", lower=0, upper=count_data-1)
tau_2 = pm.DiscreteUniform("tau_2", lower=tau_1, upper=count_data)

@pm.deterministic
def lambda_(tau_1=tau_1, tau_2=tau_2, lambda_1=lambda_1,lambda_2=lambda_2,lambda_3=lambda_3):out = np.zeros(count_data)out[:tau_1] = lambda_1out[tau_1:tau_2] = lambda_2out[tau_2:] = lambda_3return out

observation = pm.Poisson("obs",lambda_, value=sms_data, observed=True)
model = pm.Model([observation, lambda_1, lambda_2, lambda_3, tau_1, tau_2])
mcmc = pm.MCMC(model)
mcmc.sample(40000,10000)

 [-----------------100%-----------------] 40000 of 40000 complete in 17.6 sec

lambda_1_samples = mcmc.trace('lambda_1')[:]
lambda_2_samples = mcmc.trace('lambda_2')[:]
lambda_3_samples = mcmc.trace('lambda_3')[:]
tau_1_samples = mcmc.trace('tau_1')[:]
tau_2_samples = mcmc.trace('tau_2')[:]

figsize(12,10)# lambda_1
ax = plt.subplot(311)
ax.set_autoscaley_on(False)
plt.hist(lambda_1_samples, histtype='stepfilled', bins=30, label="$\lambda_1$", normed=True)
plt.legend(loc="upper left")
plt.grid(True)
plt.title(r"""Posterior distributions of the variables$\lambda_1,\;\lambda_2,\;\tau$""")
# x轴坐标范围
plt.xlim([15, 30])
plt.xlabel("$\lambda_1$ value")# lambda_2
ax = plt.subplot(312)
ax.set_autoscaley_on(False)
plt.hist(lambda_2_samples, histtype='stepfilled', bins=30,  label=" $\lambda_2$",color='#3009A6',normed=True)
plt.legend(loc="upper left")
plt.grid(True)
plt.xlim([30, 90])
plt.xlabel("$\lambda_2$ value")# lambda_3
ax = plt.subplot(313)
ax.set_autoscaley_on(False)
plt.hist(lambda_3_samples, histtype='stepfilled', bins=30,  label=" $\lambda_2$",color='#6A63A6',normed=True)
plt.legend(loc="upper left")
plt.grid(True)
plt.xlim([15, 30])
plt.xlabel("$\lambda_3$ value")

λ1,λ2,λ3\lambda_1, \lambda_2, \lambda_3λ1,λ2,λ3的后验分布

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-mCyhzodY-1601974878035)(output_20_1.png)]

τ1\tau_1τ1的后验分布

figsize(12,4)
# tau_1
w = 1.0 / tau_1_samples.shape[0] * np.ones_like(tau_1_samples)
plt.hist(tau_1_samples, bins=count_data, alpha=1,label=r"$\tau_1$",color="blue",weights=w )
plt.grid(True)
plt.legend(loc="upper left")
plt.xlabel(r"$\tau_1$ (in days)")
plt.ylabel("probability")
plt.show()

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-FMtrH7pc-1601974878036)(output_21_0.png)]

τ2\tau_2τ2的后验分布

figsize(12,4)
# tau_2
w = 1.0 / tau_2_samples.shape[0] * np.ones_like(tau_1_samples)
plt.hist(tau_2_samples, bins=count_data, alpha=1,label=r"$\tau_2$",weights=w,color="red",)
plt.xticks(np.arange(count_data))
plt.grid(True)
plt.legend(loc="upper left")
plt.ylim([0, 1.0])
plt.xlim([35, len(sms_data) - 20])
plt.xlabel(r"$\tau_2$ (in days)")
plt.ylabel("probability")

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-cherFEJu-1601974878038)(output_23_1.png)]

期望

expected_texts_per_day = np.zeros(count_data)
for day in range(0, count_data):ix_1 = day < tau_1_samplesix_2 = (day > tau_1_samples).all() and (day < tau_2_samples).all()ix_3 = day > tau_2_samplesexpected_texts_per_day[day] = (lambda_1_samples[ix_1].sum()+lambda_2_samples[ix_2].sum()+lambda_3_samples[ix_3].sum())/len(lambda_1_samples)

figsize(15,6)
plt.plot(range(count_data), expected_texts_per_day, lw=4, color="#E24A33",label="expected number ")
plt.xlim(0, count_data)
plt.xlabel("Day")
plt.ylabel("Expected # text-messages")
plt.title("Expected number of text-messages received")
plt.bar(np.arange(len(sms_data)),sms_data, color="#348ABD", label="observed texts per day")
plt.grid(True)
plt.legend(loc="best")

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-WlHgMSIQ-1601974878039)(output_25_1.png)]