對於分類問題，邏輯回歸是常用的方法，比較簡單的方法就是調sklearn的包。但引入貝葉斯思想之後，做有類似簡單的方法嗎？當然有，比如pymc3就能做。

GLM: Logistic Regression?

docs.pymc.io圖標

上面這個是官方tutorial鏈接，不過基本粗略看下是看不懂的，pymc的坑這裡面也沒仔細說。所以我們下面一步步來借鳶尾花分類來說明下pymc的邏輯，以及如何做貝葉斯邏輯回歸。

Data set:

https://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data?

archive.ics.uci.edu

問題描述：對於Fisher flower Iris數據集，根據4個特徵，train一個邏輯回歸區分是不是Virginica。

Tip: 邏輯回歸的公式 p_{logit}=frac{1}{1+e^{-x^Teta}} ，其中x是features組成的向量，beta是對應的權重。這裡x有4個features，所以需要再加一個bias項，一共5維。

  • Step 1: 預處理，train test split

預處理要把class那一列從名詞轉成0-1，是virginica取1。注意這列數據的類型，不能是「object」，相當於還是字元串。。。

split方法有很多，可以直接調sklearn的包，也可以自己在pandas裏shuffle。注意，第一個坑來了，如果想用np array作為下一步的輸入，那麼從pandas轉成np array的時候，千萬不要亂reshape，保持原樣即可。

  • Step 2: 調包，pymc3

有了數據集之後，就可以用pymc的包了。準備工作：

import pymc3 as pm
from pymc3 import MvNormal, Normal, Binomial, Model

思考下我們的模型是怎麼樣的：第一步，假設 eta 具有某個先驗分佈。第二步， eta 與x相乘。第三步，過sigmoid函數得到分類結果的概率。第四步，使用計算得到的概率和training set的準確分類結果，計算likelihood。第五步，已知先驗和likelihood，根據貝葉斯公式，求 eta 的後驗。

下面我們有3種方法：

Method #1

按照官方教程的做法，我們得到以下的模型：

with Model() as iris_classify:
priors_iris = {"Intercept":pm.Normal.dist(mu=0,sd=10),
"sepal_length":pm.Normal.dist(mu=0,sd=10),
"sepal_width":pm.Normal.dist(mu=0,sd=10),
"petal_length":pm.Normal.dist(mu=0,sd=10),
"petal_width":pm.Normal.dist(mu=0,sd=10)}
pm.glm.GLM.from_formula(classification ~ sepal_length + sepal_width + petal_length + petal_width,
train, priors=priors_iris,family=pm.glm.families.Binomial())
map_estimate = pm.find_MAP(model=iris_classify)
print(map_estimate)

其中train是訓練集的pandas dataframe。sepal_length，sepal_width，petal_length，petal_width是四列features。

第一行，先創建一個Model叫iris_classify

第二行，建立一個叫priors_iris的字典，作為先驗。字典的key是4個features+1個Bias。每個key對應的value是pm.Normal.dist，mu = mean，sd = standard deviation，這裡假設是10。如果沒有先驗，那pymc會自動加一個mu = 0，sd= 10^{12} 的先驗，考慮到默認的sd太大了，我們還是自己設置prior比較好。

第三行，pm.glm.GLM.from_formula告訴pymc我們要用公式了。輸入公式時，要用train這個dataframe的各個column名字，並且不用管 eta 和intercept。family=。。Binomial（）這句意味著這是0-1分類問題。操作完畢之後，pymc自動識別classification這列為邏輯回歸的結果，對先驗的採樣來算likelihood。

第四行，用pymc自帶的函數求MAP

Method #2

上面用formula的方法封裝的太好了，並不是太建議直接用。方法2我們用5個獨立的pm.Normal來做同樣的事情。

