[SciPy-User] help interpreting univariate spline
Elliot Hallmark
permafacture at gmail.com
Mon Apr 26 15:26:31 EDT 2010
On Fri, Apr 23, 2010 at 6:53 PM, Elliot Hallmark <permafacture at gmail.com> wrote:
> I am wanting to use a scipy interpolation routine to generate
> polynomials approximating function I have sampled.
> ...
> import numpy as np
>from scipy.interpolate import splrep, UnivariateSpline
>x = np.linspace(1, 10, 100)
>y = 1/x #approximate a hyperbola
>g = UnivariateSpline(x, y, w=None, k=3, s=.0005)
>print g.get_knots()
>print g.get_coeffs()
>and the output was:
>[ 1. 2.18181818 3.27272727 5.54545455 10. ]
>
>[ 0.98462432 0.66575572 0.423 0.25461626 0.13929847 0.11763985
0.09929961]
>That is 5 knots and 7 coefficients for a degree=3 spline. naively, I
>was expecting 4 coefficents for each interval between knots.
I found this post:
http://mail.scipy.org/pipermail/scipy-user/2009-June/021572.html
This seems very round-a-bout and ignores the piecewise nature of the
spline. I mean, really? is there no direct or documented way to get
the polynomial coefficents from a spline in python? If not, then this
code is the best I have found. For my purposes I have modified it
(see below)
"""
try_spline_coeff.py
explorations into getting polynomial coefficents from a spline.
Modified from http://mail.scipy.org/pipermail/scipy-user/2009-June/021572.html
by josef.pktd at gmail
please email offonoffoffonoff at gmail with improvements to this technique.
"""
import numpy as np
import scipy as sp
from scipy import interpolate
def cubice(x,a):
"evaluate a cubic"
return (a[:,0]+a[:,1]*x+a[:,2]*x**2+a[:,3]*x**3)
def deriv2coeff(x2, tck):
'''construct a piecewise cubic polynomial from derivatives of
spline at given knots (x2)
output is an array of knot abscissa and coefficent of poly for that section.
out = (x_knot,a,b,c,d)'''
if tck[2] != 3:
print "interpolation must be cubic!"
return np.inf
coeffr = np.vstack(interpolate.spalde(x2,tck))
coeffr[:,3] /= 6.
coeffr[:,2] = (coeffr[:,2] - 6*coeffr[:,3]*x2)/2.
coeffr[:,1] = coeffr[:,1] - 2*coeffr[:,2]*x2 - 3*coeffr[:,3]*x2**2
coeffr[:,0] = coeffr[:,0] - coeffr[:,1]*x2 - coeffr[:,2]*x2**2 -
coeffr[:,3]*x2**3
return np.concatenate((np.resize(x2,(x2.size,1)),coeffr), axis=1)
def calc_poly(xs, coeff):
'''determine which piece of the piecewise polynomial the given x
value corisponds to
and return the value of the approximation based on evaluating the
correct polynomial'''
result = np.zeros((xs.size,4))
for i,x in enumerate(xs):
a = b = c = d = 0
for x_knot,a1,b1,c1,d1 in coeff:
if x>x_knot: #x_knot is in ascending order (required by splrep also)
a=a1
b=b1
c=c1
d=d1
result[i] = a,b,c,d
#print "single: ",x, result[i]
return cubice(xs,result)
npoints = 51
x = np.linspace(1, 20, 2*npoints)
y = 1/x
tck = interpolate.splrep(x, y, k=3)
#calculate coefficents at original x points,
# would using tck knot values be more faithful to original spline?
coeffr = deriv2coeff(x, tck)
#print coeffr
#now, compare the polynomial and original spline over a different dataset
x2 = np.linspace(3, 18, 40)
y_spline = interpolate.splev(x2, tck)
y_polynomial = calc_poly(x2,coeffr)
print "max error: ", np.max(np.abs(y_spline-y_polynomial))
#using the actual knots of the spline
coeffk = deriv2coeff(tck[0][3:-3], tck)
y_polynomial = calc_poly(x2,coeffk)
print "max error using spline's knots: ", np.max(np.abs(y_spline-y_polynomial))
More information about the SciPy-User
mailing list