Hi, I have a quick question about HPDinterval (package coda). I am simulating from a trivariate multivariate normal with 3 components (A,B,C) and general (not necessarily diagonal) covariance matrix. Interest lies in describing the distribution of the function: f(A,B,C) = exp(A*B+C) The question concerns the intervals returned by HPD under two different invokations a) directly on the function f(A,B,C) b) exponentiation of the HPD interval of the function: A*B+C As long as the components of the covariance matrix of the original multivariate normal are small, the intervals returned by the 2 invokations are similar. However a simple scaling of the covariance matrix can turn the agreement to disagreement. (You can reproduce this behaviour by varying the value of X in the R code attached at the end of this email. Smaller values of x lead to greater discrepancies e.g. compare X=0.1 v.s X=5.1). >From my understanding of how HPDinterval works, the intervals returned by the >2 different invokations should be very similar. So what causes this >discrepancy? Which one of the 2 intervals should be used? Regards, Christos Argyropoulos R CODE FOLLOWS library(MASS) library(coda) ## Covariance Matrix X<-5.1 S<-matrix(c( 1.237582e-02, -5.861086e-03, -2.034300e-02 , -5.861086e-03, 1.725401e-02 , 0.234207e-01, -2.034300e-02, 0.234207e-01, 1.410518e-01),3,3)/X cov2cor(S) ## Component Means MU<-c(-0.22263475 , 0.55119463 , -0.08819141) ## Sample from MVNormal set.seed(1234) N<-100000 ran<-mvrnorm(N,MU,S) ## Interested in the following quantity ## log(Tot) = 1st*2nd MVN component + 3rd Component
logtot<-ran[,1]*ran[,2]+ran[,3] HPDinterval(as.mcmc(exp(logtot))) exp(HPDinterval(as.mcmc(logtot))) exp(quantile(logtot,c(0.025,0.975))) plot(density(logtot)) _________________________________________________________________ Hotmail: Trusted email with Microsoft’s powerful SPAM protection. ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
