拓端tecdat|R語言貝葉斯Metropolis-Hastings Gibbs 吉布斯采樣器估計變點(diǎn)指數(shù)分布分析
原文鏈接：http://tecdat.cn/?p=26578?
原文出處：拓端數(shù)據(jù)部落公眾號
最近我們被客戶要求撰寫關(guān)于吉布斯采樣器的研究報告，包括一些圖形和統(tǒng)計輸出。 指數(shù)分布是泊松過程中事件之間時間的概率分布，因此它用于預(yù)測到下一個事件的等待時間，例如，您需要在公共汽車站等待的時間，直到下一班車到了。
在本文中，我們將使用指數(shù)分布，假設(shè)它的參數(shù) λ ，即事件之間的平均時間，在某個時間點(diǎn) k 發(fā)生了變化，即：
?
我們的主要目標(biāo)是使用 Gibbs 采樣器在給定來自該分布的 n 個觀測樣本的情況下估計參數(shù) λ、α 和 k。
吉布斯Gibbs 采樣器
Gibbs 采樣器是 Metropolis-Hastings 采樣器的一個特例，通常在目標(biāo)是多元分布時使用。使用這種方法，鏈?zhǔn)峭ㄟ^從目標(biāo)分布的邊緣分布中采樣生成的，因此每個候選點(diǎn)都被接受。
Gibbs 采樣器生成馬爾可夫鏈如下：
讓?
?是 Rd 中的隨機(jī)向量，在時間 t=0 初始化 X(0)。
對于每次迭代 t=1,2,3,...重復(fù)：
設(shè)置 x1=X1(t-1)。
對于每個 j=1,...,d：
生成 X?j(t) 從?
， 其中?
?是給定 X(-j) 的 Xj的單變量條件密度。
更新?
.?
當(dāng)每個候選點(diǎn)都被接受時，設(shè)置?
.?
增加 t。
貝葉斯公式
變點(diǎn)問題的一個簡單公式假設(shè) f和 g?已知密度：
?
?其中 k 未知且 k=1,2,...,n。讓 Yi為公交車到達(dá)公交車站之間經(jīng)過的時間（以分鐘為單位）。假設(shè)變化點(diǎn)發(fā)生在第 k分鐘，即：
?
?當(dāng) Y=(Y1,Y2,...,Yn) 時，似然 L(Y|k)由下式給出：
?
假設(shè)具有獨(dú)立先驗的貝葉斯模型由下式給出：
數(shù)據(jù)和參數(shù)的聯(lián)合分布為：
其中，
正如我之前提到的，Gibbs 采樣器的實(shí)現(xiàn)需要從目標(biāo)分布的邊緣分布中采樣，因此我們需要找到 λ、α?和 k 的完整條件分布。你怎么能這樣做？簡單來說，您必須從上面介紹的連接分布中選擇僅依賴于感興趣參數(shù)的項并忽略其余項。
λ?的完整條件分布由下式給出：
?
α?的完整條件分布由下式給出：
k?的完整條件分布由下式給出：
計算方法
在這里，您將學(xué)習(xí)如何使用使用 R 的 Gibbs 采樣器來估計參數(shù) λ、α 和 k。
數(shù)據(jù)
首先，我們從具有變化點(diǎn)的下一個指數(shù)分布生成數(shù)據(jù)：
set.seed(98712)y <- c(rexp(25, rate = 2), rexp(35, rate = 10))
考慮到公交車站的情況，一開始公交車平均每2分鐘一班，但從時間i=26開始，公交車開始平均每10分鐘一班到公交車站。
Gibbs采樣器的實(shí)現(xiàn)
首先，我們需要初始化 k、λ 和 α。
n <- length(y) # 樣本的觀察值的數(shù)量lci <- 10000 # 鏈的大小aba <- alpha <- k <- numeric(lcan)k[1] <- sample(1:n,
現(xiàn)在，對于算法的每次迭代，我們需要生成 λ(t)、α(t) 和 k(t)，如下所示（記住如果 k+1>n 沒有變化點(diǎn)）：
for (i in 2:lcan){ ? ? ? ?kt <- k[i-1] ? ? ? ?# 生成lambda ? ? ? ?lambda[i] <- rgamma ? ? ? ?# 生成α ? ? ? ? ? ? ?# 產(chǎn)生k ? ? ? ? ?for (j in 1:n) { ? ? ? ? ? ? ? ?L[j] <- ((lambda[i] / alpha[i# 刪除鏈條上的前9000個值bunIn <- 9000
結(jié)果
在本節(jié)中，我們將介紹 Gibbs 采樣器生成的鏈及其參數(shù) λ、α 和 k 的分布。參數(shù)的真實(shí)值用紅線表示。
下表顯示了參數(shù)的實(shí)際值和使用 Gibbs 采樣器獲得的估計值的平均值：
res <- c(mean(k[-(1:bun)]), mean(lmba[-(1:burn)]), mean(apa[-(1:buI)]))resfil
結(jié)論
從結(jié)果中，我們可以得出結(jié)論，使用 R 中的 Gibbs 采樣器獲得的具有變點(diǎn)的指數(shù)分布對參數(shù) k、λ 和 α 的估計值的平均值接近于參數(shù)的實(shí)際值，但是我們期望更好估計。這可能是由于選擇了鏈的初始值或選擇了 λ?和 α的先驗分布。
最受歡迎的見解
1.使用R語言進(jìn)行METROPLIS-IN-GIBBS采樣和MCMC運(yùn)行
2.R語言中的Stan概率編程MCMC采樣的貝葉斯模型
3.R語言實(shí)現(xiàn)MCMC中的Metropolis–Hastings算法與吉布斯采樣
4.R語言BUGS JAGS貝葉斯分析 馬爾科夫鏈蒙特卡洛方法（MCMC）采樣
5.R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸
6.R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真分析
7.R語言用Rcpp加速M(fèi)etropolis-Hastings抽樣估計貝葉斯邏輯回歸模型的參數(shù)
8.R語言使用Metropolis- Hasting抽樣算法進(jìn)行邏輯回歸
9.R語言中基于混合數(shù)據(jù)抽樣(MIDAS)回歸的HAR-RV模型預(yù)測GDP增長
