## Normal tail precision

**I**n conjunction with the normal-Laplace comparison mentioned in the most recent post about our lack of confidence in ABC model choice, we have been working on the derivation of the exact Bayes factor and I derived an easy formula for the marginal likelihood in the Laplace case that boils down to a weighted sum of normal probabilities (with somehow obvious notations)

**I** then wrote a short R code that produced this marginal likelihood.

# ABC model comparison between Laplace and normal # L(mu,V2) versus N(mu,1) with a prior N(0,2*2) nobs=21 nsims=500 sqrtwo=sqrt(2) marnor=function(smpl){ -0.5*nobs*log(2*pi)-0.5*(nobs-1)*var(smpl)+0.5*log(1+1/(4*nobs))-0.5*mean(smpl)^2/(4+1/nobs)} marlap=function(sampl){ smpl=sort(sampl) S=sum(smpl) S=c(S,S-2*cumsum(smpl)) phi=pnorm((smpl-sqrtwo*4*(nobs-2*(1:nobs)))/2) phip=pnorm((smpl-sqrtwo*4*(nobs-2*(1:nobs)+2))/2) Dphi=log(c(phip[1],phip[-1]-phi[-nobs],1-phi[nobs])) -0.5*nobs*log(2)+log(sum(exp(-sqrtwo*S+4*(nobs-2*(0:nobs))^2+Dphi))) }

**W**hen checking it with an alternative Monte Carlo integration, Jean-Michel Marin spotted a slight but persistent numerical difference:

> test=sample(c(-1,1),nobs,rep=TRUE)*rexp(nobs,sqrt(2)) > exp(marlap(test)) [1] 3.074013e-10 > f=function(x){exp(-sqrt(2)*sum(abs(test-x)))*2^(-nobs/2)} > mean(apply(as.matrix(2*rnorm(10^6)),1,f)) [1] 3.126421e-11

**A**nd while I thought it could be due to the simulation error, he persisted in analysing the problem until he found the reason: the difference between the normal cdfs in the above marginal was replaced by zero in the tails of the sample, while it contributed in a significant manner, due to the huge weights in front of those differences! He then rewrote the marlap function so that the difference was better computed in the tails, with a much higher level of agreement!

marlap=function(test){ sigma2=4 lulu=rep(0,nobs-1) test=sort(test) for (i in 1:(nobs-1)){ cst=sqrt(2)*(nobs-2*i)*sigma2 if (test[i]<0) lulu[i]=exp(sqrt(2)*sum(test[1:i])-sqrt(2)*sum(test[(i+1):nobs])+ (nobs-2*i)^2*sigma2+pnorm((test[i+1]-cst)/sqrt(sigma2),log=TRUE)+ log(1-exp(pnorm((test[i]-cst)/sqrt(sigma2),log=TRUE)- pnorm((test[i+1]-cst)/sqrt(sigma2),log=TRUE)))) else lulu[i]=exp(sqrt(2)*sum(test[1:i])-sqrt(2)*sum(test[(i+1):nobs])+ (nobs-2*i)^2*sigma2+pnorm(-(test[i]-cst)/sqrt(sigma2),log=TRUE)+ log(1-exp(pnorm(-(test[i+1]-cst)/sqrt(sigma2),log=TRUE)- pnorm(-(test[i]-cst)/sqrt(sigma2),log=TRUE)))) if (lulu[i]==0) lulu[i]=exp(sqrt(2)*sum(test[1:i])-sqrt(2)*sum(test[(i+1):nobs])+ (nobs-2*i)^2*sigma2+log(pnorm((test[i+1]-cst)/sqrt(sigma2))- pnorm((test[i]-cst)/sqrt(sigma2)))) } lulu0=exp(-sqrt(2)*sum(test[1:nobs])+nobs^2*sigma2+ pnorm((test[1]-sqrt(2)*nobs*sigma2)/sqrt(sigma2),log=TRUE)) lulun=exp(sqrt(2)*sum(test[1:nobs])+nobs^2*sigma2+ pnorm(-(test[nobs]+sqrt(2)*nobs*sigma2)/sqrt(sigma2),log=TRUE)) 2^(-nobs/2)*sum(c(lulu0,lulu,lulun)) }

**H**ere is an example of this agreement:

> marlap(test) [1] 5.519428e-10 mean(apply(as.matrix(2*rnorm(10^6)),1,f)) [1] 5.540964e-10

February 23, 2015 at 10:50 pm

Hi Xi’an. I just started learning basic simulation in R. I’m interested in simulating multivariate correlated data and I’ve been studying the “mvtnorm” package. My question is: using the function rmvnorm, how do I generate correlated outcomes (y1 and y2) that have a skewed distribution? I know already how to generate the correlated outcomes, but I’ve been struggling how to incorporate non-normality in the data. I appreciate any help!

Thanks,

Rommel

February 24, 2015 at 6:24 am

I think you should ask this question on Stack Overflow. But first make the question more precise: do you have a specific correlated skewed distribution in mind or any skewed would do? In the later case, is there anything more you want to specify?

February 24, 2015 at 9:34 pm

I’ve actually searched on Stack Overflow but I did not quite find answers that I was looking for. Anyway, I’m looking to generate a skewed distribution that is moderately non-normal that follows a chi-square distribution with 3 df (that usually reflects a skewness = to 1.63 and a kurtosis equal=4. Thanks for the response, Xi’an.

February 25, 2015 at 9:46 am

Have you tried through copulas? You can specify both marginals and set the correlation to the level you wish.

June 28, 2011 at 2:04 am

[…] my earlier posts on the revision of Lack of confidence, here is an interesting outcome from the derivation of the exact marginal likelihood in the Laplace case. Computing the posterior probability of a normal model versus a Laplace model […]