[SciPy-User] Interpolated Univariate Spline - what are coefs, and how is integral calculated?

[Evgeni Burovski]

On Thu, Aug 11, 2016 at 4:40 AM, Andrew Nelson <andyfaff at gmail.com> wrote:
> I have a few questions about interpolate.InterpolatedUnivariateSpline (IUS).
>
> My ultimate aim, given an IUS object, a lower definite integration limit,
> and an area, is to obtain the upper definite integration limit by the
> fastest route possible.
>
> The application is to do with probability distribution functions. The IUS
> represents the PDF. I can quickly work out the CDF using IUS.integral(a, x)
> = q, where a is the lower limit of support and x is the value for which you
> want to work out the CDF.
>
> However, I would like the quickest way to calculate the percent probability
> function (PPF). I.e. given q, work out x.
>
> One way used in scipy.stats is to use optimize.brentq to find this value -
> It's known that the CDF is a monotonically increasing function. It strikes
> me that brentq might not be the fastest way of achieving this. For example,
> I could cache the integral at each of the datapoints for the IUS, which
> would immediately allow me to bracket the location of x. Given that I know
> the degree of the smoothing spline I am hoping to calculate x quickly if I
> know the polynomial coefficients of the spline in the bracketing interval.
> This would done by integrating the polynomial and using the roots of the
> integral. My question is: Is there a way of getting those polynomial
> coefficients for each interval from the output of the IUS.get_coeffs method?
>

First of all, you likely cannot beat an inverse interpolation if raw
speed is what you're after (i.e., interpolate cdf(x) vs x).

Spefically, to your question about IUS:

* note that IUS is *not* guaranteed to be monotonic even if you
tabulate what looks like a monotonic function
* IUS is defined in the B-spline basis, so .get_coeffs() returns the
coeficients of B-spline basis elements, not polynomial representation
coefficients.

If you want polynomial representation, you can use
`PPoly.from_spline()` with knots and coefficients from IUS --- check
the docstring of LSQUnivariateSpline for a possible catch with the
boundary knots though. If you go that route, you might use `splrep(x,
y, s=0)` rather than IUS though, because it's equivalent and gives
directly tck (knots, coefficients and order) which `PPoly.from_spline`
takes as a parameter.

However, if you convert to polynomials, you might as well just use the
`CubicSpline` (new in scipy 0.18.0) which constructs the spline
directly in the polynomial basis.

Unless you want to make sure you get a monotonic interpolator and can
tolerate C1 continuity, in which case you can use instead PCHIP or
Akima1DInterpolator.

For integrals specifically, one thing which might (or might not) help
is to use Bernstein polynomials (which e.g. PchipInterpolator is), and
use the fact that $\int_0^{1} b_{j, n} = \frac{1}{n+1} $, so that you
know the cumulative integrals at each breakpoint.

HTH,

Evgeni