背景知识

概率质量函数
概率质量函数（probability mass function，简写为pmf）是离散随机变量在各特定取值上的概率。
概率质量函数和概率密度函数不同之处在于：概率质量函数是对离散随机变量定义的，本身代表该值的概率；概率密度函数是对连续随机变量定义的，本身不是概率，只有对连续随机变量的概率密度函数在某区间内进行积分后才是概率。
离散

from IPython.core.pylabtools import figsize
import numpy as np
from matplotlib import pyplot as plt
import scipy.stats as statsfigsize(12.5, 4)
a = np.arange(16)
poi = stats.poisson
lambda_ = [1.5, 4.25]
colours = ["#348ABD", "#A60628"]plt.bar(a, poi.pmf(a, lambda_[0]), color=colours[0],label="$\lambda = %.1f$" % lambda_[0], alpha=0.60,edgecolor=colours[0], lw="3")plt.bar(a, poi.pmf(a, lambda_[1]), color=colours[1],label="$\lambda = %.1f$" % lambda_[1], alpha=0.60,edgecolor=colours[1], lw="3")plt.xticks(a + 0.4, a)
plt.legend()
plt.ylabel("probability of $k$")
plt.xlabel("$k$")
plt.title("Probability mass function of a Poisson random variable; differing \
$\lambda$ values");
plt.show()

连续

from IPython.core.pylabtools import figsize
import numpy as np
from matplotlib import pyplot as plt
import scipy.stats as statsfigsize(12.5, 4)poi = stats.poisson
lambda_ = [1.5, 4.25]
colours = ["#348ABD", "#A60628"]a = np.linspace(0, 4, 100)
expo = stats.expon
lambda_ = [0.5, 1]for l, c in zip(lambda_, colours):plt.plot(a, expo.pdf(a, scale=1./l), lw=3,color=c, label="$\lambda = %.1f$" % l)plt.fill_between(a, expo.pdf(a, scale=1./l), color=c, alpha=.33)plt.legend()
plt.ylabel("PDF at $z$")
plt.xlabel("$z$")
plt.ylim(0,1.2)
plt.title("Probability density function of an Exponential random variable;\differing $\lambda$")
plt.show()

实列

实例：从短信数据推断行为

from IPython.core.pylabtools import figsize
# IPython.core.pylabtools.figsize(sizex, sizey)import numpy as np
from matplotlib import pyplot as plt
import scipy.stats as stats  # 统计推断包
import pymc3 as pm
import theano.tensor as ttfigsize(12.5, 3.5)  # set figure sizecount_data = np.loadtxt("./data/txtdata.csv")
n_count_data = len(count_data)
plt.bar(np.arange(n_count_data), count_data, color="#348ABD")  # 画柱状图
plt.xlabel("Time (days)")
plt.ylabel("count of text-msgs received")
plt.title("Did the user's texting habits change over time?")
plt.xlim(0, n_count_data)  # 设置参数范围,图像所显示x轴的长度
plt.show()

