9
$\begingroup$

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.
$\endgroup$
  • $\begingroup$ In the Singh-Maddala distribution link, the pdf has two parameters - {c, k}, yet the R function dsinmad takes three parameters or am I missing something. $\endgroup$ – csgillespie Dec 7 '10 at 11:10
  • $\begingroup$ Sorry, the wikipedia link cites the wrong formula, it looked ok at first glance, when I was composing the question. I did not find a ready link, so I just put the formula in the question. $\endgroup$ – mpiktas Dec 7 '10 at 13:05
3
$\begingroup$

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.

$\endgroup$
3
$\begingroup$

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

$\endgroup$
  • $\begingroup$ I assume that you meant the steps reversed. But if this is so simple why do we need Gibbs sampling? $\endgroup$ – mpiktas Dec 21 '10 at 20:52
  • 1
    $\begingroup$ No, i meant steps 1 and 2 in the order I wrote. After all, the distribution of y is specified conditional on X, so you must generate an X before Y. As for Gibbs sampling, that is a more complicated solution meant for more complicated problems. Yours, as you describe it, is pretty straighforward, IMHO. $\endgroup$ – Mohit Dec 22 '10 at 14:15
  • 1
    $\begingroup$ You would use Gibbs sampling when you know $p(y|x)$ and $p(x|y)$, but not if you know the marginal $p(x)$ $\endgroup$ – probabilityislogic Feb 6 '11 at 1:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.