# R/normConst.R In twDEMC: parallel DEMC

#### Documented in normConstLaplacenormConstLaplaceBridge

#------ developed in scratch postnorm2.R

normConstLaplace <- function(
### Laplace Approximation of the normalizing constant
sample		##<< MCMC sample of parameters, each row of the matrix is a parameter vector
, fLogPost 	##<< function to evaluate log of unnormalized posterior density
, alpha=0.05	##<< quantile defining the multivariate normal ellipsoid for volume correction
, nd = cov.rob(sample)	##<< normal density approximation: list with entries center and cov
){
#hmu <- fPost(nd$center) #phimu <- dmvnorm(nd$center, mean=nd$center, sigma=nd$cov, log=FALSE)
#CL <- hmu/phimu
logPhimu <- dmvnorm(nd$center, mean=nd$center, sigma=nd$cov, log=TRUE) # normal approximation at mode logHmu <- fLogPost(nd$center)	# unnormalized posterior at center \hat{\theta}
logCL <- logHmu - logPhimu
# P: proportion of sample values inside normal ellipsoid defined by alpha
dist2 <- mahalanobis(sample, center=nd$center, cov=nd$cov)	# squared mahalanobis each row is a vector
distCrit2 <- qchisq(alpha, df = ncol(sample))
P <- sum(dist2 < distCrit2)/length(dist2)
#CLV <- CL * alpha/P
logCLV <- logCL + log(alpha/P)
##value<< a list with estimates of the logs of the normalizing constant of the posterior
c(logCL=logCL		##<< Laplace estimate
, logCLV=logCLV)	##<< volume corrected Laplace estimate
##end<<
}
#mtrace(normConstLaplace)
attr(normConstLaplace,"ex") <- function(){
mu = c(0,0)
sd = c(1,2)
corr = diag(nrow=2)
corr[1,2] <- corr[2,1] <- 0.49
Sigma = diag(sd, nrow=length(sd)) %*% corr %*% diag(sd,nrow=length(sd))
n <- 1000
cTrue <- 10
sample <- mvtnorm::rmvnorm(n,mean=mu,sigma=Sigma)
fPost <- function(theta){ cTrue* mvtnorm::dmvnorm(theta,mean=mu,sigma=Sigma, log=FALSE)}
fLogPost <- function(theta){ log(cTrue) + mvtnorm::dmvnorm(theta,mean=mu,sigma=Sigma, log=TRUE)}
exp(normConstLaplace( sample, fLogPost ))

#normd <- cov.rob(sample)
#(res <-  boot( sample, function(sample, i){normConstLaplace(sample[i,], fLogDen=fLogDen, nd=normd)}, 100))
#(res2 <-  boot( sample, function(sample, i){normConstLaplace(sample[i,], fLogDen=fLogDen)}, 100))
}

normConstLaplaceBridge <- function(
### Laplace bridge Approximation of the normalizing constant
sample		##<< MCMC sample of parameters, each row of the matrix is a parameter vector
, logPostSample	##<< the unnormalized posteror density of the sample
, fLogPost 		##<< function to evaluate unnormalized Log-posterior for a given parameter vector
, alpha=0.05	##<< quantile defining the multivariate normal ellipsoid for volume correction
, nd = cov.rob(sample)	##<< normal density approximation: list with entries center and cov
, M = 1e4		##<< number of additional simulations of posterior
, logCest = normConstLaplace(sample, fLogPost, alpha, nd)["logCLV"]	##<< initial estimate of the normalizing constant
){
m <- nrow(sample)
# draw a second sample from normal approximation
sample2 <- rmvnorm(M,mean=nd$center,sigma=nd$cov)
# calcalate unnormalized posterior Density for that
logPost2 <- fLogPost(sample2)
tmp.f <- function(){
postSample <- exp(logPostSample)
post2 <- exp(logPost2)
Cest <- exp(logCest)
gamma2 <- 1/(m*post2/Cest + M*dmvnorm(sample2,mean=nd$center,sigma=nd$cov))
qs1 <- dmvnorm(sample,mean=nd$center,sigma=nd$cov)
gamma1 <- 1/(m*postSample/Cest + M*dmvnorm(sample,mean=nd$center,sigma=nd$cov))
(CMW <- sum(post2*gamma2)/M / ( sum(qs1*gamma1)/m  ))
}
gamma2 <- 1/(m*exp(logPost2-logCest) + M*dmvnorm(sample2,mean=nd$center,sigma=nd$cov))
logQs1 <- dmvnorm(sample,mean=nd$center,sigma=nd$cov, log=TRUE)
gamma1 <- 1/(m*exp(logPostSample-logCest) + M*dmvnorm(sample,mean=nd$center,sigma=nd$cov))
# factor our maximum posterior density to avoid very small exponentials
maxLogPost2 <- max(logPost2)
maxLogQs1 <- max(logQs1)
logCMW <- -maxLogPost2 + log( sum(exp(logPost2+maxLogPost2)*gamma2)/M ) - (-maxLogQs1 +log( sum(exp(logQs1+maxLogQs1)*gamma1)/m ))
### log of the normalization constant
}
attr(normConstLaplaceBridge,"ex") <- function(){
mu = c(0,0)
sd = c(1,2)
corr = diag(nrow=2)
corr[1,2] <- corr[2,1] <- 0.49
Sigma = diag(sd, nrow=length(sd)) %*% corr %*% diag(sd,nrow=length(sd))
n <- 1000
cTrue <- 10
sample <- mvtnorm::rmvnorm(n,mean=mu,sigma=Sigma)
fLogPost <- function(theta){ log(cTrue) + mvtnorm::dmvnorm(theta,mean=mu,sigma=Sigma, log=TRUE)}
logPostSample <- fLogPost(sample)
(CMW <- exp(normConstLaplaceBridge( sample, logPostSample, logPostSample )) )

#normd <- cov.rob(sample)
#(res <-  boot( sample, function(sample, i){normConstLaplace(sample[i,], fLogDen=fLogDen, nd=normd)}, 100))
#(res2 <-  boot( sample, function(sample, i){normConstLaplace(sample[i,], fLogDen=fLogDen)}, 100))

#normd <- cov.rob(sample)
#(res <-  boot( sample, function(sample, i){normConstLaplaceBridge(sample[i,], fLogDen=fLogDen, logDenSample, nd=normd)}, 100))
#(res2 <-  boot( sample, function(sample, i){normConstLaplaceBridge(sample[i,], fLogDen=fLogDen, logDenSample)}, 100))
}


## Try the twDEMC package in your browser

Any scripts or data that you put into this service are public.

twDEMC documentation built on May 31, 2017, 3:44 a.m.