R Under development (unstable) (2023-10-22 r85388 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> # HEADER ####################################################
> # This is file spam/tests/jss_areal_counts.R.               #
> # It is part of the R package spam,                         #
> #  --> https://CRAN.R-project.org/package=spam              #
> #  --> https://CRAN.R-project.org/package=spam64            #
> #  --> https://git.math.uzh.ch/reinhard.furrer/spam         #
> # by Reinhard Furrer [aut, cre], Florian Gerber [aut],      #
> #    Roman Flury [aut], Daniel Gerber [ctb],                #
> #    Kaspar Moesinger [ctb]                                 #
> # HEADER END ################################################
> 
> 
> # JSS article: 
> #     "Pitfalls in the implementation of Bayesian
> #      hierarchical modeling of areal count data.
> #      An illustration using BYM and Leroux models."
> #
> # test the MCMC sampler from the paper with 30 iterations.
> 
> 
> # SETUP:
> library("spam")
Spam version 2.10-0 (2023-10-23) is loaded.
Type 'help( Spam)' or 'demo( spam)' for a short introduction 
and overview of this package.
Help for individual functions is also obtained by adding the
suffix '.spam' to the function name, e.g. 'help( chol.spam)'.

Attaching package: 'spam'

The following objects are masked from 'package:base':

    backsolve, forwardsolve

> options(spam.structurebased=TRUE)
> 
> # BYM ---------------------------------------------
> data(Oral); attach(Oral) 
> path <- system.file("demodata/germany.adjacency", package = "spam")
> A <- adjacency.landkreis(path); n <- dim(A)[1]
> 
> set.seed(2)
> hyperA <- c(1, 1); hyperB <- c(0.5, .01)  
> totalg <- 30
> 
> upost <- vpost <- array(0, c(totalg, n))  
> kpost <- array(NA, c(totalg, 2)); accept <- rep(NA, totalg)
> upost[1,] <- vpost[1,] <- rep(.001, 544); kpost[1,] <- c(10, 100) 
> 
> eta <- upost[1,] + vpost[1,]
> C <- exp(eta) * E; diagC <- diag.spam(c(rep(0, n), C))
> b <- c( rep(0, n), Y + (eta - 1) * C)
> Qu <- R <- precmat.IGMRFirreglat(A); pad(Qu) <- c(2 * n, 2 * n)
> Qv <- as.spam(rbind(cbind( diag(n), -diag(n)),
+                       cbind(-diag(n),  diag(n))))
> Q <- kpost[1,1] * Qu + kpost[1,2] * Qv + diagC
> struct <- chol(Q, memory = list(nnzcolindices = 6467))
> uRuHalf <- t(upost[1,]) %*% (R %*% upost[1,]) / 2
> vvHalf <- t(vpost[1,]) %*% vpost[1,] / 2
> postshape <- hyperA + c(n - 1, n) / 2
> 
> for (i in 2:totalg) {
+      kpost[i,] <- rgamma(2, postshape, hyperB + c(uRuHalf, vvHalf))
+  
+      etaTilde <- eta
+      for(index in 1:2){
+          C <- E * exp(etaTilde)
+          diagC <- diag.spam(c(rep(0, n), C))
+          b <- c(rep(0, 544), Y + (etaTilde - 1) * C)
+          Q <- kpost[i,1] * Qu + kpost[i,2] * Qv + diagC
+          etaTilde <- c(solve.spam(Q, b,
+                                   Rstruct = struct))[1:n + n]
+      }
+  
+      C <- exp(etaTilde) * E; diagC <- diag.spam(c(rep(0, n), C))
+      b <- c( rep(0, n), Y + (etaTilde - 1) * C)
+      Q <- kpost[i,1] * Qu + kpost[i,2] * Qv + diagC
+ 
+      x_ <- c(rmvnorm.canonical(1, b, Q, Rstruct = struct))
+      upost[i,] <- x_[1:n]; eta_ <- x_[1:n + n]; vpost[i,] <- eta_ - upost[i,]
+      uRuHalf_ <- t(upost[i,]) %*% (R %*% upost[i,]) / 2
+      vvHalf_ <- t(vpost[i,]) %*% vpost[i,] / 2
+ 
+      etaTilde_ <- eta_
+      for(index in 1:2){
+        C_ <- E * exp(etaTilde_)
+        diagC_ <- diag.spam(c(rep(0, n), C_))
+        b_ <- c(rep(0, 544), Y + (etaTilde_ - 1) * C_)
+        Q_<- kpost[i,1] * Qu + kpost[i,2] * Qv + diagC_
+        etaTilde_ <- c(solve.spam(Q_, b_,
+                                 Rstruct = struct))[1:n + n]
+      }
+      
+      C_ <- exp(etaTilde_) * E; diagC_ <- diag.spam(c(rep(0, n), C_))
+      b_ <- c( rep(0, n), Y + (etaTilde_ - 1) * C_)
+      Q_ <- kpost[i,1] * Qu + kpost[i,2] * Qv + diagC_
+ 
+      logPost_ <- sum(Y * eta_ - E * exp(eta_)) -
+          kpost[i,1] * uRuHalf_ - kpost[i, 2] * vvHalf_
+      logPost  <- sum(Y * eta - E * exp(eta)) - kpost[i,1] * uRuHalf -
+          kpost[i,2] * vvHalf
+      logApproxX_ <- - kpost[i,1] * uRuHalf_ - kpost[i,2] * vvHalf_ -
+          sum(.5 * eta_^2 * C) + sum(b * eta_)
+      logApproxX  <- - kpost[i,1] * uRuHalf  - kpost[i,2] * vvHalf -
+          sum(.5 * eta^2 * C_) + sum(b_ * eta)
+      logAlpha <- min(0, logPost_ - logPost + logApproxX - logApproxX_)
+  
+      if (log(runif(1)) < logAlpha) {    
+          uRuHalf <- uRuHalf_;  vvHalf <- vvHalf_
+          eta <- eta_; b <- b_; C <- C_; accept[i] <- 1
+      } else{                            
+        accept[i] <- 0; upost[i,] <- upost[i-1,]; vpost[i,] <- vpost[i-1,]} 
+  }
>  
> # values of  30th iteration
> head(eta)
[1] -0.45371922  0.17297575  0.02605778 -0.44984751 -0.36053283  0.01309363
> tail(b)
[1] -1.341895 -3.756725 -2.514229 -4.411210 -6.486653 -5.101487
> head(C)
[1] 15.05715 53.22751 41.15492 13.16245 19.87138 30.27948
> tail(accept)
[1] 0 1 1 1 1 0
> tail(upost[30,])
[1]  0.255439209 -0.299099415 -0.004660022 -0.392147580 -0.318792450
[6] -0.344095560
> tail(vpost[30,])
[1] -0.0152950744  0.0008418105  0.0400659014 -0.0160847172  0.0064331753
[6] -0.0154282422
> sum(accept[-1])
[1] 21
> sum(upost)
[1] -780.3814
> 
> 
> proc.time()
   user  system elapsed 
   1.53    0.35    1.81