用 Python 進行貝葉斯模型建模（1）

（點擊

上方藍字

，快速關注我們）

編譯：伯樂在線 - JLee

如有好文章投稿，請點擊 → 這裡了解詳情

本系列：

• 
  

  《

  第 0 節：導論

  》



• 
  

  《第 1 節：估計模型參數》

第1節：估計模型參數

在這一節，我們將討論貝葉斯方法是如何思考數據的，我們怎樣通過 MCMC 技術估計模型參數。

from

IPython

.

display

import

Image

import

matplotlib

.

pyplot

as

plt

import

numpy

as

np

import

pandas

as

pd

import

pymc3

as

pm

import

scipy

import

scipy

.

stats

as

stats

import

scipy

.

optimize

as

opt

import

statsmodels

.

api

as

sm

%

matplotlib

inline

plt

.

style

.

use

(

"bmh"

)

colors

=

[

"#348ABD"

,

"#A60628"

,

"#7A68A6"

,

"#467821"

,

"#D55E00"

,

"#CC79A7"

,

"#56B4E9"

,

"#009E73"

,

"#F0E442"

,

"#0072B2"

]

messages

=

pd

.

read_csv

(

"data/hangout_chat_data.csv"

)

貝葉斯方法如何思考數據？

當我開始學習如何運用貝葉斯方法的時候，我發現理解貝葉斯方法如何思考數據是很有用的。設想下述的場景：

一個好奇的男孩每天觀察從他家門前經過的汽車的數量。他每天努力地記錄汽車的總數。一個星期過去後，他的筆記本上記錄著下面的數字：12,33,20,29,20,30,18

從貝葉斯方法的角度看，這個數據是由隨機過程產生的。但是，既然數據被觀測，它便固定了並且不會改變。這個隨機過程有些模型參數被固定了。然而，貝葉斯方法用概率分布來表示這些模型參數的不確定性。

由於這個男孩調查的是計數（非負整數），一個通常的做法是用泊松分布對數據（如隨機過程）建模。泊松分布只有一個參數 μ，它既是數據的平均數，也是方差。下面是三個不同 μ 值的泊松分布。

代碼：

fig

=

plt

.

figure

(

figsize

=

(

11

,

3

))

ax

=

fig

.

add_subplot

(

111

)

x_lim

=

60

mu

=

[

5

,

20

,

40

]

for

i

in

np

.

arange

(

x_lim

)

:

plt

.

bar

(

i

,

stats

.

poisson

.

pmf

(

mu

[

0

],

i

),

color

=

colors

[

3

])

plt

.

bar

(

i

,

stats

.

poisson

.

pmf

(

mu

[

1

],

i

),

color

=

colors

[

4

])

plt

.

bar

(

i

,

stats

.

poisson

.

pmf

(

mu

[

2

],

i

),

color

=

colors

[

5

])

_

=

ax

.

set_xlim

(

0

,

x_lim

)

_

=

ax

.

set_ylim

(

0

,

0.2

)

_

=

ax

.

set_ylabel

(

"Probability mass"

)

_

=

ax

.

set_title

(

"Poisson distribution"

)

_

=

plt

.

legend

([

"$mu$ = %s"

%

mu

[

0

],

"$mu$ = %s"

%

mu

[

1

],

"$mu$ = %s"

%

mu

[

2

]])

在上一節中，我們引入我的 hangout 聊天數據集。特別地，我對我的回復時間（response_time）感興趣。鑒於 response_time 是計數數據，我們可以用泊松分布對其建模，並估計參數 μ 。我們將用頻率論方法和貝葉斯方法兩種方法來估計。

fig

=

plt

.

figure

(

figsize

=

(

11

,

3

))

_

=

plt

.

title

(

"Frequency of messages by response time"

)

_

=

plt

.

xlabel

(

"Response time (seconds)"

)

_

=

plt

.

ylabel

(

"Number of messages"

)

_

=

plt

.

hist

(

messages

[

"time_delay_seconds"

].

values

,

range

=

[

0

,

60

],

bins

=

60

,

histtype

=

"stepfilled"

)

頻率論方法估計μ

在進入貝葉斯方法之前，讓我們先看一下頻率論方法估計泊松分布參數。我們將使用優化技術使似然函數最大。

下面的函數poisson_logprob()返回在給定泊松分布模型和參數值的條件下，觀測數據總體的可能性。用方法opt.minimize_scalar找到在觀測數據基礎上參數值 μ 的最可信值（最大似然）。該方法的機理是，這個優化技術會自動迭代可能的mu值直到找到可能性最大的值。

