後驗機率分佈的抽樣

【歡迎轉發分享，轉載請註明出處】

不少介紹貝葉斯後驗機率的書中都經常使用醫學化驗爲例子，好比一種方法的檢出陽性率是多少之類的。 假設有一種檢驗方法有95%的靈敏度可以識別出一我的是否是吸血鬼vampirism,即Pr(+|vampire) = 0.95，也就是針對一個吸血鬼可以95%的可能性識別出來。固然，它也有誤診的時候，Pr(+|mortal) = 0.01，即對一個正常人，也有1%的可能被診斷爲吸血鬼。此外，在人羣中，咱們還知道，吸血鬼佔的比例是很小很小的，只有0.1%，即Pr(vampire) = 0.001。在這種狀況下，若是一個「人」被診斷爲陽性，那麼他是正常人仍是吸血鬼？機率是多少？若是咱們使用貝葉斯的方法，不難作出以下判斷：

Pr(vampire|+)=Pr(+|vampire)Pr(vampire)/Pr(+)

其中，

Pr(+)=Pr(+|vampire)Pr(vampire)+Pr(+|mortal)(1−Pr(vampire))

在R語言中，咱們能夠經過下面的過程來計算：

PrPV <- 0.95
PrPM <- 0.01
PrV <- 0.001
PrP <- PrPV*PrV + PrPM*(1-PrV)
(PrVP <- PrPV*PrV / PrP)

## [1] 0.08683729

從上面的結果能夠看出，若是一個「人」被診斷爲吸血鬼，那麼他只有8.7%的可能性是一個真正的吸血鬼。這和咱們的感受有些相悖，這麼高靈敏度(95%)的方法，只有8.7%的可能性。放在現實生活中，不少這樣的檢驗，好比HIV檢驗等。這種檢驗的一個特色是，羣體中陽性個體不多，因此即使是檢驗結果是陽性，也不能保證該個體就是患者。 不過若是我咱們換一個角度來看上面的例子，就很容易理解了，咱們不用機率問題，而是用具體數值。好比 （1）一個100,000的羣體中，有100人是吸血鬼 （2）每100個吸血鬼可以檢驗出95個陽性 （3）在剩下99,900人中，有999我的的檢驗也是陽性 這時，若是對全部100,000我的都進行檢驗，那麼有多少陽性結果是真正吸血鬼？

Pr(vampire|+)=95/(95+999) =0.087

這個時候就容易理解多了。與一些抽象的機率問題相比，人們對具體計數的問題更易於理解。 本章咱們繼續研究機率分佈的問題，不過咱們會在後驗機率中隨機抽樣，經過具體計數的方式來理解機率分佈。 爲何要用抽樣呢？首先，使用抽樣能夠避免複雜的微積分等數學問題，雖然不少科研研究人員有必定的數學基礎，可是對於微積分依然敬而遠之。而經過抽樣，就能夠把微積分問題爲題轉化爲頻率問題，接下來對頻率問題的處理對不少人就駕輕就熟了。其次，不少計算後驗機率的方法並無具體的公式，好比使用MCMC獲得的後驗機率就是經過隨機抽樣獲得的，因此沒法使用微積分來描述後驗機率分佈。

在網格估計的後驗機率分佈中抽樣

以前介紹過網格估計法。下面是經過R語言，對擲地球實驗的後驗機率的網格估計，採用1000個網格點。

p_grid <- seq(from=0 , to=1 , length.out=1000)
prior <- rep(1 , 1000 )
likelihood <- dbinom( 6 , size=9 , prob=p_grid)
posterior <- likelihood * prior
posterior <- posterior / sum(posterior)

下面咱們在上述模擬的後驗機率分佈中進行10000次隨機抽樣。

samples <- sample(p_grid, size = 10000, prob = posterior, replace = T)
plot(samples)

上圖是這10000次隨機抽樣的分佈狀況。能夠看出在機率爲0.6的時候抽樣點密度很大，在<0.25的時候，抽樣點不多。咱們換一種形式來看這個圖，改爲密度分佈圖，以下：

library(rethinking)

## rethinking (Version 1.59)

dens(samples, xlab = "Proportion water(p)")

這時這張經過抽樣估計的後驗機率分佈圖就和咱們真實獲得的後驗機率分佈圖很類似了。並且隨着抽樣樣本的增大，抽樣對後驗機率的估計會也來越精確。

抽樣結果的應用

一旦你的模型獲得了後驗機率分佈，那麼模型的任務就算完成了。可是你的任務纔算剛剛開始。你須要進一步對後驗機率分佈進行歸納總結，以及解釋。具體要根據你的研究目的來決定，通常常見的問題有：

在某一界值內的後驗機率是多少？ 在兩個界值參數之間的後驗機率是多少？ 選取什麼界值會使獲得後驗機率分佈的5%？ 等等。。。

