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

概述

最近，我們使用貝葉斯非參數(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-scale_Dirichlet過程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ī)模大時，需要維護大量未占用的成分）還是在計算負擔方面（每次迭代都需要更新大量不需要后驗推理的參數(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。

---

點擊標題查閱往期內(nèi)容

WINBUGS對隨機波動率模型進行貝葉斯估計與比較

左右滑動查看更多

01

02

03

04

下面的代碼從隨機測量中生成后驗樣本。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)的計算和存儲負擔。你可以把這種方法看作是對過程的截斷，只不過它可以被認為是*精確的截斷。事實上，在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 <- compilesamples <- 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ā)作）的風險。，每項研究都有治療和對照組。

模型的制定

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

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

、遵循共同的正態(tài)分布、和被合理地賦予非信息性先驗。參數(shù)量化了對照組和治療組之間的風險差異，而參數(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é)果表明，對照組和治療組之間存在著整體的風險差異。但是正態(tài)性假設(shè)呢？我們的結(jié)論對該假設(shè)是否穩(wěn)?。恳苍S隨機效應(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.

---