Error in log marginal likelihood estimation

[adrtod]

The problem appears when using the conjugate linear sampler ConjugateNormal_knownPrec_linearMean and a non identity observation transition.

require(rbiips)

# Kalman filter
kf = function(y, m0 = 0, v0 = 10, 
              vx = 1, vy = 1,
              fx = 1, fy = 1) {
  N = length(y)
  x_post = rep(NA, N)
  v_post = rep(NA, N)
  
  # initial prediction
  x_pred = m0
  v_pred = v0
  l = 0
  
  # iterate over time
  for (k in seq_len(N)) {
    # update log-marginal likelihood
    l = l + dnorm(y[k], 
                  mean = fy*x_pred, 
                  sd = sqrt(vy+fy^2*v_pred), 
                  log = TRUE)
    
    # update current state
    v_post[k] = 1/(1/v_pred + fy^2/vy)
    x_post[k] = x_pred + v_post[k]*fy/vy * (y[k]-fy*x_pred)
    
    if (k < N) {
      # predict next state
      x_pred = fx*x_post[k]
      v_pred = fx^2*v_post[k] + vx
    }
  }
  
  return(list(x = x_post,
              v = v_post,
              l = l))
}

# Simulate data
# ------------------------------
set.seed(123)

N = 20
x_true = y = rep(NA, N)
prec_x_init_true = prec_x_true = prec_y_true = 1
mean_x_init_true = 0
r_x_true = .85
r_y_true = 1

x_true[1] = rnorm(1, mean_x_init_true, prec_x_init_true ^ -0.5)
y[1] = rnorm(1, r_y_true * x_true[1], prec_y_true ^ -0.5)
for (t in 2:N) {
  x_true[t] = rnorm(1, r_x_true * x_true[t - 1], prec_x_true ^ -0.5)
  y[t] = rnorm(1, r_y_true * x_true[t], prec_y_true ^ -0.5)
}

data = list(
  y = y,
  prec_x_init = prec_x_init_true,
  prec_y = prec_y_true,
  mean_x_init = mean_x_init_true,
  r_x = r_x_true
)

# Check Biips vs Kalman log marginal likelihood
# -------------------------------------
model_file = tempfile()

write(
  "model {
  prec_x ~ dgamma(1,1)
  r_y ~ dgamma(1,1)
  
  x[1] ~ dnorm(mean_x_init, prec_x_init)
  y[1] ~ dnorm(r_y*x[1], prec_y)
  for (t in 2:length(y)) {
  x[t] ~ dnorm(r_x*x[t-1], prec_x)
  y[t] ~ dnorm(r_y*x[t], prec_y)
  }
  }",
 file = model_file
)

n_part = 100

n_check = 100

prec_x_seq = seq(.1, 2, len = n_check)
r_y_seq = seq(.1, 2, len = n_check)

ll_prec_x_true = ll_prec_x_biips_prior = ll_prec_x_biips_auto = rep(NA, n_check)
ll_ry_true = ll_ry_biips_prior = ll_ry_biips_auto = rep(NA, n_check)

for (i in seq_len(n_check)) {
  #  r_y fixed
  # ----------------
  prec_x = prec_x_seq[i]
  r_y = r_y_true
  
  # True (Kalman)
  out_kf = kf(y,
              mean_x_init_true,
              1 / prec_x_init_true,
              1 / prec_x,
              1 / prec_y_true,
              r_x_true,
              r_y)
  
  ll_prec_x_true[i] = out_kf$l
  
  # Biips, proposal = "prior"
  data$prec_x = prec_x
  data$r_y = r_y
  model = biips_model(model_file, data = data)
  biips_build_sampler(model, proposal = "prior")
  out_smc = biips_smc_samples(model, n_part = n_part)
  
  ll_prec_x_biips_prior[i] = out_smc$log_marg_like
  
  # Biips, proposal = "auto"
  biips_build_sampler(model, proposal = "auto")
  out_smc = biips_smc_samples(model, n_part = n_part)
  
  ll_prec_x_biips_auto[i] = out_smc$log_marg_like
  
  # prec_x fixed
  # ----------------
  prec_x = prec_x_true
  r_y = r_y_seq[i]
  
  # True (KF)
  out_kf = kf(y,
              mean_x_init_true,
              1 / prec_x_init_true,
              1 / prec_x,
              1 / prec_y_true,
              r_x_true,
              r_y)
  
  ll_ry_true[i] = out_kf$l
  
  # Biips, proposal = "prior"
  data$prec_x = prec_x
  data$r_y = r_y
  model = biips_model(model_file, data = data)
  biips_build_sampler(model, proposal = "prior")
  out_smc = biips_smc_samples(model, n_part = n_part)
  
  ll_ry_biips_prior[i] = out_smc$log_marg_like
  
  # Biips, proposal = "auto"
  biips_build_sampler(model, proposal = "auto")
  out_smc = biips_smc_samples(model, n_part = n_part)
  
  ll_ry_biips_auto[i] = out_smc$log_marg_like
}


par(mfrow = c(1, 2))

plot(
  prec_x_seq,
  ll_prec_x_true,
  ylab = "log-marginal likelihood",
  xlab = "prec_x",
  main = "r_y fixed"
)
points(prec_x_seq, ll_prec_x_biips_prior, col = 2)
points(prec_x_seq, ll_prec_x_biips_auto, col = 3)
legend(
  "bottom",
  leg = c("True (Kalman)", "Biips proposal=prior", "Biips proposal=auto"),
  text.col = 1:3,
  bty = "n"
)

plot(r_y_seq,
     ll_ry_true,
     ylab = "log-marginal likelihood",
     xlab = "r_y",
     main = "prec_x fixed")
points(r_y_seq, ll_ry_biips_prior, col = 2)
points(r_y_seq, ll_ry_biips_auto, col = 3)
legend(
  "bottom",
  leg = c("True (Kalman)", "Biips proposal=prior", "Biips proposal=auto"),
  text.col = 1:3,
  bty = "n"
)

dev.copy(png,
         file = "log_marg_like.png",
         width = 640,
         height = 480)
dev.off()