Nothing
R/sa.R
In twDEMC: parallel DEMC
twDEMCSA <- function(
### Simulated annealing DEMC
thetaPrior ##<< vector of parameters, point estimate
##, alternatively array with initial states, as returned by \code{\link{initZtwDEMCNormal}}
,covarTheta ##<< the a prior covariance of parameters, see \code{link{initZtwDEMCNormal}}
,nGen=512 ##<< number of generations, i.e. Monte Carlo steps
,nObs ##<< integer vector specifying the number of observations for each result component
## for each name in result of logDensity functions there must be an entry in nObs
,... ##<< further argument to \code{\link{twDEMCBlockInt}}
, dInfos ##<< Information on density functions (see \code{\link{twDEMCBlockInt}})
, m0 = calcM0twDEMC(length(thetaPrior),nChainPop) ##<< minimum number of samples in step for extending runs
, controlTwDEMC = list() ##<< control parameters with entry \code{thin} (see \code{\link{twDEMCBlockInt}})
, debugSequential=FALSE ##<< set to TRUE to avoid parallel execution, good for debugging
, restartFilename=NULL ##<< filename to write intermediate results to
, remoteDumpfileBasename=NULL ##<< fileBasename to write dumps to on error
#
,nChainPop=4 ##<< number of chains within population
,nPop=2 ##<< number of populations
,doIncludePrior=FALSE ##<< should the prior be part of initial population
##<< Recommendation to set to false, because if the TRUE parameter is in initial set, initial Temperature calculation is too optimistic
#
, ctrlBatch = list( ##<< list of arguments controlling batch executions, see also \code{\link{twDEMCSACont}}
##describe<<
nGen0=m0*controlTwDEMC$thin*3 ##<< number of generations for the initial batch
#,useSubspaceAdaptation=FALSE # #<< if TRUE then overall space is devided and each subspace is explored with locally adapted DEMC, see \code{\link{divideTwDEMCSACont}}
##end<<
)
, ctrlT = list( ##<< list of arguments controlling Temperature decrease, see also \code{\link{twDEMCSACont}}
##describe<<
qTempInit=0.4 ##<< quantile of logDensities used to calculate initial beginning and end temperature, with default 0.4: 40% of the space is accepted
,TBaseInit=NULL ##<< numeric scalar: initial base temperature. If given this is used for calculating initial temperature
,isVerbose=FALSE ##<< boolean scalar: set to TRUE to report stream base temperatures during batches
,TEndFixed=NULL ##<< set to a scalar end temperature, e.g. in order to decrease temperatue to a given temperature
,qBest=0.1 ##<< the proportion of best samples to base calculation of target temperature on
##end<<
)
, ctrlConvergence = list() ##<< list or arguments controlling check for convergence, see \code{\link{twDEMCSACont}}
#, ctrlSubspaces = list() # #<< list of arguments controlling splitting and merging of subspaces, see \code{\link{divideTwDEMCSACont}}
){
##detail<< \describe{\item{Continuing a previous run}{
## When supplying the first argument an object of class twDEMCPops to twDEMCSA, then this run is extended.
## All other parameters, except nGen, are then are ignored.
## In order to change parameters, modify list entry args in the twDEMCPops object.
##}}
if( inherits(thetaPrior,"twDEMCPops") && 0 != length(thetaPrior$args) ){
if( !missing(debugSequential) )
thetaPrior$args$debugSequential=debugSequential
if( !missing(restartFilename) )
thetaPrior$args$restartFilename=restartFilename
if( !missing(remoteDumpfileBasename) )
thetaPrior$args$remoteDumpfileBasename=remoteDumpfileBasename
if( !missing(controlTwDEMC) )
thetaPrior$args$controlTwDEMC <- twMergeLists( thetaPrior$args$controlTwDEMC, controlTwDEMC )
if( !missing(ctrlBatch) )
thetaPrior$args$ctrlBatch <- twMergeLists( thetaPrior$args$ctrlBatch, ctrlBatch )
if( !missing(ctrlT) )
thetaPrior$args$ctrlT <- twMergeLists( thetaPrior$args$ctrlT, ctrlT )
if( !missing(ctrlConvergence) )
thetaPrior$args$ctrlConvergence <- twMergeLists( thetaPrior$args$ctrlConvergence, ctrlConvergence )
#if( !missing(ctrlSubspaces) )
# thetaPrior$args$ctrlSubspaces <- twMergeLists( thetaPrior$args$ctrlSubspaces, ctrlSubspaces )
if( !missing(dInfos) ){
warning("twDEMCSA: continueing previous run, but with different argument dInfos")
thetaPrior$args$dInfos <- dInfos # overwrite
}
thetaPrior$failureMsg <- NULL # clear previous return and error messages
thetaPrior$finishEarly <- NULL
#ret <- if( 0!=length(thetaPrior$args$ctrlBatch) &&
# 0!=length(thetaPrior$args$ctrlBatch$useSubspaceAdaptation) &&
# thetaPrior$args$ctrlBatch$useSubspaceAdaptation
# )
# ret <- do.call(divideTwDEMCSACont,c(list(thetaPrior, nGen=nGen),thetaPrior$args))
#else
mcPops <- do.call(twDEMCSACont,c(list(thetaPrior, nGen=nGen),thetaPrior$args))
#ret$args$ctrlBatch$useSubspaceAdaptation <- ctrlBatch$useSubspaceAdaptation
return(mcPops)
}
#
#-------------- generate initial sample
# check for array before evaluating formal arguments, because m0 then changes, which is used in default values
#mtrace(initZtwDEMCNormal)
Zinit0 <- if( is.array(thetaPrior) ){
.dimZ <- dim(thetaPrior)
if( length(.dimZ) != 3 ) stop("twDEMCSA: wrong dimensions of thetaPrior")
nChain <- .dimZ[3]
if( missing(nChainPop) ) nChainPop <- nChain %/% nPop
if( missing(nPop) ) nPop <- nChain %/% nChainPop
if( .dimZ[3] != nChainPop*nPop ) stop("twDEMCSA: wrong number of chains in thetaPrior")
m0 = calcM0twDEMC( ncol(thetaPrior) ,nChainPop) ##<< minimum number of samples in step for extending runs
if( .dimZ[1] < m0 ) stop(paste("twDEMCSA: too few cases in thetaPrior (",.dimZ[1],"). Need at least m0=",m0,sep=""))
thetaPrior
} else initZtwDEMCNormal( thetaPrior, covarTheta, nChainPop=nChainPop, nPop=nPop, doIncludePrior=doIncludePrior, m0=m0)
#
#-------- fill in default argument values
if( is.null(controlTwDEMC$thin) ) controlTwDEMC$thin <- 4
thin <- controlTwDEMC$thin
# The user may provide list arguments to overwrite some default entries.
# With formals and twMergeLists, we make sure to take the other default entries in the lists
frm <- formals()
#ctrlSubspaces <- if( hasArg(ctrlSubspaces) ) twMergeLists( eval(frm[["ctrlConvergence"]]), ctrlSubspaces ) else ctrlSubspaces
ctrlBatch <- if( hasArg(ctrlBatch) ) twMergeLists( eval(frm[["ctrlBatch"]]), ctrlBatch ) else ctrlBatch
ctrlT <- if( hasArg(ctrlT) ) twMergeLists( eval(frm[["ctrlT"]]), ctrlT ) else ctrlT
#
#---------- calculate logDensity of initial components
ss <- stackChains(Zinit0)
#print("before twCalcLogDensPar"); recover()
cat("Calculating logDensity for ", length(ss)," initial states for m0=",m0,"rows per population.\n")
logDenDS0 <- (tmp <- twCalcLogDensPar(dInfos, ss, remoteDumpfileBasename=remoteDumpfileBasename))$logDenComp
boFinite <- apply(is.finite(logDenDS0), 1, all )
m0FiniteFac <- sum(boFinite) / nrow(logDenDS0)
# if too few initial states have yielded finite results of logDensity functions, add some more initial states
if( m0FiniteFac == 1){
Zinit <- Zinit0
logDenDS <- logDenDS0
}else if( m0FiniteFac < 0.05 ){
stop(paste("twDEMCSA: less than 5% finite solutions in initial exploration (",m0FiniteFac,"). Check setup."),sep="")
}else if( m0FiniteFac > 0.9 ){
sumLogDen0 <- rowSums(logDenDS0)
Zinit <- replaceZinitNonFiniteLogDens( Zinit0, sumLogDen0 )
logDenDS <- logDenDS0
logDenDS[!boFinite,] <- logDenDS[which.min(sumLogDen0),]
}else{
# generate more proposals and concatenate finite cases from both
Zinit1 <- initZtwDEMCNormal( thetaPrior, covarTheta, nChainPop=nChainPop, nPop=nPop, doIncludePrior=doIncludePrior
,m0FiniteFac=min(1,m0FiniteFac)
)
ss1 <- stackChains(Zinit1)
logDenL <- lapply( dInfos, function(dInfo){
resLogDen <- twCalcLogDenPar( dInfo$fLogDen, ss1, argsFLogDen=dInfo$argsFLogDen)$logDenComp
})
# replace missing cases
logDenDS1 <- abind( logDenL )
# logDenDS1[ 1:round(nrow(logDenDS1)/1.7),1] <- NA # for testing replacement of non-finite cases
boFinite1 <- apply(is.finite(logDenDS1), 1, all )
m0FiniteFac1 <- sum(boFinite1) / nrow(logDenDS1)
ss12 <- rbind( ss[boFinite, ,drop=FALSE], ss1[boFinite1, ,drop=FALSE] )
logDenDS12 <- rbind( logDenDS0[boFinite, ,drop=FALSE], logDenDS1[boFinite1, ,drop=FALSE])
nDiff <- nrow(ss) - nrow(ss12)
if( nDiff/nrow(ss) > 0.1 ){
stop("twDEMCSA: too many states yielding non-finite logDensities.")
}else if( nDiff > 0 ){
# duplicate missings to refill
iSample <- sample( nrow(ss12), size=nDiff )
ss12 <- rbind( ss12, ss12[ iSample, ] )
logDenDS12 <- rbind( logDenDS12, logDenDS12[iSample,] )
}else if(nDiff < 0){
iSample <- sample( nrow(ss12), size=nrow(ss))
ss12 <- ss12[ iSample, ]
logDenDS12 <- logDenDS12[iSample,]
}
## reshape to chains
ncolZ <- dim(Zinit0)[3]
tmp <- matrix(1:nrow(ss12), ncol=ncolZ)
Zinit <- abind( lapply( 1:ncolZ, function(i){ ss12[ tmp[,i], ,drop=FALSE]}), rev.along=0 )
logDenDS <- logDenDS12
} # end generating nonfinite Zinit
#
# apply overfitting correction?
# No need here because Temperature would be set to 1 for corrected and noncorrected values
#
# now the length of resComp is known (colnames(logDenDS)), so check the TFix and TMax parameters
iNoNames <- which( "" == colnames(logDenDS) )
if( length(iNoNames) ) stop("Encountered unnamed logDensitiy components. Provide names of each return component in fLogDen.")
nResComp <- ncol(logDenDS)
ctrlT$TFix <- completeResCompVec( ctrlT$TFix, colnames(logDenDS) )
iFixTemp <- which( is.finite(ctrlT$TFix) )
iNonFixTemp <- which( !is.finite(ctrlT$TFix) )
ctrlT$TMax <- tmp <- completeResCompVec( ctrlT$TMax, colnames(logDenDS) )
iMaxTemp <- which( is.finite(ctrlT$TMax) )
#
# check nObs names and set all TFix to 1
if( 0 == length(names(nObs))){
if( length(nObs) == length(iNonFixTemp) ) names(nObs) <- colnames(logDenDS)[iNonFixTemp] else
stop("Provide names with nObs, or provide one component for each temperated (i.e. not in TFix) fLogDen return component.")
}
iMissing <- which( !(colnames(logDenDS)[iNonFixTemp] %in% names(nObs)) )
if( length(iMissing) ) stop("Missing the following components in nObs: ",paste( colnames(logDenDS)[iNonFixTemp][iMissing],sep=",") )
# make nObs correspond to each data stream, set TFix Components to 1
nObsComp <- structure( rep( NA_real_, ncol(logDenDS) ), names=colnames(logDenDS) )
#nObsComp[names(nObs)] <- nObs # if there are seveval occurence of the same name, others still NA
# name <- names(nObs)[1]
for( name in names(nObs) ){ nObsComp[ names(nObsComp) == name ] <- nObs[name] }
nObs <- nObsComp
#
#if( 0 == length( TMaxInc) ) TMaxInc <- structure( rep( NA_real_, nResComp), names=colnames(logDenDS) )
#if( nResComp != length( TMaxInc) ) stop("twDEMCSA: TMaxInc must be of the same length as number of result Components.")
#iMaxIncTemp <- which( is.finite(TMaxInc) )
#
#------ intial temperatures:
#sapply( seq_along(expLogDenBest), function(i){ max(logDenDS[,i]) })
##details<< \describe{\item{Initial temperature}{
## Initial parameters are ranked according to their maximum log-Density across components.
## The parameters and logDensity results at rank position defined by argument \code{ctrlT$qTempInit} is selected.
## The stream temperatures are inferred by deviding logDensity components by the number of observations.
## From these, a common base temperatue is calculated by (see \code{\link{calcBaseTemp}}) and rescaled to stream temperatures by \code{\link{calcStreamTemp}}.
##}}
#print("twDEMCSA: before calculating initial temperature"); recover()
if( length(ctrlT$TBaseInit) ){
temp0 <- calcStreamTemp(ctrlT$TBaseInit, nObs, TFix=ctrlT$TFix, iFixTemp=iFixTemp)
} else {
iTmax <- 5
T0s <- rep(NA_real_, iTmax )
T <- nObs
T[iFixTemp] <- ctrlT$TFix[iFixTemp]
T0s[1] <- .Machine$double.xmax
relChange <- relChangeReq <- 0.05 # if below 5% change then no need to iterate further
iT <- 1
while( (iT < iTmax) && (abs(relChange) >= relChangeReq) ){
iT <- iT + 1
resLogDenT <- calcTemperatedLogDen( logDenDS, T )
iBest <- getBestModelIndices( resLogDenT, dInfos, prob=ctrlT$qTempInit )
T0 <- T0s[iT] <- calcBaseTempSk( -2*logDenDS[iBest, ,drop=FALSE], nObs, ctrlT$TFix, iFixTemp = iFixTemp, iNonFixTemp = iNonFixTemp)
T <- calcStreamTemp(T0, nObs, ctrlT$TFix, iFixTemp = iFixTemp)
relChange <- (T0 - T0s[iT-1])/T0
}
temp0 <- T
if( !all(is.finite(temp0)) ) stop("twDEMCSA: encountered non-finite Temperatures.")
}
ctrlT$TMax[iNonFixTemp] <- pmin(ifelse(is.finite(ctrlT$TMax),ctrlT$TMax, Inf), ifelse(is.finite(temp0),temp0,Inf) )[iNonFixTemp] # decrease TMax
#if( any(temp0 > 8)) stop("twDEMCSA: encountered too high temperature.")
print(paste("initial T=",paste(signif(temp0,2),collapse=",")," for ",min(nGen, ctrlBatch$nGen0)," generations"," ", date(), sep="") )
#if( sum(temp0 > 1) == 0){ dump.frames("parms/debugDump",TRUE); stop("twDEMCSA: no temperature > 0") }
#
mcPops <- res0 <- res <- twDEMCBlock( Zinit
, nGen=min(nGen, ctrlBatch$nGen0)
#,nGen=16, debugSequential=TRUE
, dInfos=dInfos
, TSpec=cbind( T0=temp0, TEnd=temp0 )
, nPop=nPop
, m0=m0, controlTwDEMC=controlTwDEMC, debugSequential=debugSequential
# XXTODO: replace lines below later on by ...
#, blocks = blocks
,...
)
print(paste("finished initial ",ctrlBatch$nGen0," out of ",nGen," gens. ", date(), sep="") )
#
#if( nGen > ctrlBatch$nGen0 ){ #depre, now need to call Continue to return proper diagnostics
#if( TRUE ){ # need to call Continue to return proper diagnostics
#nObsLocal <- nObs
#ret <- if( ctrlBatch$useSubspaceAdaptation ){
# ret <- divideTwDEMCSACont( mc=res0, nGen=nGen-ctrlBatch$nGen0, nObs=nObs
# , m0=m0, controlTwDEMC=controlTwDEMC, debugSequential=debugSequential, restartFilename=restartFilename
# , ctrlBatch=ctrlBatch, ctrlT=ctrlT, ctrlConvergence=ctrlConvergence, ctrlSubspaces=ctrlSubspaces
# , ...
# )
#}else {
mcPops <- twDEMCSACont( mc=res0, nGen=nGen-ctrlBatch$nGen0, nObs=nObs
, m0=m0, controlTwDEMC=controlTwDEMC, debugSequential=debugSequential, restartFilename=restartFilename
, ctrlBatch=ctrlBatch, ctrlT=ctrlT, ctrlConvergence=ctrlConvergence
, ...
)
#}
#}
#ret$args$ctrlBatch$useSubspaceAdaptation <- ctrlBatch$useSubspaceAdaptation # remember this argument
##value<< An object of class \code{twDEMCPops} as described in \code{\link{twDEMCBlockInt}}.
## See \code{\link{subset.twDEMCPops}} for processing further handling of this class.
mcPops
}
attr(twDEMCSA,"ex") <- function(){
# Best to have a look at the the package vignettes
#--------------- single density ----------------------------
# We will use logDenGaussian as logDensity function that compares a simple linear model with observations.
data(twLinreg1)
# collect all the arguments to the logDensity in a list (except the first argument of changing parameters)
argsFLogDen <- with( twLinreg1, 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=twLinreg1$theta0),argsFLogDen))
do.call( logDenGaussian, c(list(theta=twLinreg1$thetaTrue),argsFLogDen)) # slightly largere misfit than nObs/2=15, underestimated sdObs
.nGen=200
.nPop=2
mcPops <- twDEMCSA(
theta=twLinreg1$theta0, covarTheta=diag(twLinreg1$sdTheta^2) # to generate initial population
, nGen=.nGen
, dInfos=list(den1=list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen))
, nPop=.nPop # number of independent populations
, controlTwDEMC=list(thin=4) # see twDEMCBlockInt for all the tuning options
, ctrlConvergence=list(maxRelTChangeCrit=0.1) # ok if T changes less than 10%
, ctrlT=list(TFix=c(parms=1)) # do not use increased temperature for priors
, nObs=c(obs=length(argsFLogDen$obs)) # number of observations used in temperature calculation
)
#mcp <- twDEMCSA( mcp, nGen=2000)
mcPops <- twDEMCSA( mcPops, nGen=400) # continue run
rescoda <- as.mcmc.list(mcPops)
plot(rescoda, smooth=FALSE)
mcChains1 <- concatPops(mcPops) # array representation instead of list of pops, last dim is the chain
mcChains2 <- concatPops(stackChainsPop(mcPops)) # combining dependent chains within one population
mcChains3 <- concatPops(subsetTail(mcPops,0.5)) # take only the last part of the chains
c(getNGen(mcChains1), getNGen(mcChains3))
plot(as.mcmc.list(mcChains2), smooth=FALSE)
#--------------- multiple densities -------------------------
data(twTwoDenEx1)
thetaPrior <- twTwoDenEx1$thetaTrue
covarTheta <- diag((thetaPrior*0.3)^2)
invCovarTheta <- (thetaPrior*0.3)^2 # given as independent variances for faster calculation
thresholdCovar = 0.3 # the true value used to generate the observations
thresholdCovar = 0 # the effective model that glosses over this threshold
#str(twTwoDenEx1)
nObs <- c( obsSparse=length(twTwoDenEx1$obs$y1), obsRich=length(twTwoDenEx1$obs$y2) )
dInfos=list(
dSparse=list(fLogDen=denSparsePrior
, argsFLogDen=list(thresholdCovar=thresholdCovar, twTwoDenEx=twTwoDenEx1, theta0=thetaPrior, thetaPrior=thetaPrior, invCovarTheta=invCovarTheta)
#, maxLogDen=-1/2*nObs[c("parmsSparse","obsSparse")] # control overfitting
)
,dRich=list(fLogDen=denRichPrior
, argsFLogDen=list(thresholdCovar=thresholdCovar,twTwoDenEx=twTwoDenEx1, theta0=thetaPrior, thetaPrior=thetaPrior, invCovarTheta=invCovarTheta)
, maxLogDen=c(parmsRich=0,-1/2*nObs["obsRich"]) # control overfitting for rich datastream
)
)
blocks = list(
a=list(dInfoPos="dSparse", compPos="a")
,b=list(dInfoPos="dRich", compPos="b")
)
names(do.call( dInfos$dSparse$fLogDen, c(list(theta=twTwoDenEx1$theta0),dInfos$dSparse$argsFLogDen)))
names(do.call( dInfos$dRich$fLogDen, c(list(theta=twTwoDenEx1$theta0),dInfos$dRich$argsFLogDen)))
#trace(twDEMCSACont, recover )
#trace(twDEMCSA, recover )
res <- res0 <- twDEMCSA( thetaPrior, covarTheta, dInfos=dInfos, blocks=blocks, nObs=nObs
, nGen=3*256
, ctrlT=list( TFix=c(parmsSparse=1,parmsRich=1) ) # no increased Temperature for priors
, ctrlBatch=list( nGenBatch=256 )
, debugSequential=TRUE
, controlTwDEMC = list(
DRgamma=0.1 # use Delayed rejection
#,controlOverfittingMinNObs = 20 # use overfitting control (for obsRich), recommended on using single density
)
#, restartFilename=file.path("tmp","example_twDEMCSA.RData")
)
res <- twDEMCSA( res0, nGen=2*256 ) # extend the former run
(TCurr <- getCurrentTemp(res))
# Note that T does decrease to 1
# This accounts for structural model mismatch in addition to observation uncertinaty
mc0 <- concatPops(res)
mcE <- concatPops(subsetTail(res,0.2)) # only the last 20%
plot( as.mcmc.list(mc0) , smooth=FALSE )
matplot( mc0$temp, type="l" )
logDenT <- calcTemperatedLogDen(stackChains(mcE$resLogDen), TCurr)
iBest <- getBestModelIndex( logDenT, res$dInfos )
maxLogDenT <- logDenT[iBest, ]
ss <- stackChains(mcE$parms)
(thetaBest <- ss[iBest, ])
twTwoDenEx1$thetaTrue
(.qq <- apply(ss,2,quantile, probs=c(0.025,0.5,0.975) )) # model error really in paramter b
plot( ss[,"a"], ss[,"b"], col=rainbow(100)[twMisc::twRescale(rowSums(logDenT),c(10,100))] )
plot( ss[,"a"], ss[,"b"], col=rainbow(100)[twMisc::twRescale(logDenT[,"obsSparse"],c(10,100))] )
plot( ss[,"a"], ss[,"b"], col=rainbow(100)[twMisc::twRescale(logDenT[,"obsRich"],c(10,100))] )
plot( ss[,"a"], ss[,"b"], col=rgb(
twMisc::twRescale(logDenT[,"obsSparse"]),0, twMisc::twRescale(logDenT[,"obsRich"]) ))
apply( apply( logDenT, 2, quantile, probs=c(0.1,0.9) ),2, diff )
# density of parameters
plot( density(ss[,"a"])); abline(v=thetaPrior["a"]); abline(v=thetaBest["a"], col="blue")
plot( density(ss[,"b"])); abline(v=thetaPrior["b"]); abline(v=thetaBest["b"], col="blue")
# predictive posterior (best model only)
pred <- pred1 <- with( twTwoDenEx1, fModel(thetaBest, xSparse=xSparse, xRich=xRich) )
plot( pred$y1, twTwoDenEx1$obs$y1 ); abline(0,1)
plot( pred$y2, twTwoDenEx1$obs$y2 ); abline(0,1) # note that deviation is now in y2 - consistent with the introduced bias
} # end examples twDEMCSA
#plot predictive posterior
.tmp.f <- function(){
ss <- stackChains(d3$parms)
apply( ss, 2, mean )
apply( ss, 2, sd )
nSamp <- 40
pSamp <- ss[ sample.int(nSamp), ]
pPred <- apply( pSamp, 1, fModel, xval=xval )
statPred <- t(apply( pPred, 1, function(obsi){ c(mean(obsi), sd(obsi) )} ))
plot(statPred[,1] ~ xval, ylim=range(statPred[,1] +c(1.96,-1.96)*statPred[,2]))
lines( statPred[,1] +1.96*statPred[,2] ~ xval )
lines( statPred[,1] -1.96*statPred[,2] ~ xval )
points( obs ~ xval, col="maroon")
w <- 1/sdObs^2/sum(1/sdObs^2)
lm2 <- lm( obs ~ xval, weights=w )
w <- 1/sdObs^2 # would more correspond to real case see ?lm, but summed to 1 internally
lm2 <- lm.wfit( cbind(1,xval), obs, w=w )
summary(lm2)
obsTrue <- fModel( thetaTrue, xval)
plot(statPred[,1] ~ obsTrue, ylim=range(statPred[,1] +c(1.96,-1.96)*statPred[,2]))
lines( statPred[,1] +1.96*statPred[,2] ~ obsTrue )
lines( statPred[,1] -1.96*statPred[,2] ~ obsTrue )
lines( statPred[,1] -1.96*sdObs ~ obsTrue, col="maroon" )
lines( statPred[,1] +1.96*sdObs ~ obsTrue, col="maroon" )
abline(0,1)
points( obs ~ obsTrue, col="maroon")
}
twDEMCSACont <- function(
### continuing simulated annealing DEMC based on previous result
mc ##<< result of twDEMCBlock
,nGen=512 ##<< overall number of generations to add (within one batch provide ctrlBatch$nGenBatch)
,nObs ##<< integer vector (nResComp) specifying the number of observations for each result component
,... ##<< further argument to \code{\link{twDEMCBlockInt}}
, m0 = calcM0twDEMC(getNParms(mc),getNChainsPop(mc)) ##<< minimum number of samples in step for extending runs
, controlTwDEMC = list() ##<< list argument to \code{\link{twDEMCBlockInt}} containing entry thin
, debugSequential=FALSE ##<< set to TRUE to avoid parallel execution, good for debugging
, restartFilename=NULL ##<< filename to write intermediate results to
#
, ctrlBatch = list( ##<< list of arguments controlling batch executions
##describe<<
nSamplesBatch=m0*8 ##<< number of samples for one call to twDEMCStep (multiplied by thin to get generations), defaults to 8*m0
#nGenBatch=m0*controlTwDEMC$thin*8 ##<< number of generations for one call to twDEMCStep
#nGenBatch=m0*controlTwDEMC$thin*3 ##<< number of generations for one call to twDEMCStep
##<< default: set in a way that on average each population (assuming half are significant) is appended by 2*m0 samples
#, nSampleMin=32 ##<< minimum number of samples in each population within batch so that calculation of average density is stable
, pThinPast=0.5 ##<< in each batch thin the past to given fraction before appending new results
##end<<
)
, ctrlT = list( ##<< list of arguments controlling Temperature decrease
##describe<<
TFix=numeric(0) ##<< numeric vector (nResComp) specifying a finite value for components with fixed Temperatue, and a non-finite for others
, TMax=numeric(0) ##<< numeric vector (nResComp) specifying a maximum temperature for result components.
, TDecProp=0.9 ##<< proportion of Temperature decrease: below one to diminish risk of decreasing Temperature too fast (below what is supported by other data streams)
##end<<
)
, ctrlConvergence = list( ##<< list or arguments controlling check for convergence
##describe<<
maxRelTChangeCrit=0.05 ##<< if Temperature of the components changes less than specified value, the algorithm can finish
, minTBase0 = 1e-4 ##<< if Temperature-1 gets below this temperature, the algorithm can finish
, maxLogDenDriftCrit=0.3 ##<< if difference between mean logDensity of first and fourth quartile of the sample is less than this value, we do not need further batches because of drift in logDensity
, gelmanCrit=1.4 ##<< do not change Temperature, if variance between chains is too high, i.e. Gelman Diag is above this value
, critSpecVarRatio=20 ##<< if proprotion of spectral Density to Variation is higher than this value, signal problems and resort to subspaces
, dumpfileBasename=NULL ##<< scalar string: filename to dump stack before stopping. May set to "recover"
, maxThin=64 ##<< scalar positive integer: maximum thinning interval. If not 0 then thinning is increased on too high spectral density
##end<<
)
){
# save calling arguments to allow an continuing an interrupted run
argsF <- as.list(sys.call())[-1] # do not store the file name and the first two arguments
argsF <- argsF[ !(names(argsF) %in% c("","nGen","mc")) ] # remove positional arguments and arguments mc and nGen
argsFEval <- lapply( argsF, eval.parent ) # remember values instead of language objects, which might not be there on a repeated call
mc$args <- argsFEval
finishEarlyMsg <- NULL
#
#-- fill in default argument values
if( is.null(controlTwDEMC$thin) ){
controlTwDEMC$thin <- 4
}
thin <- controlTwDEMC$thin
frm <- formals()
ctrlConvergence <- if( hasArg(ctrlConvergence) ) twMergeLists( eval(frm[["ctrlConvergence"]]), ctrlConvergence ) else ctrlConvergence
ctrlBatch <- if( hasArg(ctrlBatch) ) twMergeLists( eval(frm[["ctrlBatch"]]), ctrlBatch ) else ctrlBatch
ctrlT <- if( hasArg(ctrlT) ) twMergeLists( eval(frm[["ctrlT"]]), ctrlT ) else ctrlT
if( debugSequential ){
if( 0==length(ctrlConvergence$dumpfileBasename) ) ctrlConvergence$dumpfileBasename <- "recover"
#if( 0==length(ctrlSubspaces$argsFSplit$debugSequential) ) ctrlSubspaces$argsFSplit$debugSequential <- TRUE
}
#
nResComp <- ncol(mc$pops[[1]]$resLogDen)
#print("saCont: before completing temperature settings."); recover()
ctrlT$TFix <- completeResCompVec( ctrlT$TFix, colnames(mc$pops[[1]]$resLogDen) )
iFixTemp <- which( is.finite(ctrlT$TFix) )
iNonFixTemp <- which( !is.finite(ctrlT$TFix) )
ctrlT$TMax <- completeResCompVec( ctrlT$TMax, colnames(mc$pops[[1]]$resLogDen) )
iMaxTemp <- which( is.finite(ctrlT$TMax) )
#
res <- mc
.tmp.f <- function(){
mc0 <- concatPops(resEnd)
#mc0 <- concatPops(res)
matplot( mc0$temp, type="l" )
matplot( mc0$pAccept[,1,], type="l" )
matplot( mc0$pAccept[,2,], type="l" )
plot( as.mcmc.list(mc0), smooth=FALSE )
tmp <- calcTemperatedLogDen(res$pops[[1]]$resLogDen[,,1], getCurrentTemp(res) )
matplot( tmp, type="l" )
plot(rowSums(tmp))
#
matplot( mc0$resLogDen[,3,], type="l" )
matplot( mc0$parms[,"a",], type="l" )
matplot( mc0$parms[1:20,"b",], type="l" )
bo <- 1:10; iPop=1
plot( mc0$parms[bo,"a",iPop], mc0$parms[bo,"b",iPop], col=rainbow(100)[twMisc::twRescale(mc0$resLogDen[bo,"parmsSparse",iPop],c(10,100))] )
}
#
# sum nObs within density - assume nObs positions correspond to resLogDenComp positions
iDens <- seq_along(mc$dInfos)
#iDen=1
iCompsNonFixDen <- lapply( iDens, function(iDen){
irc <- mc$dInfos[[iDen]]$resCompPos
irc <- irc[ !(irc %in% iFixTemp) ]
})
nObsDen <- sapply( iDens, function(iDen){ sum( nObs[iCompsNonFixDen[[iDen]] ]) })
TCurr <- getCurrentTemp(mc)
#
iBatch=1
nGenDone <- 0 # tracking the generations already done
#nBatch <- ceiling( nGen/(ctrlBatch$nSamplesBatch*controlTwDEMC$thin) )
while( nGenDone < nGen){
if((0 < length(restartFilename)) && is.character(restartFilename) && restartFilename!=""){
resRestart.twDEMCSA = res #avoid variable confusion on load by giving a longer name
resRestart.twDEMCSA$nGenDone <- nGenDone # also store updated calculation of burnin time
resRestart.twDEMCSA$args <- argsFEval # also store updated calculation of burnin time
save(resRestart.twDEMCSA, file=restartFilename)
cat(paste("Saved resRestart.twDEMCSA to ",restartFilename,"\n",sep=""))
}
# do the diagnostics on the last halv of the chain, but continue
resEnd <- thin(res, start=ceiling(getNGen(res)) * ctrlBatch$pThinPast ) # neglect the first half part
resSparse <- squeeze(res, length.out=ceiling(getNSamples(res)*ctrlBatch$pThinPast) ) # better thin the entire chain so that distant path is sparsely kept, when continuing
#----- check for high autocorrelation wihin spaces
mcEnd <- concatPops(stackChainsPop(resEnd))
#plot( as.mcmc.list(mcEnd), smooth=FALSE )
specVarRatio <- apply( mcEnd$parms, 3, function(aSamplePurePop){
spec <- spectrum0.ar(aSamplePurePop)$spec
varSample <- apply(aSamplePurePop, 2, var)
spec/varSample
})
ssc <- stackChainsPop(resEnd) # combine all chains of one population
mcl <- as.mcmc.list(ssc)
#plot( as.mcmc.list(mcl), smooth=FALSE )
#plot( res$pops[[1]]$parms[,"b",1] ~ res$pops[[1]]$parms[,"a",1], col=rainbow(255)[round(twRescale(-res$pops[[1]]$resLogDen[,"logDen1",1],c(10,200)))] )
gelmanDiagRes <- try( {tmp<-gelman.diag(mcl); if(length(tmp$mpsrf)) tmp$mpsrf else tmp$psrf[1]} ) # cholesky decomposition may throw errors
#gelmanDiagRes <- try( gelman.diag(mcl)$mpsrf ) # cholesky decomposition may throw errors
# compare intra-chain to inter-chain gelman diag
#iPop <- 4
gelmanDiagPops <- sapply( 1:getNPops(resEnd), function(iPop){
mcl <- as.mcmc.list( subPops(resEnd,iPop) )
tmp <- try(gelman.diag(mcl), silent=TRUE)
if( inherits(tmp,"try-error")) 999 else if(length(tmp$mpsrf)) tmp$mpsrf else tmp$psrf[1]
})
maxGelmanDiagPops <- max(gelmanDiagPops)
T0 <- calcBaseTemp(TCurr, nObs, TFix=ctrlT$TFix, iNonFixTemp=iNonFixTemp)-1
# calculate average resLogDen of last part
resLogDen <- stackChains(concatPops(resEnd)$resLogDen)
.resDrift <- isLogDenDrift(resLogDen, resEnd$dInfos, maxDrift=ctrlConvergence$maxLogDenDriftCrit)
lDenLastPart <- structure( attr(.resDrift,"resTTest")[2,], names=names(resEnd$dInfos) )
##details<< \describe{\item{diagnostics}{
## in resulting twDEMC an entry diagnostics is returned. Its a list with diagnostics after each batch.
## Its components have been calculated on the last halv of the chains.
##describe<<
newDiags <- c(
T0=T0 ##<< Temperature (variance inflation factor) scaled for one observation minus one
,maxGelmanDiagPops=maxGelmanDiagPops ##<< maximum of within population gelman diagnostics
,gelmanDiag=gelmanDiagRes ##<< gelman diag calculated between populations
,maxSpecVarRatio=max(specVarRatio) ##<< ratio of spectral density to sample Variance (max over all parameters). This is a measure that increases with autocorrelation.
,lDenLastPart
#,attr(.resDrift,"logDenComp")
)
# if number of densities changed, reset columns in diagnostics (else rbind fails)
if( (length(resSparse$diagnostics) > 0) && (length(newDiags) != ncol(resSparse$diagnostics)) ){
resSparse$diagnostics <- cbind( resSparse$diagnostics[,1:4]
, matrix(NA, nrow=nrow(resSparse$diagnostics), ncol=length(lDenLastPart)
, dimnames=list(NULL, names(lDenLastPart))) )
}
resSparse$diagnostics <- rbind( resSparse$diagnostics, newDiags)
rownames( resSparse$diagnostics ) <- NULL
## Last components represent the average logDensity of the last 8th of the chains.
##details<< }}
#
##details<< \describe{\item{Adaptive thinning interval (\code{maxThin})}{
## If spectral variance is larger than its critical ratio \code{ctrlConvergence$critSpecVarRatio}
## then thinning interval needs to be increased.
## If doubling the thinning interval is not larger than \code{ctrlConvergence$maxThin} then a new batch
## with higher thinning is attempted. Else the SADEMC quits with reporting an error with return values component \code{failureMsg}.
##}}
if( any(specVarRatio > ctrlConvergence$critSpecVarRatio) ){
newThin <- resSparse$thin * 2
if( length(ctrlConvergence$maxThin) && (ctrlConvergence$maxThin > 0) && (newThin <= ctrlConvergence$maxThin) ){
# increase thinning interval instead of reporting error
print(paste("twDEMCSACont: too much autocorrelation (specVarRatio=",signif(max(specVarRatio),3),"). Increasing thinning interval to ",newThin,sep=""))
resSparse <- thin(resSparse, newThin )
controlTwDEMC$thin <- newThin
}else{
res <- resSparse
res$failureMsg <- paste("twDEMCSACont: too much autocorrelation (specVarRatio=",signif(max(specVarRatio),3),"). Try using divideTwDEMC or higher thinning interval",sep=")")
print(res$failureMsg)
break
}
}
# check if can decrease temperature, must calc TEnd and TEnd0
if( !inherits(gelmanDiagRes,"try-error") &&
maxGelmanDiagPops <= ctrlConvergence$gelmanCrit && # each pop converged
(gelmanDiagRes <= 1.2*maxGelmanDiagPops) # all pops explore the same optimum
){
if( length(ctrlT$TEndFixed)==0 ){
#
# print("twDEMCSACont: before Temperature calculation."); recover()
#debug(getBestModelIndices)
resLogDenT <- calcTemperatedLogDen(resLogDen, TCurr)
#iBest <- getBestModelIndices(resLogDenT, resEnd$dInfos, 0.01)
iBest <- getBestModelIndices(resLogDenT, resEnd$dInfos, ctrlT$qBest)
SkBest <- -2*resLogDen[iBest,,drop=FALSE]
# sapply( obs, function(obsk){ sum(obsk$sd^2) })
#trace(calcBaseTempSk, recover) #untrace(calcBaseTempSk)
TEnd0Target <- calcBaseTempSk(SkBest,nObs=nObs, TFix=ctrlT$TFix, iNonFixTemp=iNonFixTemp, isVerbose=ctrlT$isVerbose) -1
TEnd <- calcStreamTemp(TEnd0Target+1, nObs, TFix=ctrlT$TFix)
TEnd[iMaxTemp] <- pmin(TEnd[iMaxTemp], ctrlT$TMax[iMaxTemp]) # do not increase T above TMax
#relTChange <- abs(TEnd - TCurr)/TEnd
# slower TDecrease to avoid Temperatues that are not supported by other datastreams
TEnd0 <- T0 - ctrlT$TDecProp*(T0-TEnd0Target)
TEnd <- TCurr - ctrlT$TDecProp*(TCurr-TEnd)
} else { # ctrlT$TEndFixed is specified
# ultimately decrease temperature to 1
TEnd0Target <- ctrlT$TEndFixed-1
TEnd0 <- T0 - ctrlT$TDecProp*(T0-(ctrlT$TEndFixed-1))
TEnd <- TCurr <- calcStreamTemp(TEnd0+1, nObs, TFix=ctrlT$TFix)
}
#relTChange <- (TEnd0Target - T0)/(max(ctrlConvergence$minTBase0,TEnd0))
relTChange <- (TEnd0 - T0)/(max(ctrlConvergence$minTBase0,TEnd0))
cat("relTChange=",signif(relTChange,2),"\n")
#if( (max(relTChange) <= maxRelTChange) )
# recover()
#trace(isLogDenDrift, recover )
# may finish early, but do at least one batch
if( relTChange > 0 ){
ctrlT$TDecProp <- ctrlT$TDecProp * 2/3
cat("decreasing ctrlT$TDecProp to ",signif(ctrlT$TDecProp*100,2),"%\n")
}
if( (nGenDone != 0) &&
((max(abs(relTChange)) <= ctrlConvergence$maxRelTChangeCrit) || T0 <= ctrlConvergence$minTBase0 ) &&
!.resDrift
){
#res <- resSparse
# keep original res
finishEarlyMsg <- paste("twDEMCSA: Maximum Temperture change only ",signif(max(relTChange)*100,2)
,"% and no drift in logDensity (lden=",paste(signif(lDenLastPart,3),collapse=","),"). Finishing early.",sep="")
print( finishEarlyMsg)
print(paste("gelmanDiagPops=",signif(gelmanDiagRes,2)
," max(gelmanDiagPop)=",signif(maxGelmanDiagPops,2)
," max(specVarRatio)=",signif(max(specVarRatio),2)
," T0=",signif(T0,2)
," logDen=",paste(signif(lDenLastPart,3),collapse=",")
#," T=",paste(signif(TCurr,2),collapse=",")
, sep=""))
break
}
}else{ # cannot decrease temperature
TEnd0 <- T0
TEnd <-TCurr
}
#TMin <- pmin(TMin, TEnd)
print(paste("gelmanDiagPops=",signif(gelmanDiagRes,2)
," max(gelmanDiagPop)=",signif(maxGelmanDiagPops,2)
," max(specVarRatio)=",signif(max(specVarRatio),2)
," T0=",signif(T0,2)
," logDen=",paste(signif(lDenLastPart,3),collapse=",")
#," T=",paste(signif(TCurr,2),collapse=",")
, sep=""))
.nGenBatch <- min( ctrlBatch$nSamplesBatch*controlTwDEMC$thin, nGen - nGenDone )
print(paste("doing ",.nGenBatch," generations to TEnd0=",signif(TEnd0,2)
,", TEnd=",paste(signif(TEnd,2),collapse=","),sep="") )
res1 <- res <- twDEMCBlock( resSparse
, nGen=.nGenBatch
, debugSequential=debugSequential
, m0=m0
, controlTwDEMC=controlTwDEMC
, TEnd=TEnd
# replace lines below later on by ...
#, blocks = blocks
,...
)
TCurr <- getCurrentTemp(res)
nGenDone <- nGenDone + .nGenBatch
print(paste("finished ",nGenDone," out of ",nGen," gens. ", date(), sep="") )
} # end while( nGenDone < nGen)
res$TGlobal <- max(getCurrentTemp(res)[unlist(iCompsNonFixDen)])
res$args <- argsFEval
res$finishEarly <- finishEarlyMsg
#---------------------- update diagnostics for last batch (replicated from start of the loop)
resEnd <- thin(res, start=ceiling(getNGen(res)) * ctrlBatch$pThinPast ) # neglect the first half part
mcEnd <- concatPops(stackChainsPop(resEnd))
#----- check for high autocorrelation wihin spaces
#plot( as.mcmc.list(mcEnd), smooth=FALSE )
specVarRatio <- apply( mcEnd$parms, 3, function(aSamplePurePop){
spec <- spectrum0.ar(aSamplePurePop)$spec
varSample <- apply(aSamplePurePop, 2, var)
spec/varSample
})
ssc <- stackChainsPop(resEnd) # combine all chains of one population
mcl <- as.mcmc.list(ssc)
#plot( as.mcmc.list(mcl), smooth=FALSE )
#plot( res$pops[[1]]$parms[,"b",1] ~ res$pops[[1]]$parms[,"a",1], col=rainbow(255)[round(twRescale(-res$pops[[1]]$resLogDen[,"logDen1",1],c(10,200)))] )
gelmanDiagRes <- try( {tmp<-gelman.diag(mcl); if(length(tmp$mpsrf)) tmp$mpsrf else tmp$psrf[1]} ) # cholesky decomposition may throw errors
#gelmanDiagRes <- try( gelman.diag(mcl)$mpsrf ) # cholesky decomposition may throw errors
# compare intra-chain to inter-chain gelman diag
#iPop <- 4
gelmanDiagPops <- sapply( 1:getNPops(resEnd), function(iPop){
mcl <- as.mcmc.list( subPops(resEnd,iPop) )
tmp <- try(gelman.diag(mcl), silent=TRUE)
if( inherits(tmp,"try-error")) 999 else if(length(tmp$mpsrf)) tmp$mpsrf else tmp$psrf[1]
})
maxGelmanDiagPops <- max(gelmanDiagPops)
T0 <- calcBaseTemp(TCurr, nObs, TFix=ctrlT$TFix, iNonFixTemp=iNonFixTemp)-1
# calculate average resLogDen of last part
resLogDen <- stackChains(concatPops(resEnd)$resLogDen)
.resDrift <- isLogDenDrift(resLogDen, resEnd$dInfos, maxDrift=ctrlConvergence$maxLogDenDriftCrit)
lDenLastPart <- structure( attr(.resDrift,"resTTest")[2,], names=names(resEnd$dInfos) )
newDiags <- c(
T0=T0 ##<< Temperature (variance inflation factor) scaled for one observation minus one
,maxGelmanDiagPops=maxGelmanDiagPops ##<< maximum of within population gelman diagnostics
,gelmanDiag=gelmanDiagRes ##<< gelman diag calculated between populations
,maxSpecVarRatio=max(specVarRatio) ##<< ratio of spectral density to sample Variance (max over all parameters). This is a measure that increases with autocorrelation.
,lDenLastPart
#,attr(.resDrift,"logDenComp")
)
res$diagnostics <- rbind( res$diagnostics, newDiags)
rownames( res$diagnostics ) <- NULL
res
}
.tmp.f <- function(){ # ggplotResults
mc0 <- concatPops(res)
.nSample <- 128
dfDen <- rbind(
cbind( data.frame( scenario="S1", {tmp <- stackChains(mc0)[,-(1:getNBlocks(mc0))]; tmp[ seq(1,nrow(tmp),length.out=.nSample),] }))
)
dfDenM <- melt(dfDen)
#str(dfDenM)
require(ggplot2)
StatDensityNonZero <- StatDensity$proto( calculate<- function(.,data, scales, dMin=0.001, ...){
res <- StatDensity$calculate(data,scales,...)
# append a zero density at the edges
if( res$density[1] > 0) res <- rbind( data.frame(x=res$x[1]-diff(res$x[1:2]),density=0,scaled=0,count=0),res)
if( res$density[nrow(res)] > 0) res <- rbind( res, data.frame(x=res$x[nrow(res)]+diff(res$x[nrow(res)-c(1,0)]),density=0,scaled=0,count=0))
bo <- (res$density < dMin) & (c(res$density[-1],0) < dMin) & (c(0,res$density[-length(res$density)]) < dMin)
res$density[ bo ] <- NA
res$count[ bo ] <- NA
res$scaled[ bo ] <- NA
res
}, required_aes=c("x"))
#with(StatDensityNonZero, mtrace(calculate))
#with(StatDensityNonZero, mtrace(calculate,F))
stat_densityNonZero <- function (
### constructs a new StatPrior statistics based on aesthetics x and parName
mapping = NULL, data = NULL, geom = "line", position = "stack",
adjust = 1, ...
){
StatDensityNonZero$new(mapping = mapping, data = data, geom = geom,
position = position, adjust = adjust, ...)
}
#pgm <- geom_ribbon( alpha=0.8, aes(ymax = ..density.., ymin = -..density..), stat = "density")
#pgm <- geom_ribbon( aes(ymax = ..density.., ymin = -..density..), stat = "density")
#pgm <- stat_densityNonZero( aes(ymax = -..density..), size=1, geom="line")
pgm <- geom_line(aes(y = +..scaled..), stat="densityNonZero", size=1)
pgm2 <- geom_line(aes(y = -..scaled..), stat="densityNonZero", size=1)
optsm <- opts(axis.title.x = theme_blank(), axis.text.y = theme_blank(), axis.ticks.y = theme_blank() )
pa <- ggplot(dfDen, aes(x = a, colour=scenario, linetype=scenario)) + pgm + pgm2 + optsm + scale_y_continuous('a')
pb <- ggplot(dfDen, aes(x = b, colour=scenario, linetype=scenario)) + pgm + pgm2 + optsm + scale_y_continuous('b')
#pb
windows(width=7, height=3)
grid.newpage()
pushViewport( viewport(layout=grid.layout(2,1)))
print(pa , vp = viewport(layout.pos.row=1,layout.pos.col=1))
print(pb + opts(legend.position = "none") , vp = viewport(layout.pos.row=2,layout.pos.col=1))
#----------- ggplot predictive posterior
#scenarios <- c("R","RS","RSw","DG","DM")
scenarios <- c("R")
nScen <- length(scenarios)
# infer quantiles of predictions
.nSample=128
y1M <- array( NA_real_, dim=c(.nSample, length(twTwoDenEx1$obs$y1),nScen), dimnames=list(sample=NULL,iObs=NULL,scenario=scenarios) )
y2M <- array( NA_real_, dim=c(.nSample, length(twTwoDenEx1$obs$y2),nScen), dimnames=list(sample=NULL,iObs=NULL,scenario=scenarios) )
resScen <- #list( res1, res2, res2b, res3a, res3 ); names(resScen) <- scenarios
resScen <- list( mc0 ); names(resScen) <- scenarios
#scen <- "R"
for( scen in scenarios ){
ss <- stackChains(resScen[[scen]])[,-(1:getNBlocks(resScen[[scen]]))]
ssThin <- ss[round(seq(1,nrow(ss),length.out=.nSample)),]
#i <- .nSample
for( i in 1:.nSample){
pred <- twTwoDenEx1$fModel(ssThin[i,], xSparse=twTwoDenEx1$xSparse, xRich=twTwoDenEx1$xRich, thresholdCovar=thresholdCovar)
y1M[i,,scen] <- pred$y1
y2M[i,,scen] <- pred$y2
}
}
predQuantilesY1 <- lapply( scenarios, function(scen){
t(apply( y1M[,,scen],2, quantile, probs=c(0.025,0.5,0.975) ))
}); names(predQuantilesY1) <- scenarios
predQuantilesY2 <- lapply( scenarios, function(scen){
t(apply( y2M[,,scen],2, quantile, probs=c(0.025,0.5,0.975) ))
}); names(predQuantilesY2) <- scenarios
str(predQuantilesY1)
rm( y1M, y2M) # free space, we only need the summary
#str(pred1)
dfPred <- rbind(
cbind( data.frame( scenario="R", variable="y1", value = pred1$y1 ), predQuantilesY1[["R"]] )
,cbind( data.frame( scenario="R", variable="y2", value = pred1$y2 ), predQuantilesY2[["R"]] )
)
dfPred$observations <- NA_real_
dfPred$observations[dfPred$variable=="y1"] <- twTwoDenEx1$obs$y1
dfPred$observations[dfPred$variable=="y2"] <- twTwoDenEx1$obs$y2
colnames(dfPred)[ match(c("2.5%","50%","97.5%"),colnames(dfPred)) ] <- c("lower","median","upper")
str(dfPred)
p1 <- ggplot(dfPred, aes(x=value, y=observations, colour=scenario) ) +
geom_errorbarh(aes(xmax = upper, xmin = lower)) +
geom_point() +
facet_wrap( ~ variable, scales="free") +
opts(axis.title.x=theme_blank() ) +
geom_abline(colour="black") +
c()
p1
}
.tmp.f <- function(){ # does not work # AllDensities
# same as example but with each parameter updated against both densities
data(twTwoDenEx1)
thetaPrior <- twTwoDenEx1$thetaTrue
covarTheta <- diag((thetaPrior*0.3)^2)
invCovarTheta <- (thetaPrior*0.3)^2 # given as independent variances for faster calculation
thresholdCovar = 0.3 # the true value
thresholdCovar = 0 # the effective model that glosses over this threshold
# for updating each in turn against both densities
dInfos=list(
dSparse=list(fLogDen=denSparsePrior, argsFLogDen=list(thresholdCovar=thresholdCovar, twTwoDenEx=twTwoDenEx1, theta0=thetaPrior, thetaPrior=thetaPrior, invCovarTheta=invCovarTheta))
,dRich=list(fLogDen=denRich, argsFLogDen=list(thresholdCovar=thresholdCovar,twTwoDenEx=twTwoDenEx1, theta0=thetaPrior))
)
blocks = list(
bSparse=list(dInfoPos="dSparse", compPos=c("a","b") )
,bRich=list(dInfoPos="dRich", compPos=c("a","b") )
)
do.call( dInfos$dSparse$fLogDen, c(list(theta=twTwoDenEx1$theta0),dInfos$dSparse$argsFLogDen))
do.call( dInfos$dRich$fLogDen, c(list(theta=twTwoDenEx1$theta0),dInfos$dRich$argsFLogDen))
#str(twTwoDenEx1)
nObs <- c( obsSparse=length(twTwoDenEx1$obs$y1), logDen1=length(twTwoDenEx1$obs$y2) )
#trace(twDEMCSA, recover )
resPops <- res <- twDEMCSA( thetaPrior, covarTheta, dInfos=dInfos, blocks=blocks, nObs=nObs, debugSequential=TRUE, ctrlBatch=list(nSamplesBatch=64) )
# continue run for 512 generations
res <- twDEMCSA( res, 512 )
(TCurr <- getCurrentTemp(resPops))
mc0 <- concatPops(res)
logDenT <- calcTemperatedLogDen(stackChains(mc0$resLogDen), TCurr)
iBest <- getBestModelIndex( logDenT, res$dInfos )
maxLogDenT <- logDenT[iBest, ]
ss <- stackChains(mc0$parms)
(thetaBest <- ss[iBest, ])
(.qq <- apply(ss,2,quantile, probs=c(0.025,0.5,0.975) ))
plot( ss[,"a"], ss[,"b"], col=rainbow(100)[twRescale(rowSums(logDenT),c(1,100))] )
plot( ss[,"a"], ss[,"b"], col=rainbow(100)[twRescale(logDenT[,"obsSparse"],c(1,100))] )
plot( ss[,"a"], ss[,"b"], col=rainbow(100)[twRescale(logDenT[,"logDen1"],c(1,100))] )
plot( ss[,"a"], ss[,"b"], col=rgb(
twRescale(logDenT[,"obsSparse"]),0, twRescale(logDenT[,"logDen1"]) ))
apply( apply( logDenT, 2, quantile, probs=c(0.1,0.9) ),2, diff )
# density of parameters
plot( density(ss[,"a"])); abline(v=thetaPrior["a"]); abline(v=thetaBest["a"], col="blue")
plot( density(ss[,"b"])); abline(v=thetaPrior["b"]); abline(v=thetaBest["b"], col="blue")
# predictive posterior (best model only)
pred <- pred1 <- with( twTwoDenEx1, fModel(thetaBest, xSparse=xSparse, xRich=xRich) )
plot( pred$y1, twTwoDenEx1$obs$y1 ); abline(0,1)
plot( pred$y2, twTwoDenEx1$obs$y2 ); abline(0,1)
}
.tmp.f <- function(){ # # AllSparse
# same as example but with b also updated against Sparse, a only against Sparse
data(twTwoDenEx1)
thetaPrior <- twTwoDenEx1$thetaTrue
covarTheta <- diag((thetaPrior*0.3)^2)
invCovarTheta <- (thetaPrior*0.3)^2 # given as independent variances for faster calculation
thresholdCovar = 0.3 # the true value
thresholdCovar = 0 # the effective model that glosses over this threshold
# for updating each in turn against both densities
dInfos=list(
dSparse=list(fLogDen=denSparsePrior, argsFLogDen=list(thresholdCovar=thresholdCovar, twTwoDenEx=twTwoDenEx1, theta0=thetaPrior, thetaPrior=thetaPrior, invCovarTheta=invCovarTheta))
,dRich=list(fLogDen=denRich, argsFLogDen=list(thresholdCovar=thresholdCovar,twTwoDenEx=twTwoDenEx1, theta0=thetaPrior))
)
blocks = list(
bSparse=list(dInfoPos="dSparse", compPos=c("a","b") )
,bRich=list(dInfoPos="dRich", compPos=c("a") )
)
do.call( dInfos$dSparse$fLogDen, c(list(theta=twTwoDenEx1$theta0),dInfos$dSparse$argsFLogDen))
do.call( dInfos$dRich$fLogDen, c(list(theta=twTwoDenEx1$theta0),dInfos$dRich$argsFLogDen))
#str(twTwoDenEx1)
nObs <- c( obsSparse=length(twTwoDenEx1$obs$y1), logDen1=length(twTwoDenEx1$obs$y2) )
#untrace(twDEMCSA )
#trace(twDEMCSA, recover )
resPops <- res <- twDEMCSA( thetaPrior, covarTheta, dInfos=dInfos, blocks=blocks, nObs=nObs, nBatch=5, debugSequential=TRUE, TFix=c(1,NA,NA) )
resPops <- res <- twDEMCSACont( res, 256, nObs=nObs, debugSequential=TRUE
, ctrlT=list( TFix=c(parmsSparse=1) )
)
(TCurr <- getCurrentTemp(resPops))
mc0 <- concatPops(res)
logDenT <- calcTemperatedLogDen(stackChains(mc0$resLogDen), TCurr)
iBest <- getBestModelIndex( logDenT, res$dInfos )
maxLogDenT <- logDenT[iBest, ]
ss <- stackChains(mc0$parms)
(thetaBest <- ss[iBest, ])
(.qq <- apply(ss,2,quantile, probs=c(0.025,0.5,0.975) ))
plot( ss[,"a"], ss[,"b"], col=rainbow(100)[twRescale(rowSums(logDenT),c(1,100))] )
plot( ss[,"a"], ss[,"b"], col=rainbow(100)[twRescale(logDenT[,"obsSparse"],c(1,100))] )
plot( ss[,"a"], ss[,"b"], col=rainbow(100)[twRescale(logDenT[,"logDen1"],c(1,100))] )
plot( ss[,"a"], ss[,"b"], col=rgb(
twRescale(logDenT[,"obsSparse"]),0, twRescale(logDenT[,"logDen1"]) ))
apply( apply( logDenT, 2, quantile, probs=c(0.1,0.9) ),2, diff )
# density of parameters
plot( density(ss[,"a"])); abline(v=thetaPrior["a"]); abline(v=thetaBest["a"], col="blue")
plot( density(ss[,"b"])); abline(v=thetaPrior["b"]); abline(v=thetaBest["b"], col="blue")
# predictive posterior (best model only)
pred <- pred1 <- with( twTwoDenEx1, fModel(thetaBest, xSparse=xSparse, xRich=xRich) )
plot( pred$y1, twTwoDenEx1$obs$y1 ); abline(0,1)
plot( pred$y2, twTwoDenEx1$obs$y2 ); abline(0,1)
}
.tmp.f <- function(){ # does not work # .tmp.2DenSwitched
# same as example but with each parameter b updated against Sparse
data(twTwoDenEx1)
thetaPrior <- twTwoDenEx1$thetaTrue
covarTheta <- diag((thetaPrior*0.3)^2)
invCovarTheta <- (thetaPrior*0.3)^2 # given as independent variances for faster calculation
thresholdCovar = 0.3 # the true value
thresholdCovar = 0 # the effective model that glosses over this threshold
# for updating each in turn against both densities
dInfos=list(
dSparse=list(fLogDen=denSparsePrior, argsFLogDen=list(thresholdCovar=thresholdCovar, twTwoDenEx=twTwoDenEx1, theta0=thetaPrior, thetaPrior=thetaPrior, invCovarTheta=invCovarTheta))
,dRich=list(fLogDen=denRich, argsFLogDen=list(thresholdCovar=thresholdCovar,twTwoDenEx=twTwoDenEx1, theta0=thetaPrior))
)
blocks = list(
bSparse=list(dInfoPos="dSparse", compPos=c("b") )
,bRich=list(dInfoPos="dRich", compPos=c("a") )
)
do.call( dInfos$dSparse$fLogDen, c(list(theta=twTwoDenEx1$theta0),dInfos$dSparse$argsFLogDen))
do.call( dInfos$dRich$fLogDen, c(list(theta=twTwoDenEx1$theta0),dInfos$dRich$argsFLogDen))
#str(twTwoDenEx1)
nObs <- c( obsSparse=length(twTwoDenEx1$obs$y1), logDen1=length(twTwoDenEx1$obs$y2) )
#untrace(twDEMCSA )
#trace(twDEMCSA, recover )
resPops <- res <- twDEMCSA( thetaPrior, covarTheta, dInfos=dInfos, blocks=blocks, nObs=nObs, nBatch=5, debugSequential=TRUE )
(TCurr <- getCurrentTemp(resPops))
mc0 <- concatPops(res)
logDenT <- calcTemperatedLogDen(stackChains(mc0$resLogDen), TCurr)
iBest <- getBestModelIndex( logDenT, res$dInfos )
maxLogDenT <- logDenT[iBest, ]
ss <- stackChains(mc0$parms)
(thetaBest <- ss[iBest, ])
(.qq <- apply(ss,2,quantile, probs=c(0.025,0.5,0.975) ))
plot( ss[,"a"], ss[,"b"], col=rainbow(100)[twRescale(rowSums(logDenT),c(1,100))] )
plot( ss[,"a"], ss[,"b"], col=rainbow(100)[twRescale(logDenT[,"obsSparse"],c(1,100))] )
plot( ss[,"a"], ss[,"b"], col=rainbow(100)[twRescale(logDenT[,"logDen1"],c(1,100))] )
plot( ss[,"a"], ss[,"b"], col=rgb(
twRescale(logDenT[,"obsSparse"]),0, twRescale(logDenT[,"logDen1"]) ))
apply( apply( logDenT, 2, quantile, probs=c(0.1,0.9) ),2, diff )
# density of parameters
plot( density(ss[,"a"])); abline(v=thetaPrior["a"]); abline(v=thetaBest["a"], col="blue")
plot( density(ss[,"b"])); abline(v=thetaPrior["b"]); abline(v=thetaBest["b"], col="blue")
# predictive posterior (best model only)
pred <- pred1 <- with( twTwoDenEx1, fModel(thetaBest, xSparse=xSparse, xRich=xRich) )
plot( pred$y1, twTwoDenEx1$obs$y1 ); abline(0,1)
plot( pred$y2, twTwoDenEx1$obs$y2 ); abline(0,1)
}
.tmp.f <- function(){ # oneDensity_andCluster
# same as example but with only one combined density for both parameters
data(twTwoDenEx1)
thetaPrior <- twTwoDenEx1$thetaTrue
covarTheta <- diag((thetaPrior*0.3)^2)
invCovarTheta <- (thetaPrior*0.3)^2 # given as independent variances for faster calculation
thresholdCovar = 0.3 # the true value
thresholdCovar = 0 # the effective model that glosses over this threshold
# for only one logDensity - Temperature of the strongest component goes to zero and other increase according to mismatch
dInfos=list( list(fLogDen=denBoth, argsFLogDen=list(thresholdCovar=thresholdCovar, twTwoDenEx=twTwoDenEx1, theta0=thetaPrior, thetaPrior=thetaPrior, invCovarTheta=invCovarTheta)) )
do.call( dInfos[[1]]$fLogDen, c(list(theta=twTwoDenEx1$theta0),dInfos[[1]]$argsFLogDen))
#str(twTwoDenEx1)
nObs <- c( y1=length(twTwoDenEx1$obs$y1), y2=length(twTwoDenEx1$obs$y2) )
#trace(twDEMCSA, recover)
argsTwDEMCSA <- list( thetaPrior=thetaPrior, covarTheta=covarTheta, dInfos=dInfos, nObs=nObs
,TFix = c(1,NA,NA)
,nGen=256
#,TMax = c(NA,1.2,NA) # do not increase Sparse observations again too much
, nBatch=3
, debugSequential=TRUE
)
resPops <- res <- do.call( twDEMCSA, argsTwDEMCSA )
.tmp.f <- function(){ #continueRun
res$pops[[2]]$spaceInd <- 4 # test spaces different from 1:nSpace
res2 <- twDEMCSA(res)
}
.tmp.f <- function(){ #byCluster
runClusterParms <- list(
fSetupCluster = function(){library(twDEMC)}
,fRun = twDEMCSA
,argsFRun = within( argsTwDEMCSA, debugSequential<-FALSE )
)
save(runClusterParms, file=file.path("..","..","projects","asom","parms","saOneDensity.RData"))
# from asom directory run:
# R CMD BATCH --vanilla '--args paramFile="parms/saOneDensity.RData"' runCluster.R Rout/runCluster_0.rout
# or from libra home directory run
# ./bsubr_i.sh runCluster.R iproc=0 nprocSinge=1 'paramFile="parms/saOneDensity.RData"'
# or
# bsub -q SLES -n 4 ./bsubr_i.sh runCluster.R nprocSingle=4 'paramFile="parms/saOneDensity.RData"'
load("parms/res_saOneDensity_1.RData")
}
(TCurr <- getCurrentTemp(res))
mc0 <- concatPops(res)
logDenT <- calcTemperatedLogDen(stackChains(mc0$resLogDen), TCurr)
iBest <- getBestModelIndex( logDenT, res$dInfos )
maxLogDenT <- logDenT[iBest, ]
ss <- stackChains(mc0$parms)
(thetaBest <- ss[iBest, ])
(.qq <- apply(ss,2,quantile, probs=c(0.025,0.5,0.975) ))
plot( ss[,"a"], ss[,"b"], col=rainbow(100)[twRescale(rowSums(logDenT),c(1,100))] )
plot( ss[,"a"], ss[,"b"], col=rainbow(100)[twRescale(logDenT[,"y1"],c(10,100))] )
plot( ss[,"a"], ss[,"b"], col=rainbow(100)[twRescale(logDenT[,"y2"],c(10,100))] )
plot( ss[,"a"], ss[,"b"], col=rgb(
twRescale(logDenT[,"y1"]),0, twRescale(logDenT[,"y2"]) ))
apply( apply( logDenT, 2, quantile, probs=c(0.2,0.8) ),2, diff )
# density of parameters
plot( density(ss[,"a"])); abline(v=thetaPrior["a"]); abline(v=thetaBest["a"], col="blue")
plot( density(ss[,"b"])); abline(v=thetaPrior["b"]); abline(v=thetaBest["b"], col="blue")
# predictive posterior (best model only)
pred <- pred1 <- with( twTwoDenEx1, fModel(thetaBest, xSparse=xSparse, xRich=xRich) )
plot( pred$y1, twTwoDenEx1$obs$y1 ); abline(0,1)
plot( pred$y2, twTwoDenEx1$obs$y2 ); abline(0,1)
}
isLogDenDrift <- function(
### check whether first quartile all the logDensities is significantly smaller than last quartile
logDenT ##<< numeric array (nStep x nResComp): logDensity (highest are best)
, dInfos ##<< list of lists with entry resCompPos (integer vector) specifying the position of result components for each density
, alpha=0.05 ##<< the significance level for a difference
, maxDrift=0.3 ##<< difference in LogDensity, below which no drift is signalled
){
##details<<
## Because of large sample sizes, very small differences may be significantly different.
## Use argument minDiff to specify below which difference a significant difference is not regarded as drift.
nr <- nrow(logDenT)
nr4 <- nr %/% 4
ds1 <- logDenT[1:nr4, ,drop=FALSE] # first quater of the dataset
ds4 <- logDenT[(nr+1-nr4):nr, ,drop=FALSE] # last quater of the dataset
#iDen <- 1
tres <- sapply( seq_along(dInfos), function(iDen){
dInfo <- dInfos[[iDen]]
rs1 <- rowSums(ds1[,dInfo$resCompPos ,drop=FALSE])
rs4 <- rowSums(ds4[,dInfo$resCompPos ,drop=FALSE])
resTTest <- try(t.test(rs1,rs4,"less"), silent=TRUE)
if( inherits(resTTest,"try-error") ){ # logDen essentially constant
#bo <- (diff(resTTest$estimate) >= maxDrift ) && (resTTest$p.value <= alpha)
c(rs1,rs4,1)
} else {
c( resTTest$estimate, p=resTTest$p.value)
}
})
#tresi <- tres[,1]
boL <- apply( tres, 2, function(tresi){
(diff(tresi[1:2]) >= maxDrift ) && (tresi["p"] <= alpha)
})
ret <- any( boL )
attr( ret, "resTTest") <- tres
attr( ret, "logDenComp") <- colMeans(ds4) # average logDensity of all components of last quarter
##value<< TRUE if any of the logDensities are a significantly greater in the fourth quantile compared to the first quantile of the samples
## , attribute \code{resTTest}: numeric maxtrix (logDenStart, logDenEnd, pTTest x nDInfo )
## , attribute \code{logDenComp}: numeric vector (nResComp): average logDenT of all components of last quater
ret
}
attr(isLogDenDrift,"ex") <- function(){
data(twdemcEx1)
logDenT <- calcTemperatedLogDen( twdemcEx1, getCurrentTemp(twdemcEx1))
#undebug(isLogDenDrift)
isLogDenDrift(logDenT, twdemcEx1$dInfos )
}
twRunDEMCSA <- function(
### wrapper to twDEMCSA compliant to runCluster.R
argsTwDEMCSA
### Arguments passed to twDEMCSA -> twDEMCBlockInt
### It is updated by \dots.
,... ##<< further arguments passed to twDEMCSA -> twDEMCBlockInt
,prevResRunCluster=NULL ##<< results of call to twRunDEMCSA, argument required to be called from runCluster.R
,restartFilename=NULL ##<< name of the file to store restart information, argument required to be called from runCluster.R
){
require(twDEMC)
#update argsDEMC to args given
argsDEMC <- argsTwDEMCSA
# update the restartFilename argument
.dots <- list(...)
argsDEMC[ names(.dots) ] <- .dots
#for( argName in names(.dots) ) argsDEMC[argName] <- .dots[argName]
if( 0 != length(restartFilename)) argsDEMC$restartFilename <- restartFilename
# do the actual call
res <- do.call( twDEMCSA, argsDEMC )
}
.tmp.f <- function(){ # testRestart
rFileName <- file.path("tmp","example_twDEMCSA.RData")
load(rFileName) # resRestart.twDEMCSA
#untrace(twDEMCSACont)
#trace(twDEMCSACont, recover )
res1 <- do.call( twDEMCSACont, c( list(mc=resRestart.twDEMCSA), resRestart.twDEMCSA$args ))
}
completeResCompVec <- function(
### given a vector with names, make a vector corresponding to resComp with components, not yet given filled by NA
x ##<< named vector to be completed
,resCompNames ##<< names of the result components
){
xAll <- x
nResComp <- length(resCompNames)
if( 0 == length( x) )
xAll <- structure( rep( NA_real_, nResComp), names=resCompNames )
#else if( nResComp != length( x) ){ # twutz 141120: can give a different names of same length, must correct
else{
iFix <- sapply( names(x), match, resCompNames )
if( any(is.na(iFix))){
misNames <- paste(names(x)[ is.na(iFix)], collapse=",")
stop(paste("completeResCompVec: x must be of the same length as resCompNames or all its names must correspond names of resCompNames. Unmatched names:",misNames))
}
xAll <- structure( rep( NA_real_, nResComp), names=resCompNames )
xfin <- x[ is.finite(x) ]
#iName <- names(xfin)[1]
for( iName in names(xfin) ){
xAll[ names(xAll)== iName ] <- xfin[iName]
}
}
xAll
### vector with names resCompNames, with corresponding values of x set, others NA
}
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.
twDEMCSA <- function(
### Simulated annealing DEMC
thetaPrior ##<< vector of parameters, point estimate
##, alternatively array with initial states, as returned by \code{\link{initZtwDEMCNormal}}
,covarTheta ##<< the a prior covariance of parameters, see \code{link{initZtwDEMCNormal}}
,nGen=512 ##<< number of generations, i.e. Monte Carlo steps
,nObs ##<< integer vector specifying the number of observations for each result component
## for each name in result of logDensity functions there must be an entry in nObs
,... ##<< further argument to \code{\link{twDEMCBlockInt}}
, dInfos ##<< Information on density functions (see \code{\link{twDEMCBlockInt}})
, m0 = calcM0twDEMC(length(thetaPrior),nChainPop) ##<< minimum number of samples in step for extending runs
, controlTwDEMC = list() ##<< control parameters with entry \code{thin} (see \code{\link{twDEMCBlockInt}})
, debugSequential=FALSE ##<< set to TRUE to avoid parallel execution, good for debugging
, restartFilename=NULL ##<< filename to write intermediate results to
, remoteDumpfileBasename=NULL ##<< fileBasename to write dumps to on error
#
,nChainPop=4 ##<< number of chains within population
,nPop=2 ##<< number of populations
,doIncludePrior=FALSE ##<< should the prior be part of initial population
##<< Recommendation to set to false, because if the TRUE parameter is in initial set, initial Temperature calculation is too optimistic
#
, ctrlBatch = list( ##<< list of arguments controlling batch executions, see also \code{\link{twDEMCSACont}}
##describe<<
nGen0=m0*controlTwDEMC$thin*3 ##<< number of generations for the initial batch
#,useSubspaceAdaptation=FALSE # #<< if TRUE then overall space is devided and each subspace is explored with locally adapted DEMC, see \code{\link{divideTwDEMCSACont}}
##end<<
)
, ctrlT = list( ##<< list of arguments controlling Temperature decrease, see also \code{\link{twDEMCSACont}}
##describe<<
qTempInit=0.4 ##<< quantile of logDensities used to calculate initial beginning and end temperature, with default 0.4: 40% of the space is accepted
,TBaseInit=NULL ##<< numeric scalar: initial base temperature. If given this is used for calculating initial temperature
,isVerbose=FALSE ##<< boolean scalar: set to TRUE to report stream base temperatures during batches
,TEndFixed=NULL ##<< set to a scalar end temperature, e.g. in order to decrease temperatue to a given temperature
,qBest=0.1 ##<< the proportion of best samples to base calculation of target temperature on
##end<<
)
, ctrlConvergence = list() ##<< list or arguments controlling check for convergence, see \code{\link{twDEMCSACont}}
#, ctrlSubspaces = list() # #<< list of arguments controlling splitting and merging of subspaces, see \code{\link{divideTwDEMCSACont}}
){
##detail<< \describe{\item{Continuing a previous run}{
## When supplying the first argument an object of class twDEMCPops to twDEMCSA, then this run is extended.
## All other parameters, except nGen, are then are ignored.
## In order to change parameters, modify list entry args in the twDEMCPops object.
##}}
if( inherits(thetaPrior,"twDEMCPops") && 0 != length(thetaPrior$args) ){
if( !missing(debugSequential) )
thetaPrior$args$debugSequential=debugSequential
if( !missing(restartFilename) )
thetaPrior$args$restartFilename=restartFilename
if( !missing(remoteDumpfileBasename) )
thetaPrior$args$remoteDumpfileBasename=remoteDumpfileBasename
if( !missing(controlTwDEMC) )
thetaPrior$args$controlTwDEMC <- twMergeLists( thetaPrior$args$controlTwDEMC, controlTwDEMC )
if( !missing(ctrlBatch) )
thetaPrior$args$ctrlBatch <- twMergeLists( thetaPrior$args$ctrlBatch, ctrlBatch )
if( !missing(ctrlT) )
thetaPrior$args$ctrlT <- twMergeLists( thetaPrior$args$ctrlT, ctrlT )
if( !missing(ctrlConvergence) )
thetaPrior$args$ctrlConvergence <- twMergeLists( thetaPrior$args$ctrlConvergence, ctrlConvergence )
#if( !missing(ctrlSubspaces) )
# thetaPrior$args$ctrlSubspaces <- twMergeLists( thetaPrior$args$ctrlSubspaces, ctrlSubspaces )
if( !missing(dInfos) ){
warning("twDEMCSA: continueing previous run, but with different argument dInfos")
thetaPrior$args$dInfos <- dInfos # overwrite
}
thetaPrior$failureMsg <- NULL # clear previous return and error messages
thetaPrior$finishEarly <- NULL
#ret <- if( 0!=length(thetaPrior$args$ctrlBatch) &&
# 0!=length(thetaPrior$args$ctrlBatch$useSubspaceAdaptation) &&
# thetaPrior$args$ctrlBatch$useSubspaceAdaptation
# )
# ret <- do.call(divideTwDEMCSACont,c(list(thetaPrior, nGen=nGen),thetaPrior$args))
#else
mcPops <- do.call(twDEMCSACont,c(list(thetaPrior, nGen=nGen),thetaPrior$args))
#ret$args$ctrlBatch$useSubspaceAdaptation <- ctrlBatch$useSubspaceAdaptation
return(mcPops)
}
#
#-------------- generate initial sample
# check for array before evaluating formal arguments, because m0 then changes, which is used in default values
#mtrace(initZtwDEMCNormal)
Zinit0 <- if( is.array(thetaPrior) ){
.dimZ <- dim(thetaPrior)
if( length(.dimZ) != 3 ) stop("twDEMCSA: wrong dimensions of thetaPrior")
nChain <- .dimZ[3]
if( missing(nChainPop) ) nChainPop <- nChain %/% nPop
if( missing(nPop) ) nPop <- nChain %/% nChainPop
if( .dimZ[3] != nChainPop*nPop ) stop("twDEMCSA: wrong number of chains in thetaPrior")
m0 = calcM0twDEMC( ncol(thetaPrior) ,nChainPop) ##<< minimum number of samples in step for extending runs
if( .dimZ[1] < m0 ) stop(paste("twDEMCSA: too few cases in thetaPrior (",.dimZ[1],"). Need at least m0=",m0,sep=""))
thetaPrior
} else initZtwDEMCNormal( thetaPrior, covarTheta, nChainPop=nChainPop, nPop=nPop, doIncludePrior=doIncludePrior, m0=m0)
#
#-------- fill in default argument values
if( is.null(controlTwDEMC$thin) ) controlTwDEMC$thin <- 4
thin <- controlTwDEMC$thin
# The user may provide list arguments to overwrite some default entries.
# With formals and twMergeLists, we make sure to take the other default entries in the lists
frm <- formals()
#ctrlSubspaces <- if( hasArg(ctrlSubspaces) ) twMergeLists( eval(frm[["ctrlConvergence"]]), ctrlSubspaces ) else ctrlSubspaces
ctrlBatch <- if( hasArg(ctrlBatch) ) twMergeLists( eval(frm[["ctrlBatch"]]), ctrlBatch ) else ctrlBatch
ctrlT <- if( hasArg(ctrlT) ) twMergeLists( eval(frm[["ctrlT"]]), ctrlT ) else ctrlT
#
#---------- calculate logDensity of initial components
ss <- stackChains(Zinit0)
#print("before twCalcLogDensPar"); recover()
cat("Calculating logDensity for ", length(ss)," initial states for m0=",m0,"rows per population.\n")
logDenDS0 <- (tmp <- twCalcLogDensPar(dInfos, ss, remoteDumpfileBasename=remoteDumpfileBasename))$logDenComp
boFinite <- apply(is.finite(logDenDS0), 1, all )
m0FiniteFac <- sum(boFinite) / nrow(logDenDS0)
# if too few initial states have yielded finite results of logDensity functions, add some more initial states
if( m0FiniteFac == 1){
Zinit <- Zinit0
logDenDS <- logDenDS0
}else if( m0FiniteFac < 0.05 ){
stop(paste("twDEMCSA: less than 5% finite solutions in initial exploration (",m0FiniteFac,"). Check setup."),sep="")
}else if( m0FiniteFac > 0.9 ){
sumLogDen0 <- rowSums(logDenDS0)
Zinit <- replaceZinitNonFiniteLogDens( Zinit0, sumLogDen0 )
logDenDS <- logDenDS0
logDenDS[!boFinite,] <- logDenDS[which.min(sumLogDen0),]
}else{
# generate more proposals and concatenate finite cases from both
Zinit1 <- initZtwDEMCNormal( thetaPrior, covarTheta, nChainPop=nChainPop, nPop=nPop, doIncludePrior=doIncludePrior
,m0FiniteFac=min(1,m0FiniteFac)
)
ss1 <- stackChains(Zinit1)
logDenL <- lapply( dInfos, function(dInfo){
resLogDen <- twCalcLogDenPar( dInfo$fLogDen, ss1, argsFLogDen=dInfo$argsFLogDen)$logDenComp
})
# replace missing cases
logDenDS1 <- abind( logDenL )
# logDenDS1[ 1:round(nrow(logDenDS1)/1.7),1] <- NA # for testing replacement of non-finite cases
boFinite1 <- apply(is.finite(logDenDS1), 1, all )
m0FiniteFac1 <- sum(boFinite1) / nrow(logDenDS1)
ss12 <- rbind( ss[boFinite, ,drop=FALSE], ss1[boFinite1, ,drop=FALSE] )
logDenDS12 <- rbind( logDenDS0[boFinite, ,drop=FALSE], logDenDS1[boFinite1, ,drop=FALSE])
nDiff <- nrow(ss) - nrow(ss12)
if( nDiff/nrow(ss) > 0.1 ){
stop("twDEMCSA: too many states yielding non-finite logDensities.")
}else if( nDiff > 0 ){
# duplicate missings to refill
iSample <- sample( nrow(ss12), size=nDiff )
ss12 <- rbind( ss12, ss12[ iSample, ] )
logDenDS12 <- rbind( logDenDS12, logDenDS12[iSample,] )
}else if(nDiff < 0){
iSample <- sample( nrow(ss12), size=nrow(ss))
ss12 <- ss12[ iSample, ]
logDenDS12 <- logDenDS12[iSample,]
}
## reshape to chains
ncolZ <- dim(Zinit0)[3]
tmp <- matrix(1:nrow(ss12), ncol=ncolZ)
Zinit <- abind( lapply( 1:ncolZ, function(i){ ss12[ tmp[,i], ,drop=FALSE]}), rev.along=0 )
logDenDS <- logDenDS12
} # end generating nonfinite Zinit
#
# apply overfitting correction?
# No need here because Temperature would be set to 1 for corrected and noncorrected values
#
# now the length of resComp is known (colnames(logDenDS)), so check the TFix and TMax parameters
iNoNames <- which( "" == colnames(logDenDS) )
if( length(iNoNames) ) stop("Encountered unnamed logDensitiy components. Provide names of each return component in fLogDen.")
nResComp <- ncol(logDenDS)
ctrlT$TFix <- completeResCompVec( ctrlT$TFix, colnames(logDenDS) )
iFixTemp <- which( is.finite(ctrlT$TFix) )
iNonFixTemp <- which( !is.finite(ctrlT$TFix) )
ctrlT$TMax <- tmp <- completeResCompVec( ctrlT$TMax, colnames(logDenDS) )
iMaxTemp <- which( is.finite(ctrlT$TMax) )
#
# check nObs names and set all TFix to 1
if( 0 == length(names(nObs))){
if( length(nObs) == length(iNonFixTemp) ) names(nObs) <- colnames(logDenDS)[iNonFixTemp] else
stop("Provide names with nObs, or provide one component for each temperated (i.e. not in TFix) fLogDen return component.")
}
iMissing <- which( !(colnames(logDenDS)[iNonFixTemp] %in% names(nObs)) )
if( length(iMissing) ) stop("Missing the following components in nObs: ",paste( colnames(logDenDS)[iNonFixTemp][iMissing],sep=",") )
# make nObs correspond to each data stream, set TFix Components to 1
nObsComp <- structure( rep( NA_real_, ncol(logDenDS) ), names=colnames(logDenDS) )
#nObsComp[names(nObs)] <- nObs # if there are seveval occurence of the same name, others still NA
# name <- names(nObs)[1]
for( name in names(nObs) ){ nObsComp[ names(nObsComp) == name ] <- nObs[name] }
nObs <- nObsComp
#
#if( 0 == length( TMaxInc) ) TMaxInc <- structure( rep( NA_real_, nResComp), names=colnames(logDenDS) )
#if( nResComp != length( TMaxInc) ) stop("twDEMCSA: TMaxInc must be of the same length as number of result Components.")
#iMaxIncTemp <- which( is.finite(TMaxInc) )
#
#------ intial temperatures:
#sapply( seq_along(expLogDenBest), function(i){ max(logDenDS[,i]) })
##details<< \describe{\item{Initial temperature}{
## Initial parameters are ranked according to their maximum log-Density across components.
## The parameters and logDensity results at rank position defined by argument \code{ctrlT$qTempInit} is selected.
## The stream temperatures are inferred by deviding logDensity components by the number of observations.
## From these, a common base temperatue is calculated by (see \code{\link{calcBaseTemp}}) and rescaled to stream temperatures by \code{\link{calcStreamTemp}}.
##}}
#print("twDEMCSA: before calculating initial temperature"); recover()
if( length(ctrlT$TBaseInit) ){
temp0 <- calcStreamTemp(ctrlT$TBaseInit, nObs, TFix=ctrlT$TFix, iFixTemp=iFixTemp)
} else {
iTmax <- 5
T0s <- rep(NA_real_, iTmax )
T <- nObs
T[iFixTemp] <- ctrlT$TFix[iFixTemp]
T0s[1] <- .Machine$double.xmax
relChange <- relChangeReq <- 0.05 # if below 5% change then no need to iterate further
iT <- 1
while( (iT < iTmax) && (abs(relChange) >= relChangeReq) ){
iT <- iT + 1
resLogDenT <- calcTemperatedLogDen( logDenDS, T )
iBest <- getBestModelIndices( resLogDenT, dInfos, prob=ctrlT$qTempInit )
T0 <- T0s[iT] <- calcBaseTempSk( -2*logDenDS[iBest, ,drop=FALSE], nObs, ctrlT$TFix, iFixTemp = iFixTemp, iNonFixTemp = iNonFixTemp)
T <- calcStreamTemp(T0, nObs, ctrlT$TFix, iFixTemp = iFixTemp)
relChange <- (T0 - T0s[iT-1])/T0
}
temp0 <- T
if( !all(is.finite(temp0)) ) stop("twDEMCSA: encountered non-finite Temperatures.")
}
ctrlT$TMax[iNonFixTemp] <- pmin(ifelse(is.finite(ctrlT$TMax),ctrlT$TMax, Inf), ifelse(is.finite(temp0),temp0,Inf) )[iNonFixTemp] # decrease TMax
#if( any(temp0 > 8)) stop("twDEMCSA: encountered too high temperature.")
print(paste("initial T=",paste(signif(temp0,2),collapse=",")," for ",min(nGen, ctrlBatch$nGen0)," generations"," ", date(), sep="") )
#if( sum(temp0 > 1) == 0){ dump.frames("parms/debugDump",TRUE); stop("twDEMCSA: no temperature > 0") }
#
mcPops <- res0 <- res <- twDEMCBlock( Zinit
, nGen=min(nGen, ctrlBatch$nGen0)
#,nGen=16, debugSequential=TRUE
, dInfos=dInfos
, TSpec=cbind( T0=temp0, TEnd=temp0 )
, nPop=nPop
, m0=m0, controlTwDEMC=controlTwDEMC, debugSequential=debugSequential
# XXTODO: replace lines below later on by ...
#, blocks = blocks
,...
)
print(paste("finished initial ",ctrlBatch$nGen0," out of ",nGen," gens. ", date(), sep="") )
#
#if( nGen > ctrlBatch$nGen0 ){ #depre, now need to call Continue to return proper diagnostics
#if( TRUE ){ # need to call Continue to return proper diagnostics
#nObsLocal <- nObs
#ret <- if( ctrlBatch$useSubspaceAdaptation ){
# ret <- divideTwDEMCSACont( mc=res0, nGen=nGen-ctrlBatch$nGen0, nObs=nObs
# , m0=m0, controlTwDEMC=controlTwDEMC, debugSequential=debugSequential, restartFilename=restartFilename
# , ctrlBatch=ctrlBatch, ctrlT=ctrlT, ctrlConvergence=ctrlConvergence, ctrlSubspaces=ctrlSubspaces
# , ...
# )
#}else {
mcPops <- twDEMCSACont( mc=res0, nGen=nGen-ctrlBatch$nGen0, nObs=nObs
, m0=m0, controlTwDEMC=controlTwDEMC, debugSequential=debugSequential, restartFilename=restartFilename
, ctrlBatch=ctrlBatch, ctrlT=ctrlT, ctrlConvergence=ctrlConvergence
, ...
)
#}
#}
#ret$args$ctrlBatch$useSubspaceAdaptation <- ctrlBatch$useSubspaceAdaptation # remember this argument
##value<< An object of class \code{twDEMCPops} as described in \code{\link{twDEMCBlockInt}}.
## See \code{\link{subset.twDEMCPops}} for processing further handling of this class.
mcPops
}
attr(twDEMCSA,"ex") <- function(){
# Best to have a look at the the package vignettes
#--------------- single density ----------------------------
# We will use logDenGaussian as logDensity function that compares a simple linear model with observations.
data(twLinreg1)
# collect all the arguments to the logDensity in a list (except the first argument of changing parameters)
argsFLogDen <- with( twLinreg1, 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=twLinreg1$theta0),argsFLogDen))
do.call( logDenGaussian, c(list(theta=twLinreg1$thetaTrue),argsFLogDen)) # slightly largere misfit than nObs/2=15, underestimated sdObs
.nGen=200
.nPop=2
mcPops <- twDEMCSA(
theta=twLinreg1$theta0, covarTheta=diag(twLinreg1$sdTheta^2) # to generate initial population
, nGen=.nGen
, dInfos=list(den1=list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen))
, nPop=.nPop # number of independent populations
, controlTwDEMC=list(thin=4) # see twDEMCBlockInt for all the tuning options
, ctrlConvergence=list(maxRelTChangeCrit=0.1) # ok if T changes less than 10%
, ctrlT=list(TFix=c(parms=1)) # do not use increased temperature for priors
, nObs=c(obs=length(argsFLogDen$obs)) # number of observations used in temperature calculation
)
#mcp <- twDEMCSA( mcp, nGen=2000)
mcPops <- twDEMCSA( mcPops, nGen=400) # continue run
rescoda <- as.mcmc.list(mcPops)
plot(rescoda, smooth=FALSE)
mcChains1 <- concatPops(mcPops) # array representation instead of list of pops, last dim is the chain
mcChains2 <- concatPops(stackChainsPop(mcPops)) # combining dependent chains within one population
mcChains3 <- concatPops(subsetTail(mcPops,0.5)) # take only the last part of the chains
c(getNGen(mcChains1), getNGen(mcChains3))
plot(as.mcmc.list(mcChains2), smooth=FALSE)
#--------------- multiple densities -------------------------
data(twTwoDenEx1)
thetaPrior <- twTwoDenEx1$thetaTrue
covarTheta <- diag((thetaPrior*0.3)^2)
invCovarTheta <- (thetaPrior*0.3)^2 # given as independent variances for faster calculation
thresholdCovar = 0.3 # the true value used to generate the observations
thresholdCovar = 0 # the effective model that glosses over this threshold
#str(twTwoDenEx1)
nObs <- c( obsSparse=length(twTwoDenEx1$obs$y1), obsRich=length(twTwoDenEx1$obs$y2) )
dInfos=list(
dSparse=list(fLogDen=denSparsePrior
, argsFLogDen=list(thresholdCovar=thresholdCovar, twTwoDenEx=twTwoDenEx1, theta0=thetaPrior, thetaPrior=thetaPrior, invCovarTheta=invCovarTheta)
#, maxLogDen=-1/2*nObs[c("parmsSparse","obsSparse")] # control overfitting
)
,dRich=list(fLogDen=denRichPrior
, argsFLogDen=list(thresholdCovar=thresholdCovar,twTwoDenEx=twTwoDenEx1, theta0=thetaPrior, thetaPrior=thetaPrior, invCovarTheta=invCovarTheta)
, maxLogDen=c(parmsRich=0,-1/2*nObs["obsRich"]) # control overfitting for rich datastream
)
)
blocks = list(
a=list(dInfoPos="dSparse", compPos="a")
,b=list(dInfoPos="dRich", compPos="b")
)
names(do.call( dInfos$dSparse$fLogDen, c(list(theta=twTwoDenEx1$theta0),dInfos$dSparse$argsFLogDen)))
names(do.call( dInfos$dRich$fLogDen, c(list(theta=twTwoDenEx1$theta0),dInfos$dRich$argsFLogDen)))
#trace(twDEMCSACont, recover )
#trace(twDEMCSA, recover )
res <- res0 <- twDEMCSA( thetaPrior, covarTheta, dInfos=dInfos, blocks=blocks, nObs=nObs
, nGen=3*256
, ctrlT=list( TFix=c(parmsSparse=1,parmsRich=1) ) # no increased Temperature for priors
, ctrlBatch=list( nGenBatch=256 )
, debugSequential=TRUE
, controlTwDEMC = list(
DRgamma=0.1 # use Delayed rejection
#,controlOverfittingMinNObs = 20 # use overfitting control (for obsRich), recommended on using single density
)
#, restartFilename=file.path("tmp","example_twDEMCSA.RData")
)
res <- twDEMCSA( res0, nGen=2*256 ) # extend the former run
(TCurr <- getCurrentTemp(res))
# Note that T does decrease to 1
# This accounts for structural model mismatch in addition to observation uncertinaty
mc0 <- concatPops(res)
mcE <- concatPops(subsetTail(res,0.2)) # only the last 20%
plot( as.mcmc.list(mc0) , smooth=FALSE )
matplot( mc0$temp, type="l" )
logDenT <- calcTemperatedLogDen(stackChains(mcE$resLogDen), TCurr)
iBest <- getBestModelIndex( logDenT, res$dInfos )
maxLogDenT <- logDenT[iBest, ]
ss <- stackChains(mcE$parms)
(thetaBest <- ss[iBest, ])
twTwoDenEx1$thetaTrue
(.qq <- apply(ss,2,quantile, probs=c(0.025,0.5,0.975) )) # model error really in paramter b
plot( ss[,"a"], ss[,"b"], col=rainbow(100)[twMisc::twRescale(rowSums(logDenT),c(10,100))] )
plot( ss[,"a"], ss[,"b"], col=rainbow(100)[twMisc::twRescale(logDenT[,"obsSparse"],c(10,100))] )
plot( ss[,"a"], ss[,"b"], col=rainbow(100)[twMisc::twRescale(logDenT[,"obsRich"],c(10,100))] )
plot( ss[,"a"], ss[,"b"], col=rgb(
twMisc::twRescale(logDenT[,"obsSparse"]),0, twMisc::twRescale(logDenT[,"obsRich"]) ))
apply( apply( logDenT, 2, quantile, probs=c(0.1,0.9) ),2, diff )
# density of parameters
plot( density(ss[,"a"])); abline(v=thetaPrior["a"]); abline(v=thetaBest["a"], col="blue")
plot( density(ss[,"b"])); abline(v=thetaPrior["b"]); abline(v=thetaBest["b"], col="blue")
# predictive posterior (best model only)
pred <- pred1 <- with( twTwoDenEx1, fModel(thetaBest, xSparse=xSparse, xRich=xRich) )
plot( pred$y1, twTwoDenEx1$obs$y1 ); abline(0,1)
plot( pred$y2, twTwoDenEx1$obs$y2 ); abline(0,1) # note that deviation is now in y2 - consistent with the introduced bias
} # end examples twDEMCSA
#plot predictive posterior
.tmp.f <- function(){
ss <- stackChains(d3$parms)
apply( ss, 2, mean )
apply( ss, 2, sd )
nSamp <- 40
pSamp <- ss[ sample.int(nSamp), ]
pPred <- apply( pSamp, 1, fModel, xval=xval )
statPred <- t(apply( pPred, 1, function(obsi){ c(mean(obsi), sd(obsi) )} ))
plot(statPred[,1] ~ xval, ylim=range(statPred[,1] +c(1.96,-1.96)*statPred[,2]))
lines( statPred[,1] +1.96*statPred[,2] ~ xval )
lines( statPred[,1] -1.96*statPred[,2] ~ xval )
points( obs ~ xval, col="maroon")
w <- 1/sdObs^2/sum(1/sdObs^2)
lm2 <- lm( obs ~ xval, weights=w )
w <- 1/sdObs^2 # would more correspond to real case see ?lm, but summed to 1 internally
lm2 <- lm.wfit( cbind(1,xval), obs, w=w )
summary(lm2)
obsTrue <- fModel( thetaTrue, xval)
plot(statPred[,1] ~ obsTrue, ylim=range(statPred[,1] +c(1.96,-1.96)*statPred[,2]))
lines( statPred[,1] +1.96*statPred[,2] ~ obsTrue )
lines( statPred[,1] -1.96*statPred[,2] ~ obsTrue )
lines( statPred[,1] -1.96*sdObs ~ obsTrue, col="maroon" )
lines( statPred[,1] +1.96*sdObs ~ obsTrue, col="maroon" )
abline(0,1)
points( obs ~ obsTrue, col="maroon")
}
twDEMCSACont <- function(
### continuing simulated annealing DEMC based on previous result
mc ##<< result of twDEMCBlock
,nGen=512 ##<< overall number of generations to add (within one batch provide ctrlBatch$nGenBatch)
,nObs ##<< integer vector (nResComp) specifying the number of observations for each result component
,... ##<< further argument to \code{\link{twDEMCBlockInt}}
, m0 = calcM0twDEMC(getNParms(mc),getNChainsPop(mc)) ##<< minimum number of samples in step for extending runs
, controlTwDEMC = list() ##<< list argument to \code{\link{twDEMCBlockInt}} containing entry thin
, debugSequential=FALSE ##<< set to TRUE to avoid parallel execution, good for debugging
, restartFilename=NULL ##<< filename to write intermediate results to
#
, ctrlBatch = list( ##<< list of arguments controlling batch executions
##describe<<
nSamplesBatch=m0*8 ##<< number of samples for one call to twDEMCStep (multiplied by thin to get generations), defaults to 8*m0
#nGenBatch=m0*controlTwDEMC$thin*8 ##<< number of generations for one call to twDEMCStep
#nGenBatch=m0*controlTwDEMC$thin*3 ##<< number of generations for one call to twDEMCStep
##<< default: set in a way that on average each population (assuming half are significant) is appended by 2*m0 samples
#, nSampleMin=32 ##<< minimum number of samples in each population within batch so that calculation of average density is stable
, pThinPast=0.5 ##<< in each batch thin the past to given fraction before appending new results
##end<<
)
, ctrlT = list( ##<< list of arguments controlling Temperature decrease
##describe<<
TFix=numeric(0) ##<< numeric vector (nResComp) specifying a finite value for components with fixed Temperatue, and a non-finite for others
, TMax=numeric(0) ##<< numeric vector (nResComp) specifying a maximum temperature for result components.
, TDecProp=0.9 ##<< proportion of Temperature decrease: below one to diminish risk of decreasing Temperature too fast (below what is supported by other data streams)
##end<<
)
, ctrlConvergence = list( ##<< list or arguments controlling check for convergence
##describe<<
maxRelTChangeCrit=0.05 ##<< if Temperature of the components changes less than specified value, the algorithm can finish
, minTBase0 = 1e-4 ##<< if Temperature-1 gets below this temperature, the algorithm can finish
, maxLogDenDriftCrit=0.3 ##<< if difference between mean logDensity of first and fourth quartile of the sample is less than this value, we do not need further batches because of drift in logDensity
, gelmanCrit=1.4 ##<< do not change Temperature, if variance between chains is too high, i.e. Gelman Diag is above this value
, critSpecVarRatio=20 ##<< if proprotion of spectral Density to Variation is higher than this value, signal problems and resort to subspaces
, dumpfileBasename=NULL ##<< scalar string: filename to dump stack before stopping. May set to "recover"
, maxThin=64 ##<< scalar positive integer: maximum thinning interval. If not 0 then thinning is increased on too high spectral density
##end<<
)
){
# save calling arguments to allow an continuing an interrupted run
argsF <- as.list(sys.call())[-1] # do not store the file name and the first two arguments
argsF <- argsF[ !(names(argsF) %in% c("","nGen","mc")) ] # remove positional arguments and arguments mc and nGen
argsFEval <- lapply( argsF, eval.parent ) # remember values instead of language objects, which might not be there on a repeated call
mc$args <- argsFEval
finishEarlyMsg <- NULL
#
#-- fill in default argument values
if( is.null(controlTwDEMC$thin) ){
controlTwDEMC$thin <- 4
}
thin <- controlTwDEMC$thin
frm <- formals()
ctrlConvergence <- if( hasArg(ctrlConvergence) ) twMergeLists( eval(frm[["ctrlConvergence"]]), ctrlConvergence ) else ctrlConvergence
ctrlBatch <- if( hasArg(ctrlBatch) ) twMergeLists( eval(frm[["ctrlBatch"]]), ctrlBatch ) else ctrlBatch
ctrlT <- if( hasArg(ctrlT) ) twMergeLists( eval(frm[["ctrlT"]]), ctrlT ) else ctrlT
if( debugSequential ){
if( 0==length(ctrlConvergence$dumpfileBasename) ) ctrlConvergence$dumpfileBasename <- "recover"
#if( 0==length(ctrlSubspaces$argsFSplit$debugSequential) ) ctrlSubspaces$argsFSplit$debugSequential <- TRUE
}
#
nResComp <- ncol(mc$pops[[1]]$resLogDen)
#print("saCont: before completing temperature settings."); recover()
ctrlT$TFix <- completeResCompVec( ctrlT$TFix, colnames(mc$pops[[1]]$resLogDen) )
iFixTemp <- which( is.finite(ctrlT$TFix) )
iNonFixTemp <- which( !is.finite(ctrlT$TFix) )
ctrlT$TMax <- completeResCompVec( ctrlT$TMax, colnames(mc$pops[[1]]$resLogDen) )
iMaxTemp <- which( is.finite(ctrlT$TMax) )
#
res <- mc
.tmp.f <- function(){
mc0 <- concatPops(resEnd)
#mc0 <- concatPops(res)
matplot( mc0$temp, type="l" )
matplot( mc0$pAccept[,1,], type="l" )
matplot( mc0$pAccept[,2,], type="l" )
plot( as.mcmc.list(mc0), smooth=FALSE )
tmp <- calcTemperatedLogDen(res$pops[[1]]$resLogDen[,,1], getCurrentTemp(res) )
matplot( tmp, type="l" )
plot(rowSums(tmp))
#
matplot( mc0$resLogDen[,3,], type="l" )
matplot( mc0$parms[,"a",], type="l" )
matplot( mc0$parms[1:20,"b",], type="l" )
bo <- 1:10; iPop=1
plot( mc0$parms[bo,"a",iPop], mc0$parms[bo,"b",iPop], col=rainbow(100)[twMisc::twRescale(mc0$resLogDen[bo,"parmsSparse",iPop],c(10,100))] )
}
#
# sum nObs within density - assume nObs positions correspond to resLogDenComp positions
iDens <- seq_along(mc$dInfos)
#iDen=1
iCompsNonFixDen <- lapply( iDens, function(iDen){
irc <- mc$dInfos[[iDen]]$resCompPos
irc <- irc[ !(irc %in% iFixTemp) ]
})
nObsDen <- sapply( iDens, function(iDen){ sum( nObs[iCompsNonFixDen[[iDen]] ]) })
TCurr <- getCurrentTemp(mc)
#
iBatch=1
nGenDone <- 0 # tracking the generations already done
#nBatch <- ceiling( nGen/(ctrlBatch$nSamplesBatch*controlTwDEMC$thin) )
while( nGenDone < nGen){
if((0 < length(restartFilename)) && is.character(restartFilename) && restartFilename!=""){
resRestart.twDEMCSA = res #avoid variable confusion on load by giving a longer name
resRestart.twDEMCSA$nGenDone <- nGenDone # also store updated calculation of burnin time
resRestart.twDEMCSA$args <- argsFEval # also store updated calculation of burnin time
save(resRestart.twDEMCSA, file=restartFilename)
cat(paste("Saved resRestart.twDEMCSA to ",restartFilename,"\n",sep=""))
}
# do the diagnostics on the last halv of the chain, but continue
resEnd <- thin(res, start=ceiling(getNGen(res)) * ctrlBatch$pThinPast ) # neglect the first half part
resSparse <- squeeze(res, length.out=ceiling(getNSamples(res)*ctrlBatch$pThinPast) ) # better thin the entire chain so that distant path is sparsely kept, when continuing
#----- check for high autocorrelation wihin spaces
mcEnd <- concatPops(stackChainsPop(resEnd))
#plot( as.mcmc.list(mcEnd), smooth=FALSE )
specVarRatio <- apply( mcEnd$parms, 3, function(aSamplePurePop){
spec <- spectrum0.ar(aSamplePurePop)$spec
varSample <- apply(aSamplePurePop, 2, var)
spec/varSample
})
ssc <- stackChainsPop(resEnd) # combine all chains of one population
mcl <- as.mcmc.list(ssc)
#plot( as.mcmc.list(mcl), smooth=FALSE )
#plot( res$pops[[1]]$parms[,"b",1] ~ res$pops[[1]]$parms[,"a",1], col=rainbow(255)[round(twRescale(-res$pops[[1]]$resLogDen[,"logDen1",1],c(10,200)))] )
gelmanDiagRes <- try( {tmp<-gelman.diag(mcl); if(length(tmp$mpsrf)) tmp$mpsrf else tmp$psrf[1]} ) # cholesky decomposition may throw errors
#gelmanDiagRes <- try( gelman.diag(mcl)$mpsrf ) # cholesky decomposition may throw errors
# compare intra-chain to inter-chain gelman diag
#iPop <- 4
gelmanDiagPops <- sapply( 1:getNPops(resEnd), function(iPop){
mcl <- as.mcmc.list( subPops(resEnd,iPop) )
tmp <- try(gelman.diag(mcl), silent=TRUE)
if( inherits(tmp,"try-error")) 999 else if(length(tmp$mpsrf)) tmp$mpsrf else tmp$psrf[1]
})
maxGelmanDiagPops <- max(gelmanDiagPops)
T0 <- calcBaseTemp(TCurr, nObs, TFix=ctrlT$TFix, iNonFixTemp=iNonFixTemp)-1
# calculate average resLogDen of last part
resLogDen <- stackChains(concatPops(resEnd)$resLogDen)
.resDrift <- isLogDenDrift(resLogDen, resEnd$dInfos, maxDrift=ctrlConvergence$maxLogDenDriftCrit)
lDenLastPart <- structure( attr(.resDrift,"resTTest")[2,], names=names(resEnd$dInfos) )
##details<< \describe{\item{diagnostics}{
## in resulting twDEMC an entry diagnostics is returned. Its a list with diagnostics after each batch.
## Its components have been calculated on the last halv of the chains.
##describe<<
newDiags <- c(
T0=T0 ##<< Temperature (variance inflation factor) scaled for one observation minus one
,maxGelmanDiagPops=maxGelmanDiagPops ##<< maximum of within population gelman diagnostics
,gelmanDiag=gelmanDiagRes ##<< gelman diag calculated between populations
,maxSpecVarRatio=max(specVarRatio) ##<< ratio of spectral density to sample Variance (max over all parameters). This is a measure that increases with autocorrelation.
,lDenLastPart
#,attr(.resDrift,"logDenComp")
)
# if number of densities changed, reset columns in diagnostics (else rbind fails)
if( (length(resSparse$diagnostics) > 0) && (length(newDiags) != ncol(resSparse$diagnostics)) ){
resSparse$diagnostics <- cbind( resSparse$diagnostics[,1:4]
, matrix(NA, nrow=nrow(resSparse$diagnostics), ncol=length(lDenLastPart)
, dimnames=list(NULL, names(lDenLastPart))) )
}
resSparse$diagnostics <- rbind( resSparse$diagnostics, newDiags)
rownames( resSparse$diagnostics ) <- NULL
## Last components represent the average logDensity of the last 8th of the chains.
##details<< }}
#
##details<< \describe{\item{Adaptive thinning interval (\code{maxThin})}{
## If spectral variance is larger than its critical ratio \code{ctrlConvergence$critSpecVarRatio}
## then thinning interval needs to be increased.
## If doubling the thinning interval is not larger than \code{ctrlConvergence$maxThin} then a new batch
## with higher thinning is attempted. Else the SADEMC quits with reporting an error with return values component \code{failureMsg}.
##}}
if( any(specVarRatio > ctrlConvergence$critSpecVarRatio) ){
newThin <- resSparse$thin * 2
if( length(ctrlConvergence$maxThin) && (ctrlConvergence$maxThin > 0) && (newThin <= ctrlConvergence$maxThin) ){
# increase thinning interval instead of reporting error
print(paste("twDEMCSACont: too much autocorrelation (specVarRatio=",signif(max(specVarRatio),3),"). Increasing thinning interval to ",newThin,sep=""))
resSparse <- thin(resSparse, newThin )
controlTwDEMC$thin <- newThin
}else{
res <- resSparse
res$failureMsg <- paste("twDEMCSACont: too much autocorrelation (specVarRatio=",signif(max(specVarRatio),3),"). Try using divideTwDEMC or higher thinning interval",sep=")")
print(res$failureMsg)
break
}
}
# check if can decrease temperature, must calc TEnd and TEnd0
if( !inherits(gelmanDiagRes,"try-error") &&
maxGelmanDiagPops <= ctrlConvergence$gelmanCrit && # each pop converged
(gelmanDiagRes <= 1.2*maxGelmanDiagPops) # all pops explore the same optimum
){
if( length(ctrlT$TEndFixed)==0 ){
#
# print("twDEMCSACont: before Temperature calculation."); recover()
#debug(getBestModelIndices)
resLogDenT <- calcTemperatedLogDen(resLogDen, TCurr)
#iBest <- getBestModelIndices(resLogDenT, resEnd$dInfos, 0.01)
iBest <- getBestModelIndices(resLogDenT, resEnd$dInfos, ctrlT$qBest)
SkBest <- -2*resLogDen[iBest,,drop=FALSE]
# sapply( obs, function(obsk){ sum(obsk$sd^2) })
#trace(calcBaseTempSk, recover) #untrace(calcBaseTempSk)
TEnd0Target <- calcBaseTempSk(SkBest,nObs=nObs, TFix=ctrlT$TFix, iNonFixTemp=iNonFixTemp, isVerbose=ctrlT$isVerbose) -1
TEnd <- calcStreamTemp(TEnd0Target+1, nObs, TFix=ctrlT$TFix)
TEnd[iMaxTemp] <- pmin(TEnd[iMaxTemp], ctrlT$TMax[iMaxTemp]) # do not increase T above TMax
#relTChange <- abs(TEnd - TCurr)/TEnd
# slower TDecrease to avoid Temperatues that are not supported by other datastreams
TEnd0 <- T0 - ctrlT$TDecProp*(T0-TEnd0Target)
TEnd <- TCurr - ctrlT$TDecProp*(TCurr-TEnd)
} else { # ctrlT$TEndFixed is specified
# ultimately decrease temperature to 1
TEnd0Target <- ctrlT$TEndFixed-1
TEnd0 <- T0 - ctrlT$TDecProp*(T0-(ctrlT$TEndFixed-1))
TEnd <- TCurr <- calcStreamTemp(TEnd0+1, nObs, TFix=ctrlT$TFix)
}
#relTChange <- (TEnd0Target - T0)/(max(ctrlConvergence$minTBase0,TEnd0))
relTChange <- (TEnd0 - T0)/(max(ctrlConvergence$minTBase0,TEnd0))
cat("relTChange=",signif(relTChange,2),"\n")
#if( (max(relTChange) <= maxRelTChange) )
# recover()
#trace(isLogDenDrift, recover )
# may finish early, but do at least one batch
if( relTChange > 0 ){
ctrlT$TDecProp <- ctrlT$TDecProp * 2/3
cat("decreasing ctrlT$TDecProp to ",signif(ctrlT$TDecProp*100,2),"%\n")
}
if( (nGenDone != 0) &&
((max(abs(relTChange)) <= ctrlConvergence$maxRelTChangeCrit) || T0 <= ctrlConvergence$minTBase0 ) &&
!.resDrift
){
#res <- resSparse
# keep original res
finishEarlyMsg <- paste("twDEMCSA: Maximum Temperture change only ",signif(max(relTChange)*100,2)
,"% and no drift in logDensity (lden=",paste(signif(lDenLastPart,3),collapse=","),"). Finishing early.",sep="")
print( finishEarlyMsg)
print(paste("gelmanDiagPops=",signif(gelmanDiagRes,2)
," max(gelmanDiagPop)=",signif(maxGelmanDiagPops,2)
," max(specVarRatio)=",signif(max(specVarRatio),2)
," T0=",signif(T0,2)
," logDen=",paste(signif(lDenLastPart,3),collapse=",")
#," T=",paste(signif(TCurr,2),collapse=",")
, sep=""))
break
}
}else{ # cannot decrease temperature
TEnd0 <- T0
TEnd <-TCurr
}
#TMin <- pmin(TMin, TEnd)
print(paste("gelmanDiagPops=",signif(gelmanDiagRes,2)
," max(gelmanDiagPop)=",signif(maxGelmanDiagPops,2)
," max(specVarRatio)=",signif(max(specVarRatio),2)
," T0=",signif(T0,2)
," logDen=",paste(signif(lDenLastPart,3),collapse=",")
#," T=",paste(signif(TCurr,2),collapse=",")
, sep=""))
.nGenBatch <- min( ctrlBatch$nSamplesBatch*controlTwDEMC$thin, nGen - nGenDone )
print(paste("doing ",.nGenBatch," generations to TEnd0=",signif(TEnd0,2)
,", TEnd=",paste(signif(TEnd,2),collapse=","),sep="") )
res1 <- res <- twDEMCBlock( resSparse
, nGen=.nGenBatch
, debugSequential=debugSequential
, m0=m0
, controlTwDEMC=controlTwDEMC
, TEnd=TEnd
# replace lines below later on by ...
#, blocks = blocks
,...
)
TCurr <- getCurrentTemp(res)
nGenDone <- nGenDone + .nGenBatch
print(paste("finished ",nGenDone," out of ",nGen," gens. ", date(), sep="") )
} # end while( nGenDone < nGen)
res$TGlobal <- max(getCurrentTemp(res)[unlist(iCompsNonFixDen)])
res$args <- argsFEval
res$finishEarly <- finishEarlyMsg
#---------------------- update diagnostics for last batch (replicated from start of the loop)
resEnd <- thin(res, start=ceiling(getNGen(res)) * ctrlBatch$pThinPast ) # neglect the first half part
mcEnd <- concatPops(stackChainsPop(resEnd))
#----- check for high autocorrelation wihin spaces
#plot( as.mcmc.list(mcEnd), smooth=FALSE )
specVarRatio <- apply( mcEnd$parms, 3, function(aSamplePurePop){
spec <- spectrum0.ar(aSamplePurePop)$spec
varSample <- apply(aSamplePurePop, 2, var)
spec/varSample
})
ssc <- stackChainsPop(resEnd) # combine all chains of one population
mcl <- as.mcmc.list(ssc)
#plot( as.mcmc.list(mcl), smooth=FALSE )
#plot( res$pops[[1]]$parms[,"b",1] ~ res$pops[[1]]$parms[,"a",1], col=rainbow(255)[round(twRescale(-res$pops[[1]]$resLogDen[,"logDen1",1],c(10,200)))] )
gelmanDiagRes <- try( {tmp<-gelman.diag(mcl); if(length(tmp$mpsrf)) tmp$mpsrf else tmp$psrf[1]} ) # cholesky decomposition may throw errors
#gelmanDiagRes <- try( gelman.diag(mcl)$mpsrf ) # cholesky decomposition may throw errors
# compare intra-chain to inter-chain gelman diag
#iPop <- 4
gelmanDiagPops <- sapply( 1:getNPops(resEnd), function(iPop){
mcl <- as.mcmc.list( subPops(resEnd,iPop) )
tmp <- try(gelman.diag(mcl), silent=TRUE)
if( inherits(tmp,"try-error")) 999 else if(length(tmp$mpsrf)) tmp$mpsrf else tmp$psrf[1]
})
maxGelmanDiagPops <- max(gelmanDiagPops)
T0 <- calcBaseTemp(TCurr, nObs, TFix=ctrlT$TFix, iNonFixTemp=iNonFixTemp)-1
# calculate average resLogDen of last part
resLogDen <- stackChains(concatPops(resEnd)$resLogDen)
.resDrift <- isLogDenDrift(resLogDen, resEnd$dInfos, maxDrift=ctrlConvergence$maxLogDenDriftCrit)
lDenLastPart <- structure( attr(.resDrift,"resTTest")[2,], names=names(resEnd$dInfos) )
newDiags <- c(
T0=T0 ##<< Temperature (variance inflation factor) scaled for one observation minus one
,maxGelmanDiagPops=maxGelmanDiagPops ##<< maximum of within population gelman diagnostics
,gelmanDiag=gelmanDiagRes ##<< gelman diag calculated between populations
,maxSpecVarRatio=max(specVarRatio) ##<< ratio of spectral density to sample Variance (max over all parameters). This is a measure that increases with autocorrelation.
,lDenLastPart
#,attr(.resDrift,"logDenComp")
)
res$diagnostics <- rbind( res$diagnostics, newDiags)
rownames( res$diagnostics ) <- NULL
res
}
.tmp.f <- function(){ # ggplotResults
mc0 <- concatPops(res)
.nSample <- 128
dfDen <- rbind(
cbind( data.frame( scenario="S1", {tmp <- stackChains(mc0)[,-(1:getNBlocks(mc0))]; tmp[ seq(1,nrow(tmp),length.out=.nSample),] }))
)
dfDenM <- melt(dfDen)
#str(dfDenM)
require(ggplot2)
StatDensityNonZero <- StatDensity$proto( calculate<- function(.,data, scales, dMin=0.001, ...){
res <- StatDensity$calculate(data,scales,...)
# append a zero density at the edges
if( res$density[1] > 0) res <- rbind( data.frame(x=res$x[1]-diff(res$x[1:2]),density=0,scaled=0,count=0),res)
if( res$density[nrow(res)] > 0) res <- rbind( res, data.frame(x=res$x[nrow(res)]+diff(res$x[nrow(res)-c(1,0)]),density=0,scaled=0,count=0))
bo <- (res$density < dMin) & (c(res$density[-1],0) < dMin) & (c(0,res$density[-length(res$density)]) < dMin)
res$density[ bo ] <- NA
res$count[ bo ] <- NA
res$scaled[ bo ] <- NA
res
}, required_aes=c("x"))
#with(StatDensityNonZero, mtrace(calculate))
#with(StatDensityNonZero, mtrace(calculate,F))
stat_densityNonZero <- function (
### constructs a new StatPrior statistics based on aesthetics x and parName
mapping = NULL, data = NULL, geom = "line", position = "stack",
adjust = 1, ...
){
StatDensityNonZero$new(mapping = mapping, data = data, geom = geom,
position = position, adjust = adjust, ...)
}
#pgm <- geom_ribbon( alpha=0.8, aes(ymax = ..density.., ymin = -..density..), stat = "density")
#pgm <- geom_ribbon( aes(ymax = ..density.., ymin = -..density..), stat = "density")
#pgm <- stat_densityNonZero( aes(ymax = -..density..), size=1, geom="line")
pgm <- geom_line(aes(y = +..scaled..), stat="densityNonZero", size=1)
pgm2 <- geom_line(aes(y = -..scaled..), stat="densityNonZero", size=1)
optsm <- opts(axis.title.x = theme_blank(), axis.text.y = theme_blank(), axis.ticks.y = theme_blank() )
pa <- ggplot(dfDen, aes(x = a, colour=scenario, linetype=scenario)) + pgm + pgm2 + optsm + scale_y_continuous('a')
pb <- ggplot(dfDen, aes(x = b, colour=scenario, linetype=scenario)) + pgm + pgm2 + optsm + scale_y_continuous('b')
#pb
windows(width=7, height=3)
grid.newpage()
pushViewport( viewport(layout=grid.layout(2,1)))
print(pa , vp = viewport(layout.pos.row=1,layout.pos.col=1))
print(pb + opts(legend.position = "none") , vp = viewport(layout.pos.row=2,layout.pos.col=1))
#----------- ggplot predictive posterior
#scenarios <- c("R","RS","RSw","DG","DM")
scenarios <- c("R")
nScen <- length(scenarios)
# infer quantiles of predictions
.nSample=128
y1M <- array( NA_real_, dim=c(.nSample, length(twTwoDenEx1$obs$y1),nScen), dimnames=list(sample=NULL,iObs=NULL,scenario=scenarios) )
y2M <- array( NA_real_, dim=c(.nSample, length(twTwoDenEx1$obs$y2),nScen), dimnames=list(sample=NULL,iObs=NULL,scenario=scenarios) )
resScen <- #list( res1, res2, res2b, res3a, res3 ); names(resScen) <- scenarios
resScen <- list( mc0 ); names(resScen) <- scenarios
#scen <- "R"
for( scen in scenarios ){
ss <- stackChains(resScen[[scen]])[,-(1:getNBlocks(resScen[[scen]]))]
ssThin <- ss[round(seq(1,nrow(ss),length.out=.nSample)),]
#i <- .nSample
for( i in 1:.nSample){
pred <- twTwoDenEx1$fModel(ssThin[i,], xSparse=twTwoDenEx1$xSparse, xRich=twTwoDenEx1$xRich, thresholdCovar=thresholdCovar)
y1M[i,,scen] <- pred$y1
y2M[i,,scen] <- pred$y2
}
}
predQuantilesY1 <- lapply( scenarios, function(scen){
t(apply( y1M[,,scen],2, quantile, probs=c(0.025,0.5,0.975) ))
}); names(predQuantilesY1) <- scenarios
predQuantilesY2 <- lapply( scenarios, function(scen){
t(apply( y2M[,,scen],2, quantile, probs=c(0.025,0.5,0.975) ))
}); names(predQuantilesY2) <- scenarios
str(predQuantilesY1)
rm( y1M, y2M) # free space, we only need the summary
#str(pred1)
dfPred <- rbind(
cbind( data.frame( scenario="R", variable="y1", value = pred1$y1 ), predQuantilesY1[["R"]] )
,cbind( data.frame( scenario="R", variable="y2", value = pred1$y2 ), predQuantilesY2[["R"]] )
)
dfPred$observations <- NA_real_
dfPred$observations[dfPred$variable=="y1"] <- twTwoDenEx1$obs$y1
dfPred$observations[dfPred$variable=="y2"] <- twTwoDenEx1$obs$y2
colnames(dfPred)[ match(c("2.5%","50%","97.5%"),colnames(dfPred)) ] <- c("lower","median","upper")
str(dfPred)
p1 <- ggplot(dfPred, aes(x=value, y=observations, colour=scenario) ) +
geom_errorbarh(aes(xmax = upper, xmin = lower)) +
geom_point() +
facet_wrap( ~ variable, scales="free") +
opts(axis.title.x=theme_blank() ) +
geom_abline(colour="black") +
c()
p1
}
.tmp.f <- function(){ # does not work # AllDensities
# same as example but with each parameter updated against both densities
data(twTwoDenEx1)
thetaPrior <- twTwoDenEx1$thetaTrue
covarTheta <- diag((thetaPrior*0.3)^2)
invCovarTheta <- (thetaPrior*0.3)^2 # given as independent variances for faster calculation
thresholdCovar = 0.3 # the true value
thresholdCovar = 0 # the effective model that glosses over this threshold
# for updating each in turn against both densities
dInfos=list(
dSparse=list(fLogDen=denSparsePrior, argsFLogDen=list(thresholdCovar=thresholdCovar, twTwoDenEx=twTwoDenEx1, theta0=thetaPrior, thetaPrior=thetaPrior, invCovarTheta=invCovarTheta))
,dRich=list(fLogDen=denRich, argsFLogDen=list(thresholdCovar=thresholdCovar,twTwoDenEx=twTwoDenEx1, theta0=thetaPrior))
)
blocks = list(
bSparse=list(dInfoPos="dSparse", compPos=c("a","b") )
,bRich=list(dInfoPos="dRich", compPos=c("a","b") )
)
do.call( dInfos$dSparse$fLogDen, c(list(theta=twTwoDenEx1$theta0),dInfos$dSparse$argsFLogDen))
do.call( dInfos$dRich$fLogDen, c(list(theta=twTwoDenEx1$theta0),dInfos$dRich$argsFLogDen))
#str(twTwoDenEx1)
nObs <- c( obsSparse=length(twTwoDenEx1$obs$y1), logDen1=length(twTwoDenEx1$obs$y2) )
#trace(twDEMCSA, recover )
resPops <- res <- twDEMCSA( thetaPrior, covarTheta, dInfos=dInfos, blocks=blocks, nObs=nObs, debugSequential=TRUE, ctrlBatch=list(nSamplesBatch=64) )
# continue run for 512 generations
res <- twDEMCSA( res, 512 )
(TCurr <- getCurrentTemp(resPops))
mc0 <- concatPops(res)
logDenT <- calcTemperatedLogDen(stackChains(mc0$resLogDen), TCurr)
iBest <- getBestModelIndex( logDenT, res$dInfos )
maxLogDenT <- logDenT[iBest, ]
ss <- stackChains(mc0$parms)
(thetaBest <- ss[iBest, ])
(.qq <- apply(ss,2,quantile, probs=c(0.025,0.5,0.975) ))
plot( ss[,"a"], ss[,"b"], col=rainbow(100)[twRescale(rowSums(logDenT),c(1,100))] )
plot( ss[,"a"], ss[,"b"], col=rainbow(100)[twRescale(logDenT[,"obsSparse"],c(1,100))] )
plot( ss[,"a"], ss[,"b"], col=rainbow(100)[twRescale(logDenT[,"logDen1"],c(1,100))] )
plot( ss[,"a"], ss[,"b"], col=rgb(
twRescale(logDenT[,"obsSparse"]),0, twRescale(logDenT[,"logDen1"]) ))
apply( apply( logDenT, 2, quantile, probs=c(0.1,0.9) ),2, diff )
# density of parameters
plot( density(ss[,"a"])); abline(v=thetaPrior["a"]); abline(v=thetaBest["a"], col="blue")
plot( density(ss[,"b"])); abline(v=thetaPrior["b"]); abline(v=thetaBest["b"], col="blue")
# predictive posterior (best model only)
pred <- pred1 <- with( twTwoDenEx1, fModel(thetaBest, xSparse=xSparse, xRich=xRich) )
plot( pred$y1, twTwoDenEx1$obs$y1 ); abline(0,1)
plot( pred$y2, twTwoDenEx1$obs$y2 ); abline(0,1)
}
.tmp.f <- function(){ # # AllSparse
# same as example but with b also updated against Sparse, a only against Sparse
data(twTwoDenEx1)
thetaPrior <- twTwoDenEx1$thetaTrue
covarTheta <- diag((thetaPrior*0.3)^2)
invCovarTheta <- (thetaPrior*0.3)^2 # given as independent variances for faster calculation
thresholdCovar = 0.3 # the true value
thresholdCovar = 0 # the effective model that glosses over this threshold
# for updating each in turn against both densities
dInfos=list(
dSparse=list(fLogDen=denSparsePrior, argsFLogDen=list(thresholdCovar=thresholdCovar, twTwoDenEx=twTwoDenEx1, theta0=thetaPrior, thetaPrior=thetaPrior, invCovarTheta=invCovarTheta))
,dRich=list(fLogDen=denRich, argsFLogDen=list(thresholdCovar=thresholdCovar,twTwoDenEx=twTwoDenEx1, theta0=thetaPrior))
)
blocks = list(
bSparse=list(dInfoPos="dSparse", compPos=c("a","b") )
,bRich=list(dInfoPos="dRich", compPos=c("a") )
)
do.call( dInfos$dSparse$fLogDen, c(list(theta=twTwoDenEx1$theta0),dInfos$dSparse$argsFLogDen))
do.call( dInfos$dRich$fLogDen, c(list(theta=twTwoDenEx1$theta0),dInfos$dRich$argsFLogDen))
#str(twTwoDenEx1)
nObs <- c( obsSparse=length(twTwoDenEx1$obs$y1), logDen1=length(twTwoDenEx1$obs$y2) )
#untrace(twDEMCSA )
#trace(twDEMCSA, recover )
resPops <- res <- twDEMCSA( thetaPrior, covarTheta, dInfos=dInfos, blocks=blocks, nObs=nObs, nBatch=5, debugSequential=TRUE, TFix=c(1,NA,NA) )
resPops <- res <- twDEMCSACont( res, 256, nObs=nObs, debugSequential=TRUE
, ctrlT=list( TFix=c(parmsSparse=1) )
)
(TCurr <- getCurrentTemp(resPops))
mc0 <- concatPops(res)
logDenT <- calcTemperatedLogDen(stackChains(mc0$resLogDen), TCurr)
iBest <- getBestModelIndex( logDenT, res$dInfos )
maxLogDenT <- logDenT[iBest, ]
ss <- stackChains(mc0$parms)
(thetaBest <- ss[iBest, ])
(.qq <- apply(ss,2,quantile, probs=c(0.025,0.5,0.975) ))
plot( ss[,"a"], ss[,"b"], col=rainbow(100)[twRescale(rowSums(logDenT),c(1,100))] )
plot( ss[,"a"], ss[,"b"], col=rainbow(100)[twRescale(logDenT[,"obsSparse"],c(1,100))] )
plot( ss[,"a"], ss[,"b"], col=rainbow(100)[twRescale(logDenT[,"logDen1"],c(1,100))] )
plot( ss[,"a"], ss[,"b"], col=rgb(
twRescale(logDenT[,"obsSparse"]),0, twRescale(logDenT[,"logDen1"]) ))
apply( apply( logDenT, 2, quantile, probs=c(0.1,0.9) ),2, diff )
# density of parameters
plot( density(ss[,"a"])); abline(v=thetaPrior["a"]); abline(v=thetaBest["a"], col="blue")
plot( density(ss[,"b"])); abline(v=thetaPrior["b"]); abline(v=thetaBest["b"], col="blue")
# predictive posterior (best model only)
pred <- pred1 <- with( twTwoDenEx1, fModel(thetaBest, xSparse=xSparse, xRich=xRich) )
plot( pred$y1, twTwoDenEx1$obs$y1 ); abline(0,1)
plot( pred$y2, twTwoDenEx1$obs$y2 ); abline(0,1)
}
.tmp.f <- function(){ # does not work # .tmp.2DenSwitched
# same as example but with each parameter b updated against Sparse
data(twTwoDenEx1)
thetaPrior <- twTwoDenEx1$thetaTrue
covarTheta <- diag((thetaPrior*0.3)^2)
invCovarTheta <- (thetaPrior*0.3)^2 # given as independent variances for faster calculation
thresholdCovar = 0.3 # the true value
thresholdCovar = 0 # the effective model that glosses over this threshold
# for updating each in turn against both densities
dInfos=list(
dSparse=list(fLogDen=denSparsePrior, argsFLogDen=list(thresholdCovar=thresholdCovar, twTwoDenEx=twTwoDenEx1, theta0=thetaPrior, thetaPrior=thetaPrior, invCovarTheta=invCovarTheta))
,dRich=list(fLogDen=denRich, argsFLogDen=list(thresholdCovar=thresholdCovar,twTwoDenEx=twTwoDenEx1, theta0=thetaPrior))
)
blocks = list(
bSparse=list(dInfoPos="dSparse", compPos=c("b") )
,bRich=list(dInfoPos="dRich", compPos=c("a") )
)
do.call( dInfos$dSparse$fLogDen, c(list(theta=twTwoDenEx1$theta0),dInfos$dSparse$argsFLogDen))
do.call( dInfos$dRich$fLogDen, c(list(theta=twTwoDenEx1$theta0),dInfos$dRich$argsFLogDen))
#str(twTwoDenEx1)
nObs <- c( obsSparse=length(twTwoDenEx1$obs$y1), logDen1=length(twTwoDenEx1$obs$y2) )
#untrace(twDEMCSA )
#trace(twDEMCSA, recover )
resPops <- res <- twDEMCSA( thetaPrior, covarTheta, dInfos=dInfos, blocks=blocks, nObs=nObs, nBatch=5, debugSequential=TRUE )
(TCurr <- getCurrentTemp(resPops))
mc0 <- concatPops(res)
logDenT <- calcTemperatedLogDen(stackChains(mc0$resLogDen), TCurr)
iBest <- getBestModelIndex( logDenT, res$dInfos )
maxLogDenT <- logDenT[iBest, ]
ss <- stackChains(mc0$parms)
(thetaBest <- ss[iBest, ])
(.qq <- apply(ss,2,quantile, probs=c(0.025,0.5,0.975) ))
plot( ss[,"a"], ss[,"b"], col=rainbow(100)[twRescale(rowSums(logDenT),c(1,100))] )
plot( ss[,"a"], ss[,"b"], col=rainbow(100)[twRescale(logDenT[,"obsSparse"],c(1,100))] )
plot( ss[,"a"], ss[,"b"], col=rainbow(100)[twRescale(logDenT[,"logDen1"],c(1,100))] )
plot( ss[,"a"], ss[,"b"], col=rgb(
twRescale(logDenT[,"obsSparse"]),0, twRescale(logDenT[,"logDen1"]) ))
apply( apply( logDenT, 2, quantile, probs=c(0.1,0.9) ),2, diff )
# density of parameters
plot( density(ss[,"a"])); abline(v=thetaPrior["a"]); abline(v=thetaBest["a"], col="blue")
plot( density(ss[,"b"])); abline(v=thetaPrior["b"]); abline(v=thetaBest["b"], col="blue")
# predictive posterior (best model only)
pred <- pred1 <- with( twTwoDenEx1, fModel(thetaBest, xSparse=xSparse, xRich=xRich) )
plot( pred$y1, twTwoDenEx1$obs$y1 ); abline(0,1)
plot( pred$y2, twTwoDenEx1$obs$y2 ); abline(0,1)
}
.tmp.f <- function(){ # oneDensity_andCluster
# same as example but with only one combined density for both parameters
data(twTwoDenEx1)
thetaPrior <- twTwoDenEx1$thetaTrue
covarTheta <- diag((thetaPrior*0.3)^2)
invCovarTheta <- (thetaPrior*0.3)^2 # given as independent variances for faster calculation
thresholdCovar = 0.3 # the true value
thresholdCovar = 0 # the effective model that glosses over this threshold
# for only one logDensity - Temperature of the strongest component goes to zero and other increase according to mismatch
dInfos=list( list(fLogDen=denBoth, argsFLogDen=list(thresholdCovar=thresholdCovar, twTwoDenEx=twTwoDenEx1, theta0=thetaPrior, thetaPrior=thetaPrior, invCovarTheta=invCovarTheta)) )
do.call( dInfos[[1]]$fLogDen, c(list(theta=twTwoDenEx1$theta0),dInfos[[1]]$argsFLogDen))
#str(twTwoDenEx1)
nObs <- c( y1=length(twTwoDenEx1$obs$y1), y2=length(twTwoDenEx1$obs$y2) )
#trace(twDEMCSA, recover)
argsTwDEMCSA <- list( thetaPrior=thetaPrior, covarTheta=covarTheta, dInfos=dInfos, nObs=nObs
,TFix = c(1,NA,NA)
,nGen=256
#,TMax = c(NA,1.2,NA) # do not increase Sparse observations again too much
, nBatch=3
, debugSequential=TRUE
)
resPops <- res <- do.call( twDEMCSA, argsTwDEMCSA )
.tmp.f <- function(){ #continueRun
res$pops[[2]]$spaceInd <- 4 # test spaces different from 1:nSpace
res2 <- twDEMCSA(res)
}
.tmp.f <- function(){ #byCluster
runClusterParms <- list(
fSetupCluster = function(){library(twDEMC)}
,fRun = twDEMCSA
,argsFRun = within( argsTwDEMCSA, debugSequential<-FALSE )
)
save(runClusterParms, file=file.path("..","..","projects","asom","parms","saOneDensity.RData"))
# from asom directory run:
# R CMD BATCH --vanilla '--args paramFile="parms/saOneDensity.RData"' runCluster.R Rout/runCluster_0.rout
# or from libra home directory run
# ./bsubr_i.sh runCluster.R iproc=0 nprocSinge=1 'paramFile="parms/saOneDensity.RData"'
# or
# bsub -q SLES -n 4 ./bsubr_i.sh runCluster.R nprocSingle=4 'paramFile="parms/saOneDensity.RData"'
load("parms/res_saOneDensity_1.RData")
}
(TCurr <- getCurrentTemp(res))
mc0 <- concatPops(res)
logDenT <- calcTemperatedLogDen(stackChains(mc0$resLogDen), TCurr)
iBest <- getBestModelIndex( logDenT, res$dInfos )
maxLogDenT <- logDenT[iBest, ]
ss <- stackChains(mc0$parms)
(thetaBest <- ss[iBest, ])
(.qq <- apply(ss,2,quantile, probs=c(0.025,0.5,0.975) ))
plot( ss[,"a"], ss[,"b"], col=rainbow(100)[twRescale(rowSums(logDenT),c(1,100))] )
plot( ss[,"a"], ss[,"b"], col=rainbow(100)[twRescale(logDenT[,"y1"],c(10,100))] )
plot( ss[,"a"], ss[,"b"], col=rainbow(100)[twRescale(logDenT[,"y2"],c(10,100))] )
plot( ss[,"a"], ss[,"b"], col=rgb(
twRescale(logDenT[,"y1"]),0, twRescale(logDenT[,"y2"]) ))
apply( apply( logDenT, 2, quantile, probs=c(0.2,0.8) ),2, diff )
# density of parameters
plot( density(ss[,"a"])); abline(v=thetaPrior["a"]); abline(v=thetaBest["a"], col="blue")
plot( density(ss[,"b"])); abline(v=thetaPrior["b"]); abline(v=thetaBest["b"], col="blue")
# predictive posterior (best model only)
pred <- pred1 <- with( twTwoDenEx1, fModel(thetaBest, xSparse=xSparse, xRich=xRich) )
plot( pred$y1, twTwoDenEx1$obs$y1 ); abline(0,1)
plot( pred$y2, twTwoDenEx1$obs$y2 ); abline(0,1)
}
isLogDenDrift <- function(
### check whether first quartile all the logDensities is significantly smaller than last quartile
logDenT ##<< numeric array (nStep x nResComp): logDensity (highest are best)
, dInfos ##<< list of lists with entry resCompPos (integer vector) specifying the position of result components for each density
, alpha=0.05 ##<< the significance level for a difference
, maxDrift=0.3 ##<< difference in LogDensity, below which no drift is signalled
){
##details<<
## Because of large sample sizes, very small differences may be significantly different.
## Use argument minDiff to specify below which difference a significant difference is not regarded as drift.
nr <- nrow(logDenT)
nr4 <- nr %/% 4
ds1 <- logDenT[1:nr4, ,drop=FALSE] # first quater of the dataset
ds4 <- logDenT[(nr+1-nr4):nr, ,drop=FALSE] # last quater of the dataset
#iDen <- 1
tres <- sapply( seq_along(dInfos), function(iDen){
dInfo <- dInfos[[iDen]]
rs1 <- rowSums(ds1[,dInfo$resCompPos ,drop=FALSE])
rs4 <- rowSums(ds4[,dInfo$resCompPos ,drop=FALSE])
resTTest <- try(t.test(rs1,rs4,"less"), silent=TRUE)
if( inherits(resTTest,"try-error") ){ # logDen essentially constant
#bo <- (diff(resTTest$estimate) >= maxDrift ) && (resTTest$p.value <= alpha)
c(rs1,rs4,1)
} else {
c( resTTest$estimate, p=resTTest$p.value)
}
})
#tresi <- tres[,1]
boL <- apply( tres, 2, function(tresi){
(diff(tresi[1:2]) >= maxDrift ) && (tresi["p"] <= alpha)
})
ret <- any( boL )
attr( ret, "resTTest") <- tres
attr( ret, "logDenComp") <- colMeans(ds4) # average logDensity of all components of last quarter
##value<< TRUE if any of the logDensities are a significantly greater in the fourth quantile compared to the first quantile of the samples
## , attribute \code{resTTest}: numeric maxtrix (logDenStart, logDenEnd, pTTest x nDInfo )
## , attribute \code{logDenComp}: numeric vector (nResComp): average logDenT of all components of last quater
ret
}
attr(isLogDenDrift,"ex") <- function(){
data(twdemcEx1)
logDenT <- calcTemperatedLogDen( twdemcEx1, getCurrentTemp(twdemcEx1))
#undebug(isLogDenDrift)
isLogDenDrift(logDenT, twdemcEx1$dInfos )
}
twRunDEMCSA <- function(
### wrapper to twDEMCSA compliant to runCluster.R
argsTwDEMCSA
### Arguments passed to twDEMCSA -> twDEMCBlockInt
### It is updated by \dots.
,... ##<< further arguments passed to twDEMCSA -> twDEMCBlockInt
,prevResRunCluster=NULL ##<< results of call to twRunDEMCSA, argument required to be called from runCluster.R
,restartFilename=NULL ##<< name of the file to store restart information, argument required to be called from runCluster.R
){
require(twDEMC)
#update argsDEMC to args given
argsDEMC <- argsTwDEMCSA
# update the restartFilename argument
.dots <- list(...)
argsDEMC[ names(.dots) ] <- .dots
#for( argName in names(.dots) ) argsDEMC[argName] <- .dots[argName]
if( 0 != length(restartFilename)) argsDEMC$restartFilename <- restartFilename
# do the actual call
res <- do.call( twDEMCSA, argsDEMC )
}
.tmp.f <- function(){ # testRestart
rFileName <- file.path("tmp","example_twDEMCSA.RData")
load(rFileName) # resRestart.twDEMCSA
#untrace(twDEMCSACont)
#trace(twDEMCSACont, recover )
res1 <- do.call( twDEMCSACont, c( list(mc=resRestart.twDEMCSA), resRestart.twDEMCSA$args ))
}
completeResCompVec <- function(
### given a vector with names, make a vector corresponding to resComp with components, not yet given filled by NA
x ##<< named vector to be completed
,resCompNames ##<< names of the result components
){
xAll <- x
nResComp <- length(resCompNames)
if( 0 == length( x) )
xAll <- structure( rep( NA_real_, nResComp), names=resCompNames )
#else if( nResComp != length( x) ){ # twutz 141120: can give a different names of same length, must correct
else{
iFix <- sapply( names(x), match, resCompNames )
if( any(is.na(iFix))){
misNames <- paste(names(x)[ is.na(iFix)], collapse=",")
stop(paste("completeResCompVec: x must be of the same length as resCompNames or all its names must correspond names of resCompNames. Unmatched names:",misNames))
}
xAll <- structure( rep( NA_real_, nResComp), names=resCompNames )
xfin <- x[ is.finite(x) ]
#iName <- names(xfin)[1]
for( iName in names(xfin) ){
xAll[ names(xAll)== iName ] <- xfin[iName]
}
}
xAll
### vector with names resCompNames, with corresponding values of x set, others NA
}
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.