I have a set of points, and on those points I have used scipy in order to calculate the interpolating polynomial. I wish to have the primitive of that function
self.p=interpolate.interp1d(self.__discreteDistribution['x'], self.__discreteDistribution['y'], kind='cubic')
I can easily use scipy to calculate the value of the integral over an interval, using
integrate.quad(self.p, 0, max)
What I want instead is to have the primitive of self.p(). I have found sympy, but I do not have an analytical version of my interpolating polynomial.
What would you do in these occasion?
Assuming that you're using a piecewise interpolator (as opposed to a global polynomial interpolation), there are several ways you can do it with scipy:
Method 1: UnivariateSpline.
In [1]: import numpy as np
In [2]: x = np.arange(8)
In [3]: y = x
In [4]: from scipy.interpolate import interp1d
In [5]: from scipy.interpolate import interp1d, UnivariateSpline
In [6]: spl = UnivariateSpline(x, y, s=0)
In [7]: spl.<TAB>
spl.antiderivative spl.get_coeffs spl.roots
spl.derivative spl.get_knots spl.set_smoothing_factor
spl.derivatives spl.get_residual
spl.ext spl.integral
In [8]: spl.integral(0, 1)
Out[8]: 0.5000000000000001
Two quirks of UnivariateSpline: first, use s=0 for interpolation (as opposed to least-square fitting). Second, watch out for extrapolation for out-of-bounds. By default, UnivariateSpline extrapolates out-of-bounds values (this can be controlled in constructor), but .integral assumes the spline is zero out of bounds.
In [9]: spl.integral(-1, 1)
Out[9]: 0.5000000000000001
Method 2: splev, splrep and splint.
In [13]: from scipy.interpolate import splev, splint, splrep
In [14]: tck = splrep(x, y, s=0)
In [15]: splint(0, 1, tck)
Out[15]: 0.5000000000000001
This is equivalent to using UnivariateSpline, only the interface is a bit different. See the docs for details.
Method 3: interp1d.
Under the hood, interp1d also uses b-splines (unless you request kind='linear' or 'nearest'), but the evaluation routines are different.
interp1d constructs a callable, which can then be fed to a general-purpose integrator.
In [18]: from scipy.interpolate import interp1d
In [19]: interp = interp1d(x, y, kind='cubic')
In [20]: from scipy.integrate import quad
In [21]: quad(interp, 0, 1)
Out[21]: (0.5000000000000024, 5.5511151231258095e-15)
Again, watch out for out-of-bounds values: The behavior of the result constructed by interp1d is not very useful (even though it's controllable to an extent).