Sampling from bivariate distribution with known density using MCMC

Question

I tried to simulate from a bivariate density $p(x,y)$ using Metropolis algorithms in R and had no luck. The density can be expressed as $p(y|x)p(x)$, where $p(x)$ is Singh-Maddala distribution

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

with parameters $a$, $q$, $b$, and $p(y|x)$ is log-normal with log-mean as a fraction of $x$, and log-sd a constant. To test whether my sample is the one I want, I looked at the marginal density of $x$, which should be $p(x)$. I tried different Metropolis algorithms from R packages MCMCpack, mcmc and dream. I discarded burn-in, used thinning, used samples with size up to million, but the resulting marginal density was never the one I supplied.

Here is the final edition of my code I used:

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

I've settled on dream package, since it samples until the convergence. I've tested whether I have the correct results in three ways. Using KS statistic, comparing quantiles, and estimating the parameters of Singh-Maddala distribution with maximum likelihood from the resulting sample:

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

When I look at the results of these comparisons, KS statistic almost always rejects the null hypothesis that sample is from Singh-Maddala distribution with supplied parameters. Maximum likelihood estimated parameters sometimes comes close to its true values, but usually too far out of comfort zone, to accept that sampling procedure was succesfull. Ditto for the quantiles, empirical quantiles are not too far, but too far away.

My question is what I am doing wrong? My own hypotheses:

1. MCMC is not appropriate for this type of sampling
2. MCMC cannot converge, due to theoretical reasons (the distribution function does not satisfy required properties, whatever they are)
3. I do not use the Metropolis algorithm correctly
4. My distribution tests are not correct, since I do not have independent sample.

Answer 1

I think the order is correct, but the labels assigned to p(x) and p(y|x) were wrong. The original problem states p(y|x) is log-normal and p(x) is Singh-Maddala. So, it's

1. Generate an X from a Singh-Maddala, and

2. generate a Y from a log-normal having a mean which is a fraction of the generated X.

Answer 2

In fact, you should not do MCMC, since your problem is so much simpler. Try this algorithm:

Step 1: Generate a X from Log Normal

Step 2: Keeping this X fixed, generate a Y from the Singh Maddala.

Voilà! Sample Ready!!!