inst/doc/SimpleSingleModel.R

In twDEMC: parallel DEMC

## ----setup, include=FALSE------------------------------------------------
library(knitr)
opts_chunk$set(out.extra='style="display:block; margin: auto"'
    #, fig.align="center"
    , fig.width=4.3, fig.height=3.2, dev.args=list(pointsize=10)
    )
knit_hooks$set(spar = function(before, options, envir) {
    if (before){
        par( las=1 )                   #also y axis labels horizontal
        par(mar=c(2.0,3.3,0,0)+0.3 )  #margins
        par(tck=0.02 )                          #axe-tick length inside plots             
        par(mgp=c(1.1,0.2,0) )  #positioning of axis title, axis labels, axis
     }
})

## ----results='hide'------------------------------------------------------
library(twDEMC)
data(twLinreg1) 
set.seed(0815)      # for reproducable results

## ------------------------------------------------------------------------
dummyTwDEMCModel 

## ------------------------------------------------------------------------
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                       ### further arguments to the model, here the vector of predictors
))
do.call( logDenGaussian, c( list(theta=c(a=10.8,b=5.2)), argsFLogDen ))

## ----eval=FALSE----------------------------------------------------------
#  dInfos=list(den1=list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen))

## ----twDEMCSA, cache=TRUE------------------------------------------------
mcPops <-  twDEMCSA( 
        dInfos=list(den1=list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen)),
        theta=twLinreg1$theta0, covarTheta=diag(twLinreg1$sdTheta^2),       # for generating an initial population
        ctrlT=list(TFix=c(parms=1)),                       # do not use increased temperature for priors
        nObs=c(obs=length(argsFLogDen$obs))               # number of records in observation data stream(s)
)

## ----plotCoda,spar=TRUE--------------------------------------------------
rescoda <- as.mcmc.list(mcPops)
plot(rescoda, smooth=FALSE)

## ----mcPopsConv, cache=TRUE, results='hide'------------------------------
mcPopsConv <- twDEMCBlock( mcPops, nGen=ceiling(256 * mcPops$thin / getNChains(mcPops)), extendRun=FALSE )

## ------------------------------------------------------------------------
ss <- stackChains(mcPopsConv)
summary(ss)

## ----mcPopsT1, cache=TRUE, results='hide'--------------------------------
 # decrease Temp to exponentially to 1
mcPopsT1 <- twDEMCBlock( mcPops, TEnd = 1, nGen=ceiling(256 * mcPops$thin / getNChains(mcPops)) )         
 # sample at this temperature
mcPopsT1Conv <- twDEMCBlock( mcPopsT1, nGen=ceiling(256 * mcPops$thin / getNChains(mcPops)), extendRun=FALSE )      
summary( stackChains(mcPopsT1Conv) )   # difference to mcPopsConv negligible