# den2dCor: den2dCor In twDEMC: parallel DEMC

## Description

Example of log-Unnormalized Density function with changing correlation and scales in variances a ~ N( theta[1], sigmaA ) b ~ N( theta[2], f(a) )

With f(a) exponentially decreasing with increasing a

## Usage

 `1` ```den2dCor(theta, muA = 0.8, sigmaA = 1, cb = 1/4, cDen = 0) ```

## Arguments

 `theta` named numeric vector with components "a" and "b" `muA` parameter for the log-normal distribution of component a `sigmaA` `cb` multiplier for density in b `cDen`

Thomas Wutzler

## Examples

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44``` ```#gridlogx <- seq(log(0.1),log(+4),length.out=91) #gridx <- exp(gridlogx) gridx <- a <- seq(-0.5,2,length.out=91) #plot( lda ~ a) gridy <- seq(-20,+40,length.out=91) gridX <- expand.grid(gridx, gridy) den2dCor(c(0.8,0.8)) luden <- apply( gridX, 1, den2dCor ) mLuden <- matrix(luden,nrow=length(gridx)) #plot( rowSums(mLuden) ~ gridx ) imax <- which( matrix(luden,nrow=length(gridx))==max(luden), arr.ind=TRUE) #c( gridx[ imax[1] ], gridy[ imax[2] ] ) image( gridx, gridy, mLuden, col = rev(heat.colors(100)), xlab="a", ylab="b" ) xyMax <- c(x=gridx[ imax[1] ], y=gridy[ imax[2] ]) image( gridx, gridy, matrix(exp(luden),nrow=length(gridx)), col = rev(heat.colors(100)), xlab="a", ylab="b" ) points( gridx[ imax[1] ], gridy[ imax[2] ] ) q20 <- quantile(luden,0.2) plot(density(luden[luden>q20])) head(sort(luden,dec=TRUE)) ### todo: normalizing constant: ##------------------ do an MCMC run (.expTheta <- c(a=0,b=0) ) (.expCovTheta <- diag(c(a=2,b=2)) ) .nPop=2 Zinit <- initZtwDEMCNormal( .expTheta, .expCovTheta, nChainPop=4, nPop=.nPop) #mtrace(twDEMCBlockInt) den2dCorTwDEMC <- twDEMCBlock(Zinit, nGen=500, dInfos=list(d1=list(fLogDen=den2dCor)), nPop=.nPop ) den2dCorTwDEMC <- twDEMCBlock(den2dCorTwDEMC, nGen=1000) plot( thinN(as.mcmc.list(den2dCorTwDEMC))) matplot( concatPops(den2dCorTwDEMC)\$pAccept[,1,], type="l" ) pps <- pps0 <- stackChains(thin(den2dCorTwDEMC,start=300)) ss <- ss0 <- pps[,-1] #plot( ss[,1], ss[,2] ) plot( ss[,1], ss[,2], ylim=c(-40,80) ) plot( density(ss[,1]) ) plot( ecdf( ss[,1] ) ) plot( ecdf( ss[,2] ) ) ```

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