Better techniques for computing CDFs from INLA?

[tbenthompson]

Currently, we use Simpson's rule in some places and Gauss-hermite quadrature in others. We haven't considered the CDF step very much. I am not sure how much error is introduced by this step. Before exploring implementation options, we should simply do a convergence test with current quadrature rules to determine how much error is being introduced at this step.

[tbenthompson]

for simpson/trapezoid, the procedure would look something like:

• generate evenly spaced grid of points for full domain.
• given a point x that we want to compute CDF for:
• identify which of the points are to the left of the point x.
• interpolate a value for the point x itself.
• use simpson’s rule for the first n_odd points
• use trapezoid rule for the remainder

[tbenthompson]

for gaussian quadrature, the procedure would look like:

• generate gaussian quadrature rule covering the full domain.
• given a point x that we want to compute CDF for:
• identify which of the points are to the left of the point x.
• say this is M points. we don’t know the quadrature weights of these points yet.
• construct the linear algebra problem relating the integrals of the first M lagrange polynomials to the unknown quadrature weights.
• solve linear system for the quadrature weights
• perform the now-reduced quadrature.

[tbenthompson]

I also would like to investigate what the "normal" way of doing this is because this CDF computation from numerical integration on a pre-computed grid of points is surely a 100% solved problem. i wouldn’t be surprised if the simpson rule version is the standard approach.

[tbenthompson]

a third approach would be:

• generate gaussian quadrature rule covering the full domain.
• given a point x that we want to compute CDF for:
• generate a new gaussian quadrature rule cover (left_end, x)
• interpolate values on the new gaussian quadrature rule grid from the points on the original grid.
• integrate.

this would involve more runtime computation because the interpolation would require a matrix vector product.