Description Usage Arguments Author(s) Examples
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
1 | den2dCor(theta, muA = 0.8, sigmaA = 1, cb = 1/4, cDen = 0)
|
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
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] ) )
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.