with Model() as iris_classify2:
beta0 = pm.Normal("beta0",0,sd=10)
beta1 = pm.Normal("beta1",0,sd=10)
beta2 = pm.Normal("beta2",0,sd=10)
beta3 = pm.Normal("beta3",0,sd=10)
beta4 = pm.Normal("beta4",0,sd=10)
likelihood = pm.Bernoulli("Y",p=pm.math.sigmoid(data_x[:,0]*beta0+data_x[:,1]*beta1+
data_x[:,2]*beta2+data_x[:,3]*beta3+data_x[:,4]*beta4),observed=data_y)

此處beta0是intercept。第七行，likelihood的定義是一個bernoulli分佈，名字叫「Y」，概率p=sigmoid函數的結果，函數自變數是beta與data_x的各列相乘。data_x是數據集，其中第0列全是1，行數為訓練集樣本數m，列為5。最重要的事情來了，此處如果沒有observed=data_y，那這個bernoulli分佈只是一個普通的分佈，就像前面beta0~4一樣。有了observerd之後，pymc知道了，這是概率模型的約束條件。其中data_y是剛才的classification那一列。注意，data_y的shape是 (m，)，而不是（m，1），否則後續包括MAP的整個結果都是錯的！！！

一個細節是，pm.math.sigmoid其實是theano的sigmoid馬甲，截止今天，pymc的各種高級運算，都是藉助theano的tensor。

Method #3

問題來了，如果我們想做的完美一點，先驗是一個聯合高斯分佈怎麼辦？或者我們不想寫那麼多行beta。

with Model() as iris_classify3:
sigma = 10
beta = pm.MvNormal("beta",mu=np.zeros(5),cov=(sigma**2)*np.eye(5),shape=5) # (5,)
likelihood = pm.Bernoulli("Y",p=pm.math.sigmoid(pm.math.dot(data_x,beta)),observed=data_y)

此處定義beta是MvNormal，一共三個輸入值，其中mu.shape = (5,1)，cov是協方差矩陣，shape要寫=5或者(5,)，但別寫（5,1）。

第四行，同樣是一個有observed的bernoulli分佈，但概率p的shape也是(5,)，使用theano的馬甲數學函數們進行運算。data_x.shape=(m,5)，bete.shape=(5,)，data_y.shape=(5,)。

  • 後續的處理，sampling，trace，Posterior distribution 和 Posterior distribution predictive

with iris_classify:
trace = pm.sample(draws=5000,chains=1,tune=1000,start=map_estimate)

trace就是一系列sample。與frequentist的觀點不同，bayesian認為一切概率都是分佈，分佈在updating。所以前面的beta雖然和固定的data_x相乘，但實際上是很多個beta。這裡採樣pymc默認是NUTS，draws是後驗樣本數，chains是trace的數目，可以多個chains對比分析採樣結果有沒有問題，tune=1000是前1000個樣本丟掉（burnin），默認=500。

pm.traceplot(trace);
pm.plots.autocorrplot(trace,max_lag=20);

第一行的命令可以看我們後驗採樣的分佈，以及整個trace

第二行的命令看trace中sample之間的相關度，如果不靠譜可以加高一點burnin，還有加點thinning，比如：

trace_thinning = trace[::10]

上面是分析Posterior distribution，針對的beta。那麼針對新的beta如何預測test set上的效果呢？

比如我們可以用後驗beta的平均值或者MAP來預測：

beta_sample = trace["beta"]
beta_mean = beta_sample.mean(axis=0)

beta_MAP = map_estimate["beta"]

那麼後驗預測的結果用beta mean和beta MAP點乘data_x，再套一層sigmoid就行了。我們可以得到後驗情況下，不同測試集裏每個sample是virginica的概率p。之後有兩種做法，一種是machine-learning-decision-theoretic，對p做切分，>0.5為virginica，<0.5就不是；另一種是繼續做一個bernoulli分佈。最後效果還行，訓練集上有~95%的準確率，測試集懶得測試了（主要是作業沒要求。。。），有結果可以分享下。

題外話：用sklearn的logistics regression可以直接得到100%的準確度，大概5行以內代碼吧。。。。


推薦閱讀：
相關文章