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

（點擊

上方藍字

，快速關注我們）

編譯：伯樂在線 - JLee

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

• 
  

  《

  第 0 節：導論

  》



• 
  

  《

  第 1 節：估計模型參數

  》



• 
  

  《第 2 節：模型檢驗》
  

第2節：模型檢驗

在這一節中，我們將用到兩種技術，旨在回答：

1. 
   

   模型和估計的參數是否對潛在數據有很好的模擬？



2. 
   

   給定兩個獨立的模型，哪個對潛在數據模擬的更好？

模型檢驗I：後驗估計檢驗

一種檢驗模型擬合的方法是後驗估計檢驗。這種方法很直觀，回顧上節中，我們通過收集 200,000 個 μ 的後驗分布樣本來對泊松分布的參數 μ進行估計，每個樣本都被認為是可信的參數值。

後驗預測檢驗需要從預測模型中產生新的數據。具體來說就是，我們已經估計了 200,000 個可信的泊松分布參數值μ，這意味著我們可以根據這些值建立 200,000 個泊松分布，然後從這些分布中隨機採樣。用公式表示為：

理論上，如果模型對潛在數據擬合的很好，那麼產生的數據應該和原始觀測數據近似。PyMC 提供一種方便的方式來從擬合模型中採樣。你可能已經注意到了上面模型應用中新的一行代碼：

y_pred = pm.Poisson("y_pred", mu=mu)

這幾乎和 y_est 一樣，只不過沒給觀測數據賦值。PyMC 把它當做隨機節點（和觀測節點相反），當 MCMC 採樣器運行時，它也從 y_est 中採集數據。

然後畫出y_pred並和觀測數據y_est比較。

import

json

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

statsmodels

.

api

as

sm

import

theano

.

tensor

as

tt

 

from

IPython

.

display

import

Image

 

%

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"

)

Applied

interval

-

transform to mu

and

added transformed mu_interval to

model

.

[

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

100

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

]

200000

of

200000

complete

in

63.5

sec

x_lim

=

60

burnin

=

50000

 

y_pred

=

trace

[

burnin

:

].

get_values

(

"y_pred"

)

mu_mean

=

trace

[

burnin

:

].

get_values

(

"mu"

).

mean

()

 

fig

=

plt

.

figure

(

figsize

=

(

10

,

6

))

fig

.

add_subplot

(

211

)

 

_

=

plt

.

hist

(

y_pred

,

range

=

[

0

,

x_lim

],

bins

=

x_lim

,

histtype

=

"stepfilled"

,

color

=

colors

[

1

])

  

_

=

plt

.

xlim

(

1

,

x_lim

)

_

=

plt

.

ylabel

(

"Frequency"

)

_

=

plt

.

title

(

"Posterior predictive distribution"

)

 

fig

.

add_subplot

(

212

)

 

_

=

plt

.

hist

(

messages

[

"time_delay_seconds"

].

values

,

range

=

[

0

,

x_lim

],

bins

=

x_lim

,

histtype

=

"stepfilled"

)

_

=

plt

.

xlabel

(

"Response time in seconds"

)

_

=

plt

.

ylabel

(

"Frequency"

)

_

=

plt

.

title

(

"Distribution of observed data"

)

 

plt

.

tight_layout

()

選擇正確的分布

我對上面的結果不是很滿意。理想情況下，我希望後驗預測分布和觀測數據的分布一定程度上近似。直觀上，如果我們正確估計了模型參數，那麼我們應該可以從模型中採樣得到類似的數據。結果很明顯不是這樣。

可能泊松分布不適合擬合這些數據。一種可選模型是負二項分布，特點和泊松分布很像，只是有兩個參數（μ 和 α），使得方差和均值獨立。回顧泊松分布只有一個參數 μ，既是均值，又是方差。

fig

=

plt

.

figure

(

figsize

=

(

10

,

5

))

fig

.

add_subplot

(

211

)

x_lim

=

70

mu

=