设Ci表示第i天收到短信的条数，Ci ~ Poi( λ ）
确定了Ci的分布类型，但是参数λ是不能确定的，因为我们并不能很直观的观察到它的取值．但是poisson分布的特点是：当λ增大的时候，取大值的概率会增加，也就是说一天内收到信息较多的概率会增加．
由图1.4.1看出，后期收到短信几率增大，可能是λ在某时刻增大了，不妨假设在τ时刻，λ突然变大。（若算出λ1=λ2，则不存在假设情况）

在贝叶斯推断下，我们需要对不同可能的λ值分配相应的先验。
λ可取任意正数—>指数分布
λ1~Exp(α）
λ2~Exp(α）
α为超参数。我们在这里不想对它赋予太多主观色彩，建议将其设定为样本中计算平均值的逆。这样比较客观，可减少超参数对模拟造成的影响。

对于参数τ，现在很难选择合适的先验，那么就假设每一天的先验估计都是一样的。

with pm.Model() as model:alpha = 1.0 / count_data.mean()  # Recall count_data is the variable that holds our txt countslambda_1 = pm.Exponential("lambda_1", alpha)lambda_2 = pm.Exponential("lambda_2", alpha)tau = pm.DiscreteUniform("tau", lower=0, upper=n_count_data)# 产生参数λ1、λ2、τ、α。其中λ1、λ2由机器随机产生，在训练过程中找到更优的τ值。with model:idx = np.arange(n_count_data)  # Indexlambda_ = pm.math.switch(tau >= idx, lambda_1, lambda_2)# The switch() function assigns lambda_1 or lambda_2 as the value of lambda_, depending on what side of tau we are on.# The values of lambda_ up until tau are lambda_1 and the values afterwards are lambda_2.# tau, lambda_1, lambda_2 are random, lambda_ will be random. we are not fixing any variables yet!# Poisson(lambda_)with model:observation = pm.Poisson("obs", lambda_, observed=count_data)#创建模型实例

MCMC
从λ1、λ2、τ的后验分布返回数千个随机变量，我们将样本（trace）收集，绘制成直方图。

    with model:step = pm.Metropolis()trace = pm.sample(10000, tune=5000, step=step)#对后验分布进行10000次采样，指定采样器pm.Metropolis()，turn要调整的迭代次数（如果适用）（默认为无）lambda_1_samples = trace['lambda_1']lambda_2_samples = trace['lambda_2']tau_samples = trace['tau']#采样值存储在trace对象中，trace对象是一个字典

这个过程会有一些慢。

未知变量的后验分布，参数可视化

figsize(12.5, 10)
# histogram of the samples:ax = plt.subplot(311)
# subplot(numRows, numCols, plotNum)
# 图表的整个绘图区域被分成 numRows 行和 numCols 列plotNum 参数指定创建的 Axes 对象所在的区域numRows,numCols和plotNum这三个数都小于10的话,可以把它们缩写为一个整数,
# 例如 subplot(323) 和 subplot(3,2,3) 是相同的.ax.set_autoscaley_on(False)
plt.hist(lambda_1_samples, histtype='stepfilled', bins=30, alpha=0.85,label="posterior of $\lambda_1$", color="#A60628", density=True)
#density=True归一化；bins=30 30个箱子，即柱子；histtype : {‘bar’, ‘barstacked’, ‘step’, ‘stepfilled’}, optional(选择展示的类型,默认为bar)plt.legend(loc="upper left")
plt.title(r"""Posterior distributions of the variables $\lambda_1,\;\lambda_2,\;\tau$""")
plt.xlim([15, 30])
plt.xlabel("$\lambda_1$ value")ax = plt.subplot(312)
ax.set_autoscaley_on(False)
plt.hist(lambda_2_samples, histtype='stepfilled', bins=30, alpha=0.85,label="posterior of $\lambda_2$", color="#7A68A6", density=True)
plt.legend(loc="upper left")
plt.xlim([15, 30])
plt.xlabel("$\lambda_2$ value")plt.subplot(313)
w = 1.0 / tau_samples.shape[0] * np.ones_like(tau_samples)  #归一化
plt.hist(tau_samples, bins=n_count_data, alpha=1, label=r"posterior of $\tau$", color="#467821", weights=w, rwidth=2.)
plt.xticks(np.arange(n_count_data))plt.legend(loc="upper left")
plt.ylim([0, .75])
plt.xlim([35, len(count_data) - 20])
plt.xlabel(r"$\tau$ (in days)")
plt.ylabel("probability")
plt.show()

λ1大概为18，λ2大概为23，这两个λ的后验分布明显是不同的，表明用户接收短信的行为确实发生了变化。
λ的后验分布并不像指数分布，事实上后验分布并不是我们从原始模型中可以辨别的任何分布，用计算机可以很好展示出来。
τ是离散变量，后验分布有所不同。在44天时，有50%把握说明用户行为改变。

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

figsize(12.5, 5)
# tau_samples, lambda_1_samples, lambda_2_samples contain
# N samples from the corresponding posterior distribution
N = tau_samples.shape[0]
expected_texts_per_day = np.zeros(n_count_data)
for day in range(0, n_count_data):# ix is a bool index of all tau samples corresponding to# the switchpoint occurring prior to value of 'day'ix = day < tau_samples# Each posterior sample corresponds to a value for tau.# for each day, that value of tau indicates whether we're "before"# (in the lambda1 "regime") or#  "after" (in the lambda2 "regime") the switchpoint.# by taking the posterior sample of lambda1/2 accordingly, we can average# over all samples to get an expected value for lambda on that day.# As explained, the "message count" random variable is Poisson distributed,# and therefore lambda (the poisson parameter) is the expected value of# "message count".expected_texts_per_day[day] = (lambda_1_samples[ix].sum()+ lambda_2_samples[~ix].sum()) / Nplt.plot(range(n_count_data), expected_texts_per_day, lw=4, color="#E24A33",label="expected number of text-messages received")
plt.xlim(0, n_count_data)
plt.xlabel("Day")
plt.ylabel("Expected # text-messages")
plt.title("Expected number of text-messages received")
plt.ylim(0, 60)
plt.bar(np.arange(len(count_data)), count_data, color="#348ABD", alpha=0.65,label="observed texts per day")plt.legend(loc="upper left")
plt.show()

观察上图中的结果发现，分析的结果很符合之前的估计，用户的行为确实发生了改变，而且变化是很突然的，所以可以推测情况产生的原因可能是：短信资费降低，或者逢年过节期间，或者天气提醒短信订阅等等．
关于scipy.stats添加链接描述
pymc官方文档添加链接描述
参考书籍《贝叶斯方法概率编程与贝叶斯推断》

PyMc01短信推断相关推荐

在真正的短信网络钓鱼攻击内部
SMS based phishing attacks (Smishing) are a real threat that we see every day. To help you spot them ...
短信验证码被盗刷了怎么办？
短信接口被恶意攻击的最常见的一种情况就是短信验证码被盗刷. 一般分两种情况: 攻击一个手机号一般是针对个人的,当攻击者通过某种途径获取被攻击者的号码后就会利用某网站或者某一个软件的短信功能对其进行短 ...
为了更有效率地偷钱，Android root木马开始试水短信扣费诈骗
本文讲的是为了更有效率地偷钱,Android root木马开始试水短信扣费诈骗, 自2006年9月以来,我们就一直在监控Google Play商店有关Ztorg木马的各种新变异版本 ,到目前为止,我 ...
Ztorg木马分析：从Android root木马演变到短信吸血鬼
阅读原文请点击摘要: 本月第二次,Google 从官方应用商店 Google Play 中移除了伪装成合法程序的恶意应用.被移除的应用都属于名叫 Ztorg 的 Android 恶意程序家族.目前为 ...
SSH2框架实现注冊发短信验证码实例
这两天開始写程序了,让用SSH2框架,曾经没有接触过Java项目更没有接触过SSH2框架,所以用注冊開始了我Java之旅.后来发现,后台代码挺easy理解的,跟.net的差点儿相同.就是层与层之间的调 ...
个人信息泄露一条短信盗走积蓄
在回复了一条短信后,一夜之间,许先生的积蓄几乎全部被人转走.银行卡.支付宝和百度钱包的钱在几个小时内被骗子轻易转出去,自己眼睁睁看着,却无能为力. 网络安全专家分析,根据许先生提供的信息,骗子在实施诈 ...
Android短信收到，语音播报
发送短信功能界面 /*** 发送短信Demo* * @description:* @author ldm* @date 2016-4-22 上午9:07:53*/ public class SmsAc ...
移动互联网（一）短信和彩信界面开发包
参与这个项目,短信和彩信功能.它被认为是该项目的一个重要组成部分,如何开发这个功能以前认为.例如,有很多订单我们永和系统,怎样让用户及时知道自己卡里的消费情况?怎样让用户心中存在安全感,试想想在你的银 ...
程序主动进行电话短信报警，自定义电话、短信、钉钉报警通知
程序主动进行电话短信报警,自定义电话.短信.钉钉报警通知一. 规则说明这里我们要利用到阿里云的云监控的手段,有一个叫做事件监控的东西,可以通过自定义事件上传来进行监控报警. 流程: 程序发现错误 ...

PyMc01短信推断

U1、从短信推断行为

背景知识

实列

PyMc01短信推断相关推荐

最新文章

热门文章