The result of integrate_1d is confusing

Modeling
,
1

I use the normal_density function provided by stan-users-guide:

real normal_density(real x, real xc, real[] theta, real[] x_r, int[] x_i) {
  real mu = theta[1];
  real sigma = theta[2];
  return 1 / (sqrt(2 * pi()) * sigma) * exp(-0.5 * ((x - mu) / sigma)^2);
}

And I compute the integrals:

print("integral [N -1.5]: ", integrate_1d(normal_density, negative_infinity(), -1.54835, {0.0, 1.0}, x_r, x_i, 1e-8))
print("integral [-1.5 P]: ", integrate_1d(normal_density, -1.54835, positive_infinity(), {0.0, 1.0}, x_r, x_i, 1e-8))
print("integral [N 1.5]: ", integrate_1d(normal_density, negative_infinity(), 1.54835, {0.0, 1.0}, x_r, x_i, 1e-8))
print("integral [1.5 P]: ", integrate_1d(normal_density, 1.54835, positive_infinity(), {0.0, 1.0}, x_r, x_i, 1e-8))

The results are:

integral [N -1.5]: 0.939231  // this one is confusing
integral [-1.5 P]: 0.939231
integral [N 1.5]: 0.939231
integral [1.5 P]: 0.060769

The integral of (negative_infinity(), -1.54835) is 0.939231. Is it a bug? I use pystan 2.19.1.1.

3 Likes
2

This is indeed weird. Looks like a bug. I just verified it also happens with recent CmdStan.

Could you please file it at https://github.com/stan-dev/stan/issues/ ? It would help if you tried if you can reproduce the issue with some even simpler integrand (say 1 / ( 1 + x^2) which should integrate to \pi over the whole real line)

Thanks for noticing this!

2 Likes
3

Done. I stick to the normal_density in case that I give a wrong density function.

3 Likes