Computing One Dimensional Integrals
Definite and indefinite one dimensional integrals can be performed in Stan using the 1D integrators, integrate_1d_double_exponential and integrate_1d_gauss_kronrod.
As an example, the normalizing constant of a left-truncated normal distribution is
\[ \int_a^\infty \frac{1}{\sqrt{2 \pi \sigma^2}} e^{-\frac{1}{2}\frac{(x - \mu)^2}{\sigma^2}}. \]
To compute this integral in Stan, the integrand must first be defined as a Stan function (see the Stan Reference Manual chapter on User-Defined Functions for more information on coding user-defined functions).
real normal_density(real x, // Function argument
real xc, // Complement of function argument
// on the domain (defined later)
real mu, // mean
real sigma) { // standard deviation
return 1 / (sqrt(2 * pi()) * sigma) * exp(-0.5 * ((x - mu) / sigma)^2);
}This function is expected to return the value of the integrand evaluated at point x. The argument xc is used in definite integrals to avoid loss of precision near the limits of integration and is set to NaN when either limit is infinite (see the section on precision loss in the chapter on Higher-Order Functions of the Stan Functions Reference for details on how to use this). For integrate_1d_gauss_kronrod, xc is always NaN.
The arguments after xc are variadic: any number of arguments of any type may follow, and they are passed unmodified from the integrator call to the integrand. They may be data or parameters; arguments that are only ever data should be marked with the data keyword for efficiency, since extra parameter arguments increase the cost of propagating derivatives through the integral.
Calling the integrator
Suppose that our model requires evaluating the lpdf of a left-truncated normal, but the truncation limit is to be estimated as a parameter. Because the truncation point is a parameter, we must include the normalization term of the truncated pdf when computing our model’s log density. Note this is just an example of how to use the 1D integrator. The more efficient way to perform the correct normalization in Stan is described in the chapter on Truncated or Censored Data of this guide.
Such a model might look like (include the function defined at the beginning of this chapter to make this code compile):
data {
int N;
array[N] real y;
}
parameters {
real mu;
real<lower=0.0> sigma;
real left_limit;
}
model {
mu ~ normal(0, 1);
sigma ~ normal(0, 1);
left_limit ~ normal(0, 1);
target += normal_lpdf(y | mu, sigma);
target += -log(integrate_1d_double_exponential(normal_density,
left_limit,
positive_infinity(),
mu, sigma));
}Choosing the integrator
integrate_1d_double_exponential is robust to algebraic or logarithmic singularities of the integrand at the integration limits, and provides the xc argument for working near those limits (see Avoiding precision loss). integrate_1d_gauss_kronrod is more accurate on smooth integrands and on integrands whose mass lies away from the limits (for example a normal density times a smooth function), but it cannot resolve endpoint singularities and always passes xc as NaN. The two functions share the same call signature; see the 1D integrators section of the Functions Reference for details.
Limits of integration
The limits of integration can be finite or infinite. The infinite limits are made available via the Stan calls negative_infinity() and positive_infinity().
If both limits are either negative_infinity() or positive_infinity(), the integral and its gradients are set to zero.
Data vs. parameters
The variadic arguments passed to the integrand may be data or parameters. Arguments that only ever involve data or transformed data variables should be marked with the data keyword, which avoids the extra cost of propagating derivatives through them.
The endpoints of integration can be data or parameters (and internally the derivatives of the integral with respect to the endpoints are handled with the Leibniz integral rule).
Integrator convergence
The integral is performed with iterative quadrature methods implemented in the Boost library (Agrawal et al. 2017). The integration succeeds when
\[ \text{error} \leq \max(\text{relative\_tolerance} \cdot |I|, \text{absolute\_tolerance}), \]
where \(|I|\) is the estimated norm of the integral. The tolerances and the maximum iteration count can be set through the _tol variants (integrate_1d_double_exponential_tol and integrate_1d_gauss_kronrod_tol), which take relative_tolerance, absolute_tolerance, and max_refinements (max_depth for the Gauss-Kronrod variant) before the variadic arguments. The non-_tol variants use relative_tolerance equal to the square root of the machine epsilon of double precision floating point numbers (about 1e-8) and absolute_tolerance equal to zero, which reduces the test to a pure relative-tolerance criterion.
A positive absolute_tolerance lets the integration succeed when \(|I|\) is so small that the relative-tolerance test is dominated by floating-point round-off, for instance when the integrand is essentially zero across the interval (as happens with nested integration in the tails of a density).
If the integrator cannot reach the requested tolerance an exception is raised with a message like “Exception: integrate: error estimate of integral 4.25366e-13 exceeds max(relative_tolerance * L1, absolute_tolerance)”. If this occurs in the transformed parameters or model block, the result has the same effect as assigning a \(-\infty\) log probability, which causes rejection of the current proposal in MCMC samplers and adjustment of search parameters in optimization. If it occurs in the generated quantities block, the returned value is NaN. In these cases, a larger tolerance can be specified.
Zero-crossing integrals
When using integrate_1d_double_exponential, integrals on the (possibly infinite) interval \((a, b)\) that cross zero are split into two integrals, one from \((a, 0)\) and one from \((0, b)\), because the double-exponential quadrature can have difficulty near zero. Each integral is separately integrated to the given tolerances. integrate_1d_gauss_kronrod does not split zero-crossing integrals.
Avoiding precision loss near limits of integration in definite integrals
The xc argument used in this section is provided by integrate_1d_double_exponential; integrate_1d_gauss_kronrod always passes xc as NaN, so these techniques apply when using the double-exponential integrator.
If care is not taken, the quadrature can suffer from numerical loss of precision near the endpoints of definite integrals.
For instance, in integrating the pdf of a beta distribution when the values of \(\alpha\) and \(\beta\) are small, most of the probability mass is lumped near zero and one.
The pdf of a beta distribution is proportional to
\[ p(x) \propto x^{\alpha - 1}(1 - x)^{\beta - 1} \]
Normalizing this distribution requires computing the integral of \(p(x)\) from zero to one. In Stan code, the integrand might look like:
real beta(real x, real xc, real alpha, real beta) {
return x^(alpha - 1.0) * (1.0 - x)^(beta - 1.0);
}The issue is that there will be numerical breakdown in the precision of 1.0 - x as x gets close to one. This is because of the limited precision of double precision floating numbers. This integral will fail to converge for values of alpha and beta much less than one.
This is where xc is useful. It is defined, for definite integrals, as a high precision version of the distance from x to the nearest endpoint — a - x or b - x for a lower endpoint a and an upper endpoint b. To make use of this for the beta integral, the integrand can be re-coded:
real beta(real x, real xc, real alpha, real beta) {
real v;
if(x > 0.5) {
v = x^(alpha - 1.0) * xc^(beta - 1.0);
} else {
v = x^(alpha - 1.0) * (1.0 - x)^(beta - 1.0);
}
return v;
}In this case, as we approach the upper limit of integration \(a = 1\), xc will take on the value of \(a - x = 1 - x\). This version of the integrand will converge for much smaller values of alpha and beta than otherwise possible.
Consider another example: let’s say we have a log-normal distribution that is both shifted away from zero by some amount \(\delta\), and truncated at some value \(b\). If we were interested in calculating the expectation of a variable \(X\) distributed in this way, we would need to calculate \[ \int_a^b xf(x)\,dx = \int_{\delta}^b xf(x)\,dx \] in the numerator, where \(f(x)\) is the probability density function for the shifted log-normal distribution. This probability density function can be coded in Stan as:
real shift_lognormal_pdf(real x,
real mu,
real sigma,
real delta) {
real p;
p = (1.0 / ((x - delta) * sigma * sqrt(2 * pi()))) *
exp(-1 * (log(x - delta) - mu)^2 / (2 * sigma^2));
return p;
}Therefore, the function that we want to integrate is:
real integrand(real x, real xc, real mu, real sigma, real delta, real b) {
real numerator;
real p;
p = shift_lognormal_pdf(x, mu, sigma, delta);
numerator = x * p;
return numerator;
}What happens here is that, given that the log-normal distribution is shifted by \(\delta\), when we then try to integrate the numerator, our x starts at values just above delta. This, in turn, causes the x - delta term to be near zero, leading to a breakdown.
We can use xc, and define the integrand as:
real integrand(real x, real xc, real mu, real sigma, real delta, real b) {
real numerator;
real p;
if (x < delta + 1) {
p = shift_lognormal_pdf(xc, mu, sigma, delta);
} else {
p = shift_lognormal_pdf(x, mu, sigma, delta);
}
numerator = x * p;
return numerator;
}Why does this work? When our values of x are less than delta + 1 (so, when they’re near delta, given that our lower bound of integration is equal to \(\delta\)), we pass xc as an argument to our shift_lognormal_pdf function. This way, instead of dealing with x - delta in shift_lognormal_pdf, we are working with xc - delta which is equal to delta - x - delta, as delta is the lower endpoint in that case. The delta terms cancel out, and we are left with a high-precision version of x. We don’t encounter the same problem at the upper limit \(b\) so we don’t adjust the code for that case.
Note, xc is only used for definite integrals. If either the left endpoint is at negative infinity or the right endpoint is at positive infinity, xc will be NaN.
For zero-crossing definite integrals (see section Zero Crossing) the integrals are broken into two pieces (\((a, 0)\) and \((0, b)\) for endpoints \(a < 0\) and \(b > 0\)) and xc is a high precision version of the distance to the limits of each of the two integrals separately. This means xc will be a high precision version of a - x, x, or b - x, depending on the value of x and the endpoints.
Example: integrating out a latent variable with Gauss-Kronrod
A common use of 1D integration is to compute a marginal likelihood by integrating out a latent variable. Consider a Poisson count y_i with a per-observation varying intercept z given a normal prior z ~ normal(0, sigmaz), and a linear predictor mu_i. The likelihood marginal over z is
\[ p(y_i \mid \mu_i, \sigma_z) = \int_{-\infty}^{\infty} \text{normal}(z \mid 0, \sigma_z)\, \text{Poisson}(y_i \mid e^{z + \mu_i})\, dz, \]
a normal density times a Poisson probability mass function. The integrand is
real integrand(real z, real notused, real sigmaz, real mu_i, data int y_i) {
real p = exp(normal_lpdf(z | 0, sigmaz)
+ poisson_log_lpmf(y_i | z + mu_i));
return (is_inf(p) || is_nan(p)) ? 0 : p;
}The notused argument is the xc slot, which integrate_1d_gauss_kronrod always passes as NaN; y_i is marked data because it is never a parameter; and the is_inf/is_nan guard returns 0 where exp over- or under-flows in the tails.
Because the normal factor has thin tails, essentially all of the integrand’s mass lies within 8 * sigmaz of zero (at the limits the normal density is about \(e^{-32}\) of its peak), so the integral can be taken over the finite interval \([-8\,\sigma_z,\ 8\,\sigma_z]\). Gauss-Kronrod places its nodes across the whole interval and resolves the peak — which sits near the conditional mode of z, away from zero when the count is large — where double-exponential quadrature over an infinite interval can undersample. In generated quantities, where mu_i is the linear predictor for observation i and sigmaz is a parameter,
log_lik[i] = log(integrate_1d_gauss_kronrod(integrand,
-8 * sigmaz,
8 * sigmaz,
sigmaz, mu_i, y[i]));The variadic arguments sigmaz, mu_i, and y[i] are passed directly to the integrand, with no packing into data and parameter arrays, and the integration limits depend on the parameter sigmaz, which is allowed.