integrate_1d with normal_density generates wrong integral

[wwguo]

Summary:

integrate_1d with normal_density generates wrong integral.

Description:

integrate_1d with the normal_density function provided in stan-users-guide gets wrong integral.

Reproducible Steps:

The normal_density function:

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);
}

Compute following 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))

Current Output:

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

Expected Output:

integral [N -1.5]: 0.060769
integral [-1.5 P]: 0.939231
integral [N 1.5]: 0.939231
integral [1.5 P]: 0.060769

Current Version:

PyStan 2.19.1.1 and recent CmdStan

[bob-carpenter]

Thanks much for the careful bug report. It's a weird bug in that it seems to be sensitive to the sign of the second argument. The error's going to be in the sign logic in stan/math/prim/functor/integrate_1d.hpp Given that there's an if (b <= 0.0) nested inside of an if (std::isinf(a)), the problem's somewhere near line 71. @charlesm93 ---can you take a look at this?

[charlesm93]

@bob-carpenter yes, I'll investigate.

[bbbales2]

@wwguo good find. Can confirm this happens in 2.25 with the provided code:

functions {
  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);
  }
}

transformed data {
  real x_r[0];
  int x_i[0];
}

parameters {
  real x;
}

model {
  x ~ normal(0.0, 1.0);

  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));
}

but it seems to be fixed in develop. I changed how the limits in the integrals work here: stan-dev/math#2183 , and it looks like there was a bug with:

auto f_wrap = [&](double x) { return f(-(x + b), NOT_A_NUMBER); };

This was for computing an integral [-inf, b] for b < 0.0, and it was done with an integral [0.0, inf).

You can see that when evaluated at 0.0, the transformed function -(x + b) puts one end of the integral at -b, which is wrong.

That pull changed how these endpoints work, and this bug seems fixed now.