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

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

最近我們被客戶要求撰寫關(guān)于貝葉斯非參數(shù)模型的研究報告，包括一些圖形和統(tǒng)計輸出。

概述

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

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

使用Dirichlet Process Mixture模型進行基本密度估計

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

，該模型的形式為

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

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

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

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

data(faithful)

觀測值?

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

.

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

模型說明

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

其中

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

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

引入輔助變量

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

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

其中

是向量

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

是第

個唯一值在

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

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

這個模型的說明是這樣的

? ?y[i] ~ dnorm(mu[i], var = s2[i]) ? ?mu[i] <- muTilde[xi[i]] ? ?s2[i] <- s2Tilde[xi[i]] ?xi[1:n] ~ dCRP(alpha, size = n) ? ?muTilde[i] ~ dnorm(0, var = s2Tilde[i]) ? ?s2Tilde[i] ~ dinvgamma(2, 1) ?alpha ~ dgamma(1, 1)

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

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

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

還需要注意的是，

的值控制著我們先驗預(yù)期的成分數(shù)量，

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

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

的默認采樣是一個自適應(yīng)的隨機游走Metropolis-Hastings算法。

運行MCMC算法

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

# 模型數(shù)據(jù)y = standlFaithful# 模型常量n = length(standlFaithful))# 參數(shù)初始化list(xi = sample(1:10, size=n, replace=TRUE),# 創(chuàng)建和編譯模型Model(code, ?data, ?inits, ?consts)##定義模型...##建立模型...##設(shè)置數(shù)據(jù)和初始值...##在模型上運行計算（隨后的任何錯誤報告可能只是反映了模型變量的缺失值）... ##檢查模型的大小和尺寸......。##模型構(gòu)建完成。## 編譯完成。#MCMC的配置、創(chuàng)建和編譯MCMC(conf)## 編譯......這可能需要一分鐘## 編譯完成。

?

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

，我們有。

# 參數(shù)的痕跡圖ts.plot(samples[ , "alpha"], xlab = "iteration", ylab = expression(alpha))

# 后驗直方圖hist(samples[ , "alpha"], xlab = expression(alpha), main = "", ylab = "Frequency")

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

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

# 參數(shù)的后驗樣本 samples[, "alpha"]# 平均值的后驗樣本 samples[, grep('muTilde', colnames(samples))] # 聚類平均數(shù)的后驗樣本。# 集群方差的后驗樣本samples[, grep('s2Tilde', colnames(samples))] # 聚類成員的后驗樣本。# 聚類成員關(guān)系的后驗樣本samples [, grep('xi', colnames(samples))] # 聚類成員的后驗樣本。hist(y, freq = FALSE, ? ? xlab = "標準化對數(shù)尺度上的等待時間")##對標準化對數(shù)網(wǎng)格的密度進行點式估計

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

standlGrid*sd(lFaithful) + mean(lFaithful) # 對數(shù)尺度上的網(wǎng)格hist(faithful$waiting, freq = FALSE ?

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

生成混合分布的樣本

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

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

中獲得后驗樣本。需要注意的是，為了從

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

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

，即truncG。后驗樣本是一個帶

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

向量的維度（在本例中為

）。

outputG <- getSamplesDPmeasure(cmcmc)

下面的代碼使用隨機測量

的后驗樣本來計算

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

? ? truncG <- outputG$trunc # G的截斷水平probY70 <- rep(0, nrow(samples)) ?# P(y.tilde>70)的后驗樣本hist(probY70 )

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

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

模型

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

其中

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

? ?y[i] ~ dgamma(shape = beta[i], scale = lambda[i]) ? ?beta[i] <- betaTilde[xi[i]] ? ?lambda[i] <- lambdaTilde[xi[i]]

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

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

運行MCMC算法

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

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

data <- list(y = waiting)Model(code, data = data)

cModel <- compile

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

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

的后驗分布。

# 參數(shù)的后驗樣本的跟蹤圖ts.plot(samples[ , 'alpha'], xlab = "iteration", ylab = expression(alpha))

# 參數(shù)的后驗樣本的直方圖 hist(samples[, 'alpha'])

從混合分布中生成樣本

和以前一樣，我們從后驗分布

中獲得樣本。

outputG <- getSamplesDPmeasure(cmcmc)

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

