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

[Christos Argyropoulos]

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.