Definite integral in RStan

Modeling
1

@bbbales2 @jonah @nhuurre @Bob_Carpenter
Hi all experts and users, I am confusing when doing definite integrals in stan. Suppose that I want to integrate a function with given endpoints as follows, This is what I understand, but I know “xc” should be defined for definite case.
So Please explain it with an example you have or help me to develop my one.

functions{
real fun_integral(real x,
real xc,
real[] theta,
real[] x_r,
int[] x_i){

  real mu = theta[1];
  real sigma2 = theta[2];
  real fun;
  fun = (0.3*exp(-square(log(x)-mu)/2*sigma2)/sqrt(2*pi()*sigma2)*x);
  return(fun);
                }

}

integrate_1d(fun_integral,50,52,{mu,sigma2},{0.0}, x_i, 1e-8)

I am wasting days for that,

Thank you.

2

Did you know that as per Stan operator precedence rules

-square(log(x)-mu)/2*sigma2

parses as

((-square((log(x) - mu)) / 2) * sigma2)

Assuming that’s a mistake and not intentionally obfuscated code the integral is the log-normal cumulative distribution function. That’s implemented in Stan math library, you don’t need a custom integrate_1d

// integrate_1d(fun_integral,50,52,{mu,sigma2},{0.0}, x_i, 1e-8)
0.3 * (lognormal_cdf(52, mu, sigma2) - lognormal_cdf(50, mu, sigma2))
1 Like
3

Thank you Niko, here I mentioned one part of my equation going to be integrated.
Complete one is

((0.3 * exp(-square(log(x) - mu)) / (2 * sigma2)) /(sqrt(2*pi()*sigma2)*x) ) * (exp(-lambda * pow(tn-x,alpha)))

4

The User Guide has an example on how to use the xc argument. I think it’s only needed when there’s a singularity at the boundary. If alpha is positive you probably don’t need to worry about it.
If you want to use xc then it’s probably something like this

real WQ_integral(real x,
      real xc,
      real[] theta,
      real[] x_r,
      int[] x_i){

  real mu = theta[1];
  real sigma2 = theta[2];
  real lambda = theta[3];
  real alpha = theta[4];
  real tn = theta[5];
  if (tn - x < 0.5)
    WQ = 0.3 * exp( -square(log(x)-mu)/(2*sigma2) - lambda*(    xc^alpha) )/(sqrt(2*pi()*sigma2)*x;
  } else {
    WQ = 0.3 * exp( -square(log(x)-mu)/(2*sigma2) - lambda*((tn-x)^alpha) )/(sqrt(2*pi()*sigma2)*x;
  }
  return(WQ);
}
3 Likes
5

It’s a very clear explanation. Many thanks. :)

1 Like