Cmdstan - 1d integration

Setup

This section contains the setup and the various utility functions used throughout.

Libraries used:

library(data.table)
library(cmdstanr)

Introduction

There are a few ways we can do integration in R and stan. R provides the integrate function, for example:

integrate(dnorm, -0.4, 0.5)

## 0.3468842 with absolute error < 3.9e-15

but you can approximate this with an expecation if you can generate random samples from the function of interest:

x <- rnorm(1e6)
mean(x >= -0.4 & x <= 0.5)

## [1] 0.346977

Stan provides a mechanism to do integration via quadrature. Basically, you code up a user defined function that returns the value of the function you want to integrate at the value and with the parameters you supply to it. So, in this case, I just code up the normal density as follows:

functions{    
  // x - function is required to return the value of the function         
  // evaluated at x    
  // xc - complement of x on the domain - set to NaN when limit is inf.    
  // theta - parameters    
  // x_r - real data - to pass args that are not a function of params    
  // x_i - integer data    
  real normal_density(real x,                 
                    real xc,                
                    array[] real theta,     
                    array[] real x_r,                        
                    array[] int x_i) {                       
    real mu = theta[1];                     
    real sigma = theta[2];    
      
    return 1 / (sqrt(2 * pi()) * sigma) * exp(-0.5 * ((x - mu) / sigma)^2);    
  }    
}    
data{    
}
transformed data{     
  array[0] real x_r;
  array[0] int x_i;  
}                       
parameters{            
}
transformed parameters{
  array[2] real theta; 
  real intgl;               
                       
  theta = { 0, 1 };        
  intgl = integrate_1d(normal_density, -0.4, 0.5, theta, x_r, x_i);
}
model {
}

And then compile the model and sample using fixed_param = T as follows.

if(basename(getwd()) == "rpubs"){
  fname <- ".//stan//integration1.stan"
} else if(basename(getwd()) == "R"){
  fname <- "..//stan//integration1.stan"
}

mod <- cmdstan_model(stan_file = fname)

f1 <- mod$sample(
  chains = 1, 
  iter_sampling = 1,
  fixed_param = TRUE,
  refresh = 0
)

## Running MCMC with 1 chain...
## 
## Chain 1 finished in 0.0 seconds.

# f1$summary()
d1 <- data.table(f1$draws(format = "df"))
d1$intgl

## [1] 0.346884

The stan manual has a slightly more complex demonstration where the integral is used in the model step to estimate the parameters of a left truncated normal distribution. It isn’t explained in much detail other than to say that the normalisation term needs to be included, but this is the URL [https://mc-stan.org/docs/2_29/stan-users-guide/calling-the-integrator.html].

All for now.