本系列：
- 《第 0 節：導論》
- 《第 1 節：估計模型引數》
- 《第 2 節：模型檢驗》
- 《第 3 節：分層模型》
- 《第 4 節：貝葉斯迴歸》
- 待續
第4節：貝葉斯迴歸
之前，我們留了個問題：“我的回覆時間受聊天的物件的影響嗎？”我們針對每個對話的人都進行了模型引數估計。但是有時候我們想了解更多的影響因素，比如：星期幾，時刻等。可以運用 GLM（通用線性模型）更好地瞭解這些因素的影響。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 |
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是 −∞ 到 ∞ 上的連續變數時，可以考慮使用線性迴歸，表示為：
![5DZL5V]9HJ_ZOUC0ABQJR]L](http://jbcdn2.b0.upaiyun.com/2016/08/28d0add7fffe66d72e62a7670cdaacb6.png)
解讀為：響應y通常服從均值為 μ，標準差為 σ 的正態分佈，μ 描述為一組解釋變數 Xβ 的線性函式，零線截距為 β0.
連線函式
在不是對 −∞ 到 ∞上連續變數進行建模的場合，你需要使用連線函式來轉換響應域。對於泊松分佈，一種權威的連線函式是 log 連線函式，可以表示為：
這被認為是一個固定效用模型，對係數 β 的估計是基於整個對話人群，而不是基於單個的人估計獨自的引數（如第三節中的合併和區域性融合模型）。
固定效用泊松迴歸
為了用PyMC3建立泊松迴歸，我們需要對 μ 使用log連線函式。PyMC3中的潛在資料模型使用了theano 庫，因此需要使用下面的 theano 張量方法 theano.tensor.exp()。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 |
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) |
|
1 |
[-----------------100%-----------------] 200000 of 200000 complete in 137.5 sec |
|
1 |
_ = 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 的不一致變化。
邊緣密度圖和二元聯合密度圖
下圖顯示了邊緣密度圖（對角線上）和二元聯合密度圖（上下窗格）。這個圖對於理解協變數的相互作用很有用，在上例中，可以看到，若參與的人數增加，則零線截距下降。
|
1 |
_ = sns.pairplot(pm.trace_to_dataframe(trace[20000:]), plot_kws={'alpha':.5}) |
混合效用泊松迴歸
對上面的模型進行小的擴充套件，引入一個隨機截距引數，這樣可以針對每個人估計一個基線引數 β0，而其他引數則針對整個群體進行估計。對於每個人i的每條訊息j，可以表示為：
通過對每個人i引入這個隨機效用引數 β0，可以使模型針對每個人建立不同的基線，同時對整個群體估計協變數的效用。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 |
# 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) |
|
1 |
[-----------------100%-----------------] 200000 of 200000 complete in 207.3 sec |
|
1 |
_ = pm.traceplot(trace[20000:]) |
上面結果的截距很有趣：
- 每個人有不同的基本回復速率（正如第3節中兩個模型展示的那樣）
- 長訊息需要較長時間編寫，通常需要較長回覆時間
- 如果你週末給我發訊息，你很可能得到較慢的回覆
- 我傾向於在群組聊天中回覆更快
經過對每個協變數的效用進行計算，模型估計出引數 β0 的值如下。
|
1 2 |
_ = plt.figure(figsize=(5, 6)) _ = pm.forestplot(trace[20000:], varnames=['intercept'], ylabels=participants) |
打賞支援我翻譯更多好文章，謝謝！
打賞譯者
打賞支援我翻譯更多好文章，謝謝！
任選一種支付方式