1. 界值區間

在擲地球的例子中，若是我要問海洋覆蓋率<0.5的後驗機率是多少？那麼，經過網格估計法獲得的後驗機率中，把小於0.5的部分的後驗機率加起來就能夠了。

sum(posterior[p_grid<0.5])

## [1] 0.1718746

也就是說，海洋覆蓋率小於50%的後驗機率是17%左右。看起來是很簡單的，可是這是咱們用網格估計法獲得的，而網格估計的缺點就是當參數特別多的時候，運算量迅速增長，因此應用範圍很窄。下面咱們使用抽樣的方法從新回答上面的問題，即只要把左右抽樣樣本中後驗機率<0.5的樣本相加，而後除以總抽樣次數(10000次)，就能夠估算出來咱們想要的結果。

sum(samples < 0.5)/10000

## [1] 0.1704

這就和咱們網格估計法的結果很接近了。(下圖左)

一樣地，若是是海洋覆蓋率在50%-75%的後驗機率，（下圖右）用抽樣估計以下：

sum(samples > 0.5 & samples < 0.75)/10000

## [1] 0.5981

2. 置信區間

好比咱們要想知道海洋覆蓋率是多少的時候，後驗機率的累積達到80%，這個時候用網格估計法就不太方便了，而使用抽樣則很簡單。下圖左

quantile(samples, 0.8)

##       80% 
## 0.7647648

一樣的，若是咱們想求中間80%的後驗機率分佈（下圖右），能夠

quantile(samples, c(0.1,0.9))

##       10%       90% 
## 0.4504505 0.8178178

上述後驗機率都是近似正態分佈的，因此經過百分區間就可以反映出整個分佈的狀況。若是後驗機率分佈不是正太分佈，那麼就有問題了。好比，咱們擲地球，前3次都擲到了海洋，這是，咱們使用網格進行後驗機率分佈，而後再估計中間50%的置信區間：

p_grid <- seq( from=0 , to=1 , length.out=1000 )
prior <- rep(1,1000)
likelihood <- dbinom( 3, size=3 , prob=p_grid )
posterior <- likelihood * prior
posterior <- posterior / sum(posterior)
samples <- sample( p_grid , size=1e4 , replace=TRUE , prob=posterior )
quantile(samples,c(0.25,0.75))

##       25%       75% 
## 0.7087087 0.9309309

中間50%的百分區間在0.70-0.93之間，可是這個區間內並無包含最高的後驗機率分佈。在本例中，最高的後驗機率分佈在p=1處。以下圖左。若是經過此區間來描述後驗機率分佈就存在嚴重的誤導。

這個時候咱們能夠用另外一個參數來描述，最高後驗機率密度區間(HPDI)，它是包含必定累積機率的最窄區間(如上圖右）。經過rethinking包中的HPDI函數能夠來求解：

HPDI(samples, prob = 0.5)

##      |0.5      0.5| 
## 0.8418418 0.9989990

這是這個區間就會包含後驗機率最大的值，此時的區間大小是最窄的，只有0.998-0.841=0.157。

百分區間和HPDI對於近似正態分佈的後驗機率分佈來說，結果是很類似的，差異就是非正態分佈的狀況。HPDI在估計的時候計算量比較大，並且對抽樣大小很敏感。

總之，在百分區間和HPDI估計差異很大的時候，那麼你最好可以把整個後驗分佈機率的狀況都畫出來，以避免產生誤導。

爲何是95%？
咱們平時最經常使用到的就是95%的置信區間，也就是咱們的參數有5%的可能性不在這個區間內。一般統計學檢驗的水平也是5%，即p<0.05。除了習慣外，恐怕沒有更科學的解釋了。這個鍋可能還要Ronald Fisher來背：「P=0.05的標準誤差是1.96，大概接近2；用這個數做爲判斷標準誤差是否有統計學差別是很方便的。」儘管不少人並不認同由於方便就使用0.05，但它彷佛成了約定俗成的一個判斷標準。在Fisher的晚年，他也很積極的反對使用0.05做爲界值。
可是，若是不用0.05，那麼要用什麼數來做爲判斷標準呢？目前彷佛尚未一致的結論。不過，你要知道0.05並非惟一的選擇。

置信區間意味着什麼？
咱們一般認爲95%的置信區間就是真實的參數有95%的可能性落在這個區間範圍以內。對於嚴格非貝葉斯者來講，這是不對的，由於他們不容許經過幾率來描述參數的不肯定性。他們認爲，若是你大量重複這個實驗，那麼你獲得區間的95%範圍會包含這個真實的參數。也就是參數是固定的。關於95%置信區間的解釋，不少人也沒有弄清楚，有些用了非貝葉斯觀點來解釋95%置信區間，有些則用了貝葉斯機率來解釋。
無論你是貝葉斯者仍是非貝葉斯者，95%的置信區間就意味着真的是95%嗎？實際上是咱們過於自信了。還記得以前提到的「小世界和大世界」嗎？這個95%僅僅是對你的「小世界」，即模型，來講的；對於真實的「大世界」，未必是你推測的95%。

3. 點估計

使用點估計來描述貝葉斯後驗機率分佈並非一個很好的選擇，一個點很難來描述整個分佈狀況。若是你就是想用一個點來描述呢？思考咱們以前提到的例子，在擲地球的實驗中，咱們前三次連續觀察到了海洋，你能夠選擇的點有最大數、平均數、中位數。

p_grid[which.max(posterior)]

## [1] 1

mean(samples)

## [1] 0.8010685

median(samples)

## [1] 0.8433433

這三個數差異很大，你應該選哪個呢？該怎麼選擇？

咱們能夠選擇一個損失函數loss function，而後基於整個分佈和咱們選擇的估計點，用該損失函數來估計選擇該估計點形成的損失大小。統計學者和博弈論者對損失函數一直很感興趣，不過在科學界用的並很少。那麼什麼是損失函數呢 舉個例子，我要和你打賭，賭海洋表面覆蓋率具體是多少，你告訴我一個你認爲最正確的地球表面海洋覆蓋率，若是你回答的正確，我就給你100美圓。若是你回答的不對，那麼我會根據你的回答和正確答案之間的差異大小，扣除相應的錢。也就是越接近正確答案，你獲得錢越多，反之越少。你由於沒有答對而形成損失大小和你的答案距離真實值的大小成正比。 若是你有後驗機率分佈，你確定要根據後驗機率分佈來選擇一個損失錢最少的答案。假如你選擇了海洋覆蓋率p=0.5做爲你的答案，那麼咱們根據後驗機率分佈計算一下此時損失多少：

sum(posterior*abs(0.5-p_grid))

## [1] 0.3128752

也就是，若是你選擇0.5做爲答案，你可能會損失3成的錢。那麼選多少做爲答案會損失最少呢？咱們能夠把因此得可能都計算一遍。

loss <- sapply(p_grid, function(d) sum(posterior*abs(d - p_grid)))
p_grid[which.min(loss)]

## [1] 0.8408408

整個損失函數以下圖：

最小值出如今0.84附近，也就是若是你回答0.84，那麼你損失的錢會是最少的。

因此爲了找一個相對準確的點估計值，你必須有一個損失函數來評估不一樣選擇形成的損失大小。可是不一樣的損失函數可能獲得不一樣的結果。上述例子中，用差值的絕對值做爲損失函數，獲得的點估計值很接近中位數，若是選擇差值的平方做爲損失函數，可能獲得的點估計值就很接近均數了。固然，若是後驗機率分佈爲近似正態分佈，那麼這兩個損失函數最後估計的結果可能差異不大。 在科學領域，人們彷佛還不太習慣使用損失函數，通常做點估計的時候一般習慣性的用均值或者中位數等。而更科學合理的方法是經過損失函數來選擇一個最合適的點。

經過抽樣進行模擬預測

1. 啞數據

在以前提到的擲地球的實驗中，咱們經過實驗，不只可以推斷出後驗機率分佈，還可以模擬出模型的預測。似然函數是能夠雙向工做的，既能夠經過觀測來推斷每個觀測發生的可能性；也能夠在給定參數下定義全部可能的觀測分佈，以此進行抽樣模擬或預測。

咱們把模擬獲得的數據稱之爲啞數據dummy data，即咱們真實觀測數據的「替代品」。對於擲地球實驗的啞數據，能夠經過下面的二項分佈似然函數來模擬。

其中w是觀測到海洋的次數，n是總的投擲數。當n=2，即咱們投擲兩次，咱們可能有3種結果：0次海洋、1次海洋和2次海洋。若是地球表面海洋覆蓋率是0.7，那麼咱們觀測到這3種結果的可能性爲：

dbinom(0:2,size = 2, prob = 0.7)

## [1] 0.09 0.42 0.49

也就是咱們有9%的機率觀測到0次海洋，42%的可能觀測到1次海洋，49%的可能觀測到2次海洋。 上面是咱們經過似然函數模擬獲得觀衆觀測的機率，若是咱們想獲得具體的模擬觀測數據，即啞數據，則可使用rbinom函數來實現：

rbinom(10, size = 2, prob = 0.7)

##  [1] 1 2 1 2 2 1 2 1 2 1

即模擬十次實驗，咱們獲得0次0觀測海洋，5次1觀測到海洋，5次2觀測海洋。 若是進行10000次模擬，每次模型投擲9次，則

dummy_w <- rbinom(10000, size = 9, prob = 0.7)
simplehist(dummy_w, xlab = "dummy water count")

經過使用似然函數進行模擬觀測，能夠看到咱們的實際觀測結果很類似。

2. 模型覈查

在上面的例子中，若是模擬觀測結果和咱們實際觀測結果差異很大，那麼就要注意了，可能在某些地方會存在錯誤。此外，還應該覈查一個模型是否充足。也就是咱們還要檢查那些沒可以被咱們模型所描述的數據，以便咱們可以對模型進行修正和改進。

沒有模型是十全十美的，總會有一些使人不滿意的地方，咱們所要判斷的就是這些地方是不是很重要的地方。沒有人但願本身建造的模型僅僅可以描述現有數據，而不能進行數據預測。模型預測出現的不完美在所不免。有不完美的地方不可怕，只要咱們知道哪兒不完美就能夠了。

再回到咱們的例子中，在上面模擬抽樣的過程當中，咱們用了抽樣機率0.7，咱們實際上使用對後驗機率進行了點估計，以此來進行抽樣的。可是，實際上，後驗機率是一個在0-1之間的分佈，而不是一個點。若是用點估計，那麼咱們就會丟失不少信息，形成「過分自信」。

實際上在抽樣過程當中，有兩個不肯定性的來源，一個是觀測的不肯定性，除非海洋覆蓋率爲0或者爲1，不然咱們沒法預測某一次投擲的結果。另外一個不肯定性來自參數不肯定性，也就是海洋覆蓋率自己，在0-1之間，咱們沒法肯定它真實的是哪個值。

若是咱們在一個參數不肯定的後驗機率分佈中進行抽樣，那麼咱們須要考慮參數的全部可能值，對於每個可能值，咱們獲得獲得一個抽樣分佈，最後把全部可能值得抽樣分佈進行平均，獲得後驗機率預測分佈。以下圖經過9個參數值的估計獲得的後驗機率預測分佈：

咱們獲得的後驗機率預測分佈是用來預測的，它考慮了參數不肯定性對預測結果的影響，因此預測會更準確一些。經過這個分佈圖，和上面咱們用0.7獲得的分佈圖比較能夠看出，這兒獲得的分佈圖更加矮胖，方差更大。而若是你要是用0.7估計獲得的高瘦的分佈圖來預測，就會出現「過分自信」的結果。

那麼咱們怎樣把參數的不肯定性考慮進去，來進行抽樣呢？很簡單，把上面的0.7換成咱們的後驗機率分佈就能夠了：

w <- rbinom(10000, size = 9, prob = samples)
simplehist(w, xlab = "number water count")

其中的「sample」是咱們在後驗機率分佈中抽樣的結果，表明了後驗機率的分佈狀況。經過上面的結果，能夠預測咱們9次抽樣最有可能看到6次結果爲「海洋」。和咱們實際觀測的異質性較高。

到這兒咱們的對模型的核查就完了嗎？尚未，還要考慮其餘可能影響模型的緣由。好比，萬一各個實驗之間不是相互獨立的怎麼辦？前一個實驗的結果會影響到後一次實驗。爲了衡量獨立性，咱們考慮兩個方面：一個是實驗中連續出現「海洋」或者「陸地」的最長長度；另外一個是相鄰實驗間，「陸地」和「海洋」之間的轉換次數。 在咱們的實際實驗中，咱們的觀測結果是依次是「W L W W W LW L W」（W-海洋，L-陸地），因此上述兩個指標在本實驗中分別是3和6。一樣咱們經過模擬，能夠獲得在相互獨立這一預期狀況下，這兩個參數的分佈狀況。以下如：

【藍色表示實際觀測值】

在相互獨立的假設下，咱們觀測的最長長度爲3和模擬的結果很吻合；而相互轉化數咱們觀測值爲6，而預期最可能的值爲4，這說明咱們的相鄰實驗結果之間可能互斥性。在這種狀況下，各次實驗之間存在的互斥性，可能使得每次實驗結果帶來的信息量大大減少，固然，咱們最終可能也能獲得比較準確地海洋覆蓋率的估計值，可是可能須要更多的試驗次數。

總結

本部份內容介紹了關於後驗機率分佈的一些操做。咱們的目的是從後驗機率分佈中進行抽樣，這樣咱們就避免了使用微積分等複雜的方式來描述後驗機率分佈了。這些在後驗機率中抽得的樣本能夠進一步用來區間估計、點估計、後驗機率預測覈查以及模擬。在點估計的時候，咱們還提到了一個重要概念，即損失函數。

後驗機率預測覈查綜合考慮了有後驗機率分佈致使的參數不肯定性和似然函數定義的抽樣自己的不肯定性，以此來覈實結果以及判斷模型是否充分。

隨着模型自己變得愈來愈複雜，後驗機率預測模擬的使用愈來愈普遍。

---