Pengambilan sampel dari distribusi bivariat dengan kerapatan yang diketahui menggunakan MCMC

Saya mencoba mensimulasikan dari kepadatan bivariat menggunakan algoritma Metropolis di R dan tidak beruntung. Kepadatan dapat dinyatakan sebagai , di mana adalah distribusi Singh-Maddala$p(x,y)$$p(y|x)p(x)$$p(x)$

$p(x)=\dfrac{aq x^{a-1}}{b^a (1 + (\frac{x}{b})^a)^{1+q}}$

dengan parameter , , , dan adalah log-normal dengan log-mean sebagai fraksi , dan log-sd adalah konstanta. Untuk menguji apakah sampel saya adalah yang saya inginkan, saya melihat kepadatan marginal , yang seharusnya . Saya mencoba berbagai algoritma Metropolis dari paket R MCMCpack, mcmc dan dream. Saya membuang burn-in, menggunakan penipisan, menggunakan sampel dengan ukuran hingga jutaan, tetapi kepadatan marginal yang dihasilkan tidak pernah yang saya berikan.$a$$q$$b$$p(y|x)$$x$$x$$p(x)$

Ini adalah edisi terakhir dari kode saya yang saya gunakan:

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)
        )

Saya sudah memilih paket mimpi, karena sampel sampai konvergensi. Saya telah menguji apakah saya memiliki hasil yang benar dalam tiga cara. Menggunakan statistik KS, membandingkan kuantil, dan memperkirakan parameter distribusi Singh-Maddala dengan kemungkinan maksimum dari sampel yang dihasilkan:

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")

Ketika saya melihat hasil perbandingan ini, statistik KS hampir selalu menolak hipotesis nol bahwa sampel berasal dari distribusi Singh-Maddala dengan parameter yang disediakan. Parameter estimasi kemungkinan maksimum terkadang mendekati nilai sebenarnya, tetapi biasanya terlalu jauh dari zona nyaman, untuk menerima bahwa prosedur pengambilan sampel berhasil. Ditto untuk kuantil, kuantil empiris tidak terlalu jauh, tetapi terlalu jauh.

Pertanyaan saya adalah apa yang saya lakukan salah? Hipotesis saya sendiri:

1. MCMC tidak sesuai untuk jenis pengambilan sampel ini
2. MCMC tidak dapat konvergen, karena alasan teoritis (fungsi distribusi tidak memenuhi properti yang diperlukan, apa pun itu)
3. Saya tidak menggunakan algoritma Metropolis dengan benar
4. Tes distribusi saya tidak benar, karena saya tidak memiliki sampel independen.

Saya pikir urutannya benar, tetapi label yang ditetapkan untuk p (x) dan p (y | x) salah. Status masalah asli p (y | x) adalah log-normal dan p (x) adalah Singh-Maddala. Jadi begitu

1. Hasilkan tanda X dari Singh-Maddala, dan

2. menghasilkan Y dari log-normal yang memiliki rata-rata yang merupakan sebagian kecil dari X yang dihasilkan.

Sebenarnya, Anda tidak boleh melakukan MCMC, karena masalah Anda jauh lebih sederhana. Coba algoritma ini:

Langkah 1: Hasilkan X dari Log Normal

Langkah 2: Menjaga X ini tetap, menghasilkan Y dari Singh Maddala.

Voa! Sampel Siap !!!