[

15

,

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

.

xlim

(

1

,

x_lim

)

_

=

plt

.

xlabel

(

"Response time in seconds"

)

_

=

plt

.

ylabel

(

"Probability mass"

)

_

=

plt

.

title

(

"Poisson distribution"

)

_

=

plt

.

legend

([

"$lambda$ = %s"

%

mu

[

0

],

                

"$lambda$ = %s"

%

mu

[

1

]])

 

# Scipy takes parameters n & p, not mu & alpha

def

get_n

(

mu

,

alpha

)

:

    

return

1.

/

alpha

*

mu

 

def

get_p

(

mu

,

alpha

)

:

    

return

get_n

(

mu

,

alpha

)

/

(

get_n

(

mu

,

alpha

)

+

mu

)

 

fig

.

add_subplot

(

212

)

 

a

=

[

2

,

4

]

 

for

i

in

np

.

arange

(

x_lim

)

:

    

plt

.

bar

(

i

,

stats

.

nbinom

.

pmf

(

i

,

n

=

get_n

(

mu

[

0

],

a

[

0

]),

p

=

get_p

(

mu

[

0

],

a

[

0

])),

color

=

colors

[

3

])

    

plt

.

bar

(

i

,

stats

.

nbinom

.

pmf

(

i

,

n

=

get_n

(

mu

[

1

],

a

[

1

]),

p

=

get_p

(

mu

[

1

],

a

[

1

])),

color

=

colors

[

4

])

 

_

=

plt

.

xlabel

(

"Response time in seconds"

)

_

=

plt

.

ylabel

(

"Probability mass"

)

_

=

plt

.

title

(

"Negative Binomial distribution"

)

_

=

plt

.

legend

([

"$mu = %s, / beta = %s$"

%

(

mu

[

0

],

a

[

0

]),

                

"$mu = %s, / beta = %s$"

%

(

mu

[

1

],

a

[

1

])])

 

plt

.

tight_layout

()

使用之前相同的數據集，繼續對負二項分布的參數進行估計。同樣地，使用均勻分布來估計 μ 和 α。模型可以表示為：

代碼：

Image("graphics/Neg Binomial Dag.png", width=400)

with

pm

.

Model

()

as

model

:

    

alpha

=

pm

.

Exponential

(

"alpha"

,

lam

=

.

2

)

    

mu

=

pm

.

Uniform

(

"mu"

,

lower

=

0

,

upper

=

100

)

    

    

y_pred

=

pm

.

NegativeBinomial

(

"y_pred"

,

mu

=

mu

,

alpha

=

alpha

)

    

y_est

=

pm

.

NegativeBinomial

(

"y_est"

,

mu

=

mu

,

alpha

=

alpha

,

observed

=

messages

[

"time_delay_seconds"

].

values

)

    

    

start

=

pm

.

find_MAP

()

    

step

=

pm

.

Metropolis

()

    

trace

=

pm

.

sample

(

200000

,

step

,

start

=

start

,

progressbar

=

True

)

Applied

log

-

transform

to

alpha

and

added transformed alpha_log

to

model

.

Applied

interval

-

transform

to

mu

and

added transformed mu_interval

to

model

.

[

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

100

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

]

200000

of

200000

complete

in

114.8

sec

_ = pm.traceplot(trace[burnin:], varnames=["alpha", "mu"])

我們看到上面的模型在平均回復時間 μ 附近有更大的不確定度：

• 
  

  泊松分布：17.5 到 18.5



• 
  

  負二項分布：16 到 21

同時，負二項分布模型有一個 1.4 到 2.4 的 α 參數，增加了估計參數 μ 的方差。讓我們看看後驗預測分布，是否和觀測數據更加近似。

x_lim

=

60

y_pred

=

trace

[

burnin

:

].

get_values

(

"y_pred"

)

 

fig

=

plt

.

figure

(

figsize

=

(

10

,

6

))

fig

.

add_subplot

(

211

)

 

fig

.

add_subplot

(

211

)

 

_

=

plt

.

hist

(

y_pred

,

range

=

[

0

,

x_lim

],

bins

=

x_lim

,

histtype

=

"stepfilled"

,

color

=

colors

[

1

])

  

_

=

plt

.

xlim

(

1

,

x_lim

)

_

=

plt

.

ylabel

(

"Frequency"

)

_

=

plt

.

title

(

"Posterior predictive distribution"

)

 

fig

.

add_subplot

(

212

)

 

_

=

plt

.

hist

(

messages

[

"time_delay_seconds"

].

values

,

range

=

[

0

,

x_lim

],

bins

=

x_lim

,

histtype

=

"stepfilled"

)

_

=

plt

.

xlabel

(

"Response time in seconds"

)

_

=

plt

.

ylabel

(

"Frequency"

)

_

=

plt

.

title

(

"Distribution of observed data"

)

 

plt

.

tight_layout

()

是的，和另一個相比，這兩個分布看起來更近似。根據後驗預測檢驗，負二項分布模型對潛在數據模擬更好。

如果你懷疑這種模型檢驗方法的嚴密性，貝葉斯方法還有一種更嚴密的分析方法。

模型檢驗II：貝葉斯因子

另一種模型檢驗的方法是計算貝葉斯因子，這種分析方法旨在兩個模型間的相互比較。

貝葉斯因子通常難以計算，因為需要整合全聯合概率分布。在低維空間中，整合是可能的，但是一旦你稍微增加建模的維度，整合全聯合概率分布的計算會變得昂貴且費時。

有一種可選的近似方法來計算貝葉斯因子。它需要把兩個模型進行比較，把它們融合成一個帶有一個模型指數（τ）的分層模型.這個指數在MCMC 過程中，會根據模型的可信度在兩個模型中切換，這樣，這個指數的軌跡就充分顯示模型 M1 相對模型 M2 的可信度。

Image("graphics/Bayes Factor DAG.png", width=540)

with

pm

.

Model

()

as

model

:

    

    

# Index to true model

    

prior_model_prob

=

0.5

    

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

    

tau

=

pm

.

Bernoulli

(

"tau"

,

prior_model_prob

)

    

    

# Poisson parameters

    

mu_p

=

pm

.

Uniform

(

"mu_p"

,

0

,

60

)

 

    

# Negative Binomial parameters

    

alpha

=

pm

.

Exponential

(

"alpha"

,

lam

=

0.2

)

    

mu_nb

=

pm

.

Uniform

(

"mu_nb"

,

lower

=

0

,

upper

=

60

)

 

    

y_like

=

pm

.

DensityDist

(

"y_like"

,

            

lambda

value

:

pm

.

switch

(

tau

,

                

pm

.

Poisson

.

dist

(

mu_p

).

logp

(

value

),

                

pm

.

NegativeBinomial

.

dist

(

mu_nb

,

alpha

).

logp

(

value

)

            

),

            

observed

=

messages

[

"time_delay_seconds"

].

values

)

    

    

start

=

pm

.

find_MAP

()

    

step1

=

pm

.

Metropolis

([

mu_p

,

alpha

,

mu_nb

])

    

step2

=

pm

.

ElemwiseCategoricalStep

(

vars

=

[

tau

],

values

=

[

0

,

1

])

    

trace

=

pm

.

sample

(

200000

,

step

=

[

step1

,

step2

],

start

=

start

)

 

_

=

pm

.

traceplot

(

trace

[

burnin

:

],

varnames

=

[

"tau"

])

Applied

interval

-

transform to mu_p

and

added transformed mu_p_interval to

model

.

Applied

log

-

transform to alpha

and

added transformed alpha_log to

model

.

Applied

interval

-

transform to mu_nb

and

added transformed mu_nb_interval to

model

.

[

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

100

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

]

200000

of

200000

complete

in

314.6

sec

根據下面的公式計算上面兩個模型的貝葉斯因子：

在上例中，我們沒有使用先驗概率，因此貝葉斯因子只是模型似然度的商。如果你發現你的 MCMC 樣本不在兩個模型間遊走，你可以引入先驗概率來幫助你更深入了解兩個模型。

# Compute the Bayes factor

prob_pois

=

trace

[

burnin

:

][

"tau"

].

mean

()

prob_nb

=

1

-

prob_pois

BF

=

(

prob_nb

/

prob_pois

)

*

(

prior_model_prob

/

(

1

-

prior_model_prob

))

print

(

"Bayes Factor: %s"

%

BF

)

Bayes Factor: 1.60114103387

貝葉斯因子大於 1 說明 M1（負二項分布）比 M2（泊松分布）更適合擬合數據。Jeffrey的貝葉斯因子優勢判別準則把貝葉斯因子為  1.60 標識為 M1 對 M2 有弱優勢。結合後驗預測檢驗和貝葉斯因子，可得出結論：負二項分布模型更優。

看完本文有收穫？請轉

發分享給更多人

關注「P

ython開發者」，提升Python技能

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

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

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