短信数量异常检测

2019-09-20 本文已影响0人 Vector_Wan

最近开始学习贝叶斯统计这门课，感觉进入到了一个新世界。今天我们来研究一个贝叶斯模型玩玩。

多长时间以来，经典的频率学派和贝叶斯学派一直争论不休，拿一个经典的抛硬币的例子来说明他们之间的不同：

我扔硬币，扔了 100 次全是正面，问下一次扔硬币还是反面的概率：

频率学派的人会说：不管前面怎样，下一次还是正面反面各百分之五十！！
贝叶斯的人会说，傻吗！！前面100次都是正面，下一次大概率正面呀！！

所以他们之间最大的差别就是贝叶斯使用了之前的经验来给出后面的分布，但是除了一些共轭分布外，计算后验分布很困难，但是没关系，我们有计算机，在这个计算机性能这么强悍的今天，计算机不用白不用嘛。我们就来结合一个实例看看贝叶斯模型是如何建立，又是如何使用计算机进行求解的。

在阅读下面的步骤之前你需要了解：概率论与数理统计，贝叶斯统计，最好学过时间序列分析，随机过程。

我们先来看看我们要解决的问题，我们获取到了一组用户一段时间内每天收到的微信信息数量的数据，我们希望根据这些数据判断用户的行为有没有发生变化，比如订阅了天气推送，开始了一段新的感情...

首先我们先导入数据并且直观的查看一下

# 加载相关包，进行相关设置
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import pymc as pm
from IPython.core.pylabtools import figsize
figsize(12.5, 4)

# 中文支持
plt.rcParams['font.sans-serif'] = ['SimHei']  
plt.rcParams['axes.unicode_minus'] = False  

# 加载数据
count_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('时间（天）')
plt.ylabel('数量（条）')
plt.title('时间序列图')
plt.show()

在建模之前我们不太好说信息的数量有没有发生变化，因为随机扰动还是很强的，但是一个直观的感受就是后期 30 条以上的天数开始变多了，所以有可能后期短信数量比前面的多。

在建模之前我们再次明确我们的问题：

从收集到的数据来看短信数量有没有明显增加（客户信息有没有异常）
如果真的有异常发生那么发生在收么时候
如果有异常情况出现出现了几次

为了方便我们建模，我们暂时先对问题进行简化，假设异常情况在我们观测的这段时间中只出现了一次，稍后我们再谈复杂的情况。

首先我们需要一个分布来描述每天的短信数量，很容易想到泊松分布（Poisson）。泊松分布的概率密度函数是：

$Poisson(N = n) = \frac{\lambda^n e^{- \lambda }}{n!}, n>0$

泊松分布很适合描述一段时间内发生某件事情的次数。参数 $\lambda$ 的含义是单位时间，或单位面积内发生某件事的次数，给定次数 k 可以给出发生 k 次的概率。设每一天的短信数量为 $C_i$ 那么 $C_i ～Poisson(\lambda_i), i=0,1,...,73$ 。

根据我们的假设异常行为只出现了一次，也即是说如果异常情况发生 $\lambda_i$ 变化了一次， $\lambda_i$ 只能取两个值 $\lambda_1$ 和 $\lambda_2$ ，设 $\lambda_i$ 在 $\tau$ 处发生了变化，那么数学语言来描述就是：

$\lambda = \left\{ \begin{array}{lr} \lambda_1, & t \leq \tau\\ \lambda_2 ,& t > \tau\\ \end{array} \right.$

假如实际上不存在异常的情况，那么我们最后就会得到 $\lambda_1 = \lambda_2$ ，并且 $0 \leq \tau \leq 73$ 。

说到这里我们先说两句题外话， $\lambda$ 到底是什么？

这个问题可以理解为统计的动机是什么。在现实生活中，我们只能感受到样本 $C_i$ 的存在，但是不知道它的参数到底是多少，比方说，我们扔一枚不均匀的硬币，这件事情服从参数为 p 的二项分布，我们只知道每次实验的结果却不知道硬币的 p 到底是多少。

确定 p 其实很困难，在少数情况下我们深入到问题中去，然后做很多的假设才可以决定，比如扔一枚均匀（假设）硬币的 p 是 0.5。但是大多数情况下我们是无法确定它到底是多少的，所以我们只能根据实验观测估计 p。对于参数 p 的估计有很多设计好的方法，比如矩估计、极大似然估计，但是因为 p 并不是一个可以真实观测到的数据，谁也不能说那种方法好。

相比于频率学派的参数估计，我觉得贝叶斯的方法更加自然。贝叶斯推断围绕着对 $p$ 的取值估计，与其不断猜测 $p$ 的精确取值，不如用一个概率分布来描述 p 的所有可能取值。

这样上面的 $\lambda_1 = \lambda_2$ 的说法就是不准确的，我们最后得到的结果应该是 $\lambda_1$ 与 $\lambda_2$ 的分布相同。

好的，继续来说我们的模型，我们首先需要确定一个 $\lambda$ 的先验分布，在这个问题中 $\lambda$ 的含义是每天收到短信的平均数量，它的分布我们应该怎么确定呢？

首先我们知道短信的数量是一个大于 0 的连续型的变量，所以我们肯定要找一个取大于零的连续型的分布，因为我们没有什么经验，所以最容易想到的就是均匀分布，但是我试了一下发现效果很不好，可以执行下面的代码测试一下。

所以我们只好考虑考虑其他的分布，后来我想到了指数分布，指数分布描述的是等待一件事情发生所需时间的分布。

概率密度函数为：

$f(x; \lambda) = 1 - e^{- \lambda x}$

期望为 $E(x) = \frac{1}{\lambda}$

所以指数分布的参数我们设置为整个观测数据的均值。

# 设置先验分布
alpha = 1.0 / count_data.mean()
# alpha = pm.Uniform('alpha', lower=0, upper=1)

lambda_1 = pm.Exponential("lambda_1", alpha)
lambda_2 = pm.Exponential("lambda_2", alpha)
# lambda_1 = pm.Uniform('lambda_1', lower=0, upper=count_data.mean())
# lambda_2 = pm.Uniform('lambda_2', lower=0, upper=count_data.mean())

tau = pm.DiscreteUniform("tau", lower=0, upper=n_count_data)

@pm.deterministic
def lambda_(tau=tau, lambda_1=lambda_1, lambda_2=lambda_2):
    out = np.zeros(n_count_data)
    out[:tau] = lambda_1
    out[tau:] = lambda_2
    return out

建立模型，使用马尔可夫蒙特卡洛计算后验。

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

D:\anaconda\lib\site-packages\pymc\MCMC.py:81: UserWarning: Instantiating a Model object directly is deprecated. We recommend passing variables directly to the Model subclass.
  warnings.warn(message)


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

lambda_1_samples = mcmc.trace('lambda_1')[:]
lambda_2_samples = mcmc.trace('lambda_2')[:]
tau_samples = mcmc.trace('tau')[:]

figsize(12.5, 10)
# histogram of the samples:

ax = plt.subplot(311)
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)
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");

