R help - Jul 2010 - HPDinterval question - nonlinear transformations/functions of parameters.

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.

Disregard the post, 
The answer to the question is actually pretty simple, had I done the math. For
completeness here what happens.

The definition of the HPD is the following:

1)? \int_{L}^{U}f(\theta|x)d\theta=1-a
2) f(U|x)=f(L|x)

A one-to-one repameterization

\theta = g(\phi) 

simultaneously transforms the integrand AND the limits of the integration. In
general, HPD regions are not symmetric and (this happens only if the density is
symmetric and unimodal). Even in this case though,? one cannot one cannot apply
the inverse of the g function to [U,L] to compute the HPD unless the?
transformed value is symmetric.
In the example I sent, the transformation (exp) does not lead to a symmetric
random variable. The skewness of f(A,B,C)=exp(A*B+C) is accentuated for larger
values of the covariance matrix of the trivariate normal (A,B,C), so that the
discrepancy increases for smaller value of X.

I apologise for generating noise in the list

C


----------------------------------------
> From: argchris at hotmail.com
> To: r-help at r-project.org
> Subject: HPDinterval question - nonlinear transformations/functions of
parameters.
> Date: Tue, 6 Jul 2010 23:01:11 +0300
>
>
>
> 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.
> https://signup.live.com/signup.aspx?id=60969
 		 	   		  
_________________________________________________________________
Hotmail: Free, trusted and rich email service.