y_obs

=

messages

[

"time_delay_seconds"

].

values

def

poisson_logprob

(

mu

,

sign

=-

1

)

:

return

np

.

sum

(

sign

*

stats

.

poisson

.

logpmf

(

y_obs

,

mu

=

mu

))

freq_results

=

opt

.

minimize_scalar

(

poisson_logprob

)

%

time

print

(

"The estimated value of mu is: %s"

%

freq_results

[

"x"

])

The estimated value of mu is: 18.0413533867

CPU times: user 63 μs, sys: 1e+03 ns, total: 64 μs

Wall time: 67.9 μs

所以，μ 的估計值是18.0413533867。優化技術沒有對不確定度進行評估，它只返回一個點，效率很高。

下圖描述的是我們優化的函數。對於每個μ值，圖線顯示給定數據和模型在μ處的似然度。優化器以登山模式工作——從曲線上隨機一點開始，不停向上攀登直到達到最高點。

x

=

np

.

linspace

(

1

,

60

)

y_min

=

np

.

min

([

poisson_logprob

(

i

,

sign

=

1

)

for

i

in

x

])

y_max

=

np

.

max

([

poisson_logprob

(

i

,

sign

=

1

)

for

i

in

x

])

fig

=

plt

.

figure

(

figsize

=

(

6

,

4

))

_

=

plt

.

plot

(

x

,

[

poisson_logprob

(

i

,

sign

=

1

)

for

i

in

x

])

_

=

plt

.

fill_between

(

x

,

[

poisson_logprob

(

i

,

sign

=

1

)

for

i

in

x

],

y_min

,

color

=

colors

[

0

],

alpha

=

0.3

)

_

=

plt

.

title

(

"Optimization of $mu$"

)

_

=

plt

.

xlabel

(

"$mu$"

)

_

=

plt

.

ylabel

(

"Log probability of $mu$ given data"

)

_

=

plt

.

vlines

(

freq_results

[

"x"

],

y_max

,

y_min

,

colors

=

"red"

,

linestyles

=

"dashed"

)

_

=

plt

.

scatter

(

freq_results

[

"x"

],

y_max

,

s

=

110

,

c

=

"red"

,

zorder

=

3

)

_

=

plt

.

ylim

(

ymin

=

y_min

,

ymax

=

0

)

_

=

plt

.

xlim

(

xmin

=

1

,

xmax

=

60

)

上述優化方法估計泊松分布的參數值為 18 .我們知道，對任意一個泊松分布，參數 μ 代表的是均值和方差。下圖描述了這個分布。

fig

=

plt

.

figure

(

figsize

=

(

11

,

3

))

ax

=

fig

.

add_subplot

(

111

)

x_lim

=

60

mu

=

np

.

int

(

freq_results

[

"x"

])

for

i

in

np

.

arange

(

x_lim

)

:

plt

.

bar

(

i

,

stats

.

poisson

.

pmf

(

mu

,

i

),

color

=

colors

[

3

])

_

=

ax

.

set_xlim

(

0

,

x_lim

)

_

=

ax

.

set_ylim

(

0

,

0.1

)

_

=

ax

.

set_xlabel

(

"Response time in seconds"

)

_

=

ax

.

set_ylabel

(

"Probability mass"

)

_

=

ax

.

set_title

(

"Estimated Poisson distribution for Hangout chat response time"

)

_

=

plt

.

legend

([

"$lambda$ = %s"

%

mu

])

上述泊松分布模型和 μ 的估計表明，觀測小於10 或大於 30 的可能性很小，絕大多數的概率分布在 10 和 30 之間。但是，這不能反映我們觀測到的介於 0 到 60 之間的數據。

貝葉斯方法估計 μ

如果你之前遇到過貝葉斯理論，下面的公式會看起來很熟悉。我從沒認同過這個框架直到我讀了John K. Kruschke的書「貝葉斯數據分析」，從他優美簡潔的貝葉斯圖表模型視角認識這個公式。

代碼：

Image("graphics/Poisson-dag.png", width=320)

上述模式可以解讀如下（從下到上）：

• 
  

  對每一個對話（i）觀測計數數據（y）



• 
  

  數據由一個隨機過程產生，我們認為可以用泊松分布模擬



• 
  

  泊松分布只有一個參數，介於 0 到 60 之間

