使用 MCMC 從已知密度的雙變量分佈中採樣

我試圖從雙變量密度進行模擬在 R 中使用 Metropolis 算法並沒有運氣。密度可以表示為 ， 在哪裡是 Singh-Maddala 分佈

帶參數,,， 和是對數正態的，對數均值是, 和 log-sd 一個常數。為了測試我的樣本是否是我想要的樣本，我查看了，應該是. 我從 R 包 MCMCpack、mcmc 和 dream 中嘗試了不同的 Metropolis 算法。我丟棄了老化，使用了減薄，使用了尺寸高達百萬的樣品，但由此產生的邊際密度從來不是我提供的。

這是我使用的代碼的最終版本：

logvrls <- function(x,el,sdlog,a,scl,q.arg) {
   if(x[2]>0) {
        dlnorm(x[1],meanlog=el*log(x[2]),sdlog=sdlog,log=TRUE)+
        dsinmad(x[2],a=a,scale=scl,q.arg=q.arg,log=TRUE)
   }
   else -Inf    
}

a <- 1.35
q <- 3.3
scale <- 10/gamma(1 + 1/a)/gamma(q - 1/a)*  gamma(q) 

Initvrls <- function(pars,nseq,meanlog,sdlog,a,scale,q) {
   cbind(rlnorm(nseq,meanlog,sdlog),rsinmad(nseq,a,scale,q))
}

library(dream)
aa <- dream(logvrls,
       func.type="logposterior.density",
       pars=list(c(0,Inf),c(0,Inf)),
       FUN.pars=list(el=0.2,sdlog=0.2,a=a,scl=scale,q.arg=q),
       INIT=Initvrls,
       INIT.pars=list(meanlog=1,sdlog=0.1,a=a,scale=scale,q=q),
       control=list(nseq=3,thin.t=10)
       )

我已經選擇了夢想包，因為它採樣直到收斂。我已經通過三種方式測試了我是否有正確的結果。使用 KS 統計量，比較分位數，並從結果樣本中以最大似然估計 Singh-Maddala 分佈的參數：

ks.test(as.numeric(aa$Seq[[2]][,2]),psinmad,a=a,scale=scale,q.arg=q)

lsinmad <- function(x,sample)
   sum(dsinmad(sample,a=x[1],scale=x[2],q.arg=x[3],log=TRUE))
optim(c(2,20,2),lsinmad,method="BFGS",sample=aa$Seq[[1]][,2])

qq <- eq(0.025,.975,by=0.025)   
tst <- cbind(qq,
             sapply(aa$Seq,function(l)round(quantile(l[,2],qq),3)),
             round(qsinmad(qq,a,scale,q),3))
colnames(tst) <- c("Quantile","S1","S2","S3","True")

library(ggplot2)
qplot(x=Quantile,y=value,
      data=melt(data.frame(tst),id=1), 
      colour=variable,group=variable,geom="line")

當我查看這些比較的結果時，KS 統計量幾乎總是拒絕零假設，即樣本來自 Singh-Maddala 分佈並提供參數。最大似然估計參數有時會接近其真實值，但通常離舒適區太遠，無法接受采樣過程是成功的。分位數也是如此，經驗分位數不是太遠，而是太遠了。

我的問題是我做錯了什麼？我自己的假設：

1. MCMC 不適用於此類採樣
2. 由於理論上的原因，MCMC 不能收斂（分佈函數不滿足所需的屬性，無論它們是什麼）
3. 我沒有正確使用 Metropolis 算法
4. 我的分佈測試不正確，因為我沒有獨立樣本。

我認為順序是正確的，但是分配給 p(x) 和 p(y|x) 的標籤是錯誤的。原始問題狀態 p(y|x) 是對數正態的，p(x) 是 Singh-Maddala。所以就是

1. 從 Singh-Maddala 生成 X，並且
2. 從具有平均值的對數正態生成 Y，該平均值是生成的 X 的一部分。

引用自：https://stats.stackexchange.com/questions/5207