[R] HPDinterval question - nonlinear transformations/functions of parameters.

[Christos Argyropoulos]

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.