Nothing
inst/unitTests/runittwDEMC.R
In twDEMC: parallel DEMC
## Test unit 'twDEMC'
# twUtestF(twDEMC)
# twUtestF(twDEMC, "test.goodStartSeqData")
# twUtestF(twDEMC, "test_int.goodStartSeqData.plot2d")
# twUtestF(twDEMC, "dump")
.setUp <- function(){
data(twLinreg1)
data(twdemcEx1)
suppressWarnings( rm( list=names(twLinreg1) ) )
attach(twLinreg1)
ex1c <- concatPops(twdemcEx1)
suppressWarnings( rm( list=names(ex1c) ) )
#attach(ex1c)
}
.tearDown <- function(){
#detach( ex1c )
detach( twLinreg1 )
}
.tmp.f <- function(){
trace( updateBlockTwDEMC, recover ) # untrace(updateBlockTwDEMC)
mtrace(logDenGaussian)
mtrace(twCalcLogDenPar)
mtrace(.doDEMCStep)
mtrace(.doDEMCSteps)
mtrace(.generateXProb)
mtrace(.generateXPropChains)
mtrace(.xStepSnooker)
mtrace(twDEMCBlockInt)
mtrace(.generateXPropThin)
mtrace(.xStepParallel)
twUtestF(twDEMC,"test.probUpDir")
twUtestF(twDEMC,"test.ZinittwDEMC") #with thinning interval
twUtestF(twDEMC,"test.goodStartSeq")
twUtestF(twDEMC,"test.goodStartPar")
twUtestF(twDEMC,"test.doAppendPrevLogDen")
twUtestF(twDEMC,"test.goodStartSeqData")
twUtestF(twDEMC)
twUtestF(twDEMCBatch)
twUtestF()
prevWd <- setwd(file.path("tests"))
dumpFileBasename="doRUnitError"
#options(error=quote(dump.frames(dumpFileBasename, TRUE)))
load(paste(dumpFileBasename,".rda"))
debugger(dumpFileBasename)
setwd(prevWd)
}
test.goodStartSeqData <- function(){
#mtrace(test_int.goodStartSeqData.plot2d)
test_int.goodStartSeqData.plot2d()
}
test_int.goodStartSeqData.plot2d <- function(){
lmDummy <- lm( obs ~ xval, weights=1/sdObs^2) # results without priors
(.expTheta <- coef(lmDummy))
(.expCovTheta <- vcov(lmDummy)) # a is very weak constrained, negative covariance between a an b
(.expSdTheta <- structure(sqrt(diag(.expCovTheta)), names=c("a","b")) )
# nice starting values
.nPop=2
argsFLogDen <- list(
fModel=dummyTwDEMCModel, ### the model function, which predicts the output based on theta
obs=obs, ### vector of data to compare with
invCovar=invCovar, ### do not constrain by data, the inverse of the Covariance of obs (its uncertainty)
#thetaPrior = thetaTrue, ### the prior estimate of the parameters
#invCovarTheta = invCovarTheta, ### the inverse of the Covariance of the prior parameter estimates
xval=xval
)
do.call( logDenGaussian, c(list(theta=theta0),argsFLogDen))
Zinit <- initZtwDEMCNormal( theta0, .expCovTheta, nChainPop=4, nPop=.nPop)
#dim(Zinit)
.nGen=100
#mtrace(logDenGaussian)
#mtrace(twDEMCBlockInt)
res <- concatPops( resBlock <- twDEMCBlock( Zinit, nGen=.nGen,
dInfos=list(list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen)),
nPop=.nPop
,controlTwDEMC=list(thin=8 )
#,debugSequential=TRUE
,doRecordProposals=TRUE
))
#str(res)
checkEquals( 8, res$thin )
checkEquals( (.nGen%/%res$thin)+1,nrow(res$logDen))
rescoda <- as.mcmc.list(res)
plot(rescoda, smooth=FALSE)
#gelman.diag(rescoda)
#summary(rescoda)
.tmp.f <- function(){
matplot( res$logDen[,1,],type="l" )
tmp <- data.frame( a=as.numeric(res$Y[,"a",]), b=as.numeric(res$Y[,"b",]), logDen=as.numeric(rowSums(res$Y[,res$YPos$resLogDen0+res$dInfos[[1]]$compPosDen,])) )
colorFun <- colorRampPalette(c("yellow","red"))
levelplot( logDen ~ a + b, tmp, col.regions = colorFun(50))
tmp2 <- data.frame( a=as.numeric(res$parms[,"a",]), b=as.numeric(res$parms[,"b",]), logDen=as.numeric(res$logDen) )
levelplot( logDen ~ a + b, tmp2, col.regions = colorFun(50))
#image( as.numeric(res$Y["a",,]), as.numeric(res$Y["b",,]), as.numeric(res$Y["logDen",,]))
}
#str(summary(rescoda))
suppressWarnings({
summary(rescoda)$statistics[,"Mean"]
summary(rescoda)$statistics[,"SD"]
thetaTrue
sdTheta
(.popmean <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"Mean"]}))
(.popsd <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"SD"]}))
})
# 1/2 orders of magnitude around prescribed sd for theta
.pop=1
for( .pop in seq(along.with=.popsd) ){
checkMagnitude( .expSdTheta, .popsd[[.pop]] )
}
# check that thetaTrue is in 95% interval
.pthetaTrue <- sapply(1:2, function(.pop){
pnorm(.expTheta, mean=.popmean[[.pop]], sd=.popsd[[.pop]])
})
checkInterval( .pthetaTrue )
# repeat with control for overfitting
resCO <- concatPops( resBlockCO <- twDEMCBlock( Zinit, nGen=.nGen,
dInfos=list( list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen, maxLogDen=c(obs=-length(obs)/2,parms=0)) ),
nPop=.nPop
,controlTwDEMC=list(thin=8)
#,debugSequential=TRUE
,doRecordProposals=TRUE
))
rescoda <- as.mcmc.list(resCO)
plot(rescoda, smooth=FALSE)
#str(summary(rescoda))
suppressWarnings({
summary(rescoda)$statistics[,"Mean"]
summary(rescoda)$statistics[,"SD"]
thetaTrue
sdTheta
(.popmean <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"Mean"]}))
(.popsd <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"SD"]}))
})
# 1/2 orders of magnitude around prescribed sd for theta
.pop=1
for( .pop in seq(along.with=.popsd) ){
checkMagnitude( .expSdTheta, .popsd[[.pop]] )
}
# check that thetaTrue is in 95% interval
.pthetaTrue <- sapply(1:2, function(.pop){
pnorm(.expTheta, mean=.popmean[[.pop]], sd=.popsd[[.pop]])
})
checkInterval( .pthetaTrue )
}
test.badStartSeqData <- function(){
lmDummy <- lm( obs ~ xval, weights=1/sdObs^2) # results without priors
(.expTheta <- coef(lmDummy))
(.expCovTheta <- vcov(lmDummy)) # a is very weak constrained, negative covariance between a an b
(.expSdTheta <- structure(sqrt(diag(.expCovTheta)), names=c("a","b")) )
.nPop=2
argsFLogDen <- list(
fModel=dummyTwDEMCModel, ### the model function, which predicts the output based on theta
obs=obs, ### vector of data to compare with
invCovar=invCovar, ### do not constrain by data, the inverse of the Covariance of obs (its uncertainty)
#thetaPrior = thetaTrue, ### the prior estimate of the parameters
#invCovarTheta = invCovarTheta, ### the inverse of the Covariance of the prior parameter estimates
xval=xval
)
do.call( logDenGaussian, c(list(theta=theta0),argsFLogDen))
#Zinit <- initZtwDEMCNormal( theta0, .expCovTheta, nChainPop=4, nPop=.nPop)
.sdThetaBad <- sdTheta*c(1/10,10)
.theta0Bad <- theta0*c(10,1/10)
Zinit <- initZtwDEMCNormal( .theta0Bad, diag(.sdThetaBad^2), nChainPop=4, nPop=.nPop)
Zinit[,,1:4] <- Zinit[,,1:4] * c(2,1)
Zinit[,,-(1:4)] <- Zinit[,,-(1:4)] * c(1,2)
#dim(Zinit)
.nGen=200
resa <- concatPops( resBlock <- twDEMCBlock( Zinit, nGen=.nGen,
dInfos=list(list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen)),
nPop=.nPop,
controlTwDEMC=list(thin=1),
debugSequential=TRUE
))
checkEquals((.nGen%/%resa$thin)+1,nrow(resa$logDen))
rescoda <- as.mcmc.list(resa)
plot(rescoda, smooth=FALSE)
#mtrace(thin.twDEMC)
res <- thin(resa, start=120)
rescoda <- as.mcmc.list(res)
plot(rescoda, smooth=FALSE)
#gelman.diag(rescoda)
#summary(rescoda)
#str(summary(rescoda))
suppressWarnings({ #glm fit in summary
summary(rescoda)$statistics[,"Mean"]
summary(rescoda)$statistics[,"SD"]
thetaTrue
sdTheta
(.popmean <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"Mean"]}))
(.popsd <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"SD"]}))
})
# 1/2 orders of magnitude around prescribed sd for theta
.pop=1
for( .pop in seq(along.with=.popsd) ){
checkMagnitude( .expSdTheta, .popsd[[.pop]] )
}
# check that thetaTrue is in 95% interval
.pthetaTrue <- sapply(1:2, function(.pop){
pnorm(.expTheta, mean=.popmean[[.pop]], sd=.popsd[[.pop]])
})
checkInterval( .pthetaTrue )
}
test.badStartSeqData1D <- function(){
# common error is dropping dimensions, if of size 1 onyle: test parameter_dimension = 1
set.seed(0815)
lmDummy <- lm( obs-theta0["a"] ~ xval-1, weights=1/sdObs^2) # results without priors
theta1d <- theta0["b"]
(.expTheta <- coef(lmDummy))
(.expCovTheta <- vcov(lmDummy))
(.expSdTheta <- structure(sqrt(diag(.expCovTheta)), names=c("b")) )
#trace( dummyTwDEMCModel, recover ) # untrace(dummyTwDEMCModel)
argsFLogDen <- list(
fModel=dummyTwDEMCModel, ### the model function, which predicts the output based on theta
obs=obs, ### vector of data to compare with
invCovar=invCovar, ### do not constrain by data, the inverse of the Covariance of obs (its uncertainty)
#thetaPrior = thetaTrue, ### the prior estimate of the parameters
invCovarTheta = invCovarTheta[2,2,drop=FALSE], ### the inverse of the Covariance of the prior parameter estimates
xval=xval
,theta0=theta0
)
do.call( logDenGaussian, c(list(theta=theta1d),argsFLogDen))
.tmp.f <- function(){
bs <- seq(2,8,length.out=30)
ll <- colSums( sapply( bs, function(b){ do.call( logDenGaussian, c(list(theta=c(b=b)),argsFLogDen)) }) )
plot( ll ~ bs )
}
#Zinit <- initZtwDEMCNormal( theta0, .expCovTheta, nChainPop=4, nPop=.nPop)
.sdThetaBad <- sdTheta["b"]*c(1/10)
.theta0Bad <- theta0["b"]*c(10)
.nPop=2
Zinit <- initZtwDEMCNormal( .theta0Bad, diag(.sdThetaBad^2,nrow = length(.sdThetaBad)), nChainPop=4, nPop=.nPop)
Zinit[,,1:4] <- Zinit[,,1:4] * c(2)
Zinit[,,-(1:4)] <- Zinit[,,-(1:4)] * c(0.5)
#dim(Zinit)
.nGen=250
#mtrace(twDEMCBlockInt)
#mtrace(.sampleStates)
#mtrace(twWhichColsEqual)
resa <- concatPops(resBlock <- twDEMCBlock( Zinit, nGen=.nGen,
dInfos=list(list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen)),
nPop=.nPop,
controlTwDEMC=list(thin=1)
,debugSequential=TRUE
))
#str(res)
plot( as.mcmc.list(resa), smooth=FALSE )
checkEquals((.nGen%/%resa$thin)+1,nrow(resa$logDen))
matplot( resa$parms[,1,],type="l")
res <- thin(resa, start=110)
matplot( res$parms[,1,],type="l")
#gelman.diag(rescoda)
#summary(rescoda)
#str(summary(rescoda))
suppressWarnings({ #glm fit in summary
(.popmean <- lapply(list(p1=1:4,p2=5:8),function(i){mean(res$parms[1,,i])}))
(.popsd <- lapply(list(p1=1:4,p2=5:8),function(i){sd(as.vector(res$parms[1,,i]))}))
})
# 1/2 orders of magnitude around prescribed sd for theta
.pop=1
for( .pop in seq(along.with=.popsd) ){
checkMagnitude( .expSdTheta, .popsd[[.pop]], msg="wrong magnitude of standard deviation" )
}
# check that thetaTrue is in 95% interval
.pthetaTrue <- sapply(1:2, function(.pop){
pnorm(.expTheta, mean=.popmean[[.pop]], sd=.popsd[[.pop]])
})
checkInterval( .pthetaTrue, 0.005, 0.995, "thetaTrue outside 0.99% confidence interval." )
}
test.goodStartPrior <- function(){
# nice starting values and prior constrains estimates
.nPop=2
argsFLogDen <- list(
fModel=dummyTwDEMCModel, ### the model function, which predicts the output based on theta
obs=obs, ### vector of data to compare with
invCovar=invCovar/1e10, ### do not constrain by data, the inverse of the Covariance of obs (its uncertainty)
thetaPrior = thetaTrue, ### the prior estimate of the parameters
invCovarTheta = invCovarTheta, ### the inverse of the Covariance of the prior parameter estimates
xval=xval
)
do.call( logDenGaussian, c(list(theta=theta0),argsFLogDen))
Zinit <- initZtwDEMCNormal( theta0, diag(sdTheta^2), nChainPop=4, nPop=.nPop)
#dim(Zinit)
.nGen=100
#mtrace(updateBlockTwDEMC)
#mtrace(twDEMCBlockInt)
#mtrace(.updateIntervalTwDEMCPar)
#mtrace(.updateBlocksTwDEMC)
resa <- res <- resSeq <- concatPops(resBlock <- twDEMCBlock( Zinit, nGen=.nGen,
dInfos=list(list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen)),
nPop=.nPop,
controlTwDEMC=list(thin=4, useMultiT=TRUE)
,debugSequential=TRUE
#,T0=10, TProp=c( obs=0.5, parms=1 )
))
#str(res)
checkEquals((.nGen%/%resa$thin)+1,nrow(resa$logDen))
#windows(record=TRUE)
rescoda <- as.mcmc.list(resa)
plot(rescoda, smooth=FALSE)
#gelman.diag(rescoda)
#summary(rescoda)
res <- thin(resa, start=30)
rescoda <- as.mcmc.list(res)
plot(rescoda, smooth=FALSE)
.tmp.f <- function(){
colnames(res$resLogDen)
matplot(res$resLogDen[,"obs",], type="l")
matplot(res$resLogDen[,"parms",], type="l")
matplot(res$logDen[,1,], type="l")
tail(res$Y[,,1])
}
#str(summary(rescoda))
suppressWarnings({ #glm fit in summary
summary(rescoda)$statistics[,"Mean"]
summary(rescoda)$statistics[,"SD"]
thetaTrue
sdTheta
(.popmean <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"Mean"]}))
(.popsd <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"SD"]}))
})
# 1/2 orders of magnitude around prescribed sd for theta
.pop=1
for( .pop in seq(along.with=.popsd) ){
checkMagnitude( sdTheta, .popsd[[.pop]] )
}
# check that thetaTrue is in 95% interval
.pthetaTrue <- sapply(1:2, function(.pop){
pnorm(thetaTrue, mean=.popmean[[.pop]], sd=.popsd[[.pop]])
})
checkInterval( .pthetaTrue ) #
#------------ again with parallel runs
resa <- resPar <- concatPops(resBlock <- twDEMCBlock( Zinit, nGen=100,
dInfos=list(list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen)),
nPop=.nPop
#fLogDenScale=-1/2
#debugSequential=TRUE
))
rescoda <- as.mcmc.list(thin(resa,start=30))
#str(summary(rescoda))
suppressWarnings({ #glm fit in summary
summary(rescoda)$statistics[,"Mean"]
summary(rescoda)$statistics[,"SD"]
thetaTrue
sdTheta
(.popmean <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"Mean"]}))
(.popsd <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"SD"]}))
})
# 1/2 orders of magnitude around prescribed sd for theta
.pop=1
for( .pop in seq(along.with=.popsd) ){
checkMagnitude( sdTheta, .popsd[[.pop]] )
}
# check that thetaTrue is in 95% interval
.pthetaTrue <- sapply(1:2, function(.pop){
pnorm(thetaTrue, mean=.popmean[[.pop]], sd=.popsd[[.pop]])
})
checkInterval( .pthetaTrue ) #
}
test.upperParBounds <- function(){
# same as goodStartPar, but giving an upper parameter bound
.nPop=2
argsFLogDen <- list(
fModel=dummyTwDEMCModel, ### the model function, which predicts the output based on theta
obs=obs, ### vector of data to compare with
invCovar=invCovar, ### the inverse of the Covariance of obs (its uncertainty)
thetaPrior = thetaTrue, ### the prior estimate of the parameters
invCovarTheta = invCovarTheta, ### the inverse of the Covariance of the prior parameter estimates
xval=xval
)
#do.call( logDenGaussian, c(list(theta=theta0),argsFLogDen))
asplit <- 10.0
Zinit0 <- initZtwDEMCNormal( theta0, diag(sdTheta^2), nChainPop=4, nPop=.nPop)
# replace all the cases with upperParBounds
bo.keep <- as.vector(Zinit0[,"a",]) <= asplit
#mtrace(replaceZinitCases)
Zinit1 <- replaceZinitCases(Zinit0, bo.keep)
Zinit2 <- replaceZinitCases(Zinit0, !bo.keep)
dim(Zinit1)
#mtrace(.doDEMCStep)
res <- res0 <- concatPops(resBlock <- twDEMCBlock( Zinit0, nGen=100,
dInfos=list(list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen)),
nPop=.nPop
#,upperParBounds=c(a=asplit)
#,debugSequential=TRUE
#,controlTwDEMC=list( DRgamma=0.1, minPCompAcceptTempDecr=0.16)
))
res <- res1 <- concatPops(resBlock <- twDEMCBlock( Zinit1, nGen=64*4,
dInfos=list(list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen)),
nPop=.nPop
,upperParBounds=c(a=asplit)
,debugSequential=TRUE
#,controlTwDEMC=list( DRgamma=0.1, minPCompAcceptTempDecr=0.16)
))
res <- res2 <- concatPops(resBlock <- twDEMCBlock( Zinit2, nGen=64*4,
dInfos=list(list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen)),
nPop=.nPop
,lowerParBounds=c(a=asplit)
,debugSequential=TRUE
))
ss0 <- stackChains(res0)
# importance sampling due to integrated probability
ss1 <- stackChains(res1)
lw1 <- twLogSumExp(ss1[,1])
ss2 <- stackChains(res2)
lw2 <- twLogSumExp(ss2[,1])
lwSum <- twLogSumExp( c(lw1,lw2) )
pSubs <- c( exp(lw1-lwSum), exp(lw2-lwSum) )
nSubs <- floor(pSubs/max(pSubs) * nrow(ss0))
ssc <- rbind( ss1[sample.int(nrow(ss1),nSubs[1]),]
, ss2[sample.int(nrow(ss2),nSubs[2]),]
)
.tmp.f <- function(){
str(res)
rescoda <- as.mcmc.list(res)
plot(rescoda, smooth=FALSE)
plot( density(ssc[,"a"])) # combined
lines(density(ss0[,"a"]), col="blue") # orig
lines(density(ss1[,"a"]), col="green") # left
lines(density(ss2[,"a"]), col="darkgreen") # right
gelman.diag(res1)
gelman.diag(res2)
}
checkTrue( all( ss1[,"a"] <= asplit ) )
checkTrue( all( ss2[,"a"] >= asplit ) )
.popmean <- colMeans( ssc[,-1])
.popsd <- apply(ssc[,-1], 2, sd )
# 1/2 orders of magnitude around prescribed sd for theta
.pop=1
for( .pop in seq(along.with=.popsd) ){
checkMagnitude( sdTheta, .popsd[[.pop]] )
}
# check that thetaTrue is in 95% interval
.pthetaTrue <- sapply(1:2, function(.pop){
pnorm(thetaTrue, mean=.popmean[[.pop]], sd=.popsd[[.pop]])
})
checkInterval( .pthetaTrue )
}
t_est.probUpDir <- function(){
#mtrace(inner.probUpDir)
innertest.probUpDir()
#twUtestF(twDEMC,"test.probUpDir")
}
i_nnertest.probUpDir <- function(){
# same as goodStartSeq, but with executing debugSequential=FALSE, i.e. parallel
.nPop=2
argsFLogDen <- list(
fModel=dummyTwDEMCModel, ### the model function, which predicts the output based on theta
obs=obs, ### vector of data to compare with
invCovar=invCovar, ### the inverse of the Covariance of obs (its uncertainty)
thetaPrior = thetaTrue, ### the prior estimate of the parameters
invCovarTheta = invCovarTheta, ### the inverse of the Covariance of the prior parameter estimates
xval=xval
)
#do.call( logDenGaussian, c(list(theta=theta0),argsFLogDen))
#
Zinit <- initZtwDEMCNormal( theta0, diag(sdTheta^2), nChainPop=4, nPop=.nPop)
dim(Zinit)
#
res <- concatPops(resBlock <- twDEMCBlock( Zinit, nGen=100,
dInfos=list(list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen)),
nPop=.nPop,
#fLogDenScale=-1/2,
controlTwDEMC = list(probUpDir=0.8)
))
str(res)
#
rescoda <- as.mcmc.list(res)
suppressWarnings({ #glm fit in summary
(.popmean <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"Mean"]}))
(.popsd <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"SD"]}))
})
#
# 1/2 orders of magnitude around prescribed sd for theta
.pop=1
for( .pop in seq(along.with=.popsd) ){
checkMagnitude( sdTheta, .popsd[[.pop]] )
}
#
# check that thetaTrue is in 95% interval
.pthetaTrue <- sapply(1:2, function(.pop){
pnorm(thetaTrue, mean=.popmean[[.pop]], sd=.popsd[[.pop]])
})
checkInterval( .pthetaTrue )
}
test.ofMulti <- function(){
# testing multiple log-Density functions
data(twTwoDenEx1)
thresholdCovar = 0.3 # the true value
#thresholdCovar = 0 # the effective model that glosses over this threshold
thetaMean <- twTwoDenEx1$thetaTrue
sdTheta <- (thetaMean*0.3)
argsFLogDen <- list(
thresholdCovar=thresholdCovar
, thetaPrior=thetaMean
, invCovarTheta=sdTheta^2
, twTwoDenEx=twTwoDenEx1
)
#
.nPop=2
.nChainPop=4
ZinitPops <- with(twTwoDenEx1, initZtwDEMCNormal( thetaMean*1.5, diag(sdTheta^2), nChainPop=.nChainPop, nPop=.nPop))
#
.nGen=128
resa <- concatPops( resBlock <- twDEMCBlock(
ZinitPops
, nGen=.nGen
,dInfos=list(
dSparse=list(fLogDen=denSparsePrior, argsFLogDen=argsFLogDen)
,dRich=list(fLogDen=denRichPrior, argsFLogDen=argsFLogDen)
)
,blocks = list(
a=list(dInfoPos="dSparse", compPos="a")
,b=list(dInfoPos="dRich", compPos="b")
)
,nPop=.nPop
,controlTwDEMC=list(thin=8)
#,debugSequential=TRUE
,doRecordProposals=TRUE
))
plot( as.mcmc.list(resa), smooth=FALSE)
rescoda <- as.mcmc.list(thin(resa,start=.nGen%/%2))
#
#i=1
.popmean <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"Mean"]})
.popsd <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"SD"]})
# 1/2 orders of magnitude around prescribed sd for theta
# no true sdTheta
#.pop=1
#for( .pop in seq(along.with=.popsd) ){
# checkMagnitude( sdTheta, .popsd[[.pop]] )
#}
#
# check that thetaTrue is in 95% interval
.pthetaTrue <- sapply(1:2, function(.pop){
pnorm(thetaMean, mean=.popmean[[.pop]], sd=.popsd[[.pop]])
})
checkInterval( .pthetaTrue )
}
test.ofMultiIntermediate <- function(){
int.ofMultiIntermediate()
}
test.ofMultiIntermediateCond <- function(){
int.ofMultiIntermediate( useConditionalProposal=TRUE )
}
test.ofMultiIntermediateControlOverfit <- function(){
int.ofMultiIntermediate( controlOverfitting=TRUE )
}
int.ofMultiIntermediate <- function(useConditionalProposal=FALSE, controlOverfitting=FALSE){
# testing multiple log-Density functions using intermediate results of each other (argument intermediate to dInfos)
data(twTwoDenEx1)
thresholdCovar = 0.3 # the true value
#thresholdCovar = 0 # the effective model that glosses over this threshold
thetaMean <- twTwoDenEx1$thetaTrue
sdTheta <- (thetaMean*0.3)
#
argsFLogDen <- list(
thresholdCovar=thresholdCovar
, thetaPrior=thetaMean
, invCovarTheta=sdTheta^2
, twTwoDenEx=twTwoDenEx1
)
.nPop=2
.nChainPop=4
ZinitPops <- with(twTwoDenEx1, initZtwDEMCNormal( thetaMean*1.5, diag(sdTheta^2), nChainPop=.nChainPop, nPop=.nPop))
#
.nGen=128
#undebug(denRichPrior)
.remoteDumpfileBasename="tmp/dump_twDEMC_ofMultiIntermediate"
resa <- concatPops( resBlock <- twDEMCBlock(
ZinitPops
, nGen=.nGen
,dInfos=list(
dSparse=list(fLogDen=denSparsePrior, argsFLogDen=argsFLogDen, maxLogDen=c("parmsSparse"=0, "obsSparse"=-10/2 ))
,dRich=list(fLogDen=denRichPrior, argsFLogDen=argsFLogDen, maxLogDen=c("parmsRich"=0, "obsRich"=-1000/2 ))
)
,blocks = list(
a=list(compPos="a", dInfoPos="dSparse", intermediateId="modTwoDenEx")
,b=list(compPos="b", dInfoPos="dRich", intermediateId="modTwoDenEx")
)
,nPop=.nPop
,controlTwDEMC=list(thin=8, useConditionalProposal=useConditionalProposal, controlOverfittingMinNObs=if(controlOverfitting) 10 else Inf )
#,debugSequential=TRUE
,remoteDumpfileBasename = .remoteDumpfileBasename
,doRecordProposals=TRUE
))
if( controlOverfitting ){
checkTrue( all(resa$resLogDen[,"obsSparse",] <= -5) ) # not across upper bound of logDen
checkTrue( all(resa$resLogDen[,"obsRich",] <= -500) ) # not across upper bound of logDen
}
# summary( resa$resLogDen[,"obsSparse",] )
#head(resa$resLogDen[,,5])
#plot( as.mcmc.list(resa), smooth=FALSE)
rescoda <- as.mcmc.list(thin(resa,start=.nGen%/%2))
#
#i=1
(.popmean <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"Mean"]}))
(.popsd <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"SD"]}))
# 1/2 orders of magnitude around prescribed sd for theta
# no true sdTheta
#.pop=1
#for( .pop in seq(along.with=.popsd) ){
# checkMagnitude( sdTheta, .popsd[[.pop]] )
#}
#
# check that thetaTrue is in 95% interval
.pthetaTrue <- sapply(1:2, function(.pop){
pnorm(thetaMean, mean=.popmean[[.pop]], sd=.popsd[[.pop]])
})
checkInterval( .pthetaTrue )
resBlock
}
test.ofMultiIntermediateConsistency <- function(){
# testing error handling on altering current state by different cost function to yield nonFininite logDen
data(twTwoDenEx1)
#
thresholdCovar = 0.3 # the true value
#thresholdCovar = 0 # the effective model that glosses over this threshold
thetaMean <- twTwoDenEx1$thetaTrue
sdTheta <- (thetaMean*0.3)
#thetaMean[2] - sdTheta[2] #1.1
denSparsePriorConst <- function(theta,...){
# usually only a updated against this density, make it fail for some values of b
ret <- denSparsePrior(theta,...)
if( theta[2] < 1.5) ret[] <- NA
ret
}
#
argsFLogDen <- list(
thresholdCovar=thresholdCovar
, thetaPrior=thetaMean
, invCovarTheta=sdTheta^2
, twTwoDenEx=twTwoDenEx1
)
#
.nPop=2
.nChainPop=4
ZinitPops <- with(twTwoDenEx1, initZtwDEMCNormal( thetaMean*1.5, diag(sdTheta^2), nChainPop=.nChainPop, nPop=.nPop))
#
.nGen=128
#undebug(denRichPrior)
.remoteDumpfileBasename="tmp/dump_twDEMC_ofMultiIntermediateConsistency"
sfExport("denSparsePrior") # used in denSparsePriorConst
checkException(
resa <- concatPops( resBlock <- twDEMCBlock(
ZinitPops
, nGen=.nGen
,dInfos=list(
dSparse=list(fLogDen=denSparsePriorConst, argsFLogDen=argsFLogDen)
,dRich=list(fLogDen=denRichPrior, argsFLogDen=argsFLogDen)
)
,blocks = list(
a=list( compPos="a",dInfoPos="dSparse", intermediateId="modTwoDenEx")
,b=list(compPos="b", dInfoPos="dRich", intermediateId="modTwoDenEx")
)
,nPop=.nPop
,controlTwDEMC=list(thin=8)
#,debugSequential=TRUE
,remoteDumpfileBasename = .remoteDumpfileBasename
,doRecordProposals=TRUE
))
)
#
denRichPriorConst <- function(theta,...){
# make denRich also return -Inf so that those states are not accepted
ret <- denRichPrior(theta,...)
if( theta[2] < 1.5) ret[] <- NA
ret
}
#
sfExport("denRichPrior") # used in denRichPriorConst
resa <- concatPops( resBlock <- twDEMCBlock(
ZinitPops
, nGen=.nGen
,dInfos=list(
dSparse=list(fLogDen=denSparsePriorConst, argsFLogDen=argsFLogDen)
,dRich=list(fLogDen=denRichPriorConst, argsFLogDen=argsFLogDen)
)
,blocks = list(
a=list(compPos="a", dInfoPos="dSparse", intermediateId="modTwoDenEx")
,b=list(compPos="b", dInfoPos="dRich", intermediateId="modTwoDenEx")
)
,nPop=.nPop
,controlTwDEMC=list(thin=8)
#,debugSequential=TRUE
,remoteDumpfileBasename = .remoteDumpfileBasename
,doRecordProposals=TRUE
))
#plot( as.mcmc.list(resa), smooth=FALSE)
rescoda <- as.mcmc.list(thin(resa,start=.nGen%/%2))
#
#i=1
(.popmean <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"Mean"]}))
(.popsd <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"SD"]}))
# 1/2 orders of magnitude around prescribed sd for theta
# no true sdTheta
#.pop=1
#for( .pop in seq(along.with=.popsd) ){
# checkMagnitude( sdTheta, .popsd[[.pop]] )
#}
#
# check that thetaTrue is in 95% interval
.pthetaTrue <- sapply(1:2, function(.pop){
pnorm(thetaMean, mean=.popmean[[.pop]], sd=.popsd[[.pop]])
})
checkInterval( .pthetaTrue )
}
.tmp.f <- function(){
.remoteDumpfile <- paste(".rda",sep="")
.remoteDumpfile <- paste(.remoteDumpfileBasename,".rda",sep="")
load(.remoteDumpfile)
debugger(get(.remoteDumpfileBasename))
}
t_est.doAppendPrevLogDen <- function(){
# same as goodStartSeq, but with extending former runs
}
.test.dump <- function(){
# test of creating a dump file on remote error
# same as goodStartPrior with changed fModel
# problem: fModel is evaluated in sequential mode before going to parallel execution, so cant test
.nPop=2
argsFLogDen <- list(
# fModel={
# iInvoc=0
# function(parms){
# iInvoc <<- iInvoc+1
# if( iInvoc %% 10 == 0 )
# stop(paste("test throwing an error, iInvoc=",iInvoc,sep=""))
# c(iInvoc=iInvoc)
# }
# }, ### throws an error every 10 invocations using a closure
fModel = function(parms){ stop(paste("test throwing an error",sep="")) },
obs=obs, ### vector of data to compare with
invCovar=invCovar, ### the inverse of the Covariance of obs (its uncertainty)
thetaPrior = thetaTrue, ### the prior estimate of the parameters
invCovarTheta = invCovarTheta ### the inverse of the Covariance of the prior parameter estimates
#xval=xval
)
#argsFLogDen$fModel()
# sfClusterApply(1, argsFLogDen$fModel ) # does not work with closures
#do.call( logDenGaussian, c(list(theta=theta0),argsFLogDen))
.nChainPop=4
Zinit <- initZtwDEMCNormal( theta0, diag(sdTheta^2), nChainPop=.nChainPop, nPop=.nPop)
#dim(Zinit)
body <- expression(
res <- concatPops(resBlock <- twDEMCBlock( Zinit, nGen=100
,dInfos=list(list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen))
,nPop=.nPop
#fLogDenScale=-1/2,
#debugSequential=TRUE
))
)
checkException( eval(body))
suppressWarnings(dir.create("tmp"))
.remoteDumpfileBasename=file.path("tmp","testDump")
.remoteDumpfile <- paste(.remoteDumpfileBasename,".rda",sep="")
# dump on initial parallel calculation, file is created in sequentail mode too
unlink(.remoteDumpfile)
checkTrue( !file.exists(.remoteDumpfile))
#mtrace(twDEMCBlockInt)
isParallel <- sfParallel()
if( !isParallel ) sfInit(TRUE,2) # remote dumpfile only created in parallel mode
body <- expression(
res <- concatPops(resBlock <- twDEMCBlock( Zinit, nGen=100
,dInfos=list(list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen))
,nPop=.nPop
,remoteDumpfileBasename=.remoteDumpfileBasename
#fLogDenScale=-1/2,
#debugSequential=TRUE
))
)
eval(body)
checkException( eval(body))
if( !isParallel ) sfStop()
checkTrue( file.exists(.remoteDumpfile))
# deprecated: dump on Metropolis step (provided initial X and logDenX)
# because fLogDen is evaluated at least once in twDEMC before invoking .doDEMCSteps
.tmp.f <- function(){
unlink(.remoteDumpfile)
checkTrue( !file.exists(.remoteDumpfile))
#mtrace(.doDEMCStep)
body <- expression(
res <- concatPops(resBlock <- twDEMCBlock( Zinit, nGen=100,
fLogDen=logDenGaussian, argsFLogDen=argsFLogDen,
nPop=.nPop,
remoteDumpfileBasename=.remoteDumpfileBasename,
#fLogDenScale=-1/2,
debugSequential=TRUE
,X = matrix(1,nrow=1,ncol=.nChains)
,logDenX = 5
,logDenCompX = character(0)
))
)
checkException( eval(body))
checkTrue( file.exists(.remoteDumpfile))
# dump on Metropolis step (provided initial X and logDenX)
unlink(.remoteDumpfile)
checkTrue( !file.exists(.remoteDumpfile))
#mtrace(.doDEMCStep)
#options(error=recover)
body <- expression({
#mtrace(twCalcLogDenPar);
res <- concatPops(resBlock <- twDEMCBlock( Zinit, nGen=100
,fLogDen=logDenGaussian, argsFLogDen=argsFLogDen
,nPop=.nPop
,remoteDumpfileBasename=.remoteDumpfileBasename
#fLogDenScale=-1/2,
#debugSequential=TRUE,
,X = matrix(1,nrow=1,ncol=.nChains)
))
})
checkException( eval(body))
checkTrue( file.exists(.remoteDumpfile))
}
.tmp.f <- function(){
load( .remoteDumpfile )
debugger(get(.remoteDumpfileBasename))
#choose last column (18)
#interactively execute commands from body of doDEMCStep
#require(debug)
mtrace(remoteFun)
do.call(remoteFun, c(remoteFunArgs, list(...)))
}
}
test.twoStepMetropolis <- function(){
# internal Metropolis step
# nice starting values
.nPop=2
argsFLogDen <- list(
fModel=dummyTwDEMCModel, ### the model function, which predicts the output based on theta
obs=obs, ### vector of data to compare with
invCovar=invCovar/1e10, ### do not constrain by data, the inverse of the Covariance of obs (its uncertainty)
thetaPrior = thetaTrue, ### the prior estimate of the parameters
invCovarTheta = invCovarTheta, ### the inverse of the Covariance of the prior parameter estimates
xval=xval
)
do.call( logDenGaussian, c(list(theta=theta0),argsFLogDen))
.nChainPop=4
Zinit <- initZtwDEMCNormal( theta0, diag(sdTheta^2), nChainPop=.nChainPop, nPop=.nPop)
#dim(Zinit)
X <- adrop(Zinit[ncol(Zinit),, ,drop=FALSE],1) #last row of Zinit
#mtrace(twCalcLogDenPar)
#mtrace(logDenGaussian)
XLogDen <- twCalcLogDenPar(fLogDen=logDenGaussian, xProp=t(X), logDenCompX=numeric(0), intResComp="parms", argsFLogDen=argsFLogDen
,debugSequential=TRUE)
.nGen=100
.thin=5
#mtrace(sfRemoteWrapper)
#mtrace(twDEMCBlockInt)
#mtrace(.doDEMCSteps )
#mtrace(logDenGaussian)
#mtrace(.doDEMCStep)
resa <- concatPops(resBlock <- twDEMCBlock( Zinit, nGen=.nGen
,dInfos=list(list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen
,intResComp="parms"))
,nPop=.nPop
,X = X, logDenCompX=t(XLogDen$logDenComp)
,controlTwDEMC = list(thin=.thin)
,debugSequential=TRUE
,doRecordProposals=TRUE
))
#str(resa)
rescoda <- as.mcmc.list(resa)
plot(rescoda, smooth=FALSE)
res <- thin(resa, start=50)
rescoda <- as.mcmc.list(res)
plot(rescoda, smooth=FALSE)
#gelman.diag(rescoda)
#summary(rescoda)
#str(summary(rescoda))
suppressWarnings({ #glm fit in summary
summary(rescoda)$statistics[,"Mean"]
summary(rescoda)$statistics[,"SD"]
thetaTrue
sdTheta
(.popmean <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"Mean"]}))
(.popsd <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"SD"]}))
})
# 1/2 orders of magnitude around prescribed sd for theta
.pop=1
for( .pop in seq(along.with=.popsd) ){
checkMagnitude( sdTheta, .popsd[[.pop]] )
}
# check that thetaTrue is in 95% interval
.pthetaTrue <- sapply(1:2, function(.pop){
pnorm(thetaTrue, mean=.popmean[[.pop]], sd=.popsd[[.pop]])
})
checkInterval( .pthetaTrue )
}
test.twoStepMetropolisTemp <- function(){
# same as goodStartPrior with internal metropolis step
.nPop=2
argsFLogDen <- list(
fModel=dummyTwDEMCModel, ### the model function, which predicts the output based on theta
obs=obs, ### vector of data to compare with
invCovar=invCovar/1e10, ### do not constrain by data, the inverse of the Covariance of obs (its uncertainty)
thetaPrior = thetaTrue, ### the prior estimate of the parameters
invCovarTheta = invCovarTheta, ### the inverse of the Covariance of the prior parameter estimates
xval=xval
)
do.call( logDenGaussian, c(list(theta=theta0),argsFLogDen))
.nChainPop=4
Zinit <- initZtwDEMCNormal( theta0, diag(sdTheta^2), nChainPop=.nChainPop, nPop=.nPop)
#dim(Zinit)
X <- adrop(Zinit[nrow(Zinit),, ,drop=FALSE],1) #last row of Zinit
#mtrace(twCalcLogDenPar)
#mtrace(logDenGaussian)
XLogDen <- twCalcLogDenPar(fLogDen=logDenGaussian, xProp=t(X), logDenCompX=numeric(0), intResComp="parms", argsFLogDen=argsFLogDen ,debugSequential=TRUE)
.nGen=100
.thin=5
#mtrace(sfRemoteWrapper)
#mtrace(twDEMCBlockInt)
#mtrace(logDenGaussian) #check metropolisStepTemp
#trace(logDenGaussian,quote(recover))
resa <- concatPops(resBlock <- twDEMCBlock( Zinit, nGen=.nGen
,dInfos=list(list(
fLogDen=logDenGaussian, argsFLogDen=argsFLogDen
,intResComp="parms" ))
,nPop=.nPop
,X = X, logDenCompX=t(XLogDen$logDenComp)
,debugSequential=TRUE
,TSpec = cbind( T0=c(obs=10,parms=1), TEnd=1/10)
,controlTwDEMC = list(thin=.thin)
))
#str(resa)
matplot(resa$temp, type="l")
checkTrue( all(resa$temp[1:10,1] > 1) )
rescoda <- as.mcmc.list(resa)
plot(rescoda, smooth=FALSE) #decreasing amplitude
res <- thin(resa, start=60)
rescoda <- as.mcmc.list(res)
plot(rescoda, smooth=FALSE)
#gelman.diag(rescoda)
#summary(rescoda)
#str(summary(rescoda))
suppressWarnings({ #glm fit in summary
summary(rescoda)$statistics[,"Mean"]
summary(rescoda)$statistics[,"SD"]
thetaTrue
sdTheta
(.popmean <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"Mean"]}))
(.popsd <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"SD"]}))
})
# 1/2 orders of magnitude around prescribed sd for theta
.pop=1
for( .pop in seq(along.with=.popsd) ){
checkMagnitude( sdTheta, .popsd[[.pop]] )
}
# check that thetaTrue is in 95% interval
.pthetaTrue <- sapply(1:2, function(.pop){
pnorm(thetaTrue, mean=.popmean[[.pop]], sd=.popsd[[.pop]])
})
checkInterval( .pthetaTrue )
}
t_est.goodStartSeqMultiTemp <- function(){
# same as goodStartSeq but with decreasing temperature and multiTemp
.nPop=2
argsFLogDen <- list(
fModel=dummyTwDEMCModel, ### the model function, which predicts the output based on theta
obs=obs, ### vector of data to compare with
invCovar=invCovar/1e10, ### do not constrain by data, the inverse of the Covariance of obs (its uncertainty)
thetaPrior = thetaTrue, ### the prior estimate of the parameters
invCovarTheta = invCovarTheta, ### the inverse of the Covariance of the prior parameter estimates
xval=xval
)
do.call( logDenGaussian, c(list(theta=theta0),argsFLogDen))
Zinit <- initZtwDEMCNormal( theta0, diag(sdTheta^2), nChainPop=4, nPop=.nPop)
#dim(Zinit)
.nGen=100
#mtrace(sfRemoteWrapper)
#mtrace(.doDEMCStep)
#mtrace(.doDEMCSteps )
#mtrace(logDenGaussian)
#mtrace(twCalcLogDenPar)
#mtrace(twDEMCBlockInt)
res <- concatPops(resBlock <- twDEMCBlock( Zinit, nGen=.nGen,
fLogDen=logDenGaussian, argsFLogDen=argsFLogDen,
nPop=.nPop,
#controlTwDEMC=list(thin=1, T0=20, TEnd=20, useMultiT=TRUE, TFix=c(parms=1)),
controlTwDEMC=list(thin=1, T0=20, TEnd=20, useMultiT=TRUE, Tprop=c(obs=0.5,parms=1), TFix=c(parms=1)),
debugSequential=TRUE
))
#str(res)
checkEquals((.nGen%/%res$thin)+1,nrow(res$logDen))
rescoda <- as.mcmc.list(res)
plot(rescoda)
#gelman.diag(rescoda)
#summary(rescoda)
#str(summary(rescoda))
suppressWarnings({ #glm fit in summary
summary(rescoda)$statistics[,"Mean"]
summary(rescoda)$statistics[,"SD"]
thetaTrue
sdTheta
(.popmean <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"Mean"]}))
(.popsd <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"SD"]}))
})
# 1/2 orders of magnitude around prescribed sd for theta
.pop=1
for( .pop in seq(along.with=.popsd) ){
checkMagnitude( sdTheta, .popsd[[.pop]] )
}
# check that thetaTrue is in 95% interval
.pthetaTrue <- sapply(1:2, function(.pop){
pnorm(thetaTrue, mean=.popmean[[.pop]], sd=.popsd[[.pop]])
})
checkInterval( .pthetaTrue )
}
.test.fLogDenScale <- function(){
# same as goodStart sequence but with logDens function returning cost instead of logDensity
.nPop=2
argsFLogDen <- list(
fModel=dummyTwDEMCModel, ### the model function, which predicts the output based on theta
obs=obs, ### vector of data to compare with
invCovar=invCovar/1e10, ### do not constrain by data, the inverse of the Covariance of obs (its uncertainty)
thetaPrior = thetaTrue, ### the prior estimate of the parameters
invCovarTheta = invCovarTheta, ### the inverse of the Covariance of the prior parameter estimates
xval=xval
,scale=1 ##<< override the -1/2
)
do.call( logDenGaussian, c(list(theta=theta0),argsFLogDen))
Zinit <- initZtwDEMCNormal( theta0, diag(sdTheta^2), nChainPop=4, nPop=.nPop)
#dim(Zinit)
.nGen=100
#mtrace(sfRemoteWrapper)
#mtrace(.doDEMCStep)
#mtrace(.doDEMCSteps )
#mtrace(logDenGaussian)
#mtrace(twCalcLogDenPar)
#mtrace(twDEMCBlockInt)
res <- concatPops(resBlock <- twDEMCBlock( Zinit, nGen=.nGen,
dInfos=list(list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen
,fLogDenScale=-1/2 # here apply the scale in twDEMCBlock
))
,nPop=.nPop
,controlTwDEMC=list(thin=1)
,debugSequential=TRUE
))
#str(res)
checkEquals((.nGen%/%res$thin)+1,nrow(res$logDen))
#windows(record=TRUE)
rescoda <- as.mcmc.list(res)
plot(rescoda)
#gelman.diag(rescoda)
#summary(rescoda)
#str(summary(rescoda))
suppressWarnings({ #glm fit in summary
summary(rescoda)$statistics[,"Mean"]
summary(rescoda)$statistics[,"SD"]
thetaTrue
sdTheta
(.popmean <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"Mean"]}))
(.popsd <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"SD"]}))
})
# 1/2 orders of magnitude around prescribed sd for theta
.pop=1
for( .pop in seq(along.with=.popsd) ){
checkMagnitude( sdTheta, .popsd[[.pop]] )
}
# check that thetaTrue is in 95% interval
.pthetaTrue <- sapply(1:2, function(.pop){
pnorm(thetaTrue, mean=.popmean[[.pop]], sd=.popsd[[.pop]])
})
checkInterval( .pthetaTrue )
}
Try the twDEMC package in your browser
Any scripts or data that you put into this service are public.
twDEMC documentation built on May 2, 2019, 5:38 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## Test unit 'twDEMC'
# twUtestF(twDEMC)
# twUtestF(twDEMC, "test.goodStartSeqData")
# twUtestF(twDEMC, "test_int.goodStartSeqData.plot2d")
# twUtestF(twDEMC, "dump")
.setUp <- function(){
data(twLinreg1)
data(twdemcEx1)
suppressWarnings( rm( list=names(twLinreg1) ) )
attach(twLinreg1)
ex1c <- concatPops(twdemcEx1)
suppressWarnings( rm( list=names(ex1c) ) )
#attach(ex1c)
}
.tearDown <- function(){
#detach( ex1c )
detach( twLinreg1 )
}
.tmp.f <- function(){
trace( updateBlockTwDEMC, recover ) # untrace(updateBlockTwDEMC)
mtrace(logDenGaussian)
mtrace(twCalcLogDenPar)
mtrace(.doDEMCStep)
mtrace(.doDEMCSteps)
mtrace(.generateXProb)
mtrace(.generateXPropChains)
mtrace(.xStepSnooker)
mtrace(twDEMCBlockInt)
mtrace(.generateXPropThin)
mtrace(.xStepParallel)
twUtestF(twDEMC,"test.probUpDir")
twUtestF(twDEMC,"test.ZinittwDEMC") #with thinning interval
twUtestF(twDEMC,"test.goodStartSeq")
twUtestF(twDEMC,"test.goodStartPar")
twUtestF(twDEMC,"test.doAppendPrevLogDen")
twUtestF(twDEMC,"test.goodStartSeqData")
twUtestF(twDEMC)
twUtestF(twDEMCBatch)
twUtestF()
prevWd <- setwd(file.path("tests"))
dumpFileBasename="doRUnitError"
#options(error=quote(dump.frames(dumpFileBasename, TRUE)))
load(paste(dumpFileBasename,".rda"))
debugger(dumpFileBasename)
setwd(prevWd)
}
test.goodStartSeqData <- function(){
#mtrace(test_int.goodStartSeqData.plot2d)
test_int.goodStartSeqData.plot2d()
}
test_int.goodStartSeqData.plot2d <- function(){
lmDummy <- lm( obs ~ xval, weights=1/sdObs^2) # results without priors
(.expTheta <- coef(lmDummy))
(.expCovTheta <- vcov(lmDummy)) # a is very weak constrained, negative covariance between a an b
(.expSdTheta <- structure(sqrt(diag(.expCovTheta)), names=c("a","b")) )
# nice starting values
.nPop=2
argsFLogDen <- list(
fModel=dummyTwDEMCModel, ### the model function, which predicts the output based on theta
obs=obs, ### vector of data to compare with
invCovar=invCovar, ### do not constrain by data, the inverse of the Covariance of obs (its uncertainty)
#thetaPrior = thetaTrue, ### the prior estimate of the parameters
#invCovarTheta = invCovarTheta, ### the inverse of the Covariance of the prior parameter estimates
xval=xval
)
do.call( logDenGaussian, c(list(theta=theta0),argsFLogDen))
Zinit <- initZtwDEMCNormal( theta0, .expCovTheta, nChainPop=4, nPop=.nPop)
#dim(Zinit)
.nGen=100
#mtrace(logDenGaussian)
#mtrace(twDEMCBlockInt)
res <- concatPops( resBlock <- twDEMCBlock( Zinit, nGen=.nGen,
dInfos=list(list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen)),
nPop=.nPop
,controlTwDEMC=list(thin=8 )
#,debugSequential=TRUE
,doRecordProposals=TRUE
))
#str(res)
checkEquals( 8, res$thin )
checkEquals( (.nGen%/%res$thin)+1,nrow(res$logDen))
rescoda <- as.mcmc.list(res)
plot(rescoda, smooth=FALSE)
#gelman.diag(rescoda)
#summary(rescoda)
.tmp.f <- function(){
matplot( res$logDen[,1,],type="l" )
tmp <- data.frame( a=as.numeric(res$Y[,"a",]), b=as.numeric(res$Y[,"b",]), logDen=as.numeric(rowSums(res$Y[,res$YPos$resLogDen0+res$dInfos[[1]]$compPosDen,])) )
colorFun <- colorRampPalette(c("yellow","red"))
levelplot( logDen ~ a + b, tmp, col.regions = colorFun(50))
tmp2 <- data.frame( a=as.numeric(res$parms[,"a",]), b=as.numeric(res$parms[,"b",]), logDen=as.numeric(res$logDen) )
levelplot( logDen ~ a + b, tmp2, col.regions = colorFun(50))
#image( as.numeric(res$Y["a",,]), as.numeric(res$Y["b",,]), as.numeric(res$Y["logDen",,]))
}
#str(summary(rescoda))
suppressWarnings({
summary(rescoda)$statistics[,"Mean"]
summary(rescoda)$statistics[,"SD"]
thetaTrue
sdTheta
(.popmean <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"Mean"]}))
(.popsd <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"SD"]}))
})
# 1/2 orders of magnitude around prescribed sd for theta
.pop=1
for( .pop in seq(along.with=.popsd) ){
checkMagnitude( .expSdTheta, .popsd[[.pop]] )
}
# check that thetaTrue is in 95% interval
.pthetaTrue <- sapply(1:2, function(.pop){
pnorm(.expTheta, mean=.popmean[[.pop]], sd=.popsd[[.pop]])
})
checkInterval( .pthetaTrue )
# repeat with control for overfitting
resCO <- concatPops( resBlockCO <- twDEMCBlock( Zinit, nGen=.nGen,
dInfos=list( list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen, maxLogDen=c(obs=-length(obs)/2,parms=0)) ),
nPop=.nPop
,controlTwDEMC=list(thin=8)
#,debugSequential=TRUE
,doRecordProposals=TRUE
))
rescoda <- as.mcmc.list(resCO)
plot(rescoda, smooth=FALSE)
#str(summary(rescoda))
suppressWarnings({
summary(rescoda)$statistics[,"Mean"]
summary(rescoda)$statistics[,"SD"]
thetaTrue
sdTheta
(.popmean <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"Mean"]}))
(.popsd <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"SD"]}))
})
# 1/2 orders of magnitude around prescribed sd for theta
.pop=1
for( .pop in seq(along.with=.popsd) ){
checkMagnitude( .expSdTheta, .popsd[[.pop]] )
}
# check that thetaTrue is in 95% interval
.pthetaTrue <- sapply(1:2, function(.pop){
pnorm(.expTheta, mean=.popmean[[.pop]], sd=.popsd[[.pop]])
})
checkInterval( .pthetaTrue )
}
test.badStartSeqData <- function(){
lmDummy <- lm( obs ~ xval, weights=1/sdObs^2) # results without priors
(.expTheta <- coef(lmDummy))
(.expCovTheta <- vcov(lmDummy)) # a is very weak constrained, negative covariance between a an b
(.expSdTheta <- structure(sqrt(diag(.expCovTheta)), names=c("a","b")) )
.nPop=2
argsFLogDen <- list(
fModel=dummyTwDEMCModel, ### the model function, which predicts the output based on theta
obs=obs, ### vector of data to compare with
invCovar=invCovar, ### do not constrain by data, the inverse of the Covariance of obs (its uncertainty)
#thetaPrior = thetaTrue, ### the prior estimate of the parameters
#invCovarTheta = invCovarTheta, ### the inverse of the Covariance of the prior parameter estimates
xval=xval
)
do.call( logDenGaussian, c(list(theta=theta0),argsFLogDen))
#Zinit <- initZtwDEMCNormal( theta0, .expCovTheta, nChainPop=4, nPop=.nPop)
.sdThetaBad <- sdTheta*c(1/10,10)
.theta0Bad <- theta0*c(10,1/10)
Zinit <- initZtwDEMCNormal( .theta0Bad, diag(.sdThetaBad^2), nChainPop=4, nPop=.nPop)
Zinit[,,1:4] <- Zinit[,,1:4] * c(2,1)
Zinit[,,-(1:4)] <- Zinit[,,-(1:4)] * c(1,2)
#dim(Zinit)
.nGen=200
resa <- concatPops( resBlock <- twDEMCBlock( Zinit, nGen=.nGen,
dInfos=list(list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen)),
nPop=.nPop,
controlTwDEMC=list(thin=1),
debugSequential=TRUE
))
checkEquals((.nGen%/%resa$thin)+1,nrow(resa$logDen))
rescoda <- as.mcmc.list(resa)
plot(rescoda, smooth=FALSE)
#mtrace(thin.twDEMC)
res <- thin(resa, start=120)
rescoda <- as.mcmc.list(res)
plot(rescoda, smooth=FALSE)
#gelman.diag(rescoda)
#summary(rescoda)
#str(summary(rescoda))
suppressWarnings({ #glm fit in summary
summary(rescoda)$statistics[,"Mean"]
summary(rescoda)$statistics[,"SD"]
thetaTrue
sdTheta
(.popmean <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"Mean"]}))
(.popsd <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"SD"]}))
})
# 1/2 orders of magnitude around prescribed sd for theta
.pop=1
for( .pop in seq(along.with=.popsd) ){
checkMagnitude( .expSdTheta, .popsd[[.pop]] )
}
# check that thetaTrue is in 95% interval
.pthetaTrue <- sapply(1:2, function(.pop){
pnorm(.expTheta, mean=.popmean[[.pop]], sd=.popsd[[.pop]])
})
checkInterval( .pthetaTrue )
}
test.badStartSeqData1D <- function(){
# common error is dropping dimensions, if of size 1 onyle: test parameter_dimension = 1
set.seed(0815)
lmDummy <- lm( obs-theta0["a"] ~ xval-1, weights=1/sdObs^2) # results without priors
theta1d <- theta0["b"]
(.expTheta <- coef(lmDummy))
(.expCovTheta <- vcov(lmDummy))
(.expSdTheta <- structure(sqrt(diag(.expCovTheta)), names=c("b")) )
#trace( dummyTwDEMCModel, recover ) # untrace(dummyTwDEMCModel)
argsFLogDen <- list(
fModel=dummyTwDEMCModel, ### the model function, which predicts the output based on theta
obs=obs, ### vector of data to compare with
invCovar=invCovar, ### do not constrain by data, the inverse of the Covariance of obs (its uncertainty)
#thetaPrior = thetaTrue, ### the prior estimate of the parameters
invCovarTheta = invCovarTheta[2,2,drop=FALSE], ### the inverse of the Covariance of the prior parameter estimates
xval=xval
,theta0=theta0
)
do.call( logDenGaussian, c(list(theta=theta1d),argsFLogDen))
.tmp.f <- function(){
bs <- seq(2,8,length.out=30)
ll <- colSums( sapply( bs, function(b){ do.call( logDenGaussian, c(list(theta=c(b=b)),argsFLogDen)) }) )
plot( ll ~ bs )
}
#Zinit <- initZtwDEMCNormal( theta0, .expCovTheta, nChainPop=4, nPop=.nPop)
.sdThetaBad <- sdTheta["b"]*c(1/10)
.theta0Bad <- theta0["b"]*c(10)
.nPop=2
Zinit <- initZtwDEMCNormal( .theta0Bad, diag(.sdThetaBad^2,nrow = length(.sdThetaBad)), nChainPop=4, nPop=.nPop)
Zinit[,,1:4] <- Zinit[,,1:4] * c(2)
Zinit[,,-(1:4)] <- Zinit[,,-(1:4)] * c(0.5)
#dim(Zinit)
.nGen=250
#mtrace(twDEMCBlockInt)
#mtrace(.sampleStates)
#mtrace(twWhichColsEqual)
resa <- concatPops(resBlock <- twDEMCBlock( Zinit, nGen=.nGen,
dInfos=list(list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen)),
nPop=.nPop,
controlTwDEMC=list(thin=1)
,debugSequential=TRUE
))
#str(res)
plot( as.mcmc.list(resa), smooth=FALSE )
checkEquals((.nGen%/%resa$thin)+1,nrow(resa$logDen))
matplot( resa$parms[,1,],type="l")
res <- thin(resa, start=110)
matplot( res$parms[,1,],type="l")
#gelman.diag(rescoda)
#summary(rescoda)
#str(summary(rescoda))
suppressWarnings({ #glm fit in summary
(.popmean <- lapply(list(p1=1:4,p2=5:8),function(i){mean(res$parms[1,,i])}))
(.popsd <- lapply(list(p1=1:4,p2=5:8),function(i){sd(as.vector(res$parms[1,,i]))}))
})
# 1/2 orders of magnitude around prescribed sd for theta
.pop=1
for( .pop in seq(along.with=.popsd) ){
checkMagnitude( .expSdTheta, .popsd[[.pop]], msg="wrong magnitude of standard deviation" )
}
# check that thetaTrue is in 95% interval
.pthetaTrue <- sapply(1:2, function(.pop){
pnorm(.expTheta, mean=.popmean[[.pop]], sd=.popsd[[.pop]])
})
checkInterval( .pthetaTrue, 0.005, 0.995, "thetaTrue outside 0.99% confidence interval." )
}
test.goodStartPrior <- function(){
# nice starting values and prior constrains estimates
.nPop=2
argsFLogDen <- list(
fModel=dummyTwDEMCModel, ### the model function, which predicts the output based on theta
obs=obs, ### vector of data to compare with
invCovar=invCovar/1e10, ### do not constrain by data, the inverse of the Covariance of obs (its uncertainty)
thetaPrior = thetaTrue, ### the prior estimate of the parameters
invCovarTheta = invCovarTheta, ### the inverse of the Covariance of the prior parameter estimates
xval=xval
)
do.call( logDenGaussian, c(list(theta=theta0),argsFLogDen))
Zinit <- initZtwDEMCNormal( theta0, diag(sdTheta^2), nChainPop=4, nPop=.nPop)
#dim(Zinit)
.nGen=100
#mtrace(updateBlockTwDEMC)
#mtrace(twDEMCBlockInt)
#mtrace(.updateIntervalTwDEMCPar)
#mtrace(.updateBlocksTwDEMC)
resa <- res <- resSeq <- concatPops(resBlock <- twDEMCBlock( Zinit, nGen=.nGen,
dInfos=list(list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen)),
nPop=.nPop,
controlTwDEMC=list(thin=4, useMultiT=TRUE)
,debugSequential=TRUE
#,T0=10, TProp=c( obs=0.5, parms=1 )
))
#str(res)
checkEquals((.nGen%/%resa$thin)+1,nrow(resa$logDen))
#windows(record=TRUE)
rescoda <- as.mcmc.list(resa)
plot(rescoda, smooth=FALSE)
#gelman.diag(rescoda)
#summary(rescoda)
res <- thin(resa, start=30)
rescoda <- as.mcmc.list(res)
plot(rescoda, smooth=FALSE)
.tmp.f <- function(){
colnames(res$resLogDen)
matplot(res$resLogDen[,"obs",], type="l")
matplot(res$resLogDen[,"parms",], type="l")
matplot(res$logDen[,1,], type="l")
tail(res$Y[,,1])
}
#str(summary(rescoda))
suppressWarnings({ #glm fit in summary
summary(rescoda)$statistics[,"Mean"]
summary(rescoda)$statistics[,"SD"]
thetaTrue
sdTheta
(.popmean <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"Mean"]}))
(.popsd <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"SD"]}))
})
# 1/2 orders of magnitude around prescribed sd for theta
.pop=1
for( .pop in seq(along.with=.popsd) ){
checkMagnitude( sdTheta, .popsd[[.pop]] )
}
# check that thetaTrue is in 95% interval
.pthetaTrue <- sapply(1:2, function(.pop){
pnorm(thetaTrue, mean=.popmean[[.pop]], sd=.popsd[[.pop]])
})
checkInterval( .pthetaTrue ) #
#------------ again with parallel runs
resa <- resPar <- concatPops(resBlock <- twDEMCBlock( Zinit, nGen=100,
dInfos=list(list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen)),
nPop=.nPop
#fLogDenScale=-1/2
#debugSequential=TRUE
))
rescoda <- as.mcmc.list(thin(resa,start=30))
#str(summary(rescoda))
suppressWarnings({ #glm fit in summary
summary(rescoda)$statistics[,"Mean"]
summary(rescoda)$statistics[,"SD"]
thetaTrue
sdTheta
(.popmean <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"Mean"]}))
(.popsd <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"SD"]}))
})
# 1/2 orders of magnitude around prescribed sd for theta
.pop=1
for( .pop in seq(along.with=.popsd) ){
checkMagnitude( sdTheta, .popsd[[.pop]] )
}
# check that thetaTrue is in 95% interval
.pthetaTrue <- sapply(1:2, function(.pop){
pnorm(thetaTrue, mean=.popmean[[.pop]], sd=.popsd[[.pop]])
})
checkInterval( .pthetaTrue ) #
}
test.upperParBounds <- function(){
# same as goodStartPar, but giving an upper parameter bound
.nPop=2
argsFLogDen <- list(
fModel=dummyTwDEMCModel, ### the model function, which predicts the output based on theta
obs=obs, ### vector of data to compare with
invCovar=invCovar, ### the inverse of the Covariance of obs (its uncertainty)
thetaPrior = thetaTrue, ### the prior estimate of the parameters
invCovarTheta = invCovarTheta, ### the inverse of the Covariance of the prior parameter estimates
xval=xval
)
#do.call( logDenGaussian, c(list(theta=theta0),argsFLogDen))
asplit <- 10.0
Zinit0 <- initZtwDEMCNormal( theta0, diag(sdTheta^2), nChainPop=4, nPop=.nPop)
# replace all the cases with upperParBounds
bo.keep <- as.vector(Zinit0[,"a",]) <= asplit
#mtrace(replaceZinitCases)
Zinit1 <- replaceZinitCases(Zinit0, bo.keep)
Zinit2 <- replaceZinitCases(Zinit0, !bo.keep)
dim(Zinit1)
#mtrace(.doDEMCStep)
res <- res0 <- concatPops(resBlock <- twDEMCBlock( Zinit0, nGen=100,
dInfos=list(list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen)),
nPop=.nPop
#,upperParBounds=c(a=asplit)
#,debugSequential=TRUE
#,controlTwDEMC=list( DRgamma=0.1, minPCompAcceptTempDecr=0.16)
))
res <- res1 <- concatPops(resBlock <- twDEMCBlock( Zinit1, nGen=64*4,
dInfos=list(list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen)),
nPop=.nPop
,upperParBounds=c(a=asplit)
,debugSequential=TRUE
#,controlTwDEMC=list( DRgamma=0.1, minPCompAcceptTempDecr=0.16)
))
res <- res2 <- concatPops(resBlock <- twDEMCBlock( Zinit2, nGen=64*4,
dInfos=list(list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen)),
nPop=.nPop
,lowerParBounds=c(a=asplit)
,debugSequential=TRUE
))
ss0 <- stackChains(res0)
# importance sampling due to integrated probability
ss1 <- stackChains(res1)
lw1 <- twLogSumExp(ss1[,1])
ss2 <- stackChains(res2)
lw2 <- twLogSumExp(ss2[,1])
lwSum <- twLogSumExp( c(lw1,lw2) )
pSubs <- c( exp(lw1-lwSum), exp(lw2-lwSum) )
nSubs <- floor(pSubs/max(pSubs) * nrow(ss0))
ssc <- rbind( ss1[sample.int(nrow(ss1),nSubs[1]),]
, ss2[sample.int(nrow(ss2),nSubs[2]),]
)
.tmp.f <- function(){
str(res)
rescoda <- as.mcmc.list(res)
plot(rescoda, smooth=FALSE)
plot( density(ssc[,"a"])) # combined
lines(density(ss0[,"a"]), col="blue") # orig
lines(density(ss1[,"a"]), col="green") # left
lines(density(ss2[,"a"]), col="darkgreen") # right
gelman.diag(res1)
gelman.diag(res2)
}
checkTrue( all( ss1[,"a"] <= asplit ) )
checkTrue( all( ss2[,"a"] >= asplit ) )
.popmean <- colMeans( ssc[,-1])
.popsd <- apply(ssc[,-1], 2, sd )
# 1/2 orders of magnitude around prescribed sd for theta
.pop=1
for( .pop in seq(along.with=.popsd) ){
checkMagnitude( sdTheta, .popsd[[.pop]] )
}
# check that thetaTrue is in 95% interval
.pthetaTrue <- sapply(1:2, function(.pop){
pnorm(thetaTrue, mean=.popmean[[.pop]], sd=.popsd[[.pop]])
})
checkInterval( .pthetaTrue )
}
t_est.probUpDir <- function(){
#mtrace(inner.probUpDir)
innertest.probUpDir()
#twUtestF(twDEMC,"test.probUpDir")
}
i_nnertest.probUpDir <- function(){
# same as goodStartSeq, but with executing debugSequential=FALSE, i.e. parallel
.nPop=2
argsFLogDen <- list(
fModel=dummyTwDEMCModel, ### the model function, which predicts the output based on theta
obs=obs, ### vector of data to compare with
invCovar=invCovar, ### the inverse of the Covariance of obs (its uncertainty)
thetaPrior = thetaTrue, ### the prior estimate of the parameters
invCovarTheta = invCovarTheta, ### the inverse of the Covariance of the prior parameter estimates
xval=xval
)
#do.call( logDenGaussian, c(list(theta=theta0),argsFLogDen))
#
Zinit <- initZtwDEMCNormal( theta0, diag(sdTheta^2), nChainPop=4, nPop=.nPop)
dim(Zinit)
#
res <- concatPops(resBlock <- twDEMCBlock( Zinit, nGen=100,
dInfos=list(list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen)),
nPop=.nPop,
#fLogDenScale=-1/2,
controlTwDEMC = list(probUpDir=0.8)
))
str(res)
#
rescoda <- as.mcmc.list(res)
suppressWarnings({ #glm fit in summary
(.popmean <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"Mean"]}))
(.popsd <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"SD"]}))
})
#
# 1/2 orders of magnitude around prescribed sd for theta
.pop=1
for( .pop in seq(along.with=.popsd) ){
checkMagnitude( sdTheta, .popsd[[.pop]] )
}
#
# check that thetaTrue is in 95% interval
.pthetaTrue <- sapply(1:2, function(.pop){
pnorm(thetaTrue, mean=.popmean[[.pop]], sd=.popsd[[.pop]])
})
checkInterval( .pthetaTrue )
}
test.ofMulti <- function(){
# testing multiple log-Density functions
data(twTwoDenEx1)
thresholdCovar = 0.3 # the true value
#thresholdCovar = 0 # the effective model that glosses over this threshold
thetaMean <- twTwoDenEx1$thetaTrue
sdTheta <- (thetaMean*0.3)
argsFLogDen <- list(
thresholdCovar=thresholdCovar
, thetaPrior=thetaMean
, invCovarTheta=sdTheta^2
, twTwoDenEx=twTwoDenEx1
)
#
.nPop=2
.nChainPop=4
ZinitPops <- with(twTwoDenEx1, initZtwDEMCNormal( thetaMean*1.5, diag(sdTheta^2), nChainPop=.nChainPop, nPop=.nPop))
#
.nGen=128
resa <- concatPops( resBlock <- twDEMCBlock(
ZinitPops
, nGen=.nGen
,dInfos=list(
dSparse=list(fLogDen=denSparsePrior, argsFLogDen=argsFLogDen)
,dRich=list(fLogDen=denRichPrior, argsFLogDen=argsFLogDen)
)
,blocks = list(
a=list(dInfoPos="dSparse", compPos="a")
,b=list(dInfoPos="dRich", compPos="b")
)
,nPop=.nPop
,controlTwDEMC=list(thin=8)
#,debugSequential=TRUE
,doRecordProposals=TRUE
))
plot( as.mcmc.list(resa), smooth=FALSE)
rescoda <- as.mcmc.list(thin(resa,start=.nGen%/%2))
#
#i=1
.popmean <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"Mean"]})
.popsd <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"SD"]})
# 1/2 orders of magnitude around prescribed sd for theta
# no true sdTheta
#.pop=1
#for( .pop in seq(along.with=.popsd) ){
# checkMagnitude( sdTheta, .popsd[[.pop]] )
#}
#
# check that thetaTrue is in 95% interval
.pthetaTrue <- sapply(1:2, function(.pop){
pnorm(thetaMean, mean=.popmean[[.pop]], sd=.popsd[[.pop]])
})
checkInterval( .pthetaTrue )
}
test.ofMultiIntermediate <- function(){
int.ofMultiIntermediate()
}
test.ofMultiIntermediateCond <- function(){
int.ofMultiIntermediate( useConditionalProposal=TRUE )
}
test.ofMultiIntermediateControlOverfit <- function(){
int.ofMultiIntermediate( controlOverfitting=TRUE )
}
int.ofMultiIntermediate <- function(useConditionalProposal=FALSE, controlOverfitting=FALSE){
# testing multiple log-Density functions using intermediate results of each other (argument intermediate to dInfos)
data(twTwoDenEx1)
thresholdCovar = 0.3 # the true value
#thresholdCovar = 0 # the effective model that glosses over this threshold
thetaMean <- twTwoDenEx1$thetaTrue
sdTheta <- (thetaMean*0.3)
#
argsFLogDen <- list(
thresholdCovar=thresholdCovar
, thetaPrior=thetaMean
, invCovarTheta=sdTheta^2
, twTwoDenEx=twTwoDenEx1
)
.nPop=2
.nChainPop=4
ZinitPops <- with(twTwoDenEx1, initZtwDEMCNormal( thetaMean*1.5, diag(sdTheta^2), nChainPop=.nChainPop, nPop=.nPop))
#
.nGen=128
#undebug(denRichPrior)
.remoteDumpfileBasename="tmp/dump_twDEMC_ofMultiIntermediate"
resa <- concatPops( resBlock <- twDEMCBlock(
ZinitPops
, nGen=.nGen
,dInfos=list(
dSparse=list(fLogDen=denSparsePrior, argsFLogDen=argsFLogDen, maxLogDen=c("parmsSparse"=0, "obsSparse"=-10/2 ))
,dRich=list(fLogDen=denRichPrior, argsFLogDen=argsFLogDen, maxLogDen=c("parmsRich"=0, "obsRich"=-1000/2 ))
)
,blocks = list(
a=list(compPos="a", dInfoPos="dSparse", intermediateId="modTwoDenEx")
,b=list(compPos="b", dInfoPos="dRich", intermediateId="modTwoDenEx")
)
,nPop=.nPop
,controlTwDEMC=list(thin=8, useConditionalProposal=useConditionalProposal, controlOverfittingMinNObs=if(controlOverfitting) 10 else Inf )
#,debugSequential=TRUE
,remoteDumpfileBasename = .remoteDumpfileBasename
,doRecordProposals=TRUE
))
if( controlOverfitting ){
checkTrue( all(resa$resLogDen[,"obsSparse",] <= -5) ) # not across upper bound of logDen
checkTrue( all(resa$resLogDen[,"obsRich",] <= -500) ) # not across upper bound of logDen
}
# summary( resa$resLogDen[,"obsSparse",] )
#head(resa$resLogDen[,,5])
#plot( as.mcmc.list(resa), smooth=FALSE)
rescoda <- as.mcmc.list(thin(resa,start=.nGen%/%2))
#
#i=1
(.popmean <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"Mean"]}))
(.popsd <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"SD"]}))
# 1/2 orders of magnitude around prescribed sd for theta
# no true sdTheta
#.pop=1
#for( .pop in seq(along.with=.popsd) ){
# checkMagnitude( sdTheta, .popsd[[.pop]] )
#}
#
# check that thetaTrue is in 95% interval
.pthetaTrue <- sapply(1:2, function(.pop){
pnorm(thetaMean, mean=.popmean[[.pop]], sd=.popsd[[.pop]])
})
checkInterval( .pthetaTrue )
resBlock
}
test.ofMultiIntermediateConsistency <- function(){
# testing error handling on altering current state by different cost function to yield nonFininite logDen
data(twTwoDenEx1)
#
thresholdCovar = 0.3 # the true value
#thresholdCovar = 0 # the effective model that glosses over this threshold
thetaMean <- twTwoDenEx1$thetaTrue
sdTheta <- (thetaMean*0.3)
#thetaMean[2] - sdTheta[2] #1.1
denSparsePriorConst <- function(theta,...){
# usually only a updated against this density, make it fail for some values of b
ret <- denSparsePrior(theta,...)
if( theta[2] < 1.5) ret[] <- NA
ret
}
#
argsFLogDen <- list(
thresholdCovar=thresholdCovar
, thetaPrior=thetaMean
, invCovarTheta=sdTheta^2
, twTwoDenEx=twTwoDenEx1
)
#
.nPop=2
.nChainPop=4
ZinitPops <- with(twTwoDenEx1, initZtwDEMCNormal( thetaMean*1.5, diag(sdTheta^2), nChainPop=.nChainPop, nPop=.nPop))
#
.nGen=128
#undebug(denRichPrior)
.remoteDumpfileBasename="tmp/dump_twDEMC_ofMultiIntermediateConsistency"
sfExport("denSparsePrior") # used in denSparsePriorConst
checkException(
resa <- concatPops( resBlock <- twDEMCBlock(
ZinitPops
, nGen=.nGen
,dInfos=list(
dSparse=list(fLogDen=denSparsePriorConst, argsFLogDen=argsFLogDen)
,dRich=list(fLogDen=denRichPrior, argsFLogDen=argsFLogDen)
)
,blocks = list(
a=list( compPos="a",dInfoPos="dSparse", intermediateId="modTwoDenEx")
,b=list(compPos="b", dInfoPos="dRich", intermediateId="modTwoDenEx")
)
,nPop=.nPop
,controlTwDEMC=list(thin=8)
#,debugSequential=TRUE
,remoteDumpfileBasename = .remoteDumpfileBasename
,doRecordProposals=TRUE
))
)
#
denRichPriorConst <- function(theta,...){
# make denRich also return -Inf so that those states are not accepted
ret <- denRichPrior(theta,...)
if( theta[2] < 1.5) ret[] <- NA
ret
}
#
sfExport("denRichPrior") # used in denRichPriorConst
resa <- concatPops( resBlock <- twDEMCBlock(
ZinitPops
, nGen=.nGen
,dInfos=list(
dSparse=list(fLogDen=denSparsePriorConst, argsFLogDen=argsFLogDen)
,dRich=list(fLogDen=denRichPriorConst, argsFLogDen=argsFLogDen)
)
,blocks = list(
a=list(compPos="a", dInfoPos="dSparse", intermediateId="modTwoDenEx")
,b=list(compPos="b", dInfoPos="dRich", intermediateId="modTwoDenEx")
)
,nPop=.nPop
,controlTwDEMC=list(thin=8)
#,debugSequential=TRUE
,remoteDumpfileBasename = .remoteDumpfileBasename
,doRecordProposals=TRUE
))
#plot( as.mcmc.list(resa), smooth=FALSE)
rescoda <- as.mcmc.list(thin(resa,start=.nGen%/%2))
#
#i=1
(.popmean <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"Mean"]}))
(.popsd <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"SD"]}))
# 1/2 orders of magnitude around prescribed sd for theta
# no true sdTheta
#.pop=1
#for( .pop in seq(along.with=.popsd) ){
# checkMagnitude( sdTheta, .popsd[[.pop]] )
#}
#
# check that thetaTrue is in 95% interval
.pthetaTrue <- sapply(1:2, function(.pop){
pnorm(thetaMean, mean=.popmean[[.pop]], sd=.popsd[[.pop]])
})
checkInterval( .pthetaTrue )
}
.tmp.f <- function(){
.remoteDumpfile <- paste(".rda",sep="")
.remoteDumpfile <- paste(.remoteDumpfileBasename,".rda",sep="")
load(.remoteDumpfile)
debugger(get(.remoteDumpfileBasename))
}
t_est.doAppendPrevLogDen <- function(){
# same as goodStartSeq, but with extending former runs
}
.test.dump <- function(){
# test of creating a dump file on remote error
# same as goodStartPrior with changed fModel
# problem: fModel is evaluated in sequential mode before going to parallel execution, so cant test
.nPop=2
argsFLogDen <- list(
# fModel={
# iInvoc=0
# function(parms){
# iInvoc <<- iInvoc+1
# if( iInvoc %% 10 == 0 )
# stop(paste("test throwing an error, iInvoc=",iInvoc,sep=""))
# c(iInvoc=iInvoc)
# }
# }, ### throws an error every 10 invocations using a closure
fModel = function(parms){ stop(paste("test throwing an error",sep="")) },
obs=obs, ### vector of data to compare with
invCovar=invCovar, ### the inverse of the Covariance of obs (its uncertainty)
thetaPrior = thetaTrue, ### the prior estimate of the parameters
invCovarTheta = invCovarTheta ### the inverse of the Covariance of the prior parameter estimates
#xval=xval
)
#argsFLogDen$fModel()
# sfClusterApply(1, argsFLogDen$fModel ) # does not work with closures
#do.call( logDenGaussian, c(list(theta=theta0),argsFLogDen))
.nChainPop=4
Zinit <- initZtwDEMCNormal( theta0, diag(sdTheta^2), nChainPop=.nChainPop, nPop=.nPop)
#dim(Zinit)
body <- expression(
res <- concatPops(resBlock <- twDEMCBlock( Zinit, nGen=100
,dInfos=list(list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen))
,nPop=.nPop
#fLogDenScale=-1/2,
#debugSequential=TRUE
))
)
checkException( eval(body))
suppressWarnings(dir.create("tmp"))
.remoteDumpfileBasename=file.path("tmp","testDump")
.remoteDumpfile <- paste(.remoteDumpfileBasename,".rda",sep="")
# dump on initial parallel calculation, file is created in sequentail mode too
unlink(.remoteDumpfile)
checkTrue( !file.exists(.remoteDumpfile))
#mtrace(twDEMCBlockInt)
isParallel <- sfParallel()
if( !isParallel ) sfInit(TRUE,2) # remote dumpfile only created in parallel mode
body <- expression(
res <- concatPops(resBlock <- twDEMCBlock( Zinit, nGen=100
,dInfos=list(list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen))
,nPop=.nPop
,remoteDumpfileBasename=.remoteDumpfileBasename
#fLogDenScale=-1/2,
#debugSequential=TRUE
))
)
eval(body)
checkException( eval(body))
if( !isParallel ) sfStop()
checkTrue( file.exists(.remoteDumpfile))
# deprecated: dump on Metropolis step (provided initial X and logDenX)
# because fLogDen is evaluated at least once in twDEMC before invoking .doDEMCSteps
.tmp.f <- function(){
unlink(.remoteDumpfile)
checkTrue( !file.exists(.remoteDumpfile))
#mtrace(.doDEMCStep)
body <- expression(
res <- concatPops(resBlock <- twDEMCBlock( Zinit, nGen=100,
fLogDen=logDenGaussian, argsFLogDen=argsFLogDen,
nPop=.nPop,
remoteDumpfileBasename=.remoteDumpfileBasename,
#fLogDenScale=-1/2,
debugSequential=TRUE
,X = matrix(1,nrow=1,ncol=.nChains)
,logDenX = 5
,logDenCompX = character(0)
))
)
checkException( eval(body))
checkTrue( file.exists(.remoteDumpfile))
# dump on Metropolis step (provided initial X and logDenX)
unlink(.remoteDumpfile)
checkTrue( !file.exists(.remoteDumpfile))
#mtrace(.doDEMCStep)
#options(error=recover)
body <- expression({
#mtrace(twCalcLogDenPar);
res <- concatPops(resBlock <- twDEMCBlock( Zinit, nGen=100
,fLogDen=logDenGaussian, argsFLogDen=argsFLogDen
,nPop=.nPop
,remoteDumpfileBasename=.remoteDumpfileBasename
#fLogDenScale=-1/2,
#debugSequential=TRUE,
,X = matrix(1,nrow=1,ncol=.nChains)
))
})
checkException( eval(body))
checkTrue( file.exists(.remoteDumpfile))
}
.tmp.f <- function(){
load( .remoteDumpfile )
debugger(get(.remoteDumpfileBasename))
#choose last column (18)
#interactively execute commands from body of doDEMCStep
#require(debug)
mtrace(remoteFun)
do.call(remoteFun, c(remoteFunArgs, list(...)))
}
}
test.twoStepMetropolis <- function(){
# internal Metropolis step
# nice starting values
.nPop=2
argsFLogDen <- list(
fModel=dummyTwDEMCModel, ### the model function, which predicts the output based on theta
obs=obs, ### vector of data to compare with
invCovar=invCovar/1e10, ### do not constrain by data, the inverse of the Covariance of obs (its uncertainty)
thetaPrior = thetaTrue, ### the prior estimate of the parameters
invCovarTheta = invCovarTheta, ### the inverse of the Covariance of the prior parameter estimates
xval=xval
)
do.call( logDenGaussian, c(list(theta=theta0),argsFLogDen))
.nChainPop=4
Zinit <- initZtwDEMCNormal( theta0, diag(sdTheta^2), nChainPop=.nChainPop, nPop=.nPop)
#dim(Zinit)
X <- adrop(Zinit[ncol(Zinit),, ,drop=FALSE],1) #last row of Zinit
#mtrace(twCalcLogDenPar)
#mtrace(logDenGaussian)
XLogDen <- twCalcLogDenPar(fLogDen=logDenGaussian, xProp=t(X), logDenCompX=numeric(0), intResComp="parms", argsFLogDen=argsFLogDen
,debugSequential=TRUE)
.nGen=100
.thin=5
#mtrace(sfRemoteWrapper)
#mtrace(twDEMCBlockInt)
#mtrace(.doDEMCSteps )
#mtrace(logDenGaussian)
#mtrace(.doDEMCStep)
resa <- concatPops(resBlock <- twDEMCBlock( Zinit, nGen=.nGen
,dInfos=list(list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen
,intResComp="parms"))
,nPop=.nPop
,X = X, logDenCompX=t(XLogDen$logDenComp)
,controlTwDEMC = list(thin=.thin)
,debugSequential=TRUE
,doRecordProposals=TRUE
))
#str(resa)
rescoda <- as.mcmc.list(resa)
plot(rescoda, smooth=FALSE)
res <- thin(resa, start=50)
rescoda <- as.mcmc.list(res)
plot(rescoda, smooth=FALSE)
#gelman.diag(rescoda)
#summary(rescoda)
#str(summary(rescoda))
suppressWarnings({ #glm fit in summary
summary(rescoda)$statistics[,"Mean"]
summary(rescoda)$statistics[,"SD"]
thetaTrue
sdTheta
(.popmean <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"Mean"]}))
(.popsd <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"SD"]}))
})
# 1/2 orders of magnitude around prescribed sd for theta
.pop=1
for( .pop in seq(along.with=.popsd) ){
checkMagnitude( sdTheta, .popsd[[.pop]] )
}
# check that thetaTrue is in 95% interval
.pthetaTrue <- sapply(1:2, function(.pop){
pnorm(thetaTrue, mean=.popmean[[.pop]], sd=.popsd[[.pop]])
})
checkInterval( .pthetaTrue )
}
test.twoStepMetropolisTemp <- function(){
# same as goodStartPrior with internal metropolis step
.nPop=2
argsFLogDen <- list(
fModel=dummyTwDEMCModel, ### the model function, which predicts the output based on theta
obs=obs, ### vector of data to compare with
invCovar=invCovar/1e10, ### do not constrain by data, the inverse of the Covariance of obs (its uncertainty)
thetaPrior = thetaTrue, ### the prior estimate of the parameters
invCovarTheta = invCovarTheta, ### the inverse of the Covariance of the prior parameter estimates
xval=xval
)
do.call( logDenGaussian, c(list(theta=theta0),argsFLogDen))
.nChainPop=4
Zinit <- initZtwDEMCNormal( theta0, diag(sdTheta^2), nChainPop=.nChainPop, nPop=.nPop)
#dim(Zinit)
X <- adrop(Zinit[nrow(Zinit),, ,drop=FALSE],1) #last row of Zinit
#mtrace(twCalcLogDenPar)
#mtrace(logDenGaussian)
XLogDen <- twCalcLogDenPar(fLogDen=logDenGaussian, xProp=t(X), logDenCompX=numeric(0), intResComp="parms", argsFLogDen=argsFLogDen ,debugSequential=TRUE)
.nGen=100
.thin=5
#mtrace(sfRemoteWrapper)
#mtrace(twDEMCBlockInt)
#mtrace(logDenGaussian) #check metropolisStepTemp
#trace(logDenGaussian,quote(recover))
resa <- concatPops(resBlock <- twDEMCBlock( Zinit, nGen=.nGen
,dInfos=list(list(
fLogDen=logDenGaussian, argsFLogDen=argsFLogDen
,intResComp="parms" ))
,nPop=.nPop
,X = X, logDenCompX=t(XLogDen$logDenComp)
,debugSequential=TRUE
,TSpec = cbind( T0=c(obs=10,parms=1), TEnd=1/10)
,controlTwDEMC = list(thin=.thin)
))
#str(resa)
matplot(resa$temp, type="l")
checkTrue( all(resa$temp[1:10,1] > 1) )
rescoda <- as.mcmc.list(resa)
plot(rescoda, smooth=FALSE) #decreasing amplitude
res <- thin(resa, start=60)
rescoda <- as.mcmc.list(res)
plot(rescoda, smooth=FALSE)
#gelman.diag(rescoda)
#summary(rescoda)
#str(summary(rescoda))
suppressWarnings({ #glm fit in summary
summary(rescoda)$statistics[,"Mean"]
summary(rescoda)$statistics[,"SD"]
thetaTrue
sdTheta
(.popmean <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"Mean"]}))
(.popsd <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"SD"]}))
})
# 1/2 orders of magnitude around prescribed sd for theta
.pop=1
for( .pop in seq(along.with=.popsd) ){
checkMagnitude( sdTheta, .popsd[[.pop]] )
}
# check that thetaTrue is in 95% interval
.pthetaTrue <- sapply(1:2, function(.pop){
pnorm(thetaTrue, mean=.popmean[[.pop]], sd=.popsd[[.pop]])
})
checkInterval( .pthetaTrue )
}
t_est.goodStartSeqMultiTemp <- function(){
# same as goodStartSeq but with decreasing temperature and multiTemp
.nPop=2
argsFLogDen <- list(
fModel=dummyTwDEMCModel, ### the model function, which predicts the output based on theta
obs=obs, ### vector of data to compare with
invCovar=invCovar/1e10, ### do not constrain by data, the inverse of the Covariance of obs (its uncertainty)
thetaPrior = thetaTrue, ### the prior estimate of the parameters
invCovarTheta = invCovarTheta, ### the inverse of the Covariance of the prior parameter estimates
xval=xval
)
do.call( logDenGaussian, c(list(theta=theta0),argsFLogDen))
Zinit <- initZtwDEMCNormal( theta0, diag(sdTheta^2), nChainPop=4, nPop=.nPop)
#dim(Zinit)
.nGen=100
#mtrace(sfRemoteWrapper)
#mtrace(.doDEMCStep)
#mtrace(.doDEMCSteps )
#mtrace(logDenGaussian)
#mtrace(twCalcLogDenPar)
#mtrace(twDEMCBlockInt)
res <- concatPops(resBlock <- twDEMCBlock( Zinit, nGen=.nGen,
fLogDen=logDenGaussian, argsFLogDen=argsFLogDen,
nPop=.nPop,
#controlTwDEMC=list(thin=1, T0=20, TEnd=20, useMultiT=TRUE, TFix=c(parms=1)),
controlTwDEMC=list(thin=1, T0=20, TEnd=20, useMultiT=TRUE, Tprop=c(obs=0.5,parms=1), TFix=c(parms=1)),
debugSequential=TRUE
))
#str(res)
checkEquals((.nGen%/%res$thin)+1,nrow(res$logDen))
rescoda <- as.mcmc.list(res)
plot(rescoda)
#gelman.diag(rescoda)
#summary(rescoda)
#str(summary(rescoda))
suppressWarnings({ #glm fit in summary
summary(rescoda)$statistics[,"Mean"]
summary(rescoda)$statistics[,"SD"]
thetaTrue
sdTheta
(.popmean <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"Mean"]}))
(.popsd <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"SD"]}))
})
# 1/2 orders of magnitude around prescribed sd for theta
.pop=1
for( .pop in seq(along.with=.popsd) ){
checkMagnitude( sdTheta, .popsd[[.pop]] )
}
# check that thetaTrue is in 95% interval
.pthetaTrue <- sapply(1:2, function(.pop){
pnorm(thetaTrue, mean=.popmean[[.pop]], sd=.popsd[[.pop]])
})
checkInterval( .pthetaTrue )
}
.test.fLogDenScale <- function(){
# same as goodStart sequence but with logDens function returning cost instead of logDensity
.nPop=2
argsFLogDen <- list(
fModel=dummyTwDEMCModel, ### the model function, which predicts the output based on theta
obs=obs, ### vector of data to compare with
invCovar=invCovar/1e10, ### do not constrain by data, the inverse of the Covariance of obs (its uncertainty)
thetaPrior = thetaTrue, ### the prior estimate of the parameters
invCovarTheta = invCovarTheta, ### the inverse of the Covariance of the prior parameter estimates
xval=xval
,scale=1 ##<< override the -1/2
)
do.call( logDenGaussian, c(list(theta=theta0),argsFLogDen))
Zinit <- initZtwDEMCNormal( theta0, diag(sdTheta^2), nChainPop=4, nPop=.nPop)
#dim(Zinit)
.nGen=100
#mtrace(sfRemoteWrapper)
#mtrace(.doDEMCStep)
#mtrace(.doDEMCSteps )
#mtrace(logDenGaussian)
#mtrace(twCalcLogDenPar)
#mtrace(twDEMCBlockInt)
res <- concatPops(resBlock <- twDEMCBlock( Zinit, nGen=.nGen,
dInfos=list(list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen
,fLogDenScale=-1/2 # here apply the scale in twDEMCBlock
))
,nPop=.nPop
,controlTwDEMC=list(thin=1)
,debugSequential=TRUE
))
#str(res)
checkEquals((.nGen%/%res$thin)+1,nrow(res$logDen))
#windows(record=TRUE)
rescoda <- as.mcmc.list(res)
plot(rescoda)
#gelman.diag(rescoda)
#summary(rescoda)
#str(summary(rescoda))
suppressWarnings({ #glm fit in summary
summary(rescoda)$statistics[,"Mean"]
summary(rescoda)$statistics[,"SD"]
thetaTrue
sdTheta
(.popmean <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"Mean"]}))
(.popsd <- lapply(list(p1=1:4,p2=5:8),function(i){summary(rescoda[i])$statistics[,"SD"]}))
})
# 1/2 orders of magnitude around prescribed sd for theta
.pop=1
for( .pop in seq(along.with=.popsd) ){
checkMagnitude( sdTheta, .popsd[[.pop]] )
}
# check that thetaTrue is in 95% interval
.pthetaTrue <- sapply(1:2, function(.pop){
pnorm(thetaTrue, mean=.popmean[[.pop]], sd=.popsd[[.pop]])
})
checkInterval( .pthetaTrue )
}
Try the twDEMC package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.