原文鏈接：http://tecdat.cn/?p=23785

原文出處：拓端數(shù)據(jù)部落公眾號

概述

最近，我們使用貝葉斯非參數(shù)（BNP）混合模型進(jìn)行馬爾科夫鏈蒙特卡洛（MCMC）推斷。

在這篇文章中，我們通過展示如何使用具有不同內(nèi)核的非參數(shù)混合模型進(jìn)行密度估計(jì)。在后面的文章中，我們將采用參數(shù)化的廣義線性混合模型，并展示如何切換到非參數(shù)化的隨機(jī)效應(yīng)表示，避免了正態(tài)分布的隨機(jī)效應(yīng)假設(shè)。

使用Dirichlet Process Mixture模型進(jìn)行基本密度估計(jì)

提供了通過Dirichlet過程混合（DPM）模型進(jìn)行非參數(shù)密度估計(jì)的機(jī)制（Ferguson, 1974; Lo, 1984; Escobar, 1994; Escobar and West, 1995）。對于一個獨(dú)立和相同分布的樣本?

，該模型的形式為

這個模型實(shí)現(xiàn)是靈活的，運(yùn)行任意核的混合。

, 可以是共軛的，也可以是不共軛的（也是任意的）基度量?

. 在共軛核/基數(shù)測量對的情況下，能夠檢測共軛的存在，并利用它來提高采樣器的性能。

為了說明這些能力，我們考慮對R中提供的Faithful火山數(shù)據(jù)集的噴發(fā)間隔時間的概率密度函數(shù)進(jìn)行估計(jì)。

data(faithful)

觀測值?

?對應(yīng)于數(shù)據(jù)框架的第二列，而?

.

使用CRP表示法擬合高斯location-scale?分布混合分布

模型說明

我們首先考慮用混合正態(tài)分布的location-scaleDirichlet過程s來擬合轉(zhuǎn)換后的數(shù)據(jù)

其中

?對應(yīng)的是正態(tài)-逆伽馬分布。這個模型可以解釋為提供一個貝葉斯版本的核密度估計(jì)?

用于使用高斯核和自適應(yīng)帶寬。在數(shù)據(jù)的原始尺度上，這可以轉(zhuǎn)化為一個自適應(yīng)的對數(shù)高斯核密度估計(jì)。

引入輔助變量

，表明混合的哪個成分產(chǎn)生了每個觀測值，并對隨機(jī)量

進(jìn)行積分，我們得到模型的CRP表示（Blackwell and MacQueen, 1973）。

其中

是向量

中唯一值的數(shù)量，

是第

個唯一值在

中出現(xiàn)的次數(shù)。這個說明清楚地表明，每個觀測值都屬于最多

正態(tài)分布聚類中的任何一個，并且CRP分布與分區(qū)結(jié)構(gòu)的先驗(yàn)分布相對應(yīng)。

這個模型的說明是這樣的

1. y[i] ~ dnorm(mu[i], var = s2[i])

2. mu[i] <- muTilde[xi[i]]

3. s2[i] <- s2Tilde[xi[i]]

4. xi[1:n] ~ dCRP(alpha, size = n)

5. muTilde[i] ~ dnorm(0, var = s2Tilde[i])

6. s2Tilde[i] ~ dinvgamma(2, 1)

7. alpha ~ dgamma(1, 1)

請注意，在模型代碼中，參數(shù)向量muTilde和s2Tilde的長度被設(shè)置為

.我們這樣做是因?yàn)槟壳暗膶?shí)現(xiàn)要求提前設(shè)置參數(shù)向量的長度，并且不允許它們的數(shù)量在迭代之間變化。因此，如果我們要確保算法總是按預(yù)期執(zhí)行，我們需要在最壞的情況下工作，即有多少個成分就有多少個觀測值的情況。但它的效率也有點(diǎn)低，無論是在內(nèi)存需求方面（當(dāng)?

規(guī)模大時，需要維護(hù)大量未占用的成分）還是在計(jì)算負(fù)擔(dān)方面（每次迭代都需要更新大量不需要后驗(yàn)推理的參數(shù)）。當(dāng)我們在下面使用伽馬分布的混合時，我們將展示一個能提高效率的計(jì)算捷徑。

還需要注意的是，

的值控制著我們先驗(yàn)預(yù)期的成分?jǐn)?shù)量，

的值越大，對應(yīng)于數(shù)據(jù)占據(jù)的成分?jǐn)?shù)量越多。因此，通過指定一個

