Description Usage Arguments Details Value Author(s) Examples
simulated annealing DEMC
1 2 3 4 5 6 7 | twDEMCSA(thetaPrior, covarTheta, nGen = 512, nObs, ..., dInfos,
m0 = calcM0twDEMC(length(thetaPrior), nChainPop), controlTwDEMC = list(),
debugSequential = FALSE, restartFilename = NULL, remoteDumpfileBasename = NULL,
nChainPop = 4, nPop = 2, doIncludePrior = FALSE, ctrlBatch = list(nGen0 = m0 *
controlTwDEMC$thin * 3), ctrlT = list(qTempInit = 0.4,
TBaseInit = NULL, isVerbose = FALSE, TEndFixed = NULL,
qBest = 0.1), ctrlConvergence = list())
|
thetaPrior |
vector of parameters, point estimate
<< , alternatively array with initial states, as returned by |
covarTheta |
the a prior covariance of parameters, see |
nGen |
number of generations |
nObs |
integer vector (names resComp) specifying the number of observations for each result component for each resCompName there must be an entry in nObs |
... |
further argument to |
dInfos |
argument to |
m0 |
minimum number of samples in step for extending runs minimum number of samples in step for extending runs |
controlTwDEMC |
list argument to |
debugSequential |
set to TRUE to avoid parallel execution, good for debugging |
restartFilename |
filename to write intermediate results to |
remoteDumpfileBasename |
fileBasename to write dumps to on error |
nChainPop |
number of chains within population |
nPop |
number of populations |
doIncludePrior |
should the prior be part of initial population << Recommendation to set to false, because if the TRUE parameter is in initial set, the Temperature is set to 1 |
ctrlBatch |
list of arguments controlling batch executions, see
|
ctrlT |
list of arguments controlling Temperature decrease, see
|
ctrlConvergence |
list or arguments controlling check for convergence, see |
Initial parameters are ranked according to their maximum log-Density across components.
The parameters and logDensity results at rank position defined by argument 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 calcBaseTemp
) and rescaled to stream temperatures by calcStreamTemp
.
An object of class twDEMCPops
as described in twDEMCBlockInt
.
See subset.twDEMCPops
for processing further handling of this class.
Thomas Wutzler
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 |
#--------------- 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
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.