# R/initZ.R In twDEMC: parallel DEMC

#### Documented in calcM0twDEMCconstrainCfStackconstrainNStackinitZtwDEMCNormalreplaceZinitCasesreplaceZinitNonFiniteLogDensreplaceZinitNonFiniteLogDensLastStep

initZtwDEMCNormal <- function(
### Generate an initial population of states for \code{\link{twDEMCBlockInt}}.
thetaPrior		##<< numeric vector (nParm) of point estimate of the mean
,covarTheta 	##<< numeric matrix (nParm x nParm) the covariance of parameters.
##<< Alternatively, can be given as vector (nParm) of the diagonal, i.e. variances, if all dimensions are independent
,nChainPop=4 	##<< number of chains to run
,nPop=2 		##<< number of independent populations among the chains
,m0=ceiling(calcM0twDEMC(length(thetaPrior),nChainPop)/(m0FiniteFac))
### number of initial states for each chain
,m0FiniteFac=1	##<< use a factor smaller than 1 to increase default m0 to account for only a portion of proposal results in finite densities
,doIncludePrior=TRUE
### If TRUE, then set last sample of chain 1 to the prior estimate,
### which might be already a kind of best estimates by an optimization.
){
# InittwDEMC

##seealso<<

##details<<
## There are several methods to establish the initial population for a twDEMC run. \itemize{
## \item{ drawing from a multivariate normal distribution: this method  }
## \item{ subsetting the result of a former twDEMC run: \code{\link{initZtwDEMCSub.twDEMC}} or sample matrix \code{\link{initZtwDEMCSub.matrix}}  }
## \item{ extending the result of a former twDEMC run to include more parameters: \code{\link{initZtwDEMCExt.twDEMC}}  }
## \item{ selecting the N closes points from a sequence of points in parameter space \code{\link{constrainNStack}}  }
## \item{ selecting the points inside a confindenc ellipsis in parameter \code{\link{constrainCfStack}}  }
## \item{ replacing cases with in initial proposals that yield non-finite density \code{\link{replaceZinitNonFiniteLogDens}}  }
## \item{ general method for replacing cases in initial proposals \code{\link{replaceZinitCases}}  }
##}

#?rmvnorm #in package mvtnorm
if( 0==length(thetaPrior))
stop("initZtwDEMCNormal: no parameters given")
if( 0==length(names(thetaPrior))){
warning("initZtwDEMCNormal: no names of parameters provided")
names(thetaPrior) <- paste("theta",1:length(thetaPrior),sep="")
}
nChains <- nChainPop*nPop
Zinit <- if( length(thetaPrior)==1 ){
abind( lapply( 1:nChains, function(i){ matrix(rnorm( m0, mean=thetaPrior, sd=covarTheta), dimnames=list(NULL,names(thetaPrior)) ) }), along=3 )
}else{
if( is.matrix(covarTheta))
abind( lapply( 1:nChains, function(i){ rmvnorm( m0, mean=thetaPrior, sigma=covarTheta ) }), along=3 )
else
# independent dimenstions, can sample each from rnorm
abind( lapply( 1:nChains, function(i){
#iComp=1
sapply( seq_along(thetaPrior), function(iComp){
rnorm(m0, mean=thetaPrior[iComp], sd=covarTheta[iComp])
})
}), along=3 )
}
# Set the last sample of chain 1 to the prior estimate, which might be already a kind of best estimates by an optimization.
if( doIncludePrior )
Zinit[m0,,1 ] <- thetaPrior
dimnames(Zinit) = list(steps=NULL,parms=names(thetaPrior),chains=NULL)
attr(Zinit,"nPop") <- nPop
Zinit
### a matrix of number of parameters by number of individuals (m0 x d x Npop), with d dimension of theta
}
attr(initZtwDEMCNormal,"ex") <- function(){
data(twLinreg1)
attach( twLinreg1 )

.nChainPop=4
.nPop=2
.nPar=length(theta0)
Zinit <- initZtwDEMCNormal( theta0, diag(sdTheta^2), nChainPop=.nChainPop, nPop=.nPop)
all.equal( c(calcM0twDEMC(.nPar,.nChainPop), .nPar, .nChainPop*.nPop), dim(Zinit) )

detach()
}

calcM0twDEMC <- function(
### Calculate appropriate number of cases for initializing twDEMC.
nPar,		##<< the number of parameters to estimate
nChainPop 	##<< the number of chains per population
){
##seealso<<

##details<< see terBraak 2006 and 2008
res <- max(4,ceiling((8*nPar)/nChainPop))
res
### length of each chain so that each population is initialized with 8*nPar cases
}

R.methodsS3::setMethodS3("initZtwDEMCSub","matrix", function(
### generates an appropriate initial sample of parameter vectors for twDEMC from subsampling a matrix
x					##<< the mcmc matrix to subsample (nCases x nParms)
,...
,vars=colnames(x)	##<< which variables to keep
,nChainPop=4		##<< number of chains per population
,nPop=1 			##<< number of populations
,m0=calcM0twDEMC(length(unique(vars)),nChainPop) ##<< number of required cases for initialization
){
# initZtwDEMCSub.matrix
##seealso<<
nChains <- nChainPop*nPop
rrc <- lapply(1:nChains, function(iChain){ sample.int( nrow(x), m0, replace=TRUE)} )
res <- abind( lapply( rrc, function(rr){ x[rr,vars ,drop=FALSE] } ),rev.along=0 )
dimnames(res) <- list( steps=NULL, parms=vars, chains=NULL )
attr(res,"nPop") <- nPop
res
### an array of dimension suitable for Zinit for twDEMCBlockInt
})
attr(initZtwDEMCSub.matrix,"ex") <- function(){
data(twLinreg1)
attach( twLinreg1 )

.nChainPop=4
.nPop=2
.nPar=length(theta0)
Zinit <- initZtwDEMCNormal( theta0, diag(sdTheta^2), nChainPop=.nChainPop, nPop=.nPop)
ss <- do.call( rbind, twMisc::twListArrDim( Zinit))	#stack the chains
ZinitSub <- initZtwDEMCSub(ss, nChainPop=.nChainPop, nPop=.nPop)
all.equal(dim(Zinit),dim(ZinitSub))

detach()
}

R.methodsS3::setMethodS3("initZtwDEMCSub","twDEMC", function(
### generates an appropriate initial sample of parameter vectors for twDEMC from subsampling a previous result
vtwdemc	##<< the twDEMC list to subsample
,...	## other parameters passed to initZtwDEMCSub.default, e.g. m0
,vars=colnames(vtwdemc$parms) ##<< which variables to keep ,nChainPop=getNChainsPop(vtwdemc) ,nPop=getNPops(vtwdemc) ){ ##seealso<< ## \code{\link{initZtwDEMCNormal}} Zinit1 <- stackChains(vtwdemc)[,-(1:getNBlocks(vtwdemc))] #one big chain initZtwDEMCSub.matrix(Zinit1,...,vars=vars,nChainPop=nChainPop,nPop=nPop) }) R.methodsS3::setMethodS3("initZtwDEMCExt","matrix", function( ### subsampling and extending twDEMC with new variables Zinit1 ##<< the matrix to subsample (nCases x nParms) ,... ,thetaPrior ##<< numeric vector: mu of multivariate gaussian prior distribtuion ,covarTheta ##<< numeric vector: sigma of multivariate gaussian prior distrbituion ,nChainPop=4 ,nPop=1 ,m0=calcM0twDEMC(length(thetaPrior), nChainPop) ### number of required cases ){ # initZtwDEMCExt ##details<< ## the variables in thetaPrior that are part of vtwdemc are subsampled ## the other variables are drawn from prior distribution ## assuming no correlations between variables present and absent in vtwdemc and if( is.null(names(thetaPrior)) ) stop("initZtwDEMCExtDEMCzsp: thetaPrior provided without names attribute") #mapping index thetaPrior to column index in pss1 pssInd <- as.numeric(lapply( 1:length(thetaPrior), function(i){ which( names(thetaPrior)[i] == colnames(Zinit1))})) mcPars <- which( !is.na(pssInd) ) extPars <- (1:length(thetaPrior))[ -mcPars ] nChain <- nChainPop*nPop Zinit <- array(NA, dim=c( m0, length(thetaPrior), nChain) ) if( length(mcPars) > 0 ){ iSteps <- sample(1:nrow(Zinit1), m0, replace=TRUE) ZinitMc <- abind( lapply( 1:nChain, function(i){ Zinit1[iSteps, pssInd[mcPars], drop=FALSE ] }), along=3 ) Zinit[,mcPars,] <- ZinitMc } if( length(extPars) > 0){ ZinitZtwDEMCExt <- initZtwDEMCNormal( thetaPrior[extPars], covarTheta[extPars,extPars,drop=FALSE], nChainPop=nChainPop, nPop=nPop,m0=m0 ) Zinit[,extPars,] <- ZinitZtwDEMCExt } dimnames(Zinit) <- list( steps=NULL, parms=names(thetaPrior), chains=NULL ) Zinit ##seealso<< ## \code{\link{initZtwDEMCNormal}} }) R.methodsS3::setMethodS3("initZtwDEMCExt","twDEMC", function( ### subsampling and extending twDEMC with new variables vtwdemc ##<< the twDEMC list to subsample ,... ##<< e.g. m0 ,thetaPrior ##<< numeric vector: mu of multivariate gaussian prior distribtuion ,covarTheta ##<< numeric vector: sigma of multivariate gaussian prior distrbituion ,nChainPop=getNChainsPop(vtwdemc) ,nPop=getNPops(vtwdemc) ){ ##seealso<< ## \code{\link{initZtwDEMCNormal}} Zinit1 <- stackChains(vtwdemc)[,-(1:getNBlocks(vtwdemc))] #one big chain initZtwDEMCExt.matrix( Zinit1, thetaPrior=thetaPrior, covarTheta=covarTheta, nChainPop=nChainPop, nPop=nPop, ... ) }) constrainNStack <- function( ### Subsetting a sequence of parameter vectors. Keeps only the n cases in pss1 that are closest to thetaPrior. pss1, ##<< numeric matrix: the stack to constrain, see details thetaPrior, ##<< the point in parameter spcae for which to select the closest values n=nrow(pss1)%/%4, ##<< the number of rows in output, defaults to 1/4 of nrow(pss1) vars=names(thetaPrior), ##<< names or indices to constrain thetaPrior and invCovarTheta invCovarTheta=if( 0<length(thetaPrior)) solve(cov(pss1[,names(thetaPrior[vars]),drop=FALSE])) else numeric(0), ### the inverse of the covaraince matrix for thetaPrior, defaults to inverse of sample covariance returnAlpha=FALSE ##<< switch to return also significance level ){ # constrainNStack ##seealso<< ## \code{\link{initZtwDEMCNormal}} ## \code{\link{constrainCfStack}} ##details<< ## pss1 is a matrix with columns corresponding to variables ## and rows corresponding to cases. ## It is typicalla the result of \code{\link{stackChains.twDEMC}} if( !(n < nrow(pss1)) ){ warning("constrainNStack: number of rows in pss1 is not greater than n, returning full pss1") return( pss1) } tmp.ind <- if( 0<length(thetaPrior)){ if( !(all(names(thetaPrior[vars]) %in% colnames(pss1))) ){ stop("constrainNStack: all components of thetaPrior[vars] must be in colnames(pss1)") } .invCovarTheta <- if( is.null(colnames(invCovarTheta)) | is.null(rownames(invCovarTheta)) ){ .lv <- length(vars) if( !identical( c(.lv,.lv), dim(invCovarTheta)) ) stop("supplied invCovarTheta without dimensions and nonmatching dimensions") invCovarTheta }else{ if( !all(vars %in% colnames(invCovarTheta)) | !all(vars %in% rownames(invCovarTheta)) ) stop("supplied invCovarTheta with dimnames not including all names of thetaPrior") invCovarTheta[vars,vars] } # calculate the distance and order by that tmp.d <- sapply(vars, function(var){ pss1[,var] - (thetaPrior[var])}) #plot( tmp.d[,1], tmp.d[,2] ) if( length(vars) > 1) tmp.mahalanobis <- apply( tmp.d, 1, function(tmp.d){ t(tmp.d) %*% .invCovarTheta %*% tmp.d }) else tmp.mahalanobis = tmp.d^2 * 1/var(pss1[,vars]) tmp.ind <- order(tmp.mahalanobis)[1:n] }else{ # no thetaPrior given: just subsample tmp.ind <- sample.int(nrow(pss1),n) } psc <- pss1[ tmp.ind, ] if( returnAlpha){ tmp.alpha <- if( 0 < length(thetaPrior)){ # calculate the corresponding alpha (p) of the confidence ellipsis # see http://www.stat.psu.edu/online/development/stat505/05_multnorm/04_multnorm_geom.html tmp.maxdist = tmp.mahalanobis[tmp.ind[n] ] tmp.alpha = pchisq( q=tmp.maxdist, df=length(vars) ) }else 1 list( res=psc, alpha=tmp.alpha) }else{ psc } ### pss1[closestValues,] ### if returnAlpha=TRUE then a list is returned with \describe{ ### \item{res}{ value as above } ### \item{alpha}{ the significance level of the corresponding confidence ellipsis }} } #mtrace(constrainNStack) #twUtestF(constrainNStack,test="test.allVars") constrainCfStack <- function( ### Subsetting a sequence of parameter vectors. Keeps only the the cases in pss1 that are inside a confidence ellipsis around thetaPrior. pss1, ##<< numeric matrix: the stack to constrain, see details thetaPrior, ##<< the point in parameter spcae for which to select the closest values alpha=0.95, ##<< the conficence range of the ellipsis vars=names(thetaPrior), ##<< names or indices to constrain thetaPrior and invCovarTheta invCovarTheta=if(0<length(thetaPrior)) solve(cov(pss1[,vars,drop=FALSE])) else numeric(0) ### the inverse of the covaraince matrix for thetaPrior, defaults to inverse of sample covariance ){ # constrainCfStack ##seealso<< ## \code{\link{initZtwDEMCNormal}} ## \code{\link{constrainNStack}} if( 0==length(thetaPrior)) return( pss1 ) # if no thetaPrior is given, just return the original matrix if( !(all(names(thetaPrior[vars]) %in% colnames(pss1))) ){ stop("constrainNStack: all components of thetaPrior[vars] must be in colnames(pss1)") } .invCovarTheta <- if( is.null(colnames(invCovarTheta)) | is.null(rownames(invCovarTheta)) ){ .lv <- length(vars) if( !identical( c(.lv,.lv), dim(invCovarTheta)) ) stop("supplied invCovarTheta without dimensions and nonmatching dimensions") invCovarTheta }else{ if( !all(vars %in% colnames(invCovarTheta)) | !all(vars %in% rownames(invCovarTheta)) ) stop("supplied invCovarTheta with dimnames not including all names of thetaPrior") invCovarTheta[vars,vars,drop=FALSE] } #first remove all values outside the 95% confidence ellipsis # see http://www.stat.psu.edu/online/development/stat505/05_multnorm/04_multnorm_geom.html tmp.d <- sapply(vars, function(var){ pss1[,var] - (thetaPrior[var])}) tmp.mahalanobis <- apply( tmp.d, 1, function(tmp.d){ t(tmp.d) %*% .invCovarTheta %*% tmp.d }) tmp.chisq = qchisq( p=alpha, df=length(vars) ) psc <- pss1[ tmp.mahalanobis < tmp.chisq, ] if( nrow(psc) < 10 ) stop(paste("constrained too strong, selected nRows=",nrow(psc))) psc ### pss1[ withinConfidenceInterval, ] } replaceZinitCases <- function( ### Replaces states of Zinit by sampling the other states that are marked good. Zinit, ##<< initial states of the form required by \code{\link{twDEMCBlockInt}} boMat ##<< boolean matrix specifying good cases with rows cases and columns chains. If it is a vector, then matrix is constructed using last dimension as the number of chains ){ ##seealso<< ## \code{\link{initZtwDEMCNormal}} ## \code{\link{replaceZinitNonFiniteLogDens}} if( is.vector(boMat)) boMat <- matrix(boMat, ncol=dim(Zinit)[3]) ##details<< ## Samples for the first half of chains are sampled from good cases of the second half of chains. ## Samples for the second half of chains are sampled from good cases of first half of chains. tmp.nChains <- dim(Zinit)[3] if( !identical(dim(Zinit)[c(1,3)], dim(boMat)) ) stop("dimensions of Zinti and boMat do not match") if( !is.logical(boMat) ) stop("boMat is not a logical matrix") goodCases <- which( boMat, arr.ind=TRUE ) goodCases1 <- goodCases[ goodCases[,2] <= tmp.nChains/2, ,drop=FALSE] #non-missing values of first half goodCases2 <- goodCases[ goodCases[,2] > tmp.nChains/2, ,drop=FALSE] #non-missing values of second half badCases <- which( !boMat, arr.ind=TRUE ) badCases1 <- badCases[ badCases[,2] <= tmp.nChains/2, ,drop=FALSE] #missing values of the first half badCases2 <- badCases[ badCases[,2] > tmp.nChains/2, ,drop=FALSE] #missing values of the second half #replace non-finite yielding states of first halv by sampling the good cases rows of the second half tmp.j = sample( 1:nrow(goodCases2), nrow(badCases1), replace=TRUE ) for( i in seq(along.with=badCases1[,1]) ) Zinit[badCases1[i,1], , badCases1[i,2] ] = Zinit[goodCases2[tmp.j[i],1], , goodCases2[tmp.j[i],2] ] #replace non-finite yielding states of second halv by sampling the good cases rows of the first half tmp.j1 = sample( 1:nrow(goodCases1), nrow(badCases2), replace=TRUE ) for( i in seq(along.with=badCases2[,1]) ) Zinit[badCases2[i,1], , badCases2[i,2] ] = Zinit[goodCases1[tmp.j1[i],1], , goodCases1[tmp.j1[i],2] ] Zinit ### Zinit, with several cols (parameter vectors) replaced by other cols } replaceZinitNonFiniteLogDens <- function( ### Replaces states of Zinit that yield non-finite rLogDen by sampling other states. Zinit, ##<< initial states see InitDEMCzsp logDen ##<< numeric matrix (nCases x nChains): calculated logDensitys for all the states in Zinit. If it is a vector then it is reshaped. ){ ##seealso<< ## \code{\link{initZtwDEMCNormal}} ## \code{\link{replaceZinitCases}} ##details<< ## In order for twDEMC to start up effectively, it is important that chains start from values, where the logDensity is finite bo <- is.finite(logDen) dim(bo) <- dim(Zinit)[c(1,3)] # steps x chains replaceZinitCases( Zinit, bo ) ### Zinit, with several cols (parameter vectors) replaced by other cols } #twUtestF(replaceZinitNonFiniteLogDens) attr(replaceZinitNonFiniteLogDens,"ex") <- function(){ data(twdemcEx1) data(twLinreg1) attach( twLinreg1 ) #mtrace(initZtwDEMCNormal) Zinit <- initZtwDEMCNormal( theta0, diag(4*sdTheta^2), nChainPop=8, nPop=2) dim(Zinit) res <- res0 <- twCalcLogDenPar(logDenGaussian, stackChains(Zinit) # chains stack to calculate in parallel ,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 )$logDen
plot(density(res))
res[res < -30] <- NA
resM <- matrix(res, ncol=dim(Zinit)[3] ) # restacking to chains work because it is the last dimension
set.seed(0815)
Zinit2 <- replaceZinitNonFiniteLogDens(Zinit, resM)
set.seed(0815)
Zinit3 <- replaceZinitNonFiniteLogDens(Zinit, res)
identical(Zinit2,Zinit3)

res2 <- twCalcLogDenPar(logDenGaussian, stackChains(Zinit2) # chains stack to calculate in parallel
,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
)$logDen d2 <- density(res2) lines( d2$x, d2$y * max(density(res0)$y)/max(d2$y), col="blue") } replaceZinitNonFiniteLogDensLastStep <- function( ### Replaces states of last step, i.e. column of Zinit that yield non-finite rLogDen by sampling other states. Zinit ##<< initial states see InitDEMCzsp ,fLogDen ##<< the logDen Function ,nPop=1 ##<< number of populations. States are only choosen from same population ,iStep=nrow(Zinit) ##<< the step for which to replace nonfinite yielding parameters. ,maxSteps=16 ##<< maximum number of iterations (when last row was replaced by another row yielding non-finite loglik) ,... ##<< arguments to \code{\link{twCalcLogDenPar}} ){ ##seealso<< ## \code{\link{initZtwDEMCNormal}} res <- twCalcLogDenPar( fLogDen, t(adrop(Zinit[iStep,,,drop=FALSE],2)), ... ) rLogDen <- res$logDen
logDenComp <- res$logDenComp iNonfinite <- which( !is.finite(rLogDen) ) # chains yielding non-finite logLik nReplace <- length(iNonfinite) if( nReplace>0 ){ # matrix finLogDen (nCase,nChain) markes the states which we do not want to sample again by FALSE dimZinit <- dim(Zinit) finLogDen <- matrix( TRUE, nrow=dimZinit[1], ncol=dimZinit[3] ) finLogDen[iStep,] <- FALSE #do not choose replacements from this step # make sure to sample a state from same population chainsPops <- matrix( 1:dimZinit[3], ncol=nPop ) # (nChainInPop x nPop) nChainPop <- nrow(chainsPops) rChain <- integer(nReplace) rStep <- integer(nReplace) i=1 while( (nReplace>0) & i<=maxSteps ){ j=1 for( j in 1:nReplace){ iPop <- ((iNonfinite[j]-1)%/%nChainPop)+1 chainsPop <- chainsPops[,iPop] rChain[j] <- sample( chainsPop, 1 ) stepsPop <- which( finLogDen[,rChain[j] ] ) rStep[j] <- sample(stepsPop , 1 ) Zinit[iStep,, iNonfinite[j] ] <- Zinit[rStep[j],,rChain[j] ] finLogDen[rStep[j],rChain[j] ] <- FALSE #do not sample those states again } res <- twCalcLogDenPar( fLogDen, t(adrop(Zinit[iStep,,iNonfinite,drop=FALSE],2)), ... ) rLogDen[iNonfinite] <- res$logDen
logDenComp[iNonfinite,] <- res\$logDenComp
iNonfinite <- which( !is.finite(rLogDen) )
nReplace <- length(iNonfinite)
i=i+1
}
if( nReplace > 0 )
warning("could not replace all non-finite states in row.")
}
##value<< list with components
list( ##describe<<
Zinit=Zinit					##<< argument \code{Zinit} with some states in last row replaced by other states from Zinit.
, rLogDen=rLogDen			##<< numeric vector: results of \code{\link{twCalcLogDenPar}} for last row for all chains
, logDenComp=logDenComp		##<< numeric matrix: results of \code{\link{twCalcLogDenPar}} for last row for all chains
) ##end<<
}
#twUtestF(replaceZinitNonFiniteLogDens)

.twDEMCStackToChains <- function(
x			##<< numeric matrix (nStep x nParm)
, nChain	##<< number of chains
){
##details<<
# if nRow(x) is not a multiple of nChain, then last rows are skipped.
nCases <- nrow(x) %/% nChain
resL <- lapply((1:nChain)-1, function(iChain0){
x[nCases*iChain0 + (1:nCases), ,drop=FALSE]
})
### list of subsets (by rows)
}
attr(.twDEMCStackToChains,"ex") <- function(){
x <- matrix( 1:16, ncol=2, dimnames=list(step=NULL,parms=c("a","b")))
tail(x)
tmp <- abind(.twDEMCStackToChains(x, 2), 3)
str(tmp)
tmp[,,2]
}


## Try the twDEMC package in your browser

Any scripts or data that you put into this service are public.

twDEMC documentation built on May 31, 2017, 3:44 a.m.