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

（點擊

上方藍字

，快速關注我們）

編譯：伯樂在線 - JLee

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

• 
  

  《

  第 0 節：導論

  》



• 
  

  《

  第 1 節：估計模型參數

  》



• 
  

  《

  第 2 節：模型檢驗

  》



• 
  

  《

  第 3 節：分層模型

  》



• 
  

  《第 4 節：貝葉斯回歸》

第4節：貝葉斯回歸

之前，我們留了個問題：「我的回復時間受聊天的對象的影響嗎？」我們針對每個對話的人都進行了模型參數估計。但是有時候我們想了解更多的影響因素，比如：星期幾，時刻等。可以運用 GLM（通用線性模型）更好地了解這些因素的影響。

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

seaborn

.

apionly

as

sns

import

statsmodels

.

api

as

sm

import

theano

.

tensor

as

tt

 

from

sklearn

import

preprocessing

 

%

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"

)

線性回歸

當響應y是 ?∞ 到  ∞ 上的連續變數時，可以考慮使用線性回歸，表示為：

解讀為：響應y通常服從均值為 μ，標準差為 σ 的正態分布，μ 描述為一組解釋變數 Xβ 的線性函數，零線截距為 β0.

連接函數

在不是對 ?∞ 到 ∞上連續變數進行建模的場合，你需要使用連接函數來轉換響應域。對於泊松分布，一種權威的連接函數是 log 連接函數，可以表示為：

這被認為是一個固定效用模型，對係數 β 的估計是基於整個對話人群，而不是基於單個的人估計獨自的參數（如第三節中的合併和局部融合模型）。

固定效用泊松回歸

為了用PyMC3建立泊松回歸，我們需要對 μ 使用log連接函數。PyMC3中的潛在數據模型使用了theano 庫，因此需要使用下面的 theano 張量方法 theano.tensor.exp()。

X

=

messages

[[

"is_weekend"

,

"day_of_week"

,

"message_length"

,

"num_participants"

]].

values

_

,

num_X

=

X

.

shape

 

with

pm

.

Model

()

as

model

:      

    

intercept

=

pm

.

Normal

(

"intercept"

,

mu

=

0

,

sd

=

100

)

    

beta_message_length

=

pm

.

Normal

(

"beta_message_length"

,

mu

=

0

,

sd

=

100

)

    

beta_is_weekend

=

pm

.

Normal

(

"beta_is_weekend"

,

mu

=

0

,

sd

=

100

)

    

beta_num_participants

=

pm

.

Normal

(

"beta_num_participants"

,

mu

=

0

,

sd

=

100

)

    

    

mu

=

tt

.

exp

(

intercept

                +

beta_message_length

*

messages

.

message_length

                +

beta_is_weekend

*

messages

.

is_weekend

                +

beta_num_participants

*

messages

.

num_participants

)

    

    

y_est

=

pm

.

Poisson

(

"y_est"

,

mu

=

mu

,

observed

=

messages

[

"time_delay_seconds"

].

values

)

    

    

start

=

pm

.

find_MAP

()

    

step

=

pm

.

Metropolis

()

    

trace

=

pm

.

sample

(

200000

,

step

,

start

=

start

,

progressbar

=

True

)

[-----------------100%-----------------] 200000 of 200000 complete in 137.5 sec

_ = pm.traceplot(trace)

從上面的結果中可以看到，零線截距 β0 估值在 2.4 到 2.8，那麼這說明什麼呢？

不幸的是，解釋泊松回歸的參數比簡單線性回歸（y=β X）更費勁，在這個線性回歸中，我們可以說 x 每增加一個單位，y 增加 β 個單位。但是，泊松回歸中我們需要考慮連接函數。following cross validated post詳細推導了下面的公式。

對泊松模型，x 變化一個單位，相應地 y 變化 y( e^β – 1)

從中得出的主要結論是，x 變化的影響取決於當前的 y 值，不像簡單線性回歸，x 同樣的變化卻引起 y 的不一致變化。

邊緣密度圖和二元聯合密度圖

下圖顯示了邊緣密度圖（對角線上）和二元聯合密度圖（上下窗格）。這個圖對於理解協變數的相互作用很有用，在上例中，可以看到，若參與的人數增加，則零線截距下降。

_ = sns.pairplot(pm.trace_to_dataframe(trace[20000:]), plot_kws={"alpha":.5})

混合效用泊松回歸

對上面的模型進行小的擴展，引入一個隨機截距參數，這樣可以針對每個人估計一個基線參數 β0，而其他參數則針對整個群體進行估計。對於每個人i的每條消息j，可以表示為：

通過對每個人i引入這個隨機效用參數 β0，可以使模型針對每個人建立不同的基線，同時對整個群體估計協變數的效用。

# Convert categorical variables to integer

le

=

preprocessing

.

LabelEncoder

()

participants_idx

=

le

.

fit_transform

(

messages

[

"prev_sender"

])

participants

=

le

.

classes_

n_participants

=

len

(

participants

)

 

with

pm

.

Model

()

as

model

:

 

    

intercept

=

pm

.

Normal

(

"intercept"

,

mu

=

0

,

sd

=

100

,

shape

=

n_participants

)

    

slope_message_length

=

pm

.

Normal

(

"slope_message_length"

,

mu

=

0

,

sd

=

100

)

    

slope_is_weekend

=

pm

.

Normal

(

"slope_is_weekend"

,

mu

=

0

,

sd

=

100

)

    

slope_num_participants

=

pm

.

Normal

(

"slope_num_participants"

,

mu

=

0

,

sd

=

100

)

    

    

mu

=

tt

.

exp

(

intercept

[

participants_idx

]

                +

slope_message_length

*

messages

.

message_length

                +

slope_is_weekend

*

messages

.

is_weekend

                +

slope_num_participants

*

messages

.

num_participants

)

    

    

y_est

=

pm

.

Poisson

(

"y_est"

,

mu

=

mu

,

observed

=

messages

[

"time_delay_seconds"

].

values

)

    

    

start

=

pm

.

find_MAP

()

    

step

=

pm

.

Metropolis

()

    

trace

=

pm

.

sample

(

200000

,

step

,

start

=

start

,

progressbar

=

True

)

[-----------------100%-----------------] 200000 of 200000 complete in 207.3 sec

_ = pm.traceplot(trace[20000:])

上面結果的截距很有趣：

• 
  

  每個人有不同的基本回復速率（正如第3節中兩個模型展示的那樣）



• 
  

  長消息需要較長時間編寫，通常需要較長回復時間



• 
  

  如果你周末給我發消息，你很可能得到較慢的回復



• 
  

  我傾向於在群組聊天中回復更快

經過對每個協變數的效用進行計算，模型估計出參數 β0 的值如下。

_

=

plt

.

figure

(

figsize

=

(

5

,

6

))

_

=

pm

.

forestplot

(

trace

[

20000

:

],

varnames

=

[

"intercept"

],

ylabels

=

participants

)

看完本文有收穫？請轉

發分享給更多人

關注「P

ython開發者」，提升Python技能

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

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

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