for(iter in seq_len)) { ?density[iter, ] <- sapply(grid, function(x) ? ?sum( weightSamples[iter, ] * dgamma)))}hist(waiting, freq = FALSE

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

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

模型

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

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

引入輔助變量，

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

其中

?是兩個獨立Gamma分布的乘積。

與

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

? ? ?y[i] ~ dgamma(shape = beta[i], scale = lambda[i]) ? ? ?beta[i] <- betaStar[z[i]] ? ? ?lambda[i] <- lambdaStar[z[i]] ? ? ?z[i] ~ dcat(w[1:Trunc]) ? ? # stick-breaking ? ? ?v[i] ~ dbeta(1, alpha) ? ?w[1:Trunc] <- stick_breaking(v[1:(Trunc-1)]) # stick-breaking 權(quán)重 ? ? ?betaStar[i] ~ dgamma(shape = 71, scale = 2)

注意，截斷水平

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

運行MCMC算法

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

data <- list(y = waiting)consts <-length(waiting)betaStar = rgamma ? ? ? ? ? ? ?lambdaStar = rgamma ? ? ? ? ? ? ?v = rbeta ? ? ? ? ? ? ?z = sample ? ? ? ? ? ? ?alpha = 1

compile(Model)

MCMC(rModel, c("w", "betaStar", "lambdaStar", 'z', 'alpha'))comp(mcmc )

MCMC(cmcmc, niter = 24000)

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

，即?

。下面的代碼使用來自

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

?densitySamples[i, ] <- sapply(grid, function(x) sum(weightSamples ?* dgamma(x, shape , ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?scale )))hist( waiting ylim=c(0,0.05),

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

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

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

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

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

模型的制定

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

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

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

、

和

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

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

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

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

? ? ? ?y[i] ~ dbin(size = m[i], prob = q[i]) # 藥物MIs ? ? ? ?x[i] ~ dbin(size = n[i], prob = p[i]) # 控制MIs ? ? ? ?q[i] <- expit(theta + gamma[i]) # 藥物的對數(shù)指數(shù) ? ? ? ?p[i] <- expit(gamma[i]) #對照組對數(shù) ? ? ? ?gamma[i] ~ dnorm(mu, var = tau2) # 研究效果 ? ?theta ~ dflat() # 藥物的影響 ? ?# 隨機效應(yīng)超參數(shù) ? ?mu ~ dnorm(0, 10) ? ?tau2 ~ dinvgamma(2, 1)

運行MCMC

讓我們來運行一個基本的MCMC。

MCMC(codeParam, ?data, inits, ? ? ? ? ? ? ? ? ? ? ?constants, monitors = c("mu", "tau2", "theta", "gamma")

par(mfrow = c(1, 4)hist(gammaMn)hist(samples[1000, gammaCols)

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

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

模型

現(xiàn)在，我們對

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

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

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

? ? ? ?y[i] ~ dbin(size = m[i], prob = q[i]) # 藥物MIs ? ? ? ?x[i] ~ dbin(size = n[i], prob = p[i]) # MIs ? ? ? ?q[i] <- expit(theta + gamma[i]) # 藥物的對數(shù)指數(shù) ? ? ? ?p[i] <- expit(gamma[i]) # 對照組對數(shù)值 ? ? ? ?gamma[i] ~ dnorm(mu[i], var = tau2[i]) # 來自混合物的隨機效應(yīng)。 ? ? ? ?mu[i]<- muTilde[xi[i]] ? ? ? ? ? ? ? ? # 來自聚類的隨機效應(yīng)的平均值 xi[i] ? ? ? ?tau2[i] <- tau2Tilde[xi[i]] ? ? ? ? ? # 來自群組xi[i]的隨機效應(yīng)變量 ? ?# 從基礎(chǔ)測量中提取混合成分參數(shù) ? ? ? ?muTilde[i] ~ dnorm(mu0, var = var0) ? ? ? ?tau2Tilde[i] ~ dinvgamma(a0, b0) ? ?# 用于將研究報告聚類為混合成分的CRP ? ?xi[1:nStudies] ~ dCRP(alpha, size = nStudies) ? ?# 超參數(shù) ? ? ?theta ~ dflat() # 藥物的影響

運行MCMC

以下代碼對模型進行了編譯，并對模型運行了一個壓縮Gibbs抽樣

inits <- list(gamma = rnorm(nStudies))MCMC(code = BNP, data = data)

hist(samplesBNP[, 'theta'], xlab = expression(theta), main = 'avandia的影響') ? ? ? ? ? ? ?main = "隨機效應(yīng)分布") ? ? ? ? ? ? ? ? ? main = "隨機效應(yīng)分布")# 推斷出了多少個混合成分？xiRes <- samplesBNP[, xiCols].

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

參考文獻

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模型實現(xiàn)

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

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

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

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

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

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

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