run_mcmc crashes with a simple linear gaussian model

[jgslazzaro]

Hi, I am trying to estimate a multivariate linear gaussian model. However, R crashes whenever I call the run_mcmc function. A toy example with simulated data is below.
The only output I get is the progress bar and after completion Rstudio tells me that the R Session aborted without further explanation.

Kalman filtering works fine.

library(bssm)

set.seed(1)

########### Simulating fake data
set.seed(1)

Tt <- 100
p<-2 #numero de observaveis
m<- 2 #numero de estados 
y <- matrix(0,Tt,p) # observaveis
Z <- matrix(0,p,m) #matriz de observacao
Tm <- matrix(0,m,m)
H <- matrix(0,2,2)
R <- matrix(0,m,m)
R[1,1] <- 0.1
R[2,2] <- 0.01
eps_y1 <- rnorm(Tt,mean = 0, sd = 0.1)
eps_y2 <- rnorm(Tt,mean = 0, sd = 0.01)
W <- matrix(0,Tt,p)
W[,1] <- eps_y1
W[,2] <- eps_y2

for (i in 1:p){
  Z[i,i] <- 1
}
Tm[1,1] <- 1
Tm[2,1] <- 0.1

for (t in 2:Tt){
  y[t,] <- Tm %*% y[t-1,] + W[t,]

}

######### Estimation

#observation with NAs
y2 <- y
y2[97,1] <- NA

update_fn <- function(theta) {

  Z <- matrix(c(1, 0, 0, 1), 2, 2)
  R <- matrix(0,2,2)
  R[1,1] <- 0.1
  R[2,2] <- 0.01
  H <- matrix(0,2,2)
  # add third dimension so we have p x n x 1, p x m x 1, p x p x 1 arrays
  dim(Z)[3] <- dim(R)[3] <- dim(H)[3] <- 1
  Tm <- matrix(0,2,2)
  Tm[1,1] <- theta
  Tm[2,1] <- 0.1
  dim(Z)[3] <- 1
  dim(R)[3] <- 1
  dim(H)[3] <- 1
  dim(Tm)[3] <- 1
  P1 <- matrix(10,m,m)
  a1 <- c(0,0)
  list(Z = Z, R = R,T = Tm, H = H,P1=P1, a1=a1)
}

init_theta <- 1

comps <- update_fn(init_theta)

prior_fn <- function(theta) {
  sum(dnorm(theta, 1, 10, log = TRUE))  # log-jacobian
}

mod <- ssm_mlg(y2,   Z = comps$Z, D = comps$D, H = comps$H, T=comps$T, R=comps$R,
               a1 = comps$a1, P1=comps$P1,state_names = c("y1","y2"),update_fn = update_fn,
               init_theta = init_theta,
               prior_fn = prior_fn)
kf <- kfilter(mod)
fit <- run_mcmc(mod, iter = 1000, thin=2, burnin = 100,seed =5,verbose = TRUE)

Session Info

```r ```

[helske]

Thanks for the reprex, I'll look into this next year; I don't see anything obviously wrong in your code, so I need to run valgrind to check for memory leaks.

[helske]

Ok, solved the issue: You are defining P1 as 2x2 matrix with every element as 10, which is not a valid covariance matrix. This leads to failure of the cholesky decomposition in state simulation, and for some reason on windows this crashes R instead of throwing an error (as it does if you run sim_smoother(mod)). Setting P1 <- diag(10, m) should fix the model and allow the MCMC to run without issues.