从上面后验的概率密度分布图中我们可以看见两个参数的分布有明显的不同， $\lambda_1$ 很有可能在 18 附近。 $\lambda_2$ 很有可能在 23 附近。这意味着均值确实发生了变化。

均值发生变化的时间很有可能是 44 或 45 天的时候。我们直观的看一下。

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()) / N


plt.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");

短信数量的均值发生变化说明了很多事情，比如在第 45 天附近的时候加入了什么群，不停的有消息进来，或者是一段新的恋情...

在之前我们对模型进行了简化，我们现在尝试将限制放开一些，我们试一下两次均值变化会发生什么。

# -*- coding: utf-8 -*-
import numpy as np
from matplotlib import pyplot as plt
import pymc as pm

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

alpha = 1.0 / count_data.mean()  # Recall count_data is the
                               # variable that holds our txt counts
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=n_count_data-1)
tau_2 = pm.DiscreteUniform("tau_2", lower=tau_1, upper=n_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(n_count_data)
    out[:tau_1] = lambda_1  # lambda before tau is lambda1
    out[tau_1:tau_2] = lambda_2  # lambda after (and including) tau is lambda2
    out[tau_2:] = lambda_3  # lambda after (and including) tau is lambda2
    return out


observation = pm.Poisson("obs", lambda_, value=count_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)

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')[:]

# histogram of the samples:

# 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_3$",color='#6A63A6',normed=True)
plt.legend(loc="upper left")
plt.grid(True)
plt.xlim([15, 30])
plt.xlabel("$\lambda_3$ value")
 

plt.show()

D:\anaconda\lib\site-packages\pymc\MCMC.py:81: UserWarning: Instantiating a Model object directly is deprecated. We recommend passing variables directly to the Model subclass.
  warnings.warn(message)


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

D:\anaconda\lib\site-packages\matplotlib\axes\_axes.py:6521: MatplotlibDeprecationWarning: 
The 'normed' kwarg was deprecated in Matplotlib 2.1 and will be removed in 3.1. Use 'density' instead.
  alternative="'density'", removal="3.1")
D:\anaconda\lib\site-packages\matplotlib\axes\_axes.py:6521: MatplotlibDeprecationWarning: 
The 'normed' kwarg was deprecated in Matplotlib 2.1 and will be removed in 3.1. Use 'density' instead.
  alternative="'density'", removal="3.1")
D:\anaconda\lib\site-packages\matplotlib\axes\_axes.py:6521: MatplotlibDeprecationWarning: 
The 'normed' kwarg was deprecated in Matplotlib 2.1 and will be removed in 3.1. Use 'density' instead.
  alternative="'density'", removal="3.1")

这两个分布还是有一定的差异性，说明在条件放开一些的情况下，短信数量还是发生了变化。

最后给出 pym3 的实现，现在 pym2 正在逐渐被弃用，目前被维护的是 pym3.

import pymc3 as pm
import theano.tensor as tt

with pm.Model() as model:
    alpha = 1.0/count_data.mean()  # Recall count_data is the
                                   # variable that holds our txt counts
    lambda_1 = pm.Exponential("lambda_1", alpha)
    lambda_2 = pm.Exponential("lambda_2", alpha)
    
    tau = pm.DiscreteUniform("tau", lower=0, upper=n_count_data - 1)
    
print(tau.random(), tau.random(), tau.random())

2 14 71

with model:
    idx = np.arange(n_count_data) # Index
    lambda_ = pm.math.switch(tau > idx, lambda_1, lambda_2)
    
with model:
    observation = pm.Poisson("obs", lambda_, observed=count_data)

with model:
    step = pm.Metropolis()
    trace = pm.sample(10000, tune=5000,step=step, cores=4)

Sampling 4 chains: 100%|██████████████████████████████████████████████████████| 60000/60000 [14:27<00:00, 69.14draws/s]
The number of effective samples is smaller than 25% for some parameters.

lambda_1_samples = trace['lambda_1']
lambda_2_samples = trace['lambda_2']
tau_samples = trace['tau']

figsize(12.5, 10)
#histogram of the samples:

ax = plt.subplot(311)
ax.set_autoscaley_on(False)

plt.hist(lambda_1_samples, histtype='stepfilled', bins=30, alpha=0.85,
         label="posterior of $\lambda_1$", color="#A60628", normed=True)
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", normed=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");

短信数量异常检测

猜你喜欢

热点阅读