biips / rbiips Goto Github PK

R package for Bayesian inference with interacting particle systems

Home Page: https://biips.github.io/rbiips/

R 73.37% C++ 25.62% sed 1.01%

r bugs-language jags biips graphical-models sequential-monte-carlo particle-filter particle-mcmc mcmc-sampler inference-engine

rbiips's Introduction

Biips

Version: 0.11.0
Last modified: 2017-01-31
Maintainer: Adrien Todeschini [email protected]
License: GPL-3
Website: https://biips.github.io

Biips is a general software for Bayesian inference with interacting particle systems, a.k.a. sequential Monte Carlo (SMC) methods. It aims at popularizing the use of these methods to non-statistician researchers and students, thanks to its automated "black box" inference engine. It borrows from the BUGS/JAGS software, widely used in Bayesian statistics, the statistical modeling with graphical models and the language associated with their descriptions.

Context

Bayesian inference consists in approximating an unknown parameter dependent conditional probability law given a set of observations. A large number of problems, e.g. non-supervised classification, filtering, etc., can be addressed having as basis the aforementioned formulation. The underlying probability law, while not calculable in analytical manner for the general case, can be approximated by using Monte Carlo Markov Chain (MCMC) methods. These methods are popular for bayesian inference thanks to BUGS software and the WinBUGS graphical interface.

Emerged as a result of recent research studies, interacting particle based algorithms – a.k.a. Sequential Monte Carlo (SMC) methods in which the most common implementation is the particle filter – proved to have superior performance when compared to classical MCMC approaches. What is more, interacting particle algorithms are well adapted for dynamic estimation problems as encountered, for example, in filtering, tracking or classification problems. They do not require burn in convergence time and include calculation of the normalizing constant.

Features

BUGS language compiler adapted from JAGS
Black-box inference engine with:
- SMC/particle algorithms for filtering and smoothing
- Static parameter estimation using particle MCMC
- Automatic choice of the proposal samplers
Core developped in C++
R interface: Rbiips
Matlab/Octave interface: Matbiips
Easy language extensions with custom R and Matlab functions
Multi-platform: Linux, Windows, Mac OSX
Free and open source (GPL)

Roadmap

Code optimization
Parallelization
Particle Gibbs algorithm

Contact

To contact us with non-public matters, please write to [email protected]. Otherwise you can use the discusssion forum.

Getting Help

The GitHub issues can be used for any question or feedback about Biips. This is the best place for getting help.

Development

The core software libraries are written in C++. The Rbiips R package allows running Biips from the R statistical software and provides posterior analysis and plotting functions. The Matbiips toolbox provides similar capabilities for Matlab/Octave.

The Biips source code is hosted on GitHub.

License

Biips is licensed under the GPL-3 License. You may freely modify and redistribute it under certain conditions (see the file COPYING.txt for details).

The source code is hosted on GitHub.

Authors

Biips code is Copyright (C) Inria, 2012-2017.

Authors:

Adrien Todeschini (main developer)
Francois Caron
Marc Fuentes

Biips code is adapted from:

JAGS, Copyright (C) Martyn Plummer, 2002-2010
SMCTC, Copyright (C) Adam M. Johansen, 2008-2009

Additional information regarding adapted open source software is included in the file NOTICES.txt.

Thanks

rbiips's People

Contributors

Stargazers

Watchers

Forkers

marcos-cascone kuonanhong

rbiips's Issues

Error in installation

Hi all,
I ried to install the package but whatever way I follow (your instructions or using git2r or manually) I receive always the same error

fatal error: Console.hpp: No such file or directory
#include <Console.hpp>

Could you help me please?
Thanks a lot
Laura Vignola

Error in log marginal likelihood estimation

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()