1

Say I have two arrays in python and I wish to get (and actually use) the cubic spline interpolation between those points. (IE: I wish to integrate the function). I would strongly prefer a way using numpy scipy.

I know about scipy.interpolate.interp1d. However that only allows me to evalute the points, for say the very simply function of:

Now I could do something simply like: 

import numpy as np
import scipy.interpolate
import matplotlib.pyplot as plt

y = np.array([0,2,3,4,8,10,12,12,12,10,9,8,7,1,0,0,1,2])
x = np.array(range(len(y)))
xvals = np.linspace(0, len(y)-1, len(y)*100, endpoint = False)
func = scipy.interpolate.interp1d(x, y, kind = "cubic")
yvals = func(xvals)
plt.plot(xvals,yvals)
plt.plot(x,y, "o")

However I wish to further process this cubic spline (ie I need to get the integration).. For manual doing things I need to get the factors, so:

a_i * x^3 + b_i * x^2 + c_i * x + d_i where i goes from 0 to n/3

(n = number of elemetns - this is just the definition of the ith cubic)

I hence expect a list of tuples (or 2d array) describing all splines. - Or a way to get the ith cubic, and really, really would love to get have a convenience "x-to-i" to find in which spline I am currently.

(Though of course this latter problem is a simple search for first value larger than reference in a sorted list - I could do that by hand easily if need be).

paul23
  • 8,799
  • 12
  • 66
  • 149

3 Answers3

1

For interpolation, you can use scipy.interpolate.UnivariateSpline(..., s=0).

It has, among other things, the integrate method.

EDIT: s=0 parameter to UnivariateSpline constructor forces the spline to pass through all the data points. The result is in the B-spline basis, you can get the knots and coefficients with get_coefs() and get_knots() methods. The format is the same as used in FITPACK, dierckx, on netlib. Note though that the tck format used internally by interp1d (which relies on splmake ATM) and UnivariateSpline (alternatively, splrep/splev) are inconsistent.

EDIT2: you can get a piecewise polynomial representation of a spline with PPoly.from_spline --- but this is not consistent with iterp1d. Use splrep(..., s=0) to get an interpolating spline, then convert the result.

ev-br
  • 24,968
  • 9
  • 65
  • 78
  • The univariatespline is not a cubic spline, instead it is a so called 'b-spline', this is a complete different thing and has very very little to do with what I wish to achieve. - for one thing it doesn't do a curve fitting through the points. A cubic spline is defined as a series of 3rd order polynomials, each fitting two points, with the additional constraint that at each end point the splines have the same derivative. – paul23 Apr 05 '15 at 18:35
  • interp1d with kind=cubic works in the b-spline basis as well. If you want a local interpolant in scipy, use BPoly.from_derivatives – ev-br Apr 06 '15 at 06:54
  • Aren't splines usually built form b-splines? – Dux Apr 13 '15 at 14:03
  • 1
    @Dux those are two completelly different mathematical fields. cubic splines are a form of interpolation. b-splines a form of approximation. cubic splines WILL go through all points, B-splines will NOT guarantee that. – paul23 Apr 13 '15 at 15:56
  • 1
    Of course you can express a cubic interpolating spline in the b-spline basis. @Dux see my edited answer. – ev-br Apr 13 '15 at 19:36
  • @paul23: following up to the last comment, b-splines form a basis of splines, therefore, you can express a spline going through all data points in a b-spline-representation of a spline. – Dux Apr 14 '15 at 08:46
  • @ev-br - well then I actually need mathematical "proof" that those two are equal. - I can't see how they are, but that might be my limited knowledge. But I certainly can't use something without knowing under what circumstances the results are true. What is the mathematics governing the B-spline, and how does it fit up to the cubic splines. Like I said: I'm less interested in looks, but more so into it's mathematical quantities. (this spline will be used for further analysis such as integration and root finding) – paul23 Apr 14 '15 at 09:10
  • @paul23 for proofs, see eg books of de Boor or Lyche or Schumaker – ev-br Apr 14 '15 at 09:17
  • @paul23 if you want "proof", use `UnivariateSpline.derivatives(x)` to find the derivatives of the spline representation at point `x`. You can then check for yourself, integrate for yourself, and compare the results. – Dux Apr 14 '15 at 09:24
  • @Dux - that's just shows things, it is in no way mathematical sound - I honestly can't use that at all in a report, and I looked into it and I can't proof the functions myself either. – paul23 Apr 14 '15 at 09:26
  • @paul23, i don't get what you want. If you want to have the cubic coefficients, you can easily calculate them from the derivatives. I can submit an answer that shows you how, if you want. If you want the spline representation, you can get them from `.get_coefs` and calculate it from the b-spline basis. I suggest you look up wikipedia on b-splines for that. – Dux Apr 14 '15 at 09:50
1

Just a thought that might help. From the docs, you can get a cubic spline interpolation in a different way which may help you:

import numpy as np
import scipy.interpolate
import matplotlib.pyplot as plt

y = np.array([0,2,3,4,8,10,12,12,12,10,9,8,7,1,0,0,1,2])
x = np.array(range(len(y)))
xvals = np.linspace(0, len(y)-1, len(y)*100, endpoint = False)
func = scipy.interpolate.splrep(x, y, s=0)
yvals = scipy.interpolate.splev(xvals, func, der=0)

# display original vs cubic spline representation for security...
plt.figure()
plt.plot(x, y, 'x', xvals, yvals, x, y, 'b')
plt.legend(['Linear', 'Cubic Spline'])
plt.axis([-0.05, 20, -2, 20])
plt.title('Cubic-spline interpolation')
plt.show()

This gives you access to coefficients should you want via

pp = scipy.interpolate.spltopp(func[0][1:-1],func[1],func[2])

#Print the coefficient arrays, one for cubed terms, one for squared etc
print(pp.coeffs)

It also gives an example on the page of how to integrate using this cubic spline representation (changed the constant to suit your situation, I hope - your mileage may vary):

def integ(x, tck, constant=0):
    x = np.atleast_1d(x)
    out = np.zeros(x.shape, dtype=x.dtype)
    for n in xrange(len(out)):
        out[n] = scipy.interpolate.splint(0, x[n], tck)
    out += constant
    return out

const_of_integration = 0
yint = integ(xvals, func, const_of_integration)
J Richard Snape
  • 20,116
  • 5
  • 51
  • 79
0

The function scipy.signal.cspline1d returns the coefficients of the b-spline terms with mirror-symmetric boundary conditions. It is useful if you need to perform additional filtering in the coefficients

Sergio
  • 677
  • 1
  • 10
  • 20