由於我們沒有任何依據對μ在這個範圍內進行預估，就用均勻分布對 μ 建模

神奇的方法MCMC

下面的動畫很好的演示了馬爾科夫鏈蒙特卡羅方法（MCMC）的過程。MCMC採樣器從先驗分布中選取參數值，計算從這些參數值決定的某個分布中得到觀測數據的可能性。

這個計算可以作為MCMC採樣器的導引。從參數的先驗分布中選取值，然後計算給定數據條件下這些參數值可能性——導引採樣器到更高概率的區域。

與上面討論的頻率論優化技術在概念上有相似之處，MCMC採樣器向可能性最高的區域運動。但是，貝葉斯方法不關心尋找絕對最大值，而是遍歷和收集概率最高區域附近的樣本。所有收集到的樣本都可認為是一個可信的參數。

Image(url="graphics/mcmc-animate.gif")

with

pm

.

Model

()

as

model

:

mu

=

pm

.

Uniform

(

"mu"

,

lower

=

0

,

upper

=

60

)

likelihood

=

pm

.

Poisson

(

"likelihood"

,

mu

=

mu

,

observed

=

messages

[

"time_delay_seconds"

].

values

)

start

=

pm

.

find_MAP

()

step

=

pm

.

Metropolis

()

trace

=

pm

.

sample

(

200000

,

step

,

start

=

start

,

progressbar

=

True

)

Applied

interval

-

transform to mu

and

added transformed mu_interval to

model

.

[

-----------------

100

%-----------------

]

200000

of

200000

complete

in

48.3

sec

上面的代碼通過遍歷 μ 的後驗分布的高概率區域，收集了 200,000 個 μ 的可信樣本。下圖（左）顯示這些值的分布，平均值與頻率論的估值（紅線）幾乎一樣。但是，我們同時也得到了不確定度，μ 值介於 17 到 19 之間都是可信的。我們稍後會看到這個不確定度很有價值。

_ = pm.traceplot(trace, varnames=["mu"], lines={"mu": freq_results["x"]})

丟棄早期樣本（壞點）

你可能疑惑上面 MCMC 代碼中pm.find.MAP()的作用。MAP 代表最大後驗估計，它幫助 MCMC 採樣器尋找合適的採樣起始點。理想情況下，模型從高概率區域開始，但有時不是這樣。因此，樣本集中的早期樣本（壞點）經常被丟棄。

fig

=

plt

.

figure

(

figsize

=

(

11

,

3

))

plt

.

subplot

(

121

)

_

=

plt

.

title

(

"Burnin trace"

)

_

=

plt

.

ylim

(

ymin

=

16.5

,

ymax

=

19.5

)

_

=

plt

.

plot

(

trace

.

get_values

(

"mu"

)[

:

1000

])

fig

=

plt

.

subplot

(

122

)

_

=

plt

.

title

(

"Full trace"

)

_

=

plt

.

ylim

(

ymin

=

16.5

,

ymax

=

19.5

)

_

=

plt

.

plot

(

trace

.

get_values

(

"mu"

))

模型的收斂性

樣本軌跡

由於上面模型對 μ 的估計不能表明估計的效果很好,可以採用一些檢驗手段。首先，觀察樣本輸出，樣本的軌跡應該輕微上下浮動，就像一條毛蟲。如果你看到軌跡上下跳躍或滯留在某點，這就是模型不收斂的標誌，通過 MCMC 採樣器得到的估計沒有意義。

自相關圖

也可采另一種測試，自相關測試（見下圖），這是評估MCMC採樣鏈中樣本前後相關關係的一種方法。相關性弱的樣本比相關性強的樣本能夠給參數估計提供更多的信息。

直觀上看，自相關圖應該快速衰減到0附近，然後在 0 附近上下波動。如果你的自相關圖沒有衰減，通常這是弱融合的標誌，你需要重新審查你的模型選擇（如似然度）和抽樣方法（如Metropolis）

_ = pm.autocorrplot(trace[:2000], varnames=["mu"])

看完本文有收穫？請轉

發分享給更多人

關注「P

ython開發者」，提升Python技能

喜歡這篇文章嗎？立刻分享出去讓更多人知道吧！

本站內容充實豐富，博大精深，小編精選每日熱門資訊，隨時更新，點擊「搶先收到最新資訊」瀏覽吧！

請您繼續閱讀更多來自 Python 的精彩文章: