Nothing
R/create.jags.syntax.R
In ctglm: ctglm
Defines functions
create.jags.syntax
make.str
moveTo.par.env
make.inverse.1
make.inverse.2
create.jags.syntax <- function ( env ) {
# get variables from env
eval( parse( text=paste0( "assign( '",ls(envir=env), "' , get('",ls(envir=env),"', envir=env ) )" ) ) )
# objects with free parameters must be overwritten or deleted
# environment to save original objects
par.env <- new.env()
# syntax matrix (1 column)
x<-matrix( "model", 1, 1 )
x<-rbind(x, "{ ")
x<-rbind(x, " ")
x<-rbind(x, " ### DATA ### ")
x<-rbind(x, " ")
x<-rbind(x, " # data dl in long format ")
x<-rbind(x, " # dl[,col.id]: person ")
x<-rbind(x, " # dl[,col.item]: item ")
x<-rbind(x, " # dl[,col.time]: time point ")
x<-rbind(x, " # dl[,col.y]: responses ")
x<-rbind(x, " ")
x<-rbind(x, " ")
x<-rbind(x, " ### MEASUREMENT MODEL ### ")
x<-rbind(x, " ")
# x<-rbind(x, " # estimated parameters ")
# x<-rbind(x, " # theta: J (person) x F (faktor) x T (time point) ")
# x<-rbind(x, " ")
x<-rbind(x, " # loop over rows in long data set dl ")
x<-rbind(x, " for (r in 1:R) { ")
x<-rbind(x, " ")
x<-rbind(x, " # distributional assumption ")
if( measurement.model$family %in% "gaussian" ) {
x<-rbind(x, " dl[r,col.y] ~ dnorm( mu.y[r], E[ dl[r,col.item], dl[r,col.item] ] ) ") }
if( measurement.model$family %in% "binomial" ) {
x<-rbind(x, " dl[r,col.y] ~ dbern( mu.y[r] ) ") }
if( measurement.model$family %in% "poisson" ) {
x<-rbind(x, " dl[r,col.y] ~ dpois( mu.y[r] ) ") }
x<-rbind(x, " ")
x<-rbind(x, " # link ")
if( measurement.model$family %in% "gaussian" ) {
x<-rbind(x, " mu.y[r] <- eta[r] ") }
if( measurement.model$link %in% "logit" ) {
x<-rbind(x, " mu.y[r] <- ilogit( eta[r] ) ") }
if( measurement.model$link %in% "log" ) {
x<-rbind(x, " log( mu.y[r] ) <- eta[r] ") }
x<-rbind(x, " ")
x<-rbind(x, " # linear predictor ")
# x<-rbind(x, " eta[r] <- sum( Lambda[ dl[r,col.item], , dl[r,col.time] ] * theta[ dl[r,col.id], , dl[r,col.time] ] ) + beta[ dl[r,col.item], 1 ] ")
x<-rbind(x, " eta[r] <- sum( Lambda[ dl[r,col.item], ] * theta[ dl[r,col.id], , dl[r,col.time] ] ) + beta[ dl[r,col.item], 1 ] ")
# x<-rbind(x, " eta[r] <- sum( Lambda[ dl[r,col.item], ] * theta[ dl[r,col.id], , dl[r,col.time] ] ) ")
x<-rbind(x, " ")
x<-rbind(x, " } ")
x<-rbind(x, " ")
if ( exists( "E" ) ) {
x<-rbind(x, " # precision of measurement errors ")
x<-rbind(x, make.str( "E" ) ); if( any ( is.parameter( E ) ) ) invisible(moveTo.par.env("E",env,par.env)) else rm("E", envir=env)
x<-rbind(x, " ") }
# if( measurement.model$family %in% "gaussian" ) {
# x<-rbind(x, " # precision of measurement errors ")
# x<-rbind(x, " for (i in 1:I) { ")
# x<-rbind(x, " E[i] ~ dgamma(1,0.005) ")
# x<-rbind(x, " } ")
# x<-rbind(x, " ") }
x<-rbind(x, " # values/priors of item easiness ")
x<-rbind(x, make.str( "beta" ) ); if( any ( is.parameter( beta ) ) ) invisible(moveTo.par.env("beta",env,par.env)) else rm("beta", envir=env)
x<-rbind(x, " ")
if ( exists( "mu.beta" ) ) {
x<-rbind(x, " # values/priors of mean of item easiness ")
x<-rbind(x, make.str( "mu.beta" ) ); if( any ( is.parameter( mu.beta ) ) ) invisible(moveTo.par.env("mu.beta",env,par.env)) else rm("mu.beta", envir=env)
x<-rbind(x, " ") }
if ( exists( "mu.beta" ) ) {
x<-rbind(x, " # values/priors of precision of item easiness ")
x<-rbind(x, make.str( "prec.beta" ) ); if( any ( is.parameter( prec.beta ) ) ) invisible(moveTo.par.env("prec.beta",env,par.env)) else rm("prec.beta", envir=env)
x<-rbind(x, " ") }
x<-rbind(x, " # values/priors of Lambda ")
x<-rbind(x, make.str( "Lambda" ) ); if( any ( is.parameter( Lambda ) ) ) invisible(moveTo.par.env("Lambda",env,par.env)) else rm("Lambda", envir=env)
x<-rbind(x, " ")
# x<-rbind(x, " # standard deviation of latent scale ")
# x<-rbind(x, " #sd.eta <- sd(eta[]) ")
# x<-rbind(x, " ")
x<-rbind(x, " ")
x<-rbind(x, " ### CONTINUOUS TIME MODEL ### ")
x<-rbind(x, " ")
# browser()
if( !exists( "stationarity" ) ) stationarity <- FALSE
if( person.var["mu.t1"] ) {
if (! stationarity ){
x<-rbind(x, " ## personal t1 means ")
x<-rbind(x, " # mu.t1 varies over persons, i.e. for each person a separate mu.t1 ")
mu.t1.j. <- make.str( "mu.t1.j" ); if( any ( is.parameter( mu.t1.j ) ) ) invisible(moveTo.par.env("mu.t1.j",env,par.env)) else rm("mu.t1.j", envir=env)
# mods to mu.t1.j.
mu.t1.j.1 <- mu.t1.j.[ !grepl( "~", mu.t1.j., fixed=TRUE ),,drop=FALSE ]
mu.t1.j.2 <- mu.t1.j.[ grepl( "~", mu.t1.j., fixed=TRUE ),,drop=FALSE ]
if( F>1 ) mu.t1.j.2[,1] <- sub( "mu.t1.j\\[\\d+", "mu.t1.j[1:F", mu.t1.j.2[,1] )
if( F==1 ) mu.t1.j.2[,1] <- sub( "mu.t1.j\\[\\d+", "mu.t1.j[1", mu.t1.j.2[,1] )
mu.t1.j.2 <- mu.t1.j.2[!duplicated(mu.t1.j.2),,drop=FALSE]
if( !identical( mu.t1.j.1, character(0) ) ) x<-rbind(x,mu.t1.j.1)
if( !identical( mu.t1.j.2, character(0) ) ) x<-rbind(x,mu.t1.j.2)
x<-rbind(x, " ")
x<-rbind(x, " # values/prior for prec.mu.t1.j ")
x<-rbind(x, ifelse( exists("prec.mu.t1.j.prior") && is.null( dim(prec.mu.t1.j.prior) ), paste0( " prec.mu.t1.j[1:F,1:F] ~ ", prec.mu.t1.j.prior, " " ), make.str( "prec.mu.t1.j" ) ) ); if( any ( is.parameter( prec.mu.t1.j ) ) ) invisible(moveTo.par.env("prec.mu.t1.j",env,par.env)) else rm("prec.mu.t1.j", envir=env)
}
if (stationarity){
mu.t1.j. <- make.str( "mu.t1.j" )
mu.t1.j.1 <- mu.t1.j.[ !grepl( "~", mu.t1.j., fixed=TRUE ),,drop=FALSE ]
mu.t1.j.2 <- mu.t1.j.[ grepl( "~", mu.t1.j., fixed=TRUE ),,drop=FALSE ]
if( F>1 ) mu.t1.j.2[,1] <- sub( "mu.t1.j\\[\\d+", "mu.t1.j[1:F", mu.t1.j.2[,1] )
if( F==1 ) mu.t1.j.2[,1] <- sub( "mu.t1.j\\[\\d+", "mu.t1.j[1", mu.t1.j.2[,1] )
mu.t1.j.2 <- mu.t1.j.2[!duplicated(mu.t1.j.2),,drop=FALSE]
mu.t1.j.3 <- sub( "^(.*)~.*$", "\\1", mu.t1.j.2 )
mu.t1.j.4 <- matrix( paste0( mu.t1.j.3, " <- -1 * A.inv[,] %*% bj[,", as.integer( sub( "^.*,(.*)\\].*$", "\\1", mu.t1.j.3 ) ) ,"]" ), ncol=1 )
if( !identical( mu.t1.j.4, character(0) ) ) x<-rbind(x,mu.t1.j.4)
# push to parameter env, because we want results for that
if( any ( is.parameter( mu.t1.j ) ) ) invisible(moveTo.par.env("mu.t1.j",env,par.env)) else rm("mu.t1.j", envir=env)
}
}
x<-rbind(x, " ")
x<-rbind(x, " # loop over persons j ")
x<-rbind(x, " for (j in 1:J) { ")
x<-rbind(x, " ")
x<-rbind(x, " # first time point ")
if( !person.var["mu.t1"] ) {
x<-rbind(x, paste0( " theta[j,1:F,1] ~ d",ifelse(F>1,"m",""),"norm( mu.t1, prec.t1 ) ") ) }
if( person.var["mu.t1"] ) {
x<-rbind(x, paste0( " theta[j,1:F,1] ~ d",ifelse(F>1,"m",""),"norm( mu.t1.j[,j], prec.t1 ) ") ) }
x<-rbind(x, " ")
x<-rbind(x, " # loop over t=2,...,Tj personal time point ")
x<-rbind(x, " for (t in 2:Tj[j]) { ")
x<-rbind(x, " ")
# browser()
if( !person.var["b"] ) {
x<-rbind(x, paste0( " theta[j,1:F,t] ~ d",ifelse(F>1,"m",""),"norm( At[,,t-1,Lpat.group[j]] %*% theta[j,,t-1] + bt[,t-1,Lpat.group[j]], Qt.prec[,,t,Lpat.group[j]] ) ") ) }
if( person.var["b"] ) {
x<-rbind(x, paste0(" theta[j,1:F,t] ~ d",ifelse(F>1,"m",""),"norm( At[,,t-1,Lpat.group[j]] %*% theta[j,,t-1] + bt[,t-1,j], Qt.prec[,,t,Lpat.group[j]] ) ") ) }
x<-rbind(x, " ")
x<-rbind(x, " } ")
x<-rbind(x, " } ")
x<-rbind(x, " ")
x<-rbind(x, " ## values/prior of parameters ")
# x<-rbind(x, " ## Note: values are commented (already set in R) ")
x<-rbind(x, " ")
x<-rbind(x, " # values/prior of mean of first time point ")
if( !person.var["mu.t1"] || !stationarity ) {
x<-rbind(x, make.str( "mu.t1" ) ); if( any ( is.parameter( mu.t1 ) ) ) invisible(moveTo.par.env("mu.t1",env,par.env)) else rm("mu.t1", envir=env)
} else {
x<-rbind(x, " for (f in 1:F) {" )
x<-rbind(x, " mu.t1[f] <- mean( mu.t1.j[f,] )" )
x<-rbind(x, " }" )
if( any ( is.parameter( mu.t1 ) ) ) invisible(moveTo.par.env("mu.t1",env,par.env)) else rm("mu.t1", envir=env)
}
x<-rbind(x, " ")
# browser()
x<-rbind(x, " # values/prior of precision of first time point ")
# x<-rbind(x, ifelse( exists("prec.t1.prior") && is.null( dim(prec.t1.prior) ), paste0( " prec.t1[1:F,1:F] ~ ", prec.t1.prior, " " ), make.str( "prec.t1" ) ) ); if( any ( is.parameter( prec.t1 ) ) ) invisible(moveTo.par.env("prec.t1",env,par.env)) else rm("prec.t1", envir=env)
x<-rbind(x, make.str( "prec.t1" ) ); if( any ( is.parameter( prec.t1 ) ) ) invisible(moveTo.par.env("prec.t1",env,par.env)) else rm("prec.t1", envir=env)
# prec.t1 is symmetric, needs definition of upper to lower elements
# nasty workaround here
# prec.t1.prior <- prec.t1
# pstr <- make.str( "prec.t1" )
# pstr <- pstr[ grepl( "<-", pstr[,1] ), ,drop=FALSE]
# if( nrow( pstr ) > 0 ) {
# x<-rbind(x,pstr)
# }
x<-rbind(x, " ")
x<-rbind(x, " # values/prior of drift matrix ")
x<-rbind(x, make.str( "A" ) ); if( any ( is.parameter( A ) ) ) invisible(moveTo.par.env("A",env,par.env)) else rm("A", envir=env)
x<-rbind(x, " ")
x<-rbind(x, " # values/prior of diffusion matrix ")
# browser()
x<-rbind(x, make.str( "Q" ) ); if( any ( is.parameter( Q ) ) ) invisible(moveTo.par.env("Q",env,par.env)) else rm("Q", envir=env)
x<-rbind(x, " ")
# browser()
if( person.var["b"] ) {
x<-rbind(x, " ## continuous time intercepts ")
x<-rbind(x, " # b varies over persons, i.e. for each person a separate b ")
# x<-rbind(x, " for( j in 1:J ) { ")
# x<-rbind(x, " bj[1:F,j] ~ dmnorm( b[,1], prec.b[,] ) ")
# x<-rbind(x, " } ") }
bj. <- make.str( "bj" ); if( any ( is.parameter( bj ) ) ) invisible(moveTo.par.env("bj",env,par.env)) else rm("bj", envir=env)
# mods to bj.
bj.1 <- bj.[ !grepl( "~", bj., fixed=TRUE ),,drop=FALSE ]
bj.2 <- bj.[ grepl( "~", bj., fixed=TRUE ),,drop=FALSE ]
if( F>1 ) bj.2[,1] <- sub( "bj\\[\\d+", "bj[1:F", bj.2[,1] )
if( F==1 ) bj.2[,1] <- sub( "bj\\[\\d+", "bj[1", bj.2[,1] )
bj.2 <- bj.2[!duplicated(bj.2),,drop=FALSE]
if( !identical( bj.1, character(0) ) ) x<-rbind(x,bj.1)
if( !identical( bj.2, character(0) ) ) x<-rbind(x,bj.2) }
x<-rbind(x, " # values/prior of continuous time intercepts ")
x<-rbind(x, make.str( "b" ) ); if( any ( is.parameter( b ) ) ) invisible(moveTo.par.env("b",env,par.env)) else rm("b", envir=env)
if( person.var["b"] ) {
x<-rbind(x, " # values/prior of precision of b ")
# x<-rbind(x, ifelse( exists("prec.b.prior") && is.null( dim(prec.b.prior) ), paste0( " prec.b[1:F,1:F] ~ ", prec.b.prior, " " ), make.str( "prec.b" ) ) ); if( any ( is.parameter( prec.b ) ) ) invisible(moveTo.par.env("prec.b",env,par.env)) else rm("prec.b", envir=env) }
# browser()
x<-rbind(x, make.str( "prec.b" ) ); if( any ( is.parameter( prec.b ) ) ) invisible(moveTo.par.env("prec.b",env,par.env)) else rm("prec.b", envir=env) }
x<-rbind(x, " ")
x<-rbind(x, " ")
x<-rbind(x, " # rowQ: stack Q row-wise into a column vector ")
x<-rbind(x, " for ( z in (1:Aw) ) { ")
x<-rbind(x, " for ( c in (1:Aw) ) { ")
x<-rbind(x, " rowQ[ c+Aw*(z-1), 1 ] <- Q[z,c] ")
x<-rbind(x, " } ")
x<-rbind(x, " } ")
x<-rbind(x, " ")
x<-rbind(x, " ## discrete/continuous time relations ")
x<-rbind(x, " # Qtv for first time point ")
x<-rbind(x, " # for technical reasons, same values are set for all lag patterns (not necessary) ")
x<-rbind(x, " for(p in 1:P) { ")
x<-rbind(x, " Qtv[1:(F*F),1,p] <- ( -1 * Ah.inv[,] ) %*% rowQ[,1] ")
x<-rbind(x, " } ")
# mexp or exp
# browser()
if( F>1 ) fexp <- "mexp" else if ( F== 1 ) fexp <- "exp"
x<-rbind(x, " # loop over p=1,...,P lag patterns ")
x<-rbind(x, " for(p in 1:P) { ")
x<-rbind(x, " # loop over t=2,...,Tp pattern-specific time point ")
x<-rbind(x, " for (t in 2:Tp[p]) { ")
x<-rbind(x, " ")
x<-rbind(x, " # autoregression matrix / drift matrix ")
x<-rbind(x, paste0( " At[1:F,1:F,t-1,p] <- ",fexp,"( A[,] * Lpat[p,t-1] ) ") )
x<-rbind(x, " ")
x<-rbind(x, " # continuous time intercepts ")
if( !person.var["b"] ) {
x<-rbind(x, paste0( " bt[1:F,t-1,p] <- ( A.inv[,] %*% ( ",fexp,"( A[,] * Lpat[p,t-1] ) - I1 ) ) %*% b[,1] ") ) }
if( person.var["b"] ) {
x<-rbind(x, " # b varies over persons, so loop over persons in specific pattern")
x<-rbind(x, " # PNj ... number of persons in pattern ")
x<-rbind(x, " # Pj (P x J) ... ragged matrix, indicators of persons in pattern")
x<-rbind(x, " for (j in 1:PNj[p]) { ")
x<-rbind(x, paste0( " bt[1:F,t-1,Pj[p,j]] <- ( A.inv[,] %*% ( ",fexp,"( A[,] * Lpat[p,t-1] ) - I1 ) ) %*% bj[,j]") )
x<-rbind(x, " } ") }
x<-rbind(x, " ")
x<-rbind(x, " # Qtv is vectorized Qt matrix ")
x<-rbind(x, paste0( " #Qtv[1:(F*F),t,p] <- ifelse( t==1, ( -1 * Ah.inv[,] ) %*% rowQ[,1] , ( Ah.inv[,] %*% ( ",fexp,"( Ah[,] * Lpat[p,t-1] ) - I2 ) ) %*% rowQ[,1] ) ") )
x<-rbind(x, paste0( " Qtv[1:(F*F),t,p] <- Ah.inv[,] %*% ( ",fexp,"( Ah[,] * Lpat[p,t-1] ) - I2 ) %*% rowQ[,1] ") )
x<-rbind(x, " } ")
x<-rbind(x, " } ")
x<-rbind(x, " ")
x<-rbind(x, " ## calculate Qt.prec ")
x<-rbind(x, " for(p in 1:P) { ")
x<-rbind(x, " for (t in 1:Tp[p]) { ")
x<-rbind(x, " ")
x<-rbind(x, " # Qt matrix (create with 'irow') ")
x<-rbind(x, " for ( z in (1:Aw) ) { ")
x<-rbind(x, " for ( c in (1:Aw) ) { ")
x<-rbind(x, " Qt[z,c,t,p] <- Qtv[ c+Aw*(z-1),t,p] ")
x<-rbind(x, " } ")
x<-rbind(x, " } ")
x<-rbind(x, " ")
# inverse of Qt; Qt has dimensions of Q; Qt.replace must exist
x<-rbind(x, make.inverse.2( "Qt", dim(Q)[1] ) )
x<-rbind(x, " } ")
x<-rbind(x, " } ")
x<-rbind(x, " ")
# inverse of A and Ah
x<-rbind(x, make.inverse.1( "A", dim(A)[1] ) )
x<-rbind(x, make.inverse.1( "Ah", dim(A)[1]^2 ) )
x<-rbind(x, " ")
x<-rbind(x, " # Kronecker of A and I1 ")
x<-rbind(x, " # https://sourceforge.net/p/mcmc-jags/discussion/610037/thread/17ddb2ce/ ")
x<-rbind(x, " for (i in 1:Aw){ # Aw = width of A ")
x<-rbind(x, " for (j in 1:Aw){ ")
x<-rbind(x, " for (k in 1:I1w){ # I1w = width of I1 ")
x<-rbind(x, " for (l in 1:I1w){ ")
x<-rbind(x, " Ah1[(i-1)*I1w+k,(j-1)*I1w+l] <- A[i,j]*I1[k,l] ")
x<-rbind(x, " } ")
x<-rbind(x, " } ")
x<-rbind(x, " } ")
x<-rbind(x, " } ")
x<-rbind(x, " # Kronecker of I1 and A ")
x<-rbind(x, " for (i in 1:I1w){ # I1w = width of I1 ")
x<-rbind(x, " for (j in 1:I1w){ ")
x<-rbind(x, " for (k in 1:Aw){ # Aw = width of A ")
x<-rbind(x, " for (l in 1:Aw){ ")
x<-rbind(x, " Ah2[(i-1)*Aw+k,(j-1)*Aw+l] <- I1[i,j]*A[k,l] ")
x<-rbind(x, " } ")
x<-rbind(x, " } ")
x<-rbind(x, " } ")
x<-rbind(x, " } ")
x<-rbind(x, " # Ah matrix ")
x<-rbind(x, " Ah[1:(I1w*Aw),1:(I1w*Aw)] <- Ah1[,] + Ah2[,] ")
x<-rbind(x, " ")
# browser()
# additional jags syntax
if( exists( "jags.add" ) && !is.null( jags.add ) ) {
if ( !is.matrix( jags.add ) ) jags.add <- matrix( jags.add, ncol=1 )
jags.add <- rbind( "## Additional/user code", jags.add )
jags.add <- rbind( jags.add, "" )
x<-rbind(x, jags.add)
}
x<-rbind(x, "} ")
# operating system
os <- Sys.info()["sysname"]
# parallelization yes/no
parall <- exists("cores") && !is.null(cores) && cores > 1 && os %in% c("Windows","Linux")
### call matrix (1 column)
y<-matrix( paste0( "### R syntax for ", model.name ), 1, 1 )
y<-rbind(y,paste0( "### engine: ", engine ) )
y<-rbind(y, "" )
y<-rbind(y, "# packages" )
y<-rbind(y, "require( 'rjags' )" )
if( parall ){
y<-rbind(y, "require( 'doParallel' ) # function foreach %dopar% ")
y<-rbind(y, "require( 'abind' ) # required by ctglm ")
y<-rbind(y, "#require( 'ctglm' ) # function mcarray.chains.combine ") }
y<-rbind(y, "" )
y<-rbind(y, "# package versions" )
y<-rbind(y, paste0( "print( installed.packages()[ installed.packages()[,1] %in% c(",paste( paste0( "'", na.omit( c('rjags' , ifelse( parall, 'ctglm', NA) ) ), "'" ) , collapse="," ),"), c(1,3) ] )" ) )
y<-rbind(y, "" )
if( !parall ){
y<-rbind(y, "# JAGS Modules" )
if( F>1 ) y<-rbind(y, "load.module('msm') # for matrix exponential, mexp()" )
y<-rbind(y, "load.module('glm') # for better glm sampler " )
y<-rbind(y, "" ) }
# create starting values
# function make.priors with mode="startingvalue" is used
if ( exists("A") && any ( is.parameter( A ) ) ) invisible( make.priors( "A", A, diag.prior = -0.5, offdiag.prior = 0, env=environment(), mode="startingvalue", verbose=FALSE ) )
# browser()
if ( exists("Q") && any ( is.parameter( Q ) ) ) invisible( make.priors( "Q", Q, diag.prior = 0.5, offdiag.prior = 0, env=environment(), mode="startingvalue", verbose=FALSE ) )
if ( exists("b") && any ( is.parameter( b ) ) ) invisible( make.priors( "b", b, prior = 0, env=environment(), mode="startingvalue", verbose=FALSE ) )
if ( exists("beta") && any ( is.parameter( beta ) ) ) invisible( make.priors( "beta", beta, prior = 0, env=environment(), mode="startingvalue", verbose=FALSE ) )
if ( exists("mu.beta") && any ( is.parameter( mu.beta ) ) ) invisible( make.priors( "mu.beta", mu.beta, prior = 0, env=environment(), mode="startingvalue", verbose=FALSE ) )
if ( exists("prec.beta") && any ( is.parameter( prec.beta ) ) ) invisible( make.priors( "prec.beta", prec.beta, diag.prior = 1, offdiag.prior = 0, env=environment(), mode="startingvalue", verbose=FALSE ) )
if ( !stationarity && exists("mu.t1") && any ( is.parameter( mu.t1 ) ) ) invisible( make.priors( "mu.t1", mu.t1, prior = 0, env=environment(), mode="startingvalue", verbose=FALSE ) )
if ( exists("prec.t1") && any ( is.parameter( prec.t1 ) ) ) {
invisible( make.priors( "prec.t1", prec.t1, diag.prior = 1, offdiag.prior = 0,env=environment(), mode="startingvalue", verbose=FALSE ) )
# !!!! because of dwish, complete starting value is needed, so assign upper with values of lower triangle
prec.t1.startingvalue[upper.tri(prec.t1.startingvalue)] <- prec.t1.startingvalue[lower.tri(prec.t1.startingvalue)]
}
if ( exists("prec.mu.t1.j") && any ( is.parameter( prec.mu.t1.j ) ) ) {
invisible( make.priors( "prec.mu.t1.j", prec.mu.t1.j, diag.prior = 1, offdiag.prior = 0,env=environment(), mode="startingvalue", verbose=FALSE ) )
# !!!! because of dwish, complete starting value is needed, so assign upper with values of lower triangle
prec.mu.t1.j.startingvalue[upper.tri(prec.mu.t1.j.startingvalue)] <- prec.mu.t1.j.startingvalue[lower.tri(prec.mu.t1.j.startingvalue)]
}
# list with starting values
startingvalues <- list()
do <- do.call( "c", sapply( c("A","Q","b","beta","mu.beta","prec.beta","mu.t1","prec.t1"), function(x) if (exists(paste0(x,".startingvalue"))) paste0( "startingvalues$", x ," <- ", x,".startingvalue" ) else NULL, simplify=FALSE ) )
if ( length( do > 0 ) ) eval( parse( text=do ) )
# push starting values to data environment
# and create code to create as many starting value lists as chains
# (chains are set in next step, run model)
if ( length(startingvalues) > 0 ){
assign( "startingvalues", startingvalues, envir=env )
y<-rbind(y, "# starting values: create as many starting value lists as chains")
y<-rbind(y, "startingvalues <- get( 'startingvalues', envir=data.env )")
y<-rbind(y, "inits <- sapply( 1:chains, function (n,l) l, startingvalues, simplify=FALSE )")
y<-rbind(y, "" )
}
# seeds
y<-rbind(y, "# seeds for each chain")
y<-rbind(y, "seeds <- parallel.seeds('base::BaseRNG', chains)")
if ( length(startingvalues) > 0 ){
y<-rbind(y, "inits <- mapply ( function ( l, k ) c(l,k), inits, seeds, SIMPLIFY=FALSE )")
} else {
y<-rbind(y, "inits <- seeds")
}
y<-rbind(y, "" )
y<-rbind(y, "# start time ")
y<-rbind(y, "start <- Sys.time() ")
y<-rbind(y, "" )
# browser()
# parallelization
indent <- ifelse( parall, " ", "" )
if( parall ){
y<-rbind(y, "## parallelization ")
y<-rbind(y, "# create cluster ")
if( os %in% "Windows" ){
y<-rbind(y, paste0( "cl <- makePSOCKcluster(",cores,") ")) }
if( os %in% "Linux" ){
y<-rbind(y, paste0( "cl <- makeForkCluster(",cores,") ")) }
y<-rbind(y, "# register cluster ")
y<-rbind(y, "registerDoParallel(cl) ")
y<-rbind(y, "" )
y<-rbind(y, "# parallel chains ")
y<-rbind(y, "res.l <- foreach(chain=1:chains, .packages='rjags') %dopar% { ")
y<-rbind(y, "" )
y<-rbind(y, paste0( indent, "# JAGS Modules (need to be loaded here in parallel mode)" ) )
if( F>1 ) y<-rbind(y, paste0( indent, "load.module('msm') # for matrix exponential, mexp()" ) )
y<-rbind(y, paste0( indent, "load.module('glm') # for better glm sampler " ) )
y<-rbind(y, "" )
}
y<-rbind(y, paste0( indent, "# initialization/adaptation ") )
# globalenv !!!
# y<-rbind(y, "eval(parse(text=paste0( name, ".ini <- jags.model ( file = mf , data=globalenv(), inits=sL, n.chains = ",chains,", n.adapt=",adapt,", quiet=FALSE )" )),envir=globalenv() )" )
# y<-rbind(y, "eval(parse(text=paste0( name, ".ini <- jags.model ( file = mf , data=globalenv(), n.chains = ",chains,", n.adapt=",adapt,", quiet=FALSE )" )),envir=globalenv() )" )
# y<-rbind(y, "eval(parse(text=paste0( model.name, .ini <- jags.model ( file = bugs.file , data=data.env, n.chains=chains, n.adapt=adapt, quiet=FALSE ) )),envir=data.env )" )
# y<-rbind(y, "eval(parse(text=paste0( model.name, .ini <- jags.model ( file = bugs.file , data=data.env, n.chains=chains, n.adapt=adapt, quiet=FALSE ) )),envir=globalenv())" )
y<-rbind(y, paste0( indent, "ini <- jags.model ( file = bugs.file , data=data.env, n.chains=",ifelse(parall,"1","chains"),", n.adapt=adapt, inits=inits",ifelse(parall,"[chain]",""),", quiet=FALSE )" ) )
y<-rbind(y, "" )
y<-rbind(y, paste0( indent, "# run ") )
# y<-rbind(y, "eval(parse(text=paste0( model.name, .res <- jags.samples ( ',model.name,'.ini , variable.names=c('A','Q','b'), n.iter=iter, thin=thin, type='trace' , progress.bar = 'none', by=20 )' )), envir=globalenv() )" )
# browser()
## add or delete person par
if( exists( "track.person.par" ) && "theta" %in% track.person.par ) {
# theta is not regularly created (makes no sense, since always estimated, never set)
# so create a structure for theta for further processing
thetas.dfr <- Reduce(function(x, y) merge(x, y, all=TRUE),list(data.frame("J"=1:J),data.frame("F"=1:F),data.frame("T"=1:T)),accumulate=FALSE )
thetas <- apply( thetas.dfr, 1, function(z) paste0( "theta_", paste( z, collapse="_" ) ) )
theta. <- array( thetas, dim=c(J,F,T) )
assign( "theta", theta., envir=par.env )
}
if( !exists( "track.person.par" ) || !"bj" %in% track.person.par ) {
# bj is regularly created, if it should not be tracked, delete it from par.env
if ( exists( "bj", envir=par.env ) ) rm( "bj", envir=par.env )
}
if( !exists( "track.person.par" ) || !"mu.t1.j" %in% track.person.par ) {
if ( exists( "mu.t1.j", envir=par.env ) ) rm( "mu.t1.j", envir=par.env )
}
# mod7: delete code in parameters
invisible( make.matrices.consistent( par.env , mods=7 ) )
# all parameters from par.env must be tracked
pars <- ls(envir=par.env)
# browser()
# mix in track.include / track.exclude
if ( exists( "track.include" ) && is.list( track.include ) && length( track.include ) > 0 ) {
# include in tracking
pars <- unique( c( pars, names( track.include ) ) )
# for automatic results retrieval structure of objects need to be known
# if not user set (i.e. NULL) it is assumed as a scalar
# in this case there is one parameter of the same name
if( any( isnull <- sapply( track.include, is.null ) ) ) track.include[ isnull ] <- names( track.include[ isnull ] )
# put into par.env
eval( parse( text=paste0( "assign( '", names(track.include), "' , track.include$'", names(track.include), "' , envir=par.env )" ) ) )
}
if ( exists( "track.exclude" ) && is.character( track.exclude ) && length( track.exclude ) > 0 ) {
# exclude from tracking
if ( any ( del <- pars %in% track.exclude ) ) pars <- pars[ !del ]
# also delete from par.env (ctglm.results depends on clean par.env)
if ( any ( del <- ls(envir=par.env) %in% track.exclude ) ) rm( list=ls(envir=par.env)[del], envir=par.env )
}
# sortieren
### ist glaub ich ziemlich sinnlos, da jags eh nicht die Reihenfolge hat dann
ord <- match( c("A","Q","b"), pars )
ord <- ord[!is.na(ord)]
if( !identical( ord, integer(0) ) ) pars <- c( pars[ord], pars[ !pars %in% pars[ord] ] )
par.string <- paste0( "c(", paste( paste0("'", pars, "'"), collapse="," ) ,")" )
y<-rbind(y, paste0( indent, "res <- jags.samples ( ini, variable.names=",par.string,", n.iter=iter, thin=thin, type='trace' , progress.bar = 'text', by=20 )" ) )
# parallelization, get results
if( parall ){
y<-rbind(y, "" )
y<-rbind(y, "}" )
y<-rbind(y, "" )
y<-rbind(y, "# stop cluster" )
y<-rbind(y, "stopCluster(cl)" )
y<-rbind(y, "" )
y<-rbind(y, "# combine results from parallel computing" )
y<-rbind(y, "res <- mcarray.chains.combine( res.l )" ) }
y<-rbind(y, "" )
y<-rbind(y, "# run time ")
y<-rbind(y, "runtime <- Sys.time() - start ")
### (over)write relevant variables to environment ###
# jags.syntax <- x
# obj <- c( "jags.syntax" )
# eval( parse ( text=paste0( "assign( '",obj, "' , get('",obj,"') , envir=env )" ) ) )
## create return object
ret <- list()
# first entry: engine
ret$engine <- engine
# second entry: model.name
ret$model.name <- model.name
# third entry: data environment
ret$data.env <- env
# fourth entry: environment with original matrices
ret$par.env <- par.env
# fifth entry: call
ret$call <- y
# sixth entry: syntax
ret$syntax <- x
# return
return( ret )
}
make.str <- function( y.name ) {
# browser()
# environment where objects are
env <- parent.frame()
# dimensions of structure, get from parent env
dim.y <- eval( parse( text=paste0( "dim(", y.name, ")" ) ), envir=env )
# if structure without dimensions (e.g. vector)
if( is.null( dim.y ) ) {
# long structure
y. <- data.frame( matrix( seq( along= eval( parse( text=paste0( "", y.name, "" ) ), envir=env ) ), ncol=1 ) )
} else {
# if structure with dimensions (e.g. matrix)
# long structure
y. <- eval(parse(text=paste0( "Reduce(function(x, y) merge(x, y, by=NULL, all=TRUE),list(", paste( paste0( "1:", dim.y ), collapse="," ), "),accumulate=FALSE )" )))
}
y.$par <- apply( y., 1, function ( z ) eval( parse( text= paste0( y.name,"[", paste(z,collapse=","), "]" ) ), envir=env ) )
# browser()
# sort out code, keep for later adding
if ( any ( is.code( y.$par ) ) ){
code <- sub( "^[^;]*;(.*)$", "\\1", y.$par[ is.code( y.$par ) ] )
add.code <- matrix( paste0( " ", gsub( "^\\s*", "", do.call( "c", strsplit(code,";") ) ) ), ncol=1 )
y.$par[ is.code( y.$par ) ] <- sub( "^([^;]*);.*$", "\\1", y.$par[ is.code( y.$par ) ] )
} else {
add.code <- NULL
}
# browser()
# tag duplicated free parameters and generate values to set duplicated
y.$dupl <- as.integer( duplicated(y.$par) & is.na(suppressWarnings(as.numeric(y.$par))) )
dupl <- unique( y.$par[ duplicated( y.$par ) & is.na(suppressWarnings(as.numeric(y.$par))) ] )
if ( length( dupl ) > 0 ) {
dl.l <- sapply( dupl, function( dl ){
y.1 <- y.[ y.$par %in% dupl, ]
y.1 <- y.1[ !duplicated( y.1$par ), ]
y.1$val <- apply( y.1, 1, function ( z ) paste0( y.name, "", "[", paste( z[-c(length(z)-1,length(z))], collapse=","), "]" ) )
return( y.1 )
}, simplify=FALSE )
dl <- do.call( "rbind", dl.l )
y.cn <- colnames( y. )
y.$order <- seq(along=rownames(y.))
y. <- merge( y., dl[,c("par","val")], by="par", sort=FALSE, all.x=TRUE )
y. <- y.[ order( y.$order), ]
y.$order <- NULL
y. <- y.[ , c( y.cn, "val" ), drop=FALSE ]
} else {
y.$val <- ""
}
# if( any( grepl( "Q", y.$par ) ) ) browser()
# loop over long structure, generate strings
make.str2 <- function( z ) {
# browser()
# if( grepl( "Q", z["par"] ) ) browser()
dupl <- z[ length(z)-1 ]
val <- z[ length(z) ]
par <- z["par"]
z <- z[ -c(length(z)-2,length(z)-1,length(z)) ]
if ( dupl %in% "0" ){
# if( grepl( "Q", par ) ) browser()
# name of parameter
nam <- paste0( y.name, "[", paste( z, collapse=","), "]" )
# value of parameter
par.val <- eval( parse( text=nam ), envir=env )
# value of prior
p.name <- paste0( y.name, ".prior" )
if ( exists( p.name, envir=env) ){
p.val <- eval( parse( text=paste0( p.name, "[", paste( z, collapse=","), "]" ) ), envir=env )
} else {
p.val <- NA
}
is.prior <- !is.na( p.val )
# determine if fixed value or freely estimated
# as.num <- suppressWarnings( as.numeric( eval( parse( text=nam ), envir=env ) ) )
# is.fixed <- !is.na( as.num )
# is fixed value or a parameter
# browser()
# operator: <- (for fixed oder parameter) or ~ (for prior)
# op <- ifelse( is.parameter( par.val ) | is.fixed( par.val ), " <- ", " ~ " )
op <- ifelse( is.prior, " ~ ", " <- " )
# if fixed, uncomment in jags syntax, because already set in R
# com <- ifelse( is.fixed, "# ", "" )
# com <- ""
# browser()
# value to be set
set.val <- "TODO"
# parameter
if( is.parameter( par.val ) ) set.val <- par
# fixed value (numeric)
if( !is.na( suppressWarnings( as.numeric( par ) ) ) ) set.val <- as.numeric( par )
# prior
if( is.prior ) set.val <- p.val
# val <- ifelse( is.fixed, as.character( as.num ), eval( parse( text=paste0( y.name, ".prior", "[", paste( z, collapse=","), "]" ) ), envir=env ) )
# val <- ifelse( is.fixed, "", eval( parse( text=paste0( y.name, ".prior", "[", paste( z, collapse=","), "]" ) ), envir=env ) )
# if( is.parameter( par.val ) ) set.val <- par
# if prior exists dann prior setzen
# get("Q",envir=env)
# get("Q.prior",envir=env)
# string
# s <- paste0( " ", com, nam, op, set.val )
s <- paste0( " ", nam, op, set.val )
} else {
s <- paste0( " ",y.name, "", "[", paste( z, collapse=","), "] <- ", val )
}
# return
return( s )
}
# browser()
# s <- apply( y.1[,-ncol(y.1),drop=FALSE], 1, make.str2 )
s <- apply( y., 1, make.str2 )
# append some white space (for optical reasons)
s <- sapply( s, function( x ) paste0( x, paste( rep( " ", 69 - nchar(x) ), collapse="" ) ) )
# as column vector
s <- matrix( s, ncol=1 )
# add.code
if ( !is.null( add.code ) ){
s <- rbind( s, add.code )
}
# return
return( s )
}
moveTo.par.env <- function( name, env, par.env ){
if( exists(name,envir=env) ){
assign( name, get( name, envir=env ), envir=par.env )
rm( list=name, envir=env )
}
TRUE
}
make.inverse.1 <- function( name, dim ){
x<-matrix( "", 1, 1 )
x<-rbind(x, paste0( " ## inverse of ",name," " ) )
if (dim %in% 1) {
x<-rbind(x, paste0( " ",name,".inv[1,1] <- (1/",name,"[1,1])" ) )
}
if (dim %in% 2) {
x<-rbind(x, paste0( " # determinant of ",name," " ) )
x<-rbind(x, paste0( " ",name,".det <- ",name,"[1,1]*",name,"[2,2] - ",name,"[2,1]*",name,"[1,2] " ) )
x<-rbind(x, paste0( " # ",name," is invertible if determinant is not 0 " ) )
x<-rbind(x, paste0( " ",name,".invertible <- ifelse( ",name,".det==0, 0, 1 ) " ) )
x<-rbind(x, paste0( " # if not invertible, modest mod of determinant " ) )
x<-rbind(x, paste0( " ",name,".detmod1 ~ dbern(0.5) " ) )
x<-rbind(x, paste0( " ",name,".detmod2 <- ifelse( ",name,".invertible==1, 0, (",name,".detmod1-0.5)/500 ) " ) )
x<-rbind(x, paste0( " ",name,".det.adj <- ",name,".det + ",name,".detmod2 " ) )
x<-rbind(x, paste0( " # inverse of ",name," " ) )
x<-rbind(x, paste0( " ",name,".help[1,1] <- ",name,"[2,2] " ) )
x<-rbind(x, paste0( " ",name,".help[2,2] <- ",name,"[1,1] " ) )
x<-rbind(x, paste0( " ",name,".help[1,2] <- -1*",name,"[1,2] " ) )
x<-rbind(x, paste0( " ",name,".help[2,1] <- -1*",name,"[2,1] " ) )
x<-rbind(x, paste0( " ",name,".inv[1:",dim,",1:",dim,"] <- (1/",name,".det.adj * ",name,".help[,]) " ) )
}
if (dim %in% 4) {
x<-rbind(x, paste0( " # http://stackoverflow.com/questions/1148309/inverting-a-4x4-matrix ") )
x<-rbind(x, paste0( " ",name,".inv.help[1,1] <- ",name,"[2,2]*",name,"[3,3]*",name,"[4,4]-",name,"[2,2]*",name,"[3,4]*",name,"[4,3]-",name,"[3,2]*",name,"[2,3]*",name,"[4,4]+",name,"[3,2]*",name,"[2,4]*",name,"[4,3]+",name,"[4,2]*",name,"[2,3]*",name,"[3,4]-",name,"[4,2]*",name,"[2,4]*",name,"[3,3] ") )
x<-rbind(x, paste0( " ",name,".inv.help[2,1] <- -",name,"[2,1]*",name,"[3,3]*",name,"[4,4]+",name,"[2,1]*",name,"[3,4]*",name,"[4,3]+",name,"[3,1]*",name,"[2,3]*",name,"[4,4]-",name,"[3,1]*",name,"[2,4]*",name,"[4,3]-",name,"[4,1]*",name,"[2,3]*",name,"[3,4]+",name,"[4,1]*",name,"[2,4]*",name,"[3,3] ") )
x<-rbind(x, paste0( " ",name,".inv.help[3,1] <- ",name,"[2,1]*",name,"[3,2]*",name,"[4,4]-",name,"[2,1]*",name,"[3,4]*",name,"[4,2]-",name,"[3,1]*",name,"[2,2]*",name,"[4,4]+",name,"[3,1]*",name,"[2,4]*",name,"[4,2]+",name,"[4,1]*",name,"[2,2]*",name,"[3,4]-",name,"[4,1]*",name,"[2,4]*",name,"[3,2] ") )
x<-rbind(x, paste0( " ",name,".inv.help[4,1] <- -",name,"[2,1]*",name,"[3,2]*",name,"[4,3]+",name,"[2,1]*",name,"[3,3]*",name,"[4,2]+",name,"[3,1]*",name,"[2,2]*",name,"[4,3]-",name,"[3,1]*",name,"[2,3]*",name,"[4,2]-",name,"[4,1]*",name,"[2,2]*",name,"[3,3]+",name,"[4,1]*",name,"[2,3]*",name,"[3,2] ") )
x<-rbind(x, paste0( " ",name,".inv.help[1,2] <- -",name,"[1,2]*",name,"[3,3]*",name,"[4,4]+",name,"[1,2]*",name,"[3,4]*",name,"[4,3]+",name,"[3,2]*",name,"[1,3]*",name,"[4,4]-",name,"[3,2]*",name,"[1,4]*",name,"[4,3]-",name,"[4,2]*",name,"[1,3]*",name,"[3,4]+",name,"[4,2]*",name,"[1,4]*",name,"[3,3] ") )
x<-rbind(x, paste0( " ",name,".inv.help[2,2] <- ",name,"[1,1]*",name,"[3,3]*",name,"[4,4]-",name,"[1,1]*",name,"[3,4]*",name,"[4,3]-",name,"[3,1]*",name,"[1,3]*",name,"[4,4]+",name,"[3,1]*",name,"[1,4]*",name,"[4,3]+",name,"[4,1]*",name,"[1,3]*",name,"[3,4]-",name,"[4,1]*",name,"[1,4]*",name,"[3,3] ") )
x<-rbind(x, paste0( " ",name,".inv.help[3,2] <- -",name,"[1,1]*",name,"[3,2]*",name,"[4,4]+",name,"[1,1]*",name,"[3,4]*",name,"[4,2]+",name,"[3,1]*",name,"[1,2]*",name,"[4,4]-",name,"[3,1]*",name,"[1,4]*",name,"[4,2]-",name,"[4,1]*",name,"[1,2]*",name,"[3,4]+",name,"[4,1]*",name,"[1,4]*",name,"[3,2] ") )
x<-rbind(x, paste0( " ",name,".inv.help[4,2] <- ",name,"[1,1]*",name,"[3,2]*",name,"[4,3]-",name,"[1,1]*",name,"[3,3]*",name,"[4,2]-",name,"[3,1]*",name,"[1,2]*",name,"[4,3]+",name,"[3,1]*",name,"[1,3]*",name,"[4,2]+",name,"[4,1]*",name,"[1,2]*",name,"[3,3]-",name,"[4,1]*",name,"[1,3]*",name,"[3,2] ") )
x<-rbind(x, paste0( " ",name,".inv.help[1,3] <- ",name,"[1,2]*",name,"[2,3]*",name,"[4,4]-",name,"[1,2]*",name,"[2,4]*",name,"[4,3]-",name,"[2,2]*",name,"[1,3]*",name,"[4,4]+",name,"[2,2]*",name,"[1,4]*",name,"[4,3]+",name,"[4,2]*",name,"[1,3]*",name,"[2,4]-",name,"[4,2]*",name,"[1,4]*",name,"[2,3] ") )
x<-rbind(x, paste0( " ",name,".inv.help[2,3] <- -",name,"[1,1]*",name,"[2,3]*",name,"[4,4]+",name,"[1,1]*",name,"[2,4]*",name,"[4,3]+",name,"[2,1]*",name,"[1,3]*",name,"[4,4]-",name,"[2,1]*",name,"[1,4]*",name,"[4,3]-",name,"[4,1]*",name,"[1,3]*",name,"[2,4]+",name,"[4,1]*",name,"[1,4]*",name,"[2,3] ") )
x<-rbind(x, paste0( " ",name,".inv.help[3,3] <- ",name,"[1,1]*",name,"[2,2]*",name,"[4,4]-",name,"[1,1]*",name,"[2,4]*",name,"[4,2]-",name,"[2,1]*",name,"[1,2]*",name,"[4,4]+",name,"[2,1]*",name,"[1,4]*",name,"[4,2]+",name,"[4,1]*",name,"[1,2]*",name,"[2,4]-",name,"[4,1]*",name,"[1,4]*",name,"[2,2] ") )
x<-rbind(x, paste0( " ",name,".inv.help[4,3] <- -",name,"[1,1]*",name,"[2,2]*",name,"[4,3]+",name,"[1,1]*",name,"[2,3]*",name,"[4,2]+",name,"[2,1]*",name,"[1,2]*",name,"[4,3]-",name,"[2,1]*",name,"[1,3]*",name,"[4,2]-",name,"[4,1]*",name,"[1,2]*",name,"[2,3]+",name,"[4,1]*",name,"[1,3]*",name,"[2,2] ") )
x<-rbind(x, paste0( " ",name,".inv.help[1,4] <- -",name,"[1,2]*",name,"[2,3]*",name,"[3,4]+",name,"[1,2]*",name,"[2,4]*",name,"[3,3]+",name,"[2,2]*",name,"[1,3]*",name,"[3,4]-",name,"[2,2]*",name,"[1,4]*",name,"[3,3]-",name,"[3,2]*",name,"[1,3]*",name,"[2,4]+",name,"[3,2]*",name,"[1,4]*",name,"[2,3] ") )
x<-rbind(x, paste0( " ",name,".inv.help[2,4] <- ",name,"[1,1]*",name,"[2,3]*",name,"[3,4]-",name,"[1,1]*",name,"[2,4]*",name,"[3,3]-",name,"[2,1]*",name,"[1,3]*",name,"[3,4]+",name,"[2,1]*",name,"[1,4]*",name,"[3,3]+",name,"[3,1]*",name,"[1,3]*",name,"[2,4]-",name,"[3,1]*",name,"[1,4]*",name,"[2,3] ") )
x<-rbind(x, paste0( " ",name,".inv.help[3,4] <- -",name,"[1,1]*",name,"[2,2]*",name,"[3,4]+",name,"[1,1]*",name,"[2,4]*",name,"[3,2]+",name,"[2,1]*",name,"[1,2]*",name,"[3,4]-",name,"[2,1]*",name,"[1,4]*",name,"[3,2]-",name,"[3,1]*",name,"[1,2]*",name,"[2,4]+",name,"[3,1]*",name,"[1,4]*",name,"[2,2] ") )
x<-rbind(x, paste0( " ",name,".inv.help[4,4] <- ",name,"[1,1]*",name,"[2,2]*",name,"[3,3]-",name,"[1,1]*",name,"[2,3]*",name,"[3,2]-",name,"[2,1]*",name,"[1,2]*",name,"[3,3]+",name,"[2,1]*",name,"[1,3]*",name,"[3,2]+",name,"[3,1]*",name,"[1,2]*",name,"[2,3]-",name,"[3,1]*",name,"[1,3]*",name,"[2,2] ") )
x<-rbind(x, paste0( " ",name,".det <- ",name,"[1,1]*",name,".inv.help[1,1]+",name,"[1,2]*",name,".inv.help[2,1]+",name,"[1,3]*",name,".inv.help[3,1]+",name,"[1,4]*",name,".inv.help[4,1] ") )
x<-rbind(x, paste0( " # ",name," is invertible if determinant is not 0 ") )
x<-rbind(x, paste0( " ",name,".invertible <- ifelse( ",name,".det==0, 0, 1 ) ") )
x<-rbind(x, paste0( " # if not invertible, modest mod of determinant ") )
x<-rbind(x, paste0( " ",name,".detmod1 ~ dbern(0.5) ") )
x<-rbind(x, paste0( " ",name,".detmod2 <- ifelse( ",name,".invertible==1, 0, (",name,".detmod1-0.5)/500 )") )
x<-rbind(x, paste0( " ",name,".det.adj <- ",name,".det + ",name,".detmod2 ") )
x<-rbind(x, paste0( " ",name,".inv <- ",name,".inv.help / ",name,".det.adj ") )
}
x<-rbind(x,"")
return( x )
}
make.inverse.2 <- function( name, dim ){
### *.prec.replace must exist (for dim>1)
x<-matrix( "", 1, 1 )
if (dim %in% 1) {
if (name %in% "Qt"){
x<-rbind(x, " # inverse of Qt ")
x<-rbind(x, " Qt.prec[1,1,t,p] <- (1/Qt[1,1,t,p]) ")
}
}
if (dim %in% 2) {
if (name %in% "Qt"){
x<-rbind(x, " ## inverse of Qt ")
x<-rbind(x, " # determinant of Qt ")
x<-rbind(x, " Qt.det[t,p] <- Qt[1,1,t,p]*Qt[2,2,t,p] - Qt[2,1,t,p]*Qt[1,2,t,p] ")
x<-rbind(x, " # Qt is invertible if determinant is not 0 ")
x<-rbind(x, " Qt.invertible[t,p] <- ifelse( Qt.det[t,p]==0, 0, 1 ) ")
x<-rbind(x, " # if not invertible, modest mod of determinant ")
x<-rbind(x, " Qt.detmod1[t,p] ~ dbern(0.5) ")
x<-rbind(x, " Qt.detmod2[t,p] <- ifelse( Qt.invertible[t,p]==1, 0, (Qt.detmod1[t,p]-0.5)/500 ) ")
x<-rbind(x, " Qt.det.adj[t,p] <- Qt.det[t,p] + Qt.detmod2[t,p] ")
x<-rbind(x, " # inverse of Qt ")
x<-rbind(x, " Qt.help[1,1,t,p] <- Qt[2,2,t,p] ")
x<-rbind(x, " Qt.help[2,2,t,p] <- Qt[1,1,t,p] ")
x<-rbind(x, " Qt.help[1,2,t,p] <- -1*Qt[1,2,t,p] ")
x<-rbind(x, " Qt.help[2,1,t,p] <- -1*Qt[2,1,t,p] ")
x<-rbind(x, " Qt.prec.orig[1:Aw,1:Aw,t,p] <- (1/Qt.det.adj[t,p] * Qt.help[,,t,p]) ")
x<-rbind(x, " ")
x<-rbind(x, " # check if Qt.prec.orig is positive definite ")
x<-rbind(x, " Qt.prec.orig.det[t,p] <- Qt.prec.orig[1,1,t,p]*Qt.prec.orig[2,2,t,p] - Qt.prec.orig[2,1,t,p]*Qt.prec.orig[1,2,t,p] ")
x<-rbind(x, " Qt.prec.orig.posdef[t,p] <- ifelse( Qt.prec.orig[1,1,t,p] > 0 && Qt.prec.orig.det[t,p] > 0, 1, 0 ) ")
x<-rbind(x, " ")
x<-rbind(x, " # if Qt.prec.orig is not pos def then replace by a pos def matrix ")
x<-rbind(x, " Qt.prec[1:2,1:2,t,p] <- ifelse( Qt.prec.orig.posdef[t,p]==1, Qt.prec.orig[,,t,p] , Qt.prec.replace[,] ) ")
}
}
x<-rbind(x,"")
return( x )
}
Try the ctglm package in your browser
Any scripts or data that you put into this service are public.
ctglm documentation built on May 31, 2017, 1:54 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions create.jags.syntax make.str moveTo.par.env make.inverse.1 make.inverse.2
create.jags.syntax <- function ( env ) {
# get variables from env
eval( parse( text=paste0( "assign( '",ls(envir=env), "' , get('",ls(envir=env),"', envir=env ) )" ) ) )
# objects with free parameters must be overwritten or deleted
# environment to save original objects
par.env <- new.env()
# syntax matrix (1 column)
x<-matrix( "model", 1, 1 )
x<-rbind(x, "{ ")
x<-rbind(x, " ")
x<-rbind(x, " ### DATA ### ")
x<-rbind(x, " ")
x<-rbind(x, " # data dl in long format ")
x<-rbind(x, " # dl[,col.id]: person ")
x<-rbind(x, " # dl[,col.item]: item ")
x<-rbind(x, " # dl[,col.time]: time point ")
x<-rbind(x, " # dl[,col.y]: responses ")
x<-rbind(x, " ")
x<-rbind(x, " ")
x<-rbind(x, " ### MEASUREMENT MODEL ### ")
x<-rbind(x, " ")
# x<-rbind(x, " # estimated parameters ")
# x<-rbind(x, " # theta: J (person) x F (faktor) x T (time point) ")
# x<-rbind(x, " ")
x<-rbind(x, " # loop over rows in long data set dl ")
x<-rbind(x, " for (r in 1:R) { ")
x<-rbind(x, " ")
x<-rbind(x, " # distributional assumption ")
if( measurement.model$family %in% "gaussian" ) {
x<-rbind(x, " dl[r,col.y] ~ dnorm( mu.y[r], E[ dl[r,col.item], dl[r,col.item] ] ) ") }
if( measurement.model$family %in% "binomial" ) {
x<-rbind(x, " dl[r,col.y] ~ dbern( mu.y[r] ) ") }
if( measurement.model$family %in% "poisson" ) {
x<-rbind(x, " dl[r,col.y] ~ dpois( mu.y[r] ) ") }
x<-rbind(x, " ")
x<-rbind(x, " # link ")
if( measurement.model$family %in% "gaussian" ) {
x<-rbind(x, " mu.y[r] <- eta[r] ") }
if( measurement.model$link %in% "logit" ) {
x<-rbind(x, " mu.y[r] <- ilogit( eta[r] ) ") }
if( measurement.model$link %in% "log" ) {
x<-rbind(x, " log( mu.y[r] ) <- eta[r] ") }
x<-rbind(x, " ")
x<-rbind(x, " # linear predictor ")
# x<-rbind(x, " eta[r] <- sum( Lambda[ dl[r,col.item], , dl[r,col.time] ] * theta[ dl[r,col.id], , dl[r,col.time] ] ) + beta[ dl[r,col.item], 1 ] ")
x<-rbind(x, " eta[r] <- sum( Lambda[ dl[r,col.item], ] * theta[ dl[r,col.id], , dl[r,col.time] ] ) + beta[ dl[r,col.item], 1 ] ")
# x<-rbind(x, " eta[r] <- sum( Lambda[ dl[r,col.item], ] * theta[ dl[r,col.id], , dl[r,col.time] ] ) ")
x<-rbind(x, " ")
x<-rbind(x, " } ")
x<-rbind(x, " ")
if ( exists( "E" ) ) {
x<-rbind(x, " # precision of measurement errors ")
x<-rbind(x, make.str( "E" ) ); if( any ( is.parameter( E ) ) ) invisible(moveTo.par.env("E",env,par.env)) else rm("E", envir=env)
x<-rbind(x, " ") }
# if( measurement.model$family %in% "gaussian" ) {
# x<-rbind(x, " # precision of measurement errors ")
# x<-rbind(x, " for (i in 1:I) { ")
# x<-rbind(x, " E[i] ~ dgamma(1,0.005) ")
# x<-rbind(x, " } ")
# x<-rbind(x, " ") }
x<-rbind(x, " # values/priors of item easiness ")
x<-rbind(x, make.str( "beta" ) ); if( any ( is.parameter( beta ) ) ) invisible(moveTo.par.env("beta",env,par.env)) else rm("beta", envir=env)
x<-rbind(x, " ")
if ( exists( "mu.beta" ) ) {
x<-rbind(x, " # values/priors of mean of item easiness ")
x<-rbind(x, make.str( "mu.beta" ) ); if( any ( is.parameter( mu.beta ) ) ) invisible(moveTo.par.env("mu.beta",env,par.env)) else rm("mu.beta", envir=env)
x<-rbind(x, " ") }
if ( exists( "mu.beta" ) ) {
x<-rbind(x, " # values/priors of precision of item easiness ")
x<-rbind(x, make.str( "prec.beta" ) ); if( any ( is.parameter( prec.beta ) ) ) invisible(moveTo.par.env("prec.beta",env,par.env)) else rm("prec.beta", envir=env)
x<-rbind(x, " ") }
x<-rbind(x, " # values/priors of Lambda ")
x<-rbind(x, make.str( "Lambda" ) ); if( any ( is.parameter( Lambda ) ) ) invisible(moveTo.par.env("Lambda",env,par.env)) else rm("Lambda", envir=env)
x<-rbind(x, " ")
# x<-rbind(x, " # standard deviation of latent scale ")
# x<-rbind(x, " #sd.eta <- sd(eta[]) ")
# x<-rbind(x, " ")
x<-rbind(x, " ")
x<-rbind(x, " ### CONTINUOUS TIME MODEL ### ")
x<-rbind(x, " ")
# browser()
if( !exists( "stationarity" ) ) stationarity <- FALSE
if( person.var["mu.t1"] ) {
if (! stationarity ){
x<-rbind(x, " ## personal t1 means ")
x<-rbind(x, " # mu.t1 varies over persons, i.e. for each person a separate mu.t1 ")
mu.t1.j. <- make.str( "mu.t1.j" ); if( any ( is.parameter( mu.t1.j ) ) ) invisible(moveTo.par.env("mu.t1.j",env,par.env)) else rm("mu.t1.j", envir=env)
# mods to mu.t1.j.
mu.t1.j.1 <- mu.t1.j.[ !grepl( "~", mu.t1.j., fixed=TRUE ),,drop=FALSE ]
mu.t1.j.2 <- mu.t1.j.[ grepl( "~", mu.t1.j., fixed=TRUE ),,drop=FALSE ]
if( F>1 ) mu.t1.j.2[,1] <- sub( "mu.t1.j\\[\\d+", "mu.t1.j[1:F", mu.t1.j.2[,1] )
if( F==1 ) mu.t1.j.2[,1] <- sub( "mu.t1.j\\[\\d+", "mu.t1.j[1", mu.t1.j.2[,1] )
mu.t1.j.2 <- mu.t1.j.2[!duplicated(mu.t1.j.2),,drop=FALSE]
if( !identical( mu.t1.j.1, character(0) ) ) x<-rbind(x,mu.t1.j.1)
if( !identical( mu.t1.j.2, character(0) ) ) x<-rbind(x,mu.t1.j.2)
x<-rbind(x, " ")
x<-rbind(x, " # values/prior for prec.mu.t1.j ")
x<-rbind(x, ifelse( exists("prec.mu.t1.j.prior") && is.null( dim(prec.mu.t1.j.prior) ), paste0( " prec.mu.t1.j[1:F,1:F] ~ ", prec.mu.t1.j.prior, " " ), make.str( "prec.mu.t1.j" ) ) ); if( any ( is.parameter( prec.mu.t1.j ) ) ) invisible(moveTo.par.env("prec.mu.t1.j",env,par.env)) else rm("prec.mu.t1.j", envir=env)
}
if (stationarity){
mu.t1.j. <- make.str( "mu.t1.j" )
mu.t1.j.1 <- mu.t1.j.[ !grepl( "~", mu.t1.j., fixed=TRUE ),,drop=FALSE ]
mu.t1.j.2 <- mu.t1.j.[ grepl( "~", mu.t1.j., fixed=TRUE ),,drop=FALSE ]
if( F>1 ) mu.t1.j.2[,1] <- sub( "mu.t1.j\\[\\d+", "mu.t1.j[1:F", mu.t1.j.2[,1] )
if( F==1 ) mu.t1.j.2[,1] <- sub( "mu.t1.j\\[\\d+", "mu.t1.j[1", mu.t1.j.2[,1] )
mu.t1.j.2 <- mu.t1.j.2[!duplicated(mu.t1.j.2),,drop=FALSE]
mu.t1.j.3 <- sub( "^(.*)~.*$", "\\1", mu.t1.j.2 )
mu.t1.j.4 <- matrix( paste0( mu.t1.j.3, " <- -1 * A.inv[,] %*% bj[,", as.integer( sub( "^.*,(.*)\\].*$", "\\1", mu.t1.j.3 ) ) ,"]" ), ncol=1 )
if( !identical( mu.t1.j.4, character(0) ) ) x<-rbind(x,mu.t1.j.4)
# push to parameter env, because we want results for that
if( any ( is.parameter( mu.t1.j ) ) ) invisible(moveTo.par.env("mu.t1.j",env,par.env)) else rm("mu.t1.j", envir=env)
}
}
x<-rbind(x, " ")
x<-rbind(x, " # loop over persons j ")
x<-rbind(x, " for (j in 1:J) { ")
x<-rbind(x, " ")
x<-rbind(x, " # first time point ")
if( !person.var["mu.t1"] ) {
x<-rbind(x, paste0( " theta[j,1:F,1] ~ d",ifelse(F>1,"m",""),"norm( mu.t1, prec.t1 ) ") ) }
if( person.var["mu.t1"] ) {
x<-rbind(x, paste0( " theta[j,1:F,1] ~ d",ifelse(F>1,"m",""),"norm( mu.t1.j[,j], prec.t1 ) ") ) }
x<-rbind(x, " ")
x<-rbind(x, " # loop over t=2,...,Tj personal time point ")
x<-rbind(x, " for (t in 2:Tj[j]) { ")
x<-rbind(x, " ")
# browser()
if( !person.var["b"] ) {
x<-rbind(x, paste0( " theta[j,1:F,t] ~ d",ifelse(F>1,"m",""),"norm( At[,,t-1,Lpat.group[j]] %*% theta[j,,t-1] + bt[,t-1,Lpat.group[j]], Qt.prec[,,t,Lpat.group[j]] ) ") ) }
if( person.var["b"] ) {
x<-rbind(x, paste0(" theta[j,1:F,t] ~ d",ifelse(F>1,"m",""),"norm( At[,,t-1,Lpat.group[j]] %*% theta[j,,t-1] + bt[,t-1,j], Qt.prec[,,t,Lpat.group[j]] ) ") ) }
x<-rbind(x, " ")
x<-rbind(x, " } ")
x<-rbind(x, " } ")
x<-rbind(x, " ")
x<-rbind(x, " ## values/prior of parameters ")
# x<-rbind(x, " ## Note: values are commented (already set in R) ")
x<-rbind(x, " ")
x<-rbind(x, " # values/prior of mean of first time point ")
if( !person.var["mu.t1"] || !stationarity ) {
x<-rbind(x, make.str( "mu.t1" ) ); if( any ( is.parameter( mu.t1 ) ) ) invisible(moveTo.par.env("mu.t1",env,par.env)) else rm("mu.t1", envir=env)
} else {
x<-rbind(x, " for (f in 1:F) {" )
x<-rbind(x, " mu.t1[f] <- mean( mu.t1.j[f,] )" )
x<-rbind(x, " }" )
if( any ( is.parameter( mu.t1 ) ) ) invisible(moveTo.par.env("mu.t1",env,par.env)) else rm("mu.t1", envir=env)
}
x<-rbind(x, " ")
# browser()
x<-rbind(x, " # values/prior of precision of first time point ")
# x<-rbind(x, ifelse( exists("prec.t1.prior") && is.null( dim(prec.t1.prior) ), paste0( " prec.t1[1:F,1:F] ~ ", prec.t1.prior, " " ), make.str( "prec.t1" ) ) ); if( any ( is.parameter( prec.t1 ) ) ) invisible(moveTo.par.env("prec.t1",env,par.env)) else rm("prec.t1", envir=env)
x<-rbind(x, make.str( "prec.t1" ) ); if( any ( is.parameter( prec.t1 ) ) ) invisible(moveTo.par.env("prec.t1",env,par.env)) else rm("prec.t1", envir=env)
# prec.t1 is symmetric, needs definition of upper to lower elements
# nasty workaround here
# prec.t1.prior <- prec.t1
# pstr <- make.str( "prec.t1" )
# pstr <- pstr[ grepl( "<-", pstr[,1] ), ,drop=FALSE]
# if( nrow( pstr ) > 0 ) {
# x<-rbind(x,pstr)
# }
x<-rbind(x, " ")
x<-rbind(x, " # values/prior of drift matrix ")
x<-rbind(x, make.str( "A" ) ); if( any ( is.parameter( A ) ) ) invisible(moveTo.par.env("A",env,par.env)) else rm("A", envir=env)
x<-rbind(x, " ")
x<-rbind(x, " # values/prior of diffusion matrix ")
# browser()
x<-rbind(x, make.str( "Q" ) ); if( any ( is.parameter( Q ) ) ) invisible(moveTo.par.env("Q",env,par.env)) else rm("Q", envir=env)
x<-rbind(x, " ")
# browser()
if( person.var["b"] ) {
x<-rbind(x, " ## continuous time intercepts ")
x<-rbind(x, " # b varies over persons, i.e. for each person a separate b ")
# x<-rbind(x, " for( j in 1:J ) { ")
# x<-rbind(x, " bj[1:F,j] ~ dmnorm( b[,1], prec.b[,] ) ")
# x<-rbind(x, " } ") }
bj. <- make.str( "bj" ); if( any ( is.parameter( bj ) ) ) invisible(moveTo.par.env("bj",env,par.env)) else rm("bj", envir=env)
# mods to bj.
bj.1 <- bj.[ !grepl( "~", bj., fixed=TRUE ),,drop=FALSE ]
bj.2 <- bj.[ grepl( "~", bj., fixed=TRUE ),,drop=FALSE ]
if( F>1 ) bj.2[,1] <- sub( "bj\\[\\d+", "bj[1:F", bj.2[,1] )
if( F==1 ) bj.2[,1] <- sub( "bj\\[\\d+", "bj[1", bj.2[,1] )
bj.2 <- bj.2[!duplicated(bj.2),,drop=FALSE]
if( !identical( bj.1, character(0) ) ) x<-rbind(x,bj.1)
if( !identical( bj.2, character(0) ) ) x<-rbind(x,bj.2) }
x<-rbind(x, " # values/prior of continuous time intercepts ")
x<-rbind(x, make.str( "b" ) ); if( any ( is.parameter( b ) ) ) invisible(moveTo.par.env("b",env,par.env)) else rm("b", envir=env)
if( person.var["b"] ) {
x<-rbind(x, " # values/prior of precision of b ")
# x<-rbind(x, ifelse( exists("prec.b.prior") && is.null( dim(prec.b.prior) ), paste0( " prec.b[1:F,1:F] ~ ", prec.b.prior, " " ), make.str( "prec.b" ) ) ); if( any ( is.parameter( prec.b ) ) ) invisible(moveTo.par.env("prec.b",env,par.env)) else rm("prec.b", envir=env) }
# browser()
x<-rbind(x, make.str( "prec.b" ) ); if( any ( is.parameter( prec.b ) ) ) invisible(moveTo.par.env("prec.b",env,par.env)) else rm("prec.b", envir=env) }
x<-rbind(x, " ")
x<-rbind(x, " ")
x<-rbind(x, " # rowQ: stack Q row-wise into a column vector ")
x<-rbind(x, " for ( z in (1:Aw) ) { ")
x<-rbind(x, " for ( c in (1:Aw) ) { ")
x<-rbind(x, " rowQ[ c+Aw*(z-1), 1 ] <- Q[z,c] ")
x<-rbind(x, " } ")
x<-rbind(x, " } ")
x<-rbind(x, " ")
x<-rbind(x, " ## discrete/continuous time relations ")
x<-rbind(x, " # Qtv for first time point ")
x<-rbind(x, " # for technical reasons, same values are set for all lag patterns (not necessary) ")
x<-rbind(x, " for(p in 1:P) { ")
x<-rbind(x, " Qtv[1:(F*F),1,p] <- ( -1 * Ah.inv[,] ) %*% rowQ[,1] ")
x<-rbind(x, " } ")
# mexp or exp
# browser()
if( F>1 ) fexp <- "mexp" else if ( F== 1 ) fexp <- "exp"
x<-rbind(x, " # loop over p=1,...,P lag patterns ")
x<-rbind(x, " for(p in 1:P) { ")
x<-rbind(x, " # loop over t=2,...,Tp pattern-specific time point ")
x<-rbind(x, " for (t in 2:Tp[p]) { ")
x<-rbind(x, " ")
x<-rbind(x, " # autoregression matrix / drift matrix ")
x<-rbind(x, paste0( " At[1:F,1:F,t-1,p] <- ",fexp,"( A[,] * Lpat[p,t-1] ) ") )
x<-rbind(x, " ")
x<-rbind(x, " # continuous time intercepts ")
if( !person.var["b"] ) {
x<-rbind(x, paste0( " bt[1:F,t-1,p] <- ( A.inv[,] %*% ( ",fexp,"( A[,] * Lpat[p,t-1] ) - I1 ) ) %*% b[,1] ") ) }
if( person.var["b"] ) {
x<-rbind(x, " # b varies over persons, so loop over persons in specific pattern")
x<-rbind(x, " # PNj ... number of persons in pattern ")
x<-rbind(x, " # Pj (P x J) ... ragged matrix, indicators of persons in pattern")
x<-rbind(x, " for (j in 1:PNj[p]) { ")
x<-rbind(x, paste0( " bt[1:F,t-1,Pj[p,j]] <- ( A.inv[,] %*% ( ",fexp,"( A[,] * Lpat[p,t-1] ) - I1 ) ) %*% bj[,j]") )
x<-rbind(x, " } ") }
x<-rbind(x, " ")
x<-rbind(x, " # Qtv is vectorized Qt matrix ")
x<-rbind(x, paste0( " #Qtv[1:(F*F),t,p] <- ifelse( t==1, ( -1 * Ah.inv[,] ) %*% rowQ[,1] , ( Ah.inv[,] %*% ( ",fexp,"( Ah[,] * Lpat[p,t-1] ) - I2 ) ) %*% rowQ[,1] ) ") )
x<-rbind(x, paste0( " Qtv[1:(F*F),t,p] <- Ah.inv[,] %*% ( ",fexp,"( Ah[,] * Lpat[p,t-1] ) - I2 ) %*% rowQ[,1] ") )
x<-rbind(x, " } ")
x<-rbind(x, " } ")
x<-rbind(x, " ")
x<-rbind(x, " ## calculate Qt.prec ")
x<-rbind(x, " for(p in 1:P) { ")
x<-rbind(x, " for (t in 1:Tp[p]) { ")
x<-rbind(x, " ")
x<-rbind(x, " # Qt matrix (create with 'irow') ")
x<-rbind(x, " for ( z in (1:Aw) ) { ")
x<-rbind(x, " for ( c in (1:Aw) ) { ")
x<-rbind(x, " Qt[z,c,t,p] <- Qtv[ c+Aw*(z-1),t,p] ")
x<-rbind(x, " } ")
x<-rbind(x, " } ")
x<-rbind(x, " ")
# inverse of Qt; Qt has dimensions of Q; Qt.replace must exist
x<-rbind(x, make.inverse.2( "Qt", dim(Q)[1] ) )
x<-rbind(x, " } ")
x<-rbind(x, " } ")
x<-rbind(x, " ")
# inverse of A and Ah
x<-rbind(x, make.inverse.1( "A", dim(A)[1] ) )
x<-rbind(x, make.inverse.1( "Ah", dim(A)[1]^2 ) )
x<-rbind(x, " ")
x<-rbind(x, " # Kronecker of A and I1 ")
x<-rbind(x, " # https://sourceforge.net/p/mcmc-jags/discussion/610037/thread/17ddb2ce/ ")
x<-rbind(x, " for (i in 1:Aw){ # Aw = width of A ")
x<-rbind(x, " for (j in 1:Aw){ ")
x<-rbind(x, " for (k in 1:I1w){ # I1w = width of I1 ")
x<-rbind(x, " for (l in 1:I1w){ ")
x<-rbind(x, " Ah1[(i-1)*I1w+k,(j-1)*I1w+l] <- A[i,j]*I1[k,l] ")
x<-rbind(x, " } ")
x<-rbind(x, " } ")
x<-rbind(x, " } ")
x<-rbind(x, " } ")
x<-rbind(x, " # Kronecker of I1 and A ")
x<-rbind(x, " for (i in 1:I1w){ # I1w = width of I1 ")
x<-rbind(x, " for (j in 1:I1w){ ")
x<-rbind(x, " for (k in 1:Aw){ # Aw = width of A ")
x<-rbind(x, " for (l in 1:Aw){ ")
x<-rbind(x, " Ah2[(i-1)*Aw+k,(j-1)*Aw+l] <- I1[i,j]*A[k,l] ")
x<-rbind(x, " } ")
x<-rbind(x, " } ")
x<-rbind(x, " } ")
x<-rbind(x, " } ")
x<-rbind(x, " # Ah matrix ")
x<-rbind(x, " Ah[1:(I1w*Aw),1:(I1w*Aw)] <- Ah1[,] + Ah2[,] ")
x<-rbind(x, " ")
# browser()
# additional jags syntax
if( exists( "jags.add" ) && !is.null( jags.add ) ) {
if ( !is.matrix( jags.add ) ) jags.add <- matrix( jags.add, ncol=1 )
jags.add <- rbind( "## Additional/user code", jags.add )
jags.add <- rbind( jags.add, "" )
x<-rbind(x, jags.add)
}
x<-rbind(x, "} ")
# operating system
os <- Sys.info()["sysname"]
# parallelization yes/no
parall <- exists("cores") && !is.null(cores) && cores > 1 && os %in% c("Windows","Linux")
### call matrix (1 column)
y<-matrix( paste0( "### R syntax for ", model.name ), 1, 1 )
y<-rbind(y,paste0( "### engine: ", engine ) )
y<-rbind(y, "" )
y<-rbind(y, "# packages" )
y<-rbind(y, "require( 'rjags' )" )
if( parall ){
y<-rbind(y, "require( 'doParallel' ) # function foreach %dopar% ")
y<-rbind(y, "require( 'abind' ) # required by ctglm ")
y<-rbind(y, "#require( 'ctglm' ) # function mcarray.chains.combine ") }
y<-rbind(y, "" )
y<-rbind(y, "# package versions" )
y<-rbind(y, paste0( "print( installed.packages()[ installed.packages()[,1] %in% c(",paste( paste0( "'", na.omit( c('rjags' , ifelse( parall, 'ctglm', NA) ) ), "'" ) , collapse="," ),"), c(1,3) ] )" ) )
y<-rbind(y, "" )
if( !parall ){
y<-rbind(y, "# JAGS Modules" )
if( F>1 ) y<-rbind(y, "load.module('msm') # for matrix exponential, mexp()" )
y<-rbind(y, "load.module('glm') # for better glm sampler " )
y<-rbind(y, "" ) }
# create starting values
# function make.priors with mode="startingvalue" is used
if ( exists("A") && any ( is.parameter( A ) ) ) invisible( make.priors( "A", A, diag.prior = -0.5, offdiag.prior = 0, env=environment(), mode="startingvalue", verbose=FALSE ) )
# browser()
if ( exists("Q") && any ( is.parameter( Q ) ) ) invisible( make.priors( "Q", Q, diag.prior = 0.5, offdiag.prior = 0, env=environment(), mode="startingvalue", verbose=FALSE ) )
if ( exists("b") && any ( is.parameter( b ) ) ) invisible( make.priors( "b", b, prior = 0, env=environment(), mode="startingvalue", verbose=FALSE ) )
if ( exists("beta") && any ( is.parameter( beta ) ) ) invisible( make.priors( "beta", beta, prior = 0, env=environment(), mode="startingvalue", verbose=FALSE ) )
if ( exists("mu.beta") && any ( is.parameter( mu.beta ) ) ) invisible( make.priors( "mu.beta", mu.beta, prior = 0, env=environment(), mode="startingvalue", verbose=FALSE ) )
if ( exists("prec.beta") && any ( is.parameter( prec.beta ) ) ) invisible( make.priors( "prec.beta", prec.beta, diag.prior = 1, offdiag.prior = 0, env=environment(), mode="startingvalue", verbose=FALSE ) )
if ( !stationarity && exists("mu.t1") && any ( is.parameter( mu.t1 ) ) ) invisible( make.priors( "mu.t1", mu.t1, prior = 0, env=environment(), mode="startingvalue", verbose=FALSE ) )
if ( exists("prec.t1") && any ( is.parameter( prec.t1 ) ) ) {
invisible( make.priors( "prec.t1", prec.t1, diag.prior = 1, offdiag.prior = 0,env=environment(), mode="startingvalue", verbose=FALSE ) )
# !!!! because of dwish, complete starting value is needed, so assign upper with values of lower triangle
prec.t1.startingvalue[upper.tri(prec.t1.startingvalue)] <- prec.t1.startingvalue[lower.tri(prec.t1.startingvalue)]
}
if ( exists("prec.mu.t1.j") && any ( is.parameter( prec.mu.t1.j ) ) ) {
invisible( make.priors( "prec.mu.t1.j", prec.mu.t1.j, diag.prior = 1, offdiag.prior = 0,env=environment(), mode="startingvalue", verbose=FALSE ) )
# !!!! because of dwish, complete starting value is needed, so assign upper with values of lower triangle
prec.mu.t1.j.startingvalue[upper.tri(prec.mu.t1.j.startingvalue)] <- prec.mu.t1.j.startingvalue[lower.tri(prec.mu.t1.j.startingvalue)]
}
# list with starting values
startingvalues <- list()
do <- do.call( "c", sapply( c("A","Q","b","beta","mu.beta","prec.beta","mu.t1","prec.t1"), function(x) if (exists(paste0(x,".startingvalue"))) paste0( "startingvalues$", x ," <- ", x,".startingvalue" ) else NULL, simplify=FALSE ) )
if ( length( do > 0 ) ) eval( parse( text=do ) )
# push starting values to data environment
# and create code to create as many starting value lists as chains
# (chains are set in next step, run model)
if ( length(startingvalues) > 0 ){
assign( "startingvalues", startingvalues, envir=env )
y<-rbind(y, "# starting values: create as many starting value lists as chains")
y<-rbind(y, "startingvalues <- get( 'startingvalues', envir=data.env )")
y<-rbind(y, "inits <- sapply( 1:chains, function (n,l) l, startingvalues, simplify=FALSE )")
y<-rbind(y, "" )
}
# seeds
y<-rbind(y, "# seeds for each chain")
y<-rbind(y, "seeds <- parallel.seeds('base::BaseRNG', chains)")
if ( length(startingvalues) > 0 ){
y<-rbind(y, "inits <- mapply ( function ( l, k ) c(l,k), inits, seeds, SIMPLIFY=FALSE )")
} else {
y<-rbind(y, "inits <- seeds")
}
y<-rbind(y, "" )
y<-rbind(y, "# start time ")
y<-rbind(y, "start <- Sys.time() ")
y<-rbind(y, "" )
# browser()
# parallelization
indent <- ifelse( parall, " ", "" )
if( parall ){
y<-rbind(y, "## parallelization ")
y<-rbind(y, "# create cluster ")
if( os %in% "Windows" ){
y<-rbind(y, paste0( "cl <- makePSOCKcluster(",cores,") ")) }
if( os %in% "Linux" ){
y<-rbind(y, paste0( "cl <- makeForkCluster(",cores,") ")) }
y<-rbind(y, "# register cluster ")
y<-rbind(y, "registerDoParallel(cl) ")
y<-rbind(y, "" )
y<-rbind(y, "# parallel chains ")
y<-rbind(y, "res.l <- foreach(chain=1:chains, .packages='rjags') %dopar% { ")
y<-rbind(y, "" )
y<-rbind(y, paste0( indent, "# JAGS Modules (need to be loaded here in parallel mode)" ) )
if( F>1 ) y<-rbind(y, paste0( indent, "load.module('msm') # for matrix exponential, mexp()" ) )
y<-rbind(y, paste0( indent, "load.module('glm') # for better glm sampler " ) )
y<-rbind(y, "" )
}
y<-rbind(y, paste0( indent, "# initialization/adaptation ") )
# globalenv !!!
# y<-rbind(y, "eval(parse(text=paste0( name, ".ini <- jags.model ( file = mf , data=globalenv(), inits=sL, n.chains = ",chains,", n.adapt=",adapt,", quiet=FALSE )" )),envir=globalenv() )" )
# y<-rbind(y, "eval(parse(text=paste0( name, ".ini <- jags.model ( file = mf , data=globalenv(), n.chains = ",chains,", n.adapt=",adapt,", quiet=FALSE )" )),envir=globalenv() )" )
# y<-rbind(y, "eval(parse(text=paste0( model.name, .ini <- jags.model ( file = bugs.file , data=data.env, n.chains=chains, n.adapt=adapt, quiet=FALSE ) )),envir=data.env )" )
# y<-rbind(y, "eval(parse(text=paste0( model.name, .ini <- jags.model ( file = bugs.file , data=data.env, n.chains=chains, n.adapt=adapt, quiet=FALSE ) )),envir=globalenv())" )
y<-rbind(y, paste0( indent, "ini <- jags.model ( file = bugs.file , data=data.env, n.chains=",ifelse(parall,"1","chains"),", n.adapt=adapt, inits=inits",ifelse(parall,"[chain]",""),", quiet=FALSE )" ) )
y<-rbind(y, "" )
y<-rbind(y, paste0( indent, "# run ") )
# y<-rbind(y, "eval(parse(text=paste0( model.name, .res <- jags.samples ( ',model.name,'.ini , variable.names=c('A','Q','b'), n.iter=iter, thin=thin, type='trace' , progress.bar = 'none', by=20 )' )), envir=globalenv() )" )
# browser()
## add or delete person par
if( exists( "track.person.par" ) && "theta" %in% track.person.par ) {
# theta is not regularly created (makes no sense, since always estimated, never set)
# so create a structure for theta for further processing
thetas.dfr <- Reduce(function(x, y) merge(x, y, all=TRUE),list(data.frame("J"=1:J),data.frame("F"=1:F),data.frame("T"=1:T)),accumulate=FALSE )
thetas <- apply( thetas.dfr, 1, function(z) paste0( "theta_", paste( z, collapse="_" ) ) )
theta. <- array( thetas, dim=c(J,F,T) )
assign( "theta", theta., envir=par.env )
}
if( !exists( "track.person.par" ) || !"bj" %in% track.person.par ) {
# bj is regularly created, if it should not be tracked, delete it from par.env
if ( exists( "bj", envir=par.env ) ) rm( "bj", envir=par.env )
}
if( !exists( "track.person.par" ) || !"mu.t1.j" %in% track.person.par ) {
if ( exists( "mu.t1.j", envir=par.env ) ) rm( "mu.t1.j", envir=par.env )
}
# mod7: delete code in parameters
invisible( make.matrices.consistent( par.env , mods=7 ) )
# all parameters from par.env must be tracked
pars <- ls(envir=par.env)
# browser()
# mix in track.include / track.exclude
if ( exists( "track.include" ) && is.list( track.include ) && length( track.include ) > 0 ) {
# include in tracking
pars <- unique( c( pars, names( track.include ) ) )
# for automatic results retrieval structure of objects need to be known
# if not user set (i.e. NULL) it is assumed as a scalar
# in this case there is one parameter of the same name
if( any( isnull <- sapply( track.include, is.null ) ) ) track.include[ isnull ] <- names( track.include[ isnull ] )
# put into par.env
eval( parse( text=paste0( "assign( '", names(track.include), "' , track.include$'", names(track.include), "' , envir=par.env )" ) ) )
}
if ( exists( "track.exclude" ) && is.character( track.exclude ) && length( track.exclude ) > 0 ) {
# exclude from tracking
if ( any ( del <- pars %in% track.exclude ) ) pars <- pars[ !del ]
# also delete from par.env (ctglm.results depends on clean par.env)
if ( any ( del <- ls(envir=par.env) %in% track.exclude ) ) rm( list=ls(envir=par.env)[del], envir=par.env )
}
# sortieren
### ist glaub ich ziemlich sinnlos, da jags eh nicht die Reihenfolge hat dann
ord <- match( c("A","Q","b"), pars )
ord <- ord[!is.na(ord)]
if( !identical( ord, integer(0) ) ) pars <- c( pars[ord], pars[ !pars %in% pars[ord] ] )
par.string <- paste0( "c(", paste( paste0("'", pars, "'"), collapse="," ) ,")" )
y<-rbind(y, paste0( indent, "res <- jags.samples ( ini, variable.names=",par.string,", n.iter=iter, thin=thin, type='trace' , progress.bar = 'text', by=20 )" ) )
# parallelization, get results
if( parall ){
y<-rbind(y, "" )
y<-rbind(y, "}" )
y<-rbind(y, "" )
y<-rbind(y, "# stop cluster" )
y<-rbind(y, "stopCluster(cl)" )
y<-rbind(y, "" )
y<-rbind(y, "# combine results from parallel computing" )
y<-rbind(y, "res <- mcarray.chains.combine( res.l )" ) }
y<-rbind(y, "" )
y<-rbind(y, "# run time ")
y<-rbind(y, "runtime <- Sys.time() - start ")
### (over)write relevant variables to environment ###
# jags.syntax <- x
# obj <- c( "jags.syntax" )
# eval( parse ( text=paste0( "assign( '",obj, "' , get('",obj,"') , envir=env )" ) ) )
## create return object
ret <- list()
# first entry: engine
ret$engine <- engine
# second entry: model.name
ret$model.name <- model.name
# third entry: data environment
ret$data.env <- env
# fourth entry: environment with original matrices
ret$par.env <- par.env
# fifth entry: call
ret$call <- y
# sixth entry: syntax
ret$syntax <- x
# return
return( ret )
}
make.str <- function( y.name ) {
# browser()
# environment where objects are
env <- parent.frame()
# dimensions of structure, get from parent env
dim.y <- eval( parse( text=paste0( "dim(", y.name, ")" ) ), envir=env )
# if structure without dimensions (e.g. vector)
if( is.null( dim.y ) ) {
# long structure
y. <- data.frame( matrix( seq( along= eval( parse( text=paste0( "", y.name, "" ) ), envir=env ) ), ncol=1 ) )
} else {
# if structure with dimensions (e.g. matrix)
# long structure
y. <- eval(parse(text=paste0( "Reduce(function(x, y) merge(x, y, by=NULL, all=TRUE),list(", paste( paste0( "1:", dim.y ), collapse="," ), "),accumulate=FALSE )" )))
}
y.$par <- apply( y., 1, function ( z ) eval( parse( text= paste0( y.name,"[", paste(z,collapse=","), "]" ) ), envir=env ) )
# browser()
# sort out code, keep for later adding
if ( any ( is.code( y.$par ) ) ){
code <- sub( "^[^;]*;(.*)$", "\\1", y.$par[ is.code( y.$par ) ] )
add.code <- matrix( paste0( " ", gsub( "^\\s*", "", do.call( "c", strsplit(code,";") ) ) ), ncol=1 )
y.$par[ is.code( y.$par ) ] <- sub( "^([^;]*);.*$", "\\1", y.$par[ is.code( y.$par ) ] )
} else {
add.code <- NULL
}
# browser()
# tag duplicated free parameters and generate values to set duplicated
y.$dupl <- as.integer( duplicated(y.$par) & is.na(suppressWarnings(as.numeric(y.$par))) )
dupl <- unique( y.$par[ duplicated( y.$par ) & is.na(suppressWarnings(as.numeric(y.$par))) ] )
if ( length( dupl ) > 0 ) {
dl.l <- sapply( dupl, function( dl ){
y.1 <- y.[ y.$par %in% dupl, ]
y.1 <- y.1[ !duplicated( y.1$par ), ]
y.1$val <- apply( y.1, 1, function ( z ) paste0( y.name, "", "[", paste( z[-c(length(z)-1,length(z))], collapse=","), "]" ) )
return( y.1 )
}, simplify=FALSE )
dl <- do.call( "rbind", dl.l )
y.cn <- colnames( y. )
y.$order <- seq(along=rownames(y.))
y. <- merge( y., dl[,c("par","val")], by="par", sort=FALSE, all.x=TRUE )
y. <- y.[ order( y.$order), ]
y.$order <- NULL
y. <- y.[ , c( y.cn, "val" ), drop=FALSE ]
} else {
y.$val <- ""
}
# if( any( grepl( "Q", y.$par ) ) ) browser()
# loop over long structure, generate strings
make.str2 <- function( z ) {
# browser()
# if( grepl( "Q", z["par"] ) ) browser()
dupl <- z[ length(z)-1 ]
val <- z[ length(z) ]
par <- z["par"]
z <- z[ -c(length(z)-2,length(z)-1,length(z)) ]
if ( dupl %in% "0" ){
# if( grepl( "Q", par ) ) browser()
# name of parameter
nam <- paste0( y.name, "[", paste( z, collapse=","), "]" )
# value of parameter
par.val <- eval( parse( text=nam ), envir=env )
# value of prior
p.name <- paste0( y.name, ".prior" )
if ( exists( p.name, envir=env) ){
p.val <- eval( parse( text=paste0( p.name, "[", paste( z, collapse=","), "]" ) ), envir=env )
} else {
p.val <- NA
}
is.prior <- !is.na( p.val )
# determine if fixed value or freely estimated
# as.num <- suppressWarnings( as.numeric( eval( parse( text=nam ), envir=env ) ) )
# is.fixed <- !is.na( as.num )
# is fixed value or a parameter
# browser()
# operator: <- (for fixed oder parameter) or ~ (for prior)
# op <- ifelse( is.parameter( par.val ) | is.fixed( par.val ), " <- ", " ~ " )
op <- ifelse( is.prior, " ~ ", " <- " )
# if fixed, uncomment in jags syntax, because already set in R
# com <- ifelse( is.fixed, "# ", "" )
# com <- ""
# browser()
# value to be set
set.val <- "TODO"
# parameter
if( is.parameter( par.val ) ) set.val <- par
# fixed value (numeric)
if( !is.na( suppressWarnings( as.numeric( par ) ) ) ) set.val <- as.numeric( par )
# prior
if( is.prior ) set.val <- p.val
# val <- ifelse( is.fixed, as.character( as.num ), eval( parse( text=paste0( y.name, ".prior", "[", paste( z, collapse=","), "]" ) ), envir=env ) )
# val <- ifelse( is.fixed, "", eval( parse( text=paste0( y.name, ".prior", "[", paste( z, collapse=","), "]" ) ), envir=env ) )
# if( is.parameter( par.val ) ) set.val <- par
# if prior exists dann prior setzen
# get("Q",envir=env)
# get("Q.prior",envir=env)
# string
# s <- paste0( " ", com, nam, op, set.val )
s <- paste0( " ", nam, op, set.val )
} else {
s <- paste0( " ",y.name, "", "[", paste( z, collapse=","), "] <- ", val )
}
# return
return( s )
}
# browser()
# s <- apply( y.1[,-ncol(y.1),drop=FALSE], 1, make.str2 )
s <- apply( y., 1, make.str2 )
# append some white space (for optical reasons)
s <- sapply( s, function( x ) paste0( x, paste( rep( " ", 69 - nchar(x) ), collapse="" ) ) )
# as column vector
s <- matrix( s, ncol=1 )
# add.code
if ( !is.null( add.code ) ){
s <- rbind( s, add.code )
}
# return
return( s )
}
moveTo.par.env <- function( name, env, par.env ){
if( exists(name,envir=env) ){
assign( name, get( name, envir=env ), envir=par.env )
rm( list=name, envir=env )
}
TRUE
}
make.inverse.1 <- function( name, dim ){
x<-matrix( "", 1, 1 )
x<-rbind(x, paste0( " ## inverse of ",name," " ) )
if (dim %in% 1) {
x<-rbind(x, paste0( " ",name,".inv[1,1] <- (1/",name,"[1,1])" ) )
}
if (dim %in% 2) {
x<-rbind(x, paste0( " # determinant of ",name," " ) )
x<-rbind(x, paste0( " ",name,".det <- ",name,"[1,1]*",name,"[2,2] - ",name,"[2,1]*",name,"[1,2] " ) )
x<-rbind(x, paste0( " # ",name," is invertible if determinant is not 0 " ) )
x<-rbind(x, paste0( " ",name,".invertible <- ifelse( ",name,".det==0, 0, 1 ) " ) )
x<-rbind(x, paste0( " # if not invertible, modest mod of determinant " ) )
x<-rbind(x, paste0( " ",name,".detmod1 ~ dbern(0.5) " ) )
x<-rbind(x, paste0( " ",name,".detmod2 <- ifelse( ",name,".invertible==1, 0, (",name,".detmod1-0.5)/500 ) " ) )
x<-rbind(x, paste0( " ",name,".det.adj <- ",name,".det + ",name,".detmod2 " ) )
x<-rbind(x, paste0( " # inverse of ",name," " ) )
x<-rbind(x, paste0( " ",name,".help[1,1] <- ",name,"[2,2] " ) )
x<-rbind(x, paste0( " ",name,".help[2,2] <- ",name,"[1,1] " ) )
x<-rbind(x, paste0( " ",name,".help[1,2] <- -1*",name,"[1,2] " ) )
x<-rbind(x, paste0( " ",name,".help[2,1] <- -1*",name,"[2,1] " ) )
x<-rbind(x, paste0( " ",name,".inv[1:",dim,",1:",dim,"] <- (1/",name,".det.adj * ",name,".help[,]) " ) )
}
if (dim %in% 4) {
x<-rbind(x, paste0( " # http://stackoverflow.com/questions/1148309/inverting-a-4x4-matrix ") )
x<-rbind(x, paste0( " ",name,".inv.help[1,1] <- ",name,"[2,2]*",name,"[3,3]*",name,"[4,4]-",name,"[2,2]*",name,"[3,4]*",name,"[4,3]-",name,"[3,2]*",name,"[2,3]*",name,"[4,4]+",name,"[3,2]*",name,"[2,4]*",name,"[4,3]+",name,"[4,2]*",name,"[2,3]*",name,"[3,4]-",name,"[4,2]*",name,"[2,4]*",name,"[3,3] ") )
x<-rbind(x, paste0( " ",name,".inv.help[2,1] <- -",name,"[2,1]*",name,"[3,3]*",name,"[4,4]+",name,"[2,1]*",name,"[3,4]*",name,"[4,3]+",name,"[3,1]*",name,"[2,3]*",name,"[4,4]-",name,"[3,1]*",name,"[2,4]*",name,"[4,3]-",name,"[4,1]*",name,"[2,3]*",name,"[3,4]+",name,"[4,1]*",name,"[2,4]*",name,"[3,3] ") )
x<-rbind(x, paste0( " ",name,".inv.help[3,1] <- ",name,"[2,1]*",name,"[3,2]*",name,"[4,4]-",name,"[2,1]*",name,"[3,4]*",name,"[4,2]-",name,"[3,1]*",name,"[2,2]*",name,"[4,4]+",name,"[3,1]*",name,"[2,4]*",name,"[4,2]+",name,"[4,1]*",name,"[2,2]*",name,"[3,4]-",name,"[4,1]*",name,"[2,4]*",name,"[3,2] ") )
x<-rbind(x, paste0( " ",name,".inv.help[4,1] <- -",name,"[2,1]*",name,"[3,2]*",name,"[4,3]+",name,"[2,1]*",name,"[3,3]*",name,"[4,2]+",name,"[3,1]*",name,"[2,2]*",name,"[4,3]-",name,"[3,1]*",name,"[2,3]*",name,"[4,2]-",name,"[4,1]*",name,"[2,2]*",name,"[3,3]+",name,"[4,1]*",name,"[2,3]*",name,"[3,2] ") )
x<-rbind(x, paste0( " ",name,".inv.help[1,2] <- -",name,"[1,2]*",name,"[3,3]*",name,"[4,4]+",name,"[1,2]*",name,"[3,4]*",name,"[4,3]+",name,"[3,2]*",name,"[1,3]*",name,"[4,4]-",name,"[3,2]*",name,"[1,4]*",name,"[4,3]-",name,"[4,2]*",name,"[1,3]*",name,"[3,4]+",name,"[4,2]*",name,"[1,4]*",name,"[3,3] ") )
x<-rbind(x, paste0( " ",name,".inv.help[2,2] <- ",name,"[1,1]*",name,"[3,3]*",name,"[4,4]-",name,"[1,1]*",name,"[3,4]*",name,"[4,3]-",name,"[3,1]*",name,"[1,3]*",name,"[4,4]+",name,"[3,1]*",name,"[1,4]*",name,"[4,3]+",name,"[4,1]*",name,"[1,3]*",name,"[3,4]-",name,"[4,1]*",name,"[1,4]*",name,"[3,3] ") )
x<-rbind(x, paste0( " ",name,".inv.help[3,2] <- -",name,"[1,1]*",name,"[3,2]*",name,"[4,4]+",name,"[1,1]*",name,"[3,4]*",name,"[4,2]+",name,"[3,1]*",name,"[1,2]*",name,"[4,4]-",name,"[3,1]*",name,"[1,4]*",name,"[4,2]-",name,"[4,1]*",name,"[1,2]*",name,"[3,4]+",name,"[4,1]*",name,"[1,4]*",name,"[3,2] ") )
x<-rbind(x, paste0( " ",name,".inv.help[4,2] <- ",name,"[1,1]*",name,"[3,2]*",name,"[4,3]-",name,"[1,1]*",name,"[3,3]*",name,"[4,2]-",name,"[3,1]*",name,"[1,2]*",name,"[4,3]+",name,"[3,1]*",name,"[1,3]*",name,"[4,2]+",name,"[4,1]*",name,"[1,2]*",name,"[3,3]-",name,"[4,1]*",name,"[1,3]*",name,"[3,2] ") )
x<-rbind(x, paste0( " ",name,".inv.help[1,3] <- ",name,"[1,2]*",name,"[2,3]*",name,"[4,4]-",name,"[1,2]*",name,"[2,4]*",name,"[4,3]-",name,"[2,2]*",name,"[1,3]*",name,"[4,4]+",name,"[2,2]*",name,"[1,4]*",name,"[4,3]+",name,"[4,2]*",name,"[1,3]*",name,"[2,4]-",name,"[4,2]*",name,"[1,4]*",name,"[2,3] ") )
x<-rbind(x, paste0( " ",name,".inv.help[2,3] <- -",name,"[1,1]*",name,"[2,3]*",name,"[4,4]+",name,"[1,1]*",name,"[2,4]*",name,"[4,3]+",name,"[2,1]*",name,"[1,3]*",name,"[4,4]-",name,"[2,1]*",name,"[1,4]*",name,"[4,3]-",name,"[4,1]*",name,"[1,3]*",name,"[2,4]+",name,"[4,1]*",name,"[1,4]*",name,"[2,3] ") )
x<-rbind(x, paste0( " ",name,".inv.help[3,3] <- ",name,"[1,1]*",name,"[2,2]*",name,"[4,4]-",name,"[1,1]*",name,"[2,4]*",name,"[4,2]-",name,"[2,1]*",name,"[1,2]*",name,"[4,4]+",name,"[2,1]*",name,"[1,4]*",name,"[4,2]+",name,"[4,1]*",name,"[1,2]*",name,"[2,4]-",name,"[4,1]*",name,"[1,4]*",name,"[2,2] ") )
x<-rbind(x, paste0( " ",name,".inv.help[4,3] <- -",name,"[1,1]*",name,"[2,2]*",name,"[4,3]+",name,"[1,1]*",name,"[2,3]*",name,"[4,2]+",name,"[2,1]*",name,"[1,2]*",name,"[4,3]-",name,"[2,1]*",name,"[1,3]*",name,"[4,2]-",name,"[4,1]*",name,"[1,2]*",name,"[2,3]+",name,"[4,1]*",name,"[1,3]*",name,"[2,2] ") )
x<-rbind(x, paste0( " ",name,".inv.help[1,4] <- -",name,"[1,2]*",name,"[2,3]*",name,"[3,4]+",name,"[1,2]*",name,"[2,4]*",name,"[3,3]+",name,"[2,2]*",name,"[1,3]*",name,"[3,4]-",name,"[2,2]*",name,"[1,4]*",name,"[3,3]-",name,"[3,2]*",name,"[1,3]*",name,"[2,4]+",name,"[3,2]*",name,"[1,4]*",name,"[2,3] ") )
x<-rbind(x, paste0( " ",name,".inv.help[2,4] <- ",name,"[1,1]*",name,"[2,3]*",name,"[3,4]-",name,"[1,1]*",name,"[2,4]*",name,"[3,3]-",name,"[2,1]*",name,"[1,3]*",name,"[3,4]+",name,"[2,1]*",name,"[1,4]*",name,"[3,3]+",name,"[3,1]*",name,"[1,3]*",name,"[2,4]-",name,"[3,1]*",name,"[1,4]*",name,"[2,3] ") )
x<-rbind(x, paste0( " ",name,".inv.help[3,4] <- -",name,"[1,1]*",name,"[2,2]*",name,"[3,4]+",name,"[1,1]*",name,"[2,4]*",name,"[3,2]+",name,"[2,1]*",name,"[1,2]*",name,"[3,4]-",name,"[2,1]*",name,"[1,4]*",name,"[3,2]-",name,"[3,1]*",name,"[1,2]*",name,"[2,4]+",name,"[3,1]*",name,"[1,4]*",name,"[2,2] ") )
x<-rbind(x, paste0( " ",name,".inv.help[4,4] <- ",name,"[1,1]*",name,"[2,2]*",name,"[3,3]-",name,"[1,1]*",name,"[2,3]*",name,"[3,2]-",name,"[2,1]*",name,"[1,2]*",name,"[3,3]+",name,"[2,1]*",name,"[1,3]*",name,"[3,2]+",name,"[3,1]*",name,"[1,2]*",name,"[2,3]-",name,"[3,1]*",name,"[1,3]*",name,"[2,2] ") )
x<-rbind(x, paste0( " ",name,".det <- ",name,"[1,1]*",name,".inv.help[1,1]+",name,"[1,2]*",name,".inv.help[2,1]+",name,"[1,3]*",name,".inv.help[3,1]+",name,"[1,4]*",name,".inv.help[4,1] ") )
x<-rbind(x, paste0( " # ",name," is invertible if determinant is not 0 ") )
x<-rbind(x, paste0( " ",name,".invertible <- ifelse( ",name,".det==0, 0, 1 ) ") )
x<-rbind(x, paste0( " # if not invertible, modest mod of determinant ") )
x<-rbind(x, paste0( " ",name,".detmod1 ~ dbern(0.5) ") )
x<-rbind(x, paste0( " ",name,".detmod2 <- ifelse( ",name,".invertible==1, 0, (",name,".detmod1-0.5)/500 )") )
x<-rbind(x, paste0( " ",name,".det.adj <- ",name,".det + ",name,".detmod2 ") )
x<-rbind(x, paste0( " ",name,".inv <- ",name,".inv.help / ",name,".det.adj ") )
}
x<-rbind(x,"")
return( x )
}
make.inverse.2 <- function( name, dim ){
### *.prec.replace must exist (for dim>1)
x<-matrix( "", 1, 1 )
if (dim %in% 1) {
if (name %in% "Qt"){
x<-rbind(x, " # inverse of Qt ")
x<-rbind(x, " Qt.prec[1,1,t,p] <- (1/Qt[1,1,t,p]) ")
}
}
if (dim %in% 2) {
if (name %in% "Qt"){
x<-rbind(x, " ## inverse of Qt ")
x<-rbind(x, " # determinant of Qt ")
x<-rbind(x, " Qt.det[t,p] <- Qt[1,1,t,p]*Qt[2,2,t,p] - Qt[2,1,t,p]*Qt[1,2,t,p] ")
x<-rbind(x, " # Qt is invertible if determinant is not 0 ")
x<-rbind(x, " Qt.invertible[t,p] <- ifelse( Qt.det[t,p]==0, 0, 1 ) ")
x<-rbind(x, " # if not invertible, modest mod of determinant ")
x<-rbind(x, " Qt.detmod1[t,p] ~ dbern(0.5) ")
x<-rbind(x, " Qt.detmod2[t,p] <- ifelse( Qt.invertible[t,p]==1, 0, (Qt.detmod1[t,p]-0.5)/500 ) ")
x<-rbind(x, " Qt.det.adj[t,p] <- Qt.det[t,p] + Qt.detmod2[t,p] ")
x<-rbind(x, " # inverse of Qt ")
x<-rbind(x, " Qt.help[1,1,t,p] <- Qt[2,2,t,p] ")
x<-rbind(x, " Qt.help[2,2,t,p] <- Qt[1,1,t,p] ")
x<-rbind(x, " Qt.help[1,2,t,p] <- -1*Qt[1,2,t,p] ")
x<-rbind(x, " Qt.help[2,1,t,p] <- -1*Qt[2,1,t,p] ")
x<-rbind(x, " Qt.prec.orig[1:Aw,1:Aw,t,p] <- (1/Qt.det.adj[t,p] * Qt.help[,,t,p]) ")
x<-rbind(x, " ")
x<-rbind(x, " # check if Qt.prec.orig is positive definite ")
x<-rbind(x, " Qt.prec.orig.det[t,p] <- Qt.prec.orig[1,1,t,p]*Qt.prec.orig[2,2,t,p] - Qt.prec.orig[2,1,t,p]*Qt.prec.orig[1,2,t,p] ")
x<-rbind(x, " Qt.prec.orig.posdef[t,p] <- ifelse( Qt.prec.orig[1,1,t,p] > 0 && Qt.prec.orig.det[t,p] > 0, 1, 0 ) ")
x<-rbind(x, " ")
x<-rbind(x, " # if Qt.prec.orig is not pos def then replace by a pos def matrix ")
x<-rbind(x, " Qt.prec[1:2,1:2,t,p] <- ifelse( Qt.prec.orig.posdef[t,p]==1, Qt.prec.orig[,,t,p] , Qt.prec.replace[,] ) ")
}
}
x<-rbind(x,"")
return( x )
}
Try the ctglm package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.