inst/doc/SamplingObsUncertainty.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 reproducible results

## ------------------------------------------------------------------------
denSigmaEx1Obs <- function(
    theta                ##<< named numeric vector a,b,logSigma2
    ,twLinreg=twLinreg1  ##<< list with components fModel, xval, and obs 
){
    pred <- twLinreg$fModel(theta[1:2],twLinreg$xval)
    resid <- pred - twLinreg$obs
    varObs <- exp(theta["logSigma2"])
    structure(-1/2 * sum(resid^2)/varObs, names="obs")          
}

## ----eval=FALSE----------------------------------------------------------
#      dInfos=list(
#          dObs=list(fLogDen=denSigmaEx1Obs, argsFLogDen=list(twLinreg=twLinreg1), compPosDen=c("a","b","logSigma2"))
#      )

## ----eval=FALSE----------------------------------------------------------
#      bObs=list(compPos=c("a","b"), dInfoPos="dObs" )

## ------------------------------------------------------------------------
fResidDummyModel <- function(theta, twLinreg=twLinreg1){
    pred <- twLinreg$fModel(theta,twLinreg$xval)
    resid <- pred - twLinreg$obs
}

## ----eval=FALSE, tidy=FALSE----------------------------------------------
#      bSigma=list( compPos=c("logSigma2"), fUpdateBlock=updateSigmaByGibbsSamplingInvchisq,
#                argsFUpdate=list(fResid=fResidDummyModel, twLinreg=twLinreg1))

## ----sa1, cache=TRUE, results='hide', tidy=FALSE-------------------------
logSigma2 <- log( mean(twLinreg1$sdObs^2) )     # expected sigma2
varLogSigma2=0.8*logSigma2                      # variance for initialization

resBlock <- twDEMCSA( 
    theta=c(twLinreg1$theta0,logSigma2=logSigma2), 
    covarTheta=diag(c(twLinreg1$sdTheta^2,varLogSigma2)),
    dInfos=list(
        dObs=list(fLogDen=denSigmaEx1Obs, 
                argsFLogDen=list(twLinreg=twLinreg1), 
                compPosDen=c("a","b","logSigma2"))
    ),
    blocks=list(
        bObs=list(compPos=c("a","b"), dInfoPos="dObs" ),
        bSigma=list( compPos=c("logSigma2"), fUpdateBlock=updateSigmaByGibbsSamplingInvchisq,
                argsFUpdate=list(fResid=fResidDummyModel, twLinreg=twLinreg1))
    ),
    nObs=c(obs=length(twLinreg1$obs)),
    ctrlT=list(TBaseInit=1, TEndFixed=1)
)

## ----,spar=TRUE----------------------------------------------------------
plot( as.mcmc.list(resBlock), smooth=FALSE )

## ------------------------------------------------------------------------
denSigmaEx1ObsIntermediate <- function(
        ### unnormalized log density for observations for given parameters
        theta                ##<< named numeric vector a,b,logSigma2
        ,intermediate = list()  
        ,twLinreg=twLinreg1  ##<< list with components xval and obs 
){
    # the predictions by the forward model are an intermediate retuls that can be reused
    if( !length(intermediate) ){
        intermediate <- twLinreg$fModel(theta,twLinreg$xval)
    } 
    pred <- intermediate
    resid <- pred - twLinreg$obs
    varObs <- exp(theta[3])
    structure(-1/2 * sum(resid^2)/varObs, names="obs", intermediate=intermediate) 
}

fResidDummyModelIntermediate <- function(theta, intermediate=list(), twLinreg){
    if( !length(intermediate) ){
        intermediate <- twLinreg$fModel(theta,twLinreg$xval)
    } 
    pred <- intermediate
    resid <- structure(pred - twLinreg$obs, intermediate=intermediate)
}

## ----SAIntermediate, results='hide', tidy=FALSE--------------------------
resBlock <- twDEMCSA( 
    theta=c(twLinreg1$theta0,logSigma2=logSigma2), 
    covarTheta=diag(c(twLinreg1$sdTheta^2,varLogSigma2)),
    dInfos=list(
        dObs=list(fLogDen=denSigmaEx1ObsIntermediate,
            argsFLogDen=list(twLinreg=twLinreg1),
            compPosDen=c("a","b","logSigma2")
        )
    ),
    blocks=list(
        bObs=list(compPos=c("a","b"), dInfoPos="dObs", 
            intermediateId="fModel"
        ),
        bSigma=list(compPos=c("logSigma2"), fUpdateBlock=updateSigmaByGibbsSamplingInvchisq,
            argsFUpdate=list(fResid=fResidDummyModelIntermediate, twLinreg=twLinreg1), 
            intermediateId="fModel"
        )
    ),
    nObs=c(obs=length(twLinreg1$obs)),
    ctrlT=list(TBaseInit=1, TEndFixed=1),
)

## ----,spar=TRUE----------------------------------------------------------
plot( as.mcmc.list(resBlock), smooth=FALSE )