# Sampling from bivariate distribution with known density using MCMC

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)+
}
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) {
}

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,