短信數(shù)量異常檢測(cè)

最近開(kāi)始學(xué)習(xí)貝葉斯統(tǒng)計(jì)這門(mén)課，感覺(jué)進(jìn)入到了一個(gè)新世界。今天我們來(lái)研究一個(gè)貝葉斯模型玩玩。

多長(zhǎng)時(shí)間以來(lái)，經(jīng)典的頻率學(xué)派和貝葉斯學(xué)派一直爭(zhēng)論不休，拿一個(gè)經(jīng)典的拋硬幣的例子來(lái)說(shuō)明他們之間的不同：

我扔硬幣，扔了 100 次全是正面，問(wèn)下一次扔硬幣還是反面的概率：

  • 頻率學(xué)派的人會(huì)說(shuō)：不管前面怎樣，下一次還是正面反面各百分之五十！！
  • 貝葉斯的人會(huì)說(shuō)，傻嗎?。∏懊?00次都是正面，下一次大概率正面呀?。?/li>

所以他們之間最大的差別就是貝葉斯使用了之前的經(jīng)驗(yàn)來(lái)給出后面的分布，但是除了一些共軛分布外，計(jì)算后驗(yàn)分布很困難，但是沒(méi)關(guān)系，我們有計(jì)算機(jī)，在這個(gè)計(jì)算機(jī)性能這么強(qiáng)悍的今天，計(jì)算機(jī)不用白不用嘛。我們就來(lái)結(jié)合一個(gè)實(shí)例看看貝葉斯模型是如何建立，又是如何使用計(jì)算機(jī)進(jìn)行求解的。

在閱讀下面的步驟之前你需要了解：概率論與數(shù)理統(tǒng)計(jì)，貝葉斯統(tǒng)計(jì)，最好學(xué)過(guò)時(shí)間序列分析，隨機(jī)過(guò)程。

我們先來(lái)看看我們要解決的問(wèn)題，我們獲取到了一組用戶(hù)一段時(shí)間內(nèi)每天收到的微信信息數(shù)量的數(shù)據(jù)，我們希望根據(jù)這些數(shù)據(jù)判斷用戶(hù)的行為有沒(méi)有發(fā)生變化，比如訂閱了天氣推送，開(kāi)始了一段新的感情...

首先我們先導(dǎo)入數(shù)據(jù)并且直觀的查看一下

# 加載相關(guān)包，進(jìn)行相關(guān)設(shè)置
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  

# 加載數(shù)據(jù)
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('時(shí)間（天）')
plt.ylabel('數(shù)量（條）')
plt.title('時(shí)間序列圖')
plt.show()

在建模之前我們不太好說(shuō)信息的數(shù)量有沒(méi)有發(fā)生變化，因?yàn)殡S機(jī)擾動(dòng)還是很強(qiáng)的，但是一個(gè)直觀的感受就是后期 30 條以上的天數(shù)開(kāi)始變多了，所以有可能后期短信數(shù)量比前面的多。

在建模之前我們?cè)俅蚊鞔_我們的問(wèn)題：

  1. 從收集到的數(shù)據(jù)來(lái)看短信數(shù)量有沒(méi)有明顯增加（客戶(hù)信息有沒(méi)有異常）
  2. 如果真的有異常發(fā)生那么發(fā)生在收么時(shí)候
  3. 如果有異常情況出現(xiàn)出現(xiàn)了幾次

為了方便我們建模，我們暫時(shí)先對(duì)問(wèn)題進(jìn)行簡(jiǎn)化，假設(shè)異常情況在我們觀測(cè)的這段時(shí)間中只出現(xiàn)了一次，稍后我們?cè)僬剰?fù)雜的情況。

首先我們需要一個(gè)分布來(lái)描述每天的短信數(shù)量，很容易想到泊松分布（Poisson）。泊松分布的概率密度函數(shù)是：

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

泊松分布很適合描述一段時(shí)間內(nèi)發(fā)生某件事情的次數(shù)。參數(shù) \lambda 的含義是單位時(shí)間，或單位面積內(nèi)發(fā)生某件事的次數(shù)，給定次數(shù) k 可以給出發(fā)生 k 次的概率。設(shè)每一天的短信數(shù)量為 C_i 那么 C_i ～Poisson(\lambda_i), i=0,1,...,73。

根據(jù)我們的假設(shè)異常行為只出現(xiàn)了一次，也即是說(shuō)如果異常情況發(fā)生 \lambda_i 變化了一次， \lambda_i 只能取兩個(gè)值\lambda_1 和 \lambda_2 ，設(shè) \lambda_i 在 \tau 處發(fā)生了變化，那么數(shù)學(xué)語(yǔ)言來(lái)描述就是：

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

假如實(shí)際上不存在異常的情況，那么我們最后就會(huì)得到 \lambda_1 = \lambda_2，并且 0 \leq \tau \leq 73。

說(shuō)到這里我們先說(shuō)兩句題外話(huà)，\lambda 到底是什么？

這個(gè)問(wèn)題可以理解為統(tǒng)計(jì)的動(dòng)機(jī)是什么。在現(xiàn)實(shí)生活中，我們只能感受到樣本 C_i 的存在，但是不知道它的參數(shù)到底是多少，比方說(shuō)，我們?nèi)右幻恫痪鶆虻挠矌?，這件事情服從參數(shù)為 p 的二項(xiàng)分布，我們只知道每次實(shí)驗(yàn)的結(jié)果卻不知道硬幣的 p 到底是多少。

確定 p 其實(shí)很困難，在少數(shù)情況下我們深入到問(wèn)題中去，然后做很多的假設(shè)才可以決定，比如扔一枚均勻（假設(shè)）硬幣的 p 是 0.5。但是大多數(shù)情況下我們是無(wú)法確定它到底是多少的，所以我們只能根據(jù)實(shí)驗(yàn)觀測(cè)估計(jì) p。對(duì)于參數(shù) p 的估計(jì)有很多設(shè)計(jì)好的方法，比如矩估計(jì)、極大似然估計(jì)，但是因?yàn)?p 并不是一個(gè)可以真實(shí)觀測(cè)到的數(shù)據(jù)，誰(shuí)也不能說(shuō)那種方法好。

相比于頻率學(xué)派的參數(shù)估計(jì)，我覺(jué)得貝葉斯的方法更加自然。貝葉斯推斷圍繞著對(duì) p 的取值估計(jì)，與其不斷猜測(cè) p 的精確取值，不如用一個(gè)概率分布來(lái)描述 p 的所有可能取值。

這樣上面的 \lambda_1 = \lambda_2 的說(shuō)法就是不準(zhǔn)確的，我們最后得到的結(jié)果應(yīng)該是 \lambda_1 與 \lambda_2 的分布相同。

好的，繼續(xù)來(lái)說(shuō)我們的模型，我們首先需要確定一個(gè) \lambda 的先驗(yàn)分布，在這個(gè)問(wèn)題中 \lambda 的含義是每天收到短信的平均數(shù)量，它的分布我們應(yīng)該怎么確定呢？

首先我們知道短信的數(shù)量是一個(gè)大于 0 的連續(xù)型的變量，所以我們肯定要找一個(gè)取大于零的連續(xù)型的分布，因?yàn)槲覀儧](méi)有什么經(jīng)驗(yàn)，所以最容易想到的就是均勻分布，但是我試了一下發(fā)現(xiàn)效果很不好，可以執(zhí)行下面的代碼測(cè)試一下。

所以我們只好考慮考慮其他的分布，后來(lái)我想到了指數(shù)分布，指數(shù)分布描述的是等待一件事情發(fā)生所需時(shí)間的分布。

概率密度函數(shù)為：

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

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

所以指數(shù)分布的參數(shù)我們?cè)O(shè)置為整個(gè)觀測(cè)數(shù)據(jù)的均值。

# 設(shè)置先驗(yàn)分布
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

建立模型，使用馬爾可夫蒙特卡洛計(jì)算后驗(yàn)。

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

從上面后驗(yàn)的概率密度分布圖中我們可以看見(jiàn)兩個(gè)參數(shù)的分布有明顯的不同，\lambda_1 很有可能在 18 附近。 \lambda_2 很有可能在 23 附近。這意味著均值確實(shí)發(fā)生了變化。

均值發(fā)生變化的時(shí)間很有可能是 44 或 45 天的時(shí)候。我們直觀的看一下。

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

短信數(shù)量的均值發(fā)生變化說(shuō)明了很多事情，比如在第 45 天附近的時(shí)候加入了什么群，不停的有消息進(jìn)來(lái)，或者是一段新的戀情...

在之前我們對(duì)模型進(jìn)行了簡(jiǎn)化，我們現(xiàn)在嘗試將限制放開(kāi)一些，我們?cè)囈幌聝纱尉底兓瘯?huì)發(fā)生什么。

# -*- 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軸坐標(biāo)范圍
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")

這兩個(gè)分布還是有一定的差異性，說(shuō)明在條件放開(kāi)一些的情況下，短信數(shù)量還是發(fā)生了變化。

最后給出 pym3 的實(shí)現(xiàn)，現(xiàn)在 pym2 正在逐漸被棄用，目前被維護(hù)的是 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");
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書(shū)系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容