先驗(yàn)值，我們?yōu)槟Ｐ驮黾恿遂`活性。對Gamma先驗(yàn)的特殊選擇允許使用數(shù)據(jù)增強(qiáng)方案從相應(yīng)的全條件分布中有效取樣。也可以選擇其他的先驗(yàn)，在這種情況下，這個參數(shù)

的默認(rèn)采樣是一個自適應(yīng)的隨機(jī)游走M(jìn)etropolis-Hastings算法。

運(yùn)行MCMC算法

下面的代碼設(shè)置了數(shù)據(jù)和常數(shù)，初始化了參數(shù)，定義了模型對象，并建立和運(yùn)行了MCMC算法。默認(rèn)采樣器是一個折疊的吉布斯采樣器（Neal, 2000）。

2. # 模型數(shù)據(jù)

3. y = standlFaithful

4. # 模型常量

5. n = length(standlFaithful))

6. # 參數(shù)初始化

7. list(xi = sample(1:10, size=n, replace=TRUE),

8. # 創(chuàng)建和編譯模型

9. Model(code, ?data, ?inits, ?consts)

10. ##定義模型...

11. ##建立模型...

12. ##設(shè)置數(shù)據(jù)和初始值...

13. ##在模型上運(yùn)行計(jì)算（隨后的任何錯誤報告可能只是反映了模型變量的缺失值）...

14. ##檢查模型的大小和尺寸......。

15. ##模型構(gòu)建完成。

16. ## 編譯完成。

17. #MCMC的配置、創(chuàng)建和編譯

18. MCMC(conf)

19. ## 編譯......這可能需要一分鐘

20. ## 編譯完成。

21. 


?

我們可以從參數(shù)的后驗(yàn)分布中提取樣本，并創(chuàng)建痕跡圖、直方圖和任何其他感興趣的總結(jié)。例如，對于參數(shù)

，我們有。
?

1. # 參數(shù)的痕跡圖

2. ts.plot(samples[ , "alpha"], xlab = "iteration", ylab = expression(alpha))

1. 


2. # 后驗(yàn)直方圖

3. hist(samples[ , "alpha"], xlab = expression(alpha), main = "", ylab = "Frequency")

4. 


在這個模型下，對于一個新的觀察

，后驗(yàn)預(yù)測分布是最佳密度估計(jì)（在平方誤差損失下）。這個估計(jì)的樣本可以很容易地從我們的MCMC產(chǎn)生的樣本中計(jì)算出來。

1. 


2. # 參數(shù)的后驗(yàn)樣本

3. samples[, "alpha"]

4. # 平均值的后驗(yàn)樣本

5. samples[, grep('muTilde', colnames(samples))] # 聚類平均數(shù)的后驗(yàn)樣本。

6. # 集群方差的后驗(yàn)樣本

7. samples[, grep('s2Tilde', colnames(samples))] # 聚類成員的后驗(yàn)樣本。

8. # 聚類成員關(guān)系的后驗(yàn)樣本

9. samples [, grep('xi', colnames(samples))] # 聚類成員的后驗(yàn)樣本。

10. 


11. hist(y, freq = FALSE,

12. xlab = "標(biāo)準(zhǔn)化對數(shù)尺度上的等待時間")

13. ##對標(biāo)準(zhǔn)化對數(shù)網(wǎng)格的密度進(jìn)行點(diǎn)式估計(jì)

14. 


然而，回顧一下，這是對等待時間的對數(shù)的密度估計(jì)。為了獲得原始尺度上的密度，我們需要對內(nèi)核進(jìn)行適當(dāng)?shù)霓D(zhuǎn)換。?

1. standlGrid*sd(lFaithful) + mean(lFaithful) # 對數(shù)尺度上的網(wǎng)格

2. 


3. hist(faithful$waiting, freq = FALSE

4. 


無論是哪種情況，都有明顯的證據(jù)表明，數(shù)據(jù)中的等待時間有兩個組成部分。

生成混合分布的樣本

雖然混合分布的線性函數(shù)的后驗(yàn)分布的樣本（比如上面的預(yù)測分布）可以直接從折疊采樣器的實(shí)現(xiàn)中計(jì)算出來，但是對于非線性函數(shù)

的推斷需要我們首先從混合分布中生成樣本。我們可以從隨機(jī)度量

中獲得后驗(yàn)樣本。需要注意的是，為了從

?，得到后驗(yàn)樣本，我們需要監(jiān)控所有參與其計(jì)算的隨機(jī)變量，即成員變量xi，聚類參數(shù)muTilde和s2Tilde，以及濃度參數(shù)alpha。

下面的代碼從隨機(jī)測量中生成后驗(yàn)樣本。cMCMC對象包括模型和參數(shù)的后驗(yàn)樣本。函數(shù)估計(jì)了一個截斷水平

，即truncG。后驗(yàn)樣本是一個帶

列的矩陣，其中參數(shù)分布

向量的維度（在本例中為

）。

outputG <- getSamplesDPmeasure(cmcmc)

下面的代碼使用隨機(jī)測量

的后驗(yàn)樣本來計(jì)算

的后驗(yàn)樣本。請注意，這些樣本是基于轉(zhuǎn)換后的模型計(jì)算的，大于70的值對應(yīng)于上述定義的網(wǎng)格上大于0.035的值。

1. 


2. truncG <- outputG$trunc # G的截斷水平

3. 


4. 


5. 


6. probY70 <- rep(0, nrow(samples)) ?# P(y.tilde>70)的后驗(yàn)樣本

7. 


8. hist(probY70 )

9. 


使用CRP表示法擬合伽馬混合分布

不限于在DPM模型中使用高斯核。就Old Faithful數(shù)據(jù)而言，除了我們在上一節(jié)中介紹的對數(shù)尺度上的高斯核的混合分布外，還有一種選擇是數(shù)據(jù)原始尺度上的伽馬混合分布。

模型

在這種情況下，模型的形式為

其中

對應(yīng)于兩個獨(dú)立Gamma分布的乘積。下面的代碼提供了該模型。

1. y[i] ~ dgamma(shape = beta[i], scale = lambda[i])

2. beta[i] <- betaTilde[xi[i]]

3. lambda[i] <- lambdaTilde[xi[i]]

請注意，在這種情況下，向量beta和lambda的長度為?

。這樣做是為了減少與采樣算法有關(guān)的計(jì)算和存儲負(fù)擔(dān)。你可以把這種方法看作是對過程的截斷，只不過它可以被認(rèn)為是*精確的截斷。事實(shí)上，在CRP表示法下，只要采樣器的成分?jǐn)?shù)嚴(yán)格低于采樣器每次迭代的參數(shù)向量的長度，使用長度短于樣本中觀察值的參數(shù)向量就會生成一個合適的算法。

運(yùn)行MCMC算法

下面的代碼設(shè)置了模型數(shù)據(jù)和常數(shù)，初始化了參數(shù)，定義了模型對象，并建立和運(yùn)行了Gamma混合分布的MCMC算法。請注意，在構(gòu)建MCMC時，會產(chǎn)生一個關(guān)于聚類參數(shù)數(shù)量的警告信息。這是因?yàn)閎etaTilde和lambdaTilde的長度小于

。另外，請注意，在執(zhí)行過程中沒有產(chǎn)生錯誤信息，這表明所需的集群數(shù)量未超過50個的上限。

1. data <- list(y = waiting)

2. Model(code, data = data)

cModel <- compilesamples <- runMCMC(cmcmc, niter = 7000, nburnin = 2000, setSeed = TRUE)

在這種情況下，我們使用參數(shù)的后驗(yàn)樣本來構(gòu)建一個軌跡圖，并估計(jì)

的后驗(yàn)分布。
?

1. # 參數(shù)的后驗(yàn)樣本的跟蹤圖

2. ts.plot(samples[ , 'alpha'], xlab = "iteration", ylab = expression(alpha))

1. 


2. # 參數(shù)的后驗(yàn)樣本的直方圖

3. hist(samples[, 'alpha'])

4. 


從混合分布中生成樣本

和以前一樣，我們從后驗(yàn)分布

中獲得樣本。
?

outputG <- getSamplesDPmeasure(cmcmc)

我們使用這些樣本來創(chuàng)建一個數(shù)據(jù)密度的估計(jì)值，以及一個95%置信帶。?

1. 


2. for(iter in seq_len)) {

3. density[iter, ] <- sapply(grid, function(x)

4. sum( weightSamples[iter, ] * dgamma)))

5. }

8. hist(waiting, freq = FALSE

我們再次看到，數(shù)據(jù)的密度是雙峰的，看起來與我們之前得到的數(shù)據(jù)非常相似。

使用stick-breaking?表示法擬合伽馬DP混合分布

模型

Dirichlet過程混合物的另一種表示方法是使用隨機(jī)分布

的stick-breaking表示（Sethuraman, 1994）。

引入輔助變量，

表明哪個成分產(chǎn)生了每個觀測值，上一節(jié)討論的Gamma密度的混合物的相應(yīng)模型的形式為

其中

?是兩個獨(dú)立Gamma分布的乘積。

與

. 下面的代碼提供了該模型說明。?

1. 


2. y[i] ~ dgamma(shape = beta[i], scale = lambda[i])

3. beta[i] <- betaStar[z[i]]

4. lambda[i] <- lambdaStar[z[i]]

5. z[i] ~ dcat(w[1:Trunc])

6. # stick-breaking

7. v[i] ~ dbeta(1, alpha)

8. w[1:Trunc] <- stick_breaking(v[1:(Trunc-1)]) # stick-breaking 權(quán)重

9. betaStar[i] ~ dgamma(shape = 71, scale = 2)

注意，截斷水平

已被設(shè)置為Trunc值，該值將在函數(shù)的常數(shù)參數(shù)中定義。

運(yùn)行MCMC算法

下面的代碼設(shè)置了模型數(shù)據(jù)和常量，初始化了參數(shù)，定義了模型對象，并建立和運(yùn)行了Gamma混合分布的MCMC算法。當(dāng)使用stick-breaking表示時，會指定一個分塊Gibbs抽樣器（Ishwaran, 2001; Ishwaran and James, 2002）。

1. data <- list(y = waiting)

2. consts <-length(waiting)

3. betaStar = rgamma

4. lambdaStar = rgamma

5. v = rbeta

6. z = sample

7. alpha = 1

8. 


compile(Model)

1. MCMC(rModel, c("w", "betaStar", "lambdaStar", 'z', 'alpha'))

2. comp(mcmc )

MCMC(cmcmc, niter = 24000)

使用stick-breaking近似法會自動提供隨機(jī)分布的近似值

，即?

。下面的代碼使用來自

樣本對象的后驗(yàn)樣本計(jì)算后驗(yàn)樣本，并從中計(jì)算出數(shù)據(jù)的密度估計(jì)。

1. 


2. densitySamples[i, ] <- sapply(grid, function(x) sum(weightSamples ?* dgamma(x, shape ,

3. scale )))

4. 


5. hist( waiting ylim=c(0,0.05),

6. 


正如預(yù)期的那樣，這個估計(jì)值看起來與我們通過CRP表示的過程獲得的估計(jì)值相同。

貝葉斯非參數(shù)化：非參數(shù)化隨機(jī)效應(yīng)

我們將采用一個參數(shù)化的廣義線性混合模型，并展示如何切換到非參數(shù)化的隨機(jī)效應(yīng)表示，避免了正態(tài)分布的隨機(jī)效應(yīng)假設(shè)。

心肌梗死（MIs）的參數(shù)化meta分析

我們將在對以前非常流行的糖尿病藥物 "Avandia "的副作用進(jìn)行meta分析的背景下，說明使用非參數(shù)混合模型對隨機(jī)效應(yīng)分布進(jìn)行建模。我們分析的數(shù)據(jù)在引起對這種藥物的安全性的嚴(yán)重質(zhì)疑方面發(fā)揮了作用。問題是使用"Avandia "是否會增加心肌梗死（心臟病發(fā)作）的風(fēng)險。，每項(xiàng)研究都有治療和對照組。

模型的制定

我們首先進(jìn)行基于標(biāo)準(zhǔn)的廣義線性混合模型（GLMM）的meta分析。向量n和x分別包含對照組的患者總數(shù)和每項(xiàng)研究中對照組的心肌梗死患者人數(shù)。同樣，向量m和y包含接受藥物的病人的類似信息。該模型的形式為

其中，隨機(jī)效應(yīng)

、遵循共同的正態(tài)分布

、

和

被合理地賦予非信息性先驗(yàn)。參數(shù)

量化了對照組和治療組之間的風(fēng)險差異，而參數(shù)

則量化了研究的具體變化。

這個模型可以用以下代碼指定。

1. 


2. y[i] ~ dbin(size = m[i], prob = q[i]) # 藥物MIs

3. x[i] ~ dbin(size = n[i], prob = p[i]) # 控制MIs

4. q[i] <- expit(theta + gamma[i]) # 藥物的對數(shù)指數(shù)

5. p[i] <- expit(gamma[i]) #對照組對數(shù)

6. gamma[i] ~ dnorm(mu, var = tau2) # 研究效果

7. theta ~ dflat() # 藥物的影響

8. # 隨機(jī)效應(yīng)超參數(shù)

9. mu ~ dnorm(0, 10)

10. tau2 ~ dinvgamma(2, 1)

11. 


運(yùn)行MCMC

讓我們來運(yùn)行一個基本的MCMC。

1. MCMC(codeParam, ?data, inits,

2. constants, monitors = c("mu", "tau2", "theta", "gamma")

1. par(mfrow = c(1, 4)

2. 
   

3. hist(gammaMn)

4. hist(samples[1000, gammaCols)

結(jié)果表明，對照組和治療組之間存在著整體的風(fēng)險差異。但是正態(tài)性假設(shè)呢？我們的結(jié)論對該假設(shè)是否穩(wěn)健？也許隨機(jī)效應(yīng)的分布是偏斜的。

用于meta分析的基于DP的隨機(jī)效應(yīng)模型

模型

現(xiàn)在，我們對

使用非參數(shù)分布。更具體地說，我們假設(shè)每個

都是由位置尺度的正態(tài)混合分布產(chǎn)生的。

這種模型引起了隨機(jī)效應(yīng)之間的聚類。與密度估計(jì)問題的情況一樣，DP先驗(yàn)允許數(shù)據(jù)決定分量的數(shù)量，從最少的一個分量（即簡化為參數(shù)模型）到最多的分量，即每個觀測值有一個分量。如果數(shù)據(jù)支持這種行為，這允許隨機(jī)效應(yīng)的分布是多模態(tài)的，大大增加了其靈活性。這個模型可以用以下代碼指定。

1. 


2. y[i] ~ dbin(size = m[i], prob = q[i]) # 藥物MIs

3. x[i] ~ dbin(size = n[i], prob = p[i]) # MIs

4. q[i] <- expit(theta + gamma[i]) # 藥物的對數(shù)指數(shù)

5. p[i] <- expit(gamma[i]) # 對照組對數(shù)值

6. gamma[i] ~ dnorm(mu[i], var = tau2[i]) # 來自混合物的隨機(jī)效應(yīng)。

7. mu[i]<- muTilde[xi[i]] ? ? ? ? ? ? ? ? # 來自聚類的隨機(jī)效應(yīng)的平均值 xi[i]

8. tau2[i] <- tau2Tilde[xi[i]] ? ? ? ? ? # 來自群組xi[i]的隨機(jī)效應(yīng)變量

9. # 從基礎(chǔ)測量中提取混合成分參數(shù)

10. muTilde[i] ~ dnorm(mu0, var = var0)

11. tau2Tilde[i] ~ dinvgamma(a0, b0)

12. # 用于將研究報告聚類為混合成分的CRP

13. xi[1:nStudies] ~ dCRP(alpha, size = nStudies)

14. # 超參數(shù)

15. 


16. theta ~ dflat() # 藥物的影響

17. 


運(yùn)行MCMC

以下代碼對模型進(jìn)行了編譯，并對模型運(yùn)行了一個壓縮Gibbs抽樣

1. inits <- list(gamma = rnorm(nStudies))

2. 


3. MCMC(code = BNP, data = data)

1. 
   

2. hist(samplesBNP[, 'theta'], xlab = expression(theta), main = 'avandia的影響')

3. main = "隨機(jī)效應(yīng)分布")

4. main = "隨機(jī)效應(yīng)分布")

5. 
   

6. # 推斷出了多少個混合成分？

7. xiRes <- samplesBNP[, xiCols].

8. 
   

主要推論似乎對原始的參數(shù)化假設(shè)很穩(wěn)健。這可能是由于沒有太多證據(jù)表明隨機(jī)效應(yīng)分布中缺乏正態(tài)性。

參考文獻(xiàn)

Blackwell, D. and MacQueen, J. 1973. Ferguson distributions via Polya urn schemes. The Annals of Statistics 1:353-355.

Ferguson, T.S. 1974. Prior distribution on the spaces of probability measures. Annals of Statistics 2:615-629.

Lo, A.Y. 1984. On a class of Bayesian nonparametric estimates I: Density estimates. The Annals of Statistics 12:351-357.

最受歡迎的見解

1.matlab使用貝葉斯優(yōu)化的深度學(xué)習(xí)

2.matlab貝葉斯隱馬爾可夫hmm模型實(shí)現(xiàn)

3.R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真

4.R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸

5.R語言中的Stan概率編程MCMC采樣的貝葉斯模型

6.Python用PyMC3實(shí)現(xiàn)貝葉斯線性回歸模型

7.R語言使用貝葉斯 層次模型進(jìn)行空間數(shù)據(jù)分析

8.R語言隨機(jī)搜索變量選擇SSVS估計(jì)貝葉斯向量自回歸（BVAR）模型

9.matlab貝葉斯隱馬爾可夫hmm模型實(shí)現(xiàn)