inst/doc/TwoDen.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(twTwoDenEx1)
set.seed(0815)      # for reproducable results

## ----eval=FALSE----------------------------------------------------------
#  twTwoDenEx1$fModel = function(
#          theta           ##<< model parameters a and b
#          , xSparse       ##<< numeric vector of Sparse input/output relationship
#          , xRich         ##<< numeric vector of rich input/output relationship
#          , thresholdCovar=0      ##<< model structural deficiency
#  ){
#      list(
#              y1 = as.numeric(theta[1]*xSparse + theta[2]*mean(xRich)/10)
#              ,y2 = as.numeric(theta[1]*xSparse[1] + theta[2]*pmax(0,xRich-thresholdCovar) )
#      )
#  }

## ----results='hide'------------------------------------------------------
 #thresholdCovar = 0.3    # the true value
thresholdCovar = 0      # the effective biased model that glosses over this threshold

## ------------------------------------------------------------------------
denSparse <- function( theta, twTwoDenEx=twTwoDenEx1, ... ){
    pred <- twTwoDenEx$fModel(theta, xSparse=twTwoDenEx$xSparse, xRich=twTwoDenEx$xRich, ...) 
    misfit <- twTwoDenEx$obs$y1 - pred$y1
    structure( -1/2 * sum((misfit/twTwoDenEx$sdObs$y1)^2), names='y1' )
}

denRich <- function( theta, twTwoDenEx=twTwoDenEx1, ... ){
    pred <- twTwoDenEx$fModel(theta, xSparse=twTwoDenEx$xSparse, xRich=twTwoDenEx$xRich, ...) 
    misfit <- twTwoDenEx$obs$y2 - pred$y2
    structure( -1/2 * sum((misfit/twTwoDenEx$sdObs$y2)^2), names='y2' )
}

## ----eval=FALSE, tidy=FALSE----------------------------------------------
#      dInfos=list(
#          dSparse=list(fLogDen=denSparse,
#               argsFLogDen=list(twTwoDenEx=twTwoDenEx1,
#               thresholdCovar=thresholdCovar)),
#          dRich=list(fLogDen=denRich,
#               argsFLogDen=list(twTwoDenEx=twTwoDenEx1,
#               thresholdCovar=thresholdCovar))
#      )

## ----eval=FALSE, tidy=FALSE----------------------------------------------
#      blocks = list(
#          a=list(compPos="a", dInfoPos="dSparse"),
#          b=list(compPos="b", dInfoPos="dRich")
#      )

## ----sa1, cache=TRUE, results='hide', tidy=FALSE-------------------------
resBlock <- resBlock0 <- twDEMCSA(
        theta=twTwoDenEx1$thetaTrue, 
        covarTheta=diag((twTwoDenEx1$thetaTrue*0.3)^2),
        dInfos=list(
            dSparse=list(fLogDen=denSparse,
                argsFLogDen=list(twTwoDenEx=twTwoDenEx1,
                thresholdCovar=thresholdCovar)),
            dRich=list(fLogDen=denRich, 
                argsFLogDen=list(twTwoDenEx=twTwoDenEx1,
                thresholdCovar=thresholdCovar))
        ),
        blocks = list(
            a=list(compPos="a", dInfoPos="dSparse"),
            b=list(compPos="b", dInfoPos="dRich")
        ),
        nObs = c( y1=length(twTwoDenEx1$obs$y1),
                  y2=length(twTwoDenEx1$obs$y2) ),
        nGen=2048
    )
    
resConv <- twDEMCBlock( resBlock0, nGen=ceiling(256 * resBlock$thin / getNChains(resBlock)), extendRun=FALSE )    

## ----plotCoda,spar=TRUE--------------------------------------------------
plot(as.mcmc.list(resBlock), smooth=FALSE)
 # plot(as.mcmc.list(resConv), smooth=FALSE)
ss <- stackChains(resConv)
summary(ss)

## ----plotRanks,spar=TRUE, fig.width=7.5, fig.height=1.7------------------
oldPar <- par(mfrow=c(1,5), par(mar=c(2.0,3.3,2.0,0)+0.3 ) )
colh <- heat.colors(nrow(ss))
plot( b ~ a, data=as.data.frame(ss[ rev(order(ss[,"dSparse"])),])
    , col=colh, main="rank dSparse" )
plot( b ~ a, data=as.data.frame(ss[ rev(order(ss[,"dRich"])),])
    , col=colh, main="rank dRich" )
plot( b ~ a, data=as.data.frame(ss[ rev(order(rowSums(ss[,c("dSparse","dRich")]))),])
    , col=colh, main="rank sum" )
plot( b ~ a, data=as.data.frame(ss[ rev(order(rowSums(apply(ss[,c("dSparse","dRich")],2,rank)) )),])
    , col=colh, main="sum ranks" )
plot( b ~ a, data=as.data.frame(ss[ orderLogDen(ss),])
    , col=colh, main="min ranks" )

## ------------------------------------------------------------------------
denSparseInt <- function( theta, intermediate=list(), twTwoDenEx=twTwoDenEx1, ... ){
    if( !length(intermediate) ){
        intermediate <- twTwoDenEx$fModel(theta, xSparse=twTwoDenEx$xSparse, xRich=twTwoDenEx$xRich, ...) 
    }
    pred <- intermediate  
    misfit <- twTwoDenEx$obs$y1 - pred$y1
    structure( -1/2 * sum((misfit/twTwoDenEx$sdObs$y1)^2), names='y1', intermediate=intermediate )
}

denRichInt <- function( theta, intermediate=list(), twTwoDenEx=twTwoDenEx1, ... ){
    if( !length(intermediate) ){
        intermediate <- twTwoDenEx$fModel(theta, xSparse=twTwoDenEx$xSparse, xRich=twTwoDenEx$xRich, ...) 
    }
    pred <- intermediate  
    misfit <- twTwoDenEx$obs$y2 - pred$y2
    structure( -1/2 * sum((misfit/twTwoDenEx$sdObs$y2)^2), names='y2', intermediate=intermediate )
}

## ----sa1Int, cache=TRUE, results='hide', tidy=FALSE----------------------
resBlock <- resBlockInt <- twDEMCSA(
        theta=twTwoDenEx1$thetaTrue, 
        covarTheta=diag((twTwoDenEx1$thetaTrue*0.3)^2),
        dInfos=list(
            dSparse=list(fLogDen=denSparseInt,
                argsFLogDen=list(twTwoDenEx=twTwoDenEx1,
                thresholdCovar=thresholdCovar)),
            dRich=list(fLogDen=denRichInt, 
                argsFLogDen=list(twTwoDenEx=twTwoDenEx1,
                thresholdCovar=thresholdCovar))
        ),
        blocks = list(
            a=list(compPos="a", dInfoPos="dSparse", intermediateId="resFModel"),
            b=list(compPos="b", dInfoPos="dRich", intermediateId="resFModel")
        ),
        nObs = c( y1=length(twTwoDenEx1$obs$y1),
                  y2=length(twTwoDenEx1$obs$y2) ),
        nGen=2048
    )
    
resConvInt <- twDEMCBlock( resBlockInt, extendRun=FALSE, 
                nGen=ceiling(256 * resBlock$thin / getNChains(resBlock)))    

## ------------------------------------------------------------------------
ssInt <- stackChains(resConvInt)
summary(ssInt)