不想去健身房的我，最后被貝葉斯分析說服了...

作者：大數(shù)據(jù)文摘 2018-07-27 13:55:12

作者從小被人說體型瘦小，于是他用了貝葉斯分析最終得出了自己的體重確實(shí)低于本國(guó)平均水平的結(jié)論。這篇文章將敘述一個(gè)在線性回歸理論中應(yīng)用貝葉斯分析的有趣試驗(yàn)。

大數(shù)據(jù)文摘出品

編譯：Zhifu、JonyKai、湯圓、夏雅薇

可能經(jīng)常你會(huì)聽到一些很主觀的評(píng)價(jià)比如“你太瘦了”或者“你怎么那么高”，但這里瘦或者高都是基于評(píng)價(jià)者的主觀判斷和視覺記憶做出的評(píng)述，并沒有嚴(yán)格的參照。

作者從小被人說體型瘦小，于是他用了貝葉斯分析最終得出了自己的體重確實(shí)低于本國(guó)平均水平的結(jié)論。沒有比直白的數(shù)字更有說服力了，想要說服自己健身的小伙伴不妨也試試作者的統(tǒng)計(jì)學(xué)暴擊法!

這篇文章將敘述一個(gè)在線性回歸理論中應(yīng)用貝葉斯分析的有趣試驗(yàn)(一個(gè)小秘密：我在這篇文章中使用了公制計(jì)量單位呦)。

如文章標(biāo)題所述，我將會(huì)對(duì)自己的體格進(jìn)行一番科學(xué)的研究。

在開始之前，我希望大家多了解我一些：我在越南出生，在新加坡上高中，現(xiàn)在在美國(guó)上大學(xué)。我經(jīng)常因?yàn)樯硇问菪”蝗藗內(nèi)⌒Γf我應(yīng)該多去健身房鍛煉增肌，擁有更強(qiáng)壯的體魄。這些評(píng)價(jià)我一般都一笑置之。對(duì)于一個(gè)身高169厘米(5尺6寸)，體重58公斤(127磅)的人來說，我的BMI指數(shù)(20.3)幾近***。

[[238155]]

仔細(xì)一想，大家可能沒說錯(cuò)：我比一般的越南男性要高，但卻只有平均體重(維基百科里越南男性的平均體重是58公斤,平均身高是162厘米)，“看起來”可能是要稍微瘦一些。

這里“看起來”是關(guān)鍵：背后的邏輯很清楚，不是嗎?如果你把自己抻長(zhǎng)一些，體重不變，那確實(shí)應(yīng)該看起來苗條一些。我把這個(gè)看作是嚴(yán)肅的科學(xué)問題，并準(zhǔn)備深入研究。

對(duì)于一個(gè)169厘米高的越南男性來說，我到底輕了多少呢?

我們需要一種有理有據(jù)的方式來研究這個(gè)問題。有個(gè)好方法是盡可能多地找到越南男子身高和體重的數(shù)據(jù)，來判斷我在這個(gè)樣本中的位置。

越南人口數(shù)據(jù)

在瀏覽各種網(wǎng)頁(yè)后，我找到了一份調(diào)查研究數(shù)據(jù)，包含超過10,000名越南人的人口統(tǒng)計(jì)信息。將抽樣條件設(shè)置為年齡18-29歲的越南男性，從而得到數(shù)量為383的樣本，這個(gè)樣本足以用來進(jìn)行接下來的分析啦。

首先，通過人口體重的直方圖，看看我在年輕越南男性中的體重分布位置。

紅線表示樣本的中位數(shù)，而橙色線表示平均值

這張直方圖表明我的體重略低于383名年輕越南男子的體重平均值和中位數(shù)。看起來是與我們要研究的相關(guān)呀!然而問題并不在于我的體重與這個(gè)樣本本身的比較，是假設(shè)這383個(gè)人可以代表越南男性，在身高169厘米的情況下，我的體重與整個(gè)越南人口相比處于一個(gè)什么位置。為此，我們需要進(jìn)行回歸分析。

首先繪制一個(gè)身高和體重的二維散點(diǎn)圖

好吧，看起來我處在平均水平。但是如果我們只看身高169厘米的數(shù)據(jù)(想象一條垂直x軸于169厘米這個(gè)刻度并穿過紅點(diǎn)的直線)，我的體重在他們之中處于下游。

另一個(gè)重要的發(fā)現(xiàn)是越南男性身高和體重呈正相關(guān)。我們將進(jìn)行定量分析來進(jìn)一步了解這種關(guān)系。

首先，讓我們快速添加“普通最小二乘”線。我稍后會(huì)詳細(xì)介紹這一點(diǎn)，現(xiàn)在先在圖上展示出來。

最小二乘線可以表示為y = -86.32 + 0.889x，這表明通常情況下，我這個(gè)年齡的越南男性，每增加1厘米的身高，體重會(huì)增加0.88千克。

但是，這并沒有解決我們的問題;身高169厘米，體重58公斤到底是太沉，太輕還是剛剛好呢?要以定量的方式進(jìn)一步解釋這個(gè)問題，如果有所有身高1米68的人的體重分布，那么我的體重排在前25%，50%或75%的幾率是多少?要弄清這一點(diǎn)，我們需要深入學(xué)習(xí)并理解回歸背后的理論。

線性回歸理論

在線性回歸模型中，Y變量(在我們的例子中，是人的體重)是x(身高)的線性函數(shù)。在這個(gè)線性關(guān)系中截距和斜率分別為β0和β1;也就是說，假設(shè)E(Y | X = x)=β0+β1x。我們不知道β0和β1是多少，所以將它們視為未知參數(shù)。

在大多數(shù)標(biāo)準(zhǔn)線性回歸模型中，我們進(jìn)一步假設(shè)給定X = x的情況下，Y的條件分布是正態(tài)分布的。

這就是基本的線性回歸模型：

可以被改寫成：

注意，在許多模型中，我們可以用精度參數(shù)τ替換方差參數(shù)σ，其中τ= 1 /σ。

總結(jié)：因變量Y遵循由平均數(shù)μi和精度參數(shù)τ決定的正態(tài)分布。μi是由β0和β1決定的X的線性函數(shù)。

***，我們還需假設(shè)未知方差不依賴于x;這種假設(shè)稱為同方差性。

涉及的內(nèi)容很多，都涵括在下面這張圖里啦。

圖像來源：Joseph Chang(耶魯大學(xué))

在實(shí)際的數(shù)據(jù)分析問題中，我們掌握的只是黑點(diǎn) - 數(shù)據(jù)。使用這些數(shù)據(jù)，我們的目標(biāo)是推斷不知道的事情，包括β0，β1(在圖片中的藍(lán)色虛線)和σ(它決定了在給定一個(gè)y值的時(shí)候，紅色正態(tài)分布密度的寬度)。注意，每個(gè)黑點(diǎn)周圍的正態(tài)分布看起來完全相同。這是同方差性的本質(zhì)。

估算參數(shù)

現(xiàn)在，有幾種方法可以估算β0和β1。如果你使用普通最小二乘估計(jì)這樣的模型，你不必?fù)?dān)心概率公式，因?yàn)槟阏趯ふ?beta;0和β1的***值，而這是通過最小化擬合值與預(yù)測(cè)值的平方誤差得到的。

另一方面，你可以使用***似然估計(jì)(MLE)來估計(jì)此類模型：通過***化似然函數(shù)來尋找參數(shù)的***值。

注意：一個(gè)有趣的結(jié)果(我沒有放上具體的數(shù)學(xué)證明)是，如果我們假設(shè)誤差項(xiàng)也屬于正態(tài)分布的話，那么最小二乘估計(jì)的結(jié)果也是***似然估計(jì)的結(jié)果。

貝葉斯方法的線性回歸

不同于***化似然函數(shù)，貝葉斯方法會(huì)假設(shè)參數(shù)服從一個(gè)先驗(yàn)分布。根據(jù)貝葉斯公式計(jì)算出參數(shù)后驗(yàn)概率：

貝葉斯方法的似然函數(shù)同上面的類似，不同之處在于加入了對(duì)估計(jì)參數(shù)β0，β1，τ的先驗(yàn)分布：

等等，什么是先驗(yàn)分布，為什么這會(huì)讓公式看上去更加復(fù)雜?

請(qǐng)相信我，先驗(yàn)分布雖然看上去很奇怪，但實(shí)際上很直觀。事實(shí)上，我們對(duì)未知的參數(shù)(例子中的β0，β1，τ)分配一個(gè)看上去很隨意的先驗(yàn)分布，這里存在著很強(qiáng)的哲學(xué)緣由。

先驗(yàn)分布能夠反映出在沒看見數(shù)據(jù)之前我們對(duì)數(shù)據(jù)的假設(shè)理解。在觀察過一些數(shù)據(jù)之后，我們應(yīng)用貝葉斯公式，就得到了同時(shí)考慮到了先驗(yàn)和數(shù)據(jù)的未知參數(shù)后驗(yàn)分布。根據(jù)后驗(yàn)分布，我們就能預(yù)測(cè)出未來的數(shù)據(jù)的分布。

最終的參數(shù)估計(jì)雖然取決于數(shù)據(jù)和先驗(yàn)分布，但是如果數(shù)據(jù)中包含的信息越多，那先驗(yàn)的影響就越小。

那么我該如何選擇先驗(yàn)分布

這是個(gè)好問題，因?yàn)檫@里存在著無(wú)數(shù)種可能。理論上只有存在一個(gè)正確的先驗(yàn)?zāi)軌蚰軌蚍磻?yīng)出你的先驗(yàn)假設(shè)。但是在實(shí)際中，先驗(yàn)分布的選取是相當(dāng)?shù)闹饔^，甚至有時(shí)可以是任意的。

我們可以選擇一個(gè)有著較大標(biāo)準(zhǔn)方差(意味著精度很低)的正態(tài)分布先驗(yàn)。比如，我們假設(shè)參數(shù)β0和β1是服從均值為0標(biāo)準(zhǔn)方差為10000的正態(tài)分布。這種分布是毫無(wú)信息的分布，因?yàn)榉植际制教?這意味著，參數(shù)在任意區(qū)間的取值概率幾乎相同)。

如果選取了這種類型的先驗(yàn)分布，那么我們就不用考慮在這類分布中哪種分布更好，因?yàn)榉植紟缀醵己芷教梗诿總€(gè)地方的概率都可以忽略不計(jì)。此外，后驗(yàn)分布不會(huì)受這種分布的影響。

同樣，對(duì)于精度τ，因?yàn)槠浔仨毷欠秦?fù)的，所以需要選取一個(gè)取值限定在非負(fù)范圍的分布。比如，在這里可以選取一個(gè)帶有較小形狀和尺度參數(shù)的伽馬分布。

另外一種很有用的不附帶信息的分布是均勻分布。如果對(duì) σ 或 τ選擇了均勻分布，那么你最終將會(huì)得到這樣的模型。如下圖所示，由John K. Kruschke繪制。

John K. Kruschke繪制的模型

用R語(yǔ)言和JAGS模擬數(shù)據(jù)

到目前為止，我們?nèi)灾煌Ａ粼诶碚撾A段。大多數(shù)情況下。后驗(yàn)分布并不能直接得到(想想正態(tài)分布和伽馬分布有多復(fù)雜，然后還要再將他們乘起來)。Markov Chain Monte Carlo方法常常用來估計(jì)模型的參數(shù)。利用JAGS就能幫我們完成。

“等一下!!!這個(gè)工具真能夠幫我們解決這些復(fù)雜的公式么?”

JAGS模型是一個(gè)基于Markov Chain Monte Carlo(MCMC)的仿真過程，它能夠生成出參數(shù)空間θ=(β0;β1;τ)的許多次迭代。由參數(shù)空間中每個(gè)參數(shù)生成的樣本分布會(huì)估計(jì)出參數(shù)最有可能的分布

為什們是這樣呢?這個(gè)解釋起來十分復(fù)雜，已經(jīng)超出了本篇介紹的范圍。直接來說就是：MCMC通過構(gòu)建一個(gè)馬爾可夫鏈產(chǎn)生了服從后驗(yàn)分布的樣本，這個(gè)馬爾可夫鏈有著同樣的目標(biāo)后驗(yàn)分布!?

這個(gè)過程很沒意思。最快的方法就是：不去解等式(2)，因?yàn)橥ǔ２豢赡艿玫浇馕鼋狻Ｎ覀兡茏龅木褪锹斆鞯牟蓸樱跀?shù)學(xué)上證明了這些樣本確實(shí)服從以β0，β1，τ為參數(shù)的分布。

那么揮別了數(shù)學(xué)之后，我該怎么使用JAGS?

我們按下面的步驟在R語(yǔ)言中運(yùn)行JAGS

首先以文本的形式寫下模型

然后，我們讓JAGS執(zhí)行仿真模擬。這里我使用JAGs對(duì)參數(shù)空間θ進(jìn)行10000次模擬

抽樣之后，我們就得到了θ=(β0;β1;τ)的抽樣數(shù)據(jù)，如下表所示：