Muestreo de distribución bivariante con conocidos de la densidad de uso de MCMC

Question

Traté de simular de una densidad bivariante $p(x,y)$ el uso de Metropolis algoritmos en R y no tuvo suerte. La densidad puede ser expresado como $p(y|x)p(x)$ donde $p(x)$ es Singh-Maddala distribución

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

con los parámetros de $a$, $q$, $b$, y $p(y|x)$ es log-normal log-media como una fracción de $x$, y log-sd una constante. Para probar si mi muestra es la que yo quiero, me miró a los marginales de la densidad de $x$, que debe ser $p(x)$. He probado diferentes Metrópolis algoritmos de R paquetes MCMCpack, mcmc y el sueño. Yo lo descarté de burn-in, usa adelgazamiento, utilizado muestras con un tamaño de hasta millones, pero el resultado marginal de la densidad fue nunca el que me suministra.

Aquí está la edición final de mi código que he usado:

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

Me he asentado en paquete sueño, ya que las muestras hasta la convergencia. He probado ya tengo los resultados correctos en tres maneras. El uso de KS estadística, la comparación de cuantiles, y la estimación de los parámetros de Singh-Maddala de distribución con la máxima probabilidad de la muestra resultante:

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

Cuando veo los resultados de estas comparaciones, KS estadística casi siempre se rechaza la hipótesis nula de que la muestra es de Singh-Maddala de distribución con los parámetros proporcionados. La máxima probabilidad de parámetros estimados a veces viene cerca de sus verdaderos valores, pero por lo general demasiado lejos de la zona de confort, a aceptar que el procedimiento de muestreo fue exitoso. Lo mismo para los cuantiles empíricos de cuantiles no están muy lejos, pero muy lejos.

Mi pregunta es, ¿qué estoy haciendo mal? Mi propia hipótesis:

1. MCMC no es adecuada para este tipo de muestreo
2. MCMC no convergen, debido a razones teóricas (la función de distribución no cumple con las propiedades requeridas, sean las que sean)
3. Yo no uso el algoritmo de Metropolis correctamente
4. Mi distribución de las pruebas no son correctos, ya que no tengo de muestras independientes.

Answer 1

De hecho, usted no debe hacer MCMC, ya que su problema es mucho más sencillo. Pruebe este algoritmo:

Paso 1: Generar un X de Registro Normal

Paso 2: Mantener esta X fijo, generar Y de la Singh Maddala.

Voila! Ejemplo De Lista!!!

Answer 2

Creo que el orden es correcto, pero las etiquetas asignadas a p(x) y p(y|x) estaban equivocados. El problema original de los estados p(y|x) es log-normal y p(x) es Singh-Maddala. Así, es

1. Generar una X de una Singh-Maddala, y

2. generar Y de una log-normal con una media que es una fracción de la que genera X.