Skip to content

Commit

Permalink
Fix errors with INLA_24.06.27
Browse files Browse the repository at this point in the history
  • Loading branch information
aritz-adin committed Aug 13, 2024
1 parent 6e64a09 commit 0d1d92f
Show file tree
Hide file tree
Showing 5 changed files with 112 additions and 112 deletions.
2 changes: 1 addition & 1 deletion DESCRIPTION
Original file line number Diff line number Diff line change
Expand Up @@ -31,5 +31,5 @@ License: GPL-3
Encoding: UTF-8
LazyData: true
LazyDataCompression: xz
RoxygenNote: 7.3.1
RoxygenNote: 7.3.2
Config/testthat/edition: 3
48 changes: 24 additions & 24 deletions R/Mmodel_icar.R
Original file line number Diff line number Diff line change
Expand Up @@ -20,102 +20,102 @@
#' @param theta Vector of hyperparameters.
#'
#' @return This is used internally by the \code{INLA::inla.rgeneric.define()} function.
#'
#'
#' @import Matrix
#'
#'
#'
#' @export
########################################################################
## Mmodels - iCAR - FE (BARTLETT DECOMPOSITION)
########################################################################
Mmodel_icar <- function(cmd=c("graph","Q","mu","initial","log.norm.const","log.prior","quit"), theta=NULL){

envir <- parent.env(environment())
if(!exists("cache.done", envir=envir)){
BlockIW <- kronecker(Diagonal(J),Diagonal(x=colSums(W))-W)
BlockIW <- kronecker(Matrix::Diagonal(J),Matrix::Diagonal(x=colSums(W))-W)

assign("BlockIW", BlockIW, envir=envir)
assign("cache.done", TRUE, envir=envir)
}

########################################################################
## theta
########################################################################
interpret.theta <- function(){
diag.N <- sapply(theta[as.integer(1:J)], function(x) { exp(x) })
no.diag.N <- theta[as.integer(J+1:(J*(J-1)/2))]
N <- diag(diag.N,J)

N <- diag(diag.N,J)
N[lower.tri(N, diag=FALSE)] <- no.diag.N

Covar <- N %*% t(N)

e <- eigen(Covar)
M <- t(e$vectors %*% diag(sqrt(e$values)))
# S <- svd(Covar)
# M <- t(S$u %*% diag(sqrt(S$d)))

return(list(Covar=Covar, M=M))
}

########################################################################
## Graph of precision function; i.e., a 0/1 representation of precision matrix
########################################################################
graph <- function(){ return(Q()) }

########################################################################
## Precision matrix
########################################################################
Q <- function(){
param <- interpret.theta()
M.inv <- solve(param$M)
MI <- kronecker(M.inv, Diagonal(nrow(W)))
Q <- (MI %*% BlockIW) %*% t(MI)
MI <- kronecker(M.inv, Matrix::Diagonal(nrow(W)))
Q <- (MI %*% BlockIW) %*% Matrix::t(MI)
return(Q)
}

########################################################################
## Mean of model
########################################################################
mu <- function(){ return(numeric(0)) }

########################################################################
## log.norm.const
########################################################################
log.norm.const <- function(){
val <- numeric(0)
return(val)
}

########################################################################
## log.prior: return the log-prior for the hyperparameters
########################################################################
log.prior <- function(){
param <- interpret.theta()

## n^2_jj ~ chisq(J-j+1) ##
val <- J*log(2) + 2 *sum(theta[1:J]) + sum(dchisq(exp(2*theta[1:J]), df=(J+2)-1:J+1, log=TRUE))

## n_ji ~ N(0,1) ##
val <- val + sum(dnorm(theta[as.integer(J+1:(J*(J-1)/2))], mean=0, sd=1, log=TRUE))

return(val)
}

########################################################################
## initial: return initial values
########################################################################
initial <- function(){ return(as.vector(initial.values)) }

########################################################################
########################################################################
quit <- function(){ return(invisible()) }

if(!length(theta)) theta <- initial()
val <- do.call(match.arg(cmd), args=list())

return(val)
}

utils::globalVariables(c("J","W","interpret.theta","Q","log.prior","initial",
"dchisq","dnorm","initial.values"))
"dchisq","dnorm","initial.values"))
50 changes: 25 additions & 25 deletions R/Mmodel_iid.R
Original file line number Diff line number Diff line change
Expand Up @@ -20,101 +20,101 @@
#' @param theta Vector of hyperparameters.
#'
#' @return This is used internally by the \code{INLA::inla.rgeneric.define()} function.
#'
#'
#' @import Matrix
#'
#'
#'
#' @export
########################################################################
## Mmodels - idd - FE (BARTLETT DECOMPOSITION)
########################################################################
Mmodel_iid <- function(cmd=c("graph","Q","mu","initial","log.norm.const","log.prior","quit"), theta=NULL){

envir <- parent.env(environment())
if(!exists("cache.done", envir=envir)){
BlockIW <- Diagonal(nrow(W)*J)
BlockIW <- Matrix::Diagonal(nrow(W)*J)

assign("BlockIW", BlockIW, envir=envir)
assign("cache.done", TRUE, envir=envir)
}

########################################################################
## theta
########################################################################
interpret.theta <- function(){
diag.N <- sapply(theta[as.integer(1:J)], function(x) { exp(x) })
no.diag.N <- theta[as.integer(J+1:(J*(J-1)/2))]
N <- diag(diag.N,J)

N <- diag(diag.N,J)
N[lower.tri(N, diag=FALSE)] <- no.diag.N

Covar <- N %*% t(N)

e <- eigen(Covar)
M <- t(e$vectors %*% diag(sqrt(e$values)))
# S <- svd(Covar)
# M <- t(S$u %*% diag(sqrt(S$d)))

return(list(Covar=Covar, M=M))
}

########################################################################
## Graph of precision function; i.e., a 0/1 representation of precision matrix
########################################################################
graph <- function(){ return(Q()) }

########################################################################
## Precision matrix
########################################################################
Q <- function(){
param <- interpret.theta()
M.inv <- solve(param$M)
MI <- kronecker(M.inv, Diagonal(nrow(W)))
Q <- (MI %*% BlockIW) %*% t(MI)
MI <- kronecker(M.inv, Matrix::Diagonal(nrow(W)))
Q <- (MI %*% BlockIW) %*% Matrix::t(MI)
return(Q)
}

########################################################################
## Mean of model
########################################################################
mu <- function(){ return(numeric(0)) }

########################################################################
## log.norm.const
########################################################################
log.norm.const <- function(){
val <- numeric(0)
return(val)
}

########################################################################
## log.prior: return the log-prior for the hyperparameters
########################################################################
log.prior <- function(){
param <- interpret.theta()

## n^2_jj ~ chisq(J-j+1) ##
val <- J*log(2) + 2 *sum(theta[1:J]) + sum(dchisq(exp(2*theta[1:J]), df=(J+2)-1:J+1, log=TRUE))

## n_ji ~ N(0,1) ##
val <- val + sum(dnorm(theta[as.integer(J+1:(J*(J-1)/2))], mean=0, sd=1, log=TRUE))

return(val)
}

########################################################################
## initial: return initial values
########################################################################
initial <- function(){ return(as.vector(initial.values)) }

########################################################################
########################################################################
quit <- function(){ return(invisible()) }

if(!length(theta)) theta <- initial()
val <- do.call(match.arg(cmd), args=list())

return(val)
}

utils::globalVariables(c("dchisq","dnorm","initial.values"))
utils::globalVariables(c("dchisq","dnorm","initial.values"))
Loading

0 comments on commit 0d1d92f

Please sign in to comment.