lagrange - Python interp1d vs. UnivariateSpline

signal interpolation python

Question

Estoy tratando de portar algún código de MatLab a Scipy, y he probado dos funciones diferentes de scipy.interpolate, interp1d y UnivariateSpline . Los resultados interp1d coinciden con la función InterLaT MatLab, pero los números Univariados se ven diferentes y, en algunos casos, muy diferentes.

f = interp1d(row1,row2,kind=''cubic'',bounds_error=False,fill_value=numpy.max(row2))
return  f(interp)

f = UnivariateSpline(row1,row2,k=3,s=0)
return f(interp)

¿Alguien podría ofrecer alguna idea? Mis vals x no están igualmente espaciados, aunque no estoy seguro de por qué eso sería importante.

Answer 1

UnivariateSpline: un contenedor más reciente de las rutinas de FITPACK.

esto podría explicar los valores ligeramente diferentes? (También experimenté que UnivariateSpline es mucho más rápido que interp1d).

Answer 2

Funciona para mi,

from scipy import allclose, linspace
from scipy.interpolate import interp1d, UnivariateSpline

from numpy.random import normal

from pylab import plot, show

n = 2**5

x = linspace(0,3,n)
y = (2*x**2 + 3*x + 1) + normal(0.0,2.0,n)

i = interp1d(x,y,kind=3)
u = UnivariateSpline(x,y,k=3,s=0)

m = 2**4

t = linspace(1,2,m)

plot(x,y,''r,'')
plot(t,i(t),''b'')
plot(t,u(t),''g'')

print allclose(i(t),u(t)) # evaluates to True

show()

Esto me da,

Answer 3

El motivo por el cual los resultados son diferentes (pero ambos probablemente sean correctos) es que las rutinas de interpolación utilizadas por UnivariateSpline e interp1d son diferentes.

• interp1d construye una B-spline suave utilizando los puntos x que le asignó como nudos

• UnivariateSpline se basa en FITPACK, que también construye una B-spline suave. Sin embargo, FITPACK intenta elegir nuevos nudos para la spline, para ajustar mejor los datos (probablemente para minimizar chi ^ 2 más alguna penalización por curvatura, o algo similar). Puede averiguar qué puntos de nudo utilizó a través de g.get_knots() .

Entonces, la razón por la que obtienes diferentes resultados es que el algoritmo de interpolación es diferente. Si desea B-splines con nudos en puntos de datos, use interp1d o splmake . Si quiere lo que hace FITPACK, use UnivariateSpline . En el límite de datos densos, ambos métodos dan los mismos resultados, pero cuando los datos son escasos, puede obtener resultados diferentes.

(¿Cómo sé todo esto? Leí el código :-)

Answer 4

Me encontré con el mismo problema.

Respuesta corta

Use InterpolatedUnivariateSpline en su lugar:

f = InterpolatedUnivariateSpline(row1, row2)
return f(interp)

Respuesta larga

UnivariateSpline es una "spline de suavizado unidimensional que se ajusta a un conjunto dado de puntos de datos", mientras que InterpolatedUnivariateSpline es una "spline de interpolación unidimensional para un conjunto determinado de puntos de datos". El primero suaviza los datos mientras que el último es un método de interpolación más convencional y reproduce los resultados esperados de interp1d . La figura a continuación ilustra la diferencia.

El código para reproducir la figura se muestra a continuación.

import scipy.interpolate as ip

#Define independent variable
sparse = linspace(0, 2 * pi, num = 20)
dense = linspace(0, 2 * pi, num = 200)

#Define function and calculate dependent variable
f = lambda x: sin(x) + 2
fsparse = f(sparse)
fdense = f(dense)

ax = subplot(2, 1, 1)

#Plot the sparse samples and the true function
plot(sparse, fsparse, label = ''Sparse samples'', linestyle = ''None'', marker = ''o'')
plot(dense, fdense, label = ''True function'')

#Plot the different interpolation results
interpolate = ip.InterpolatedUnivariateSpline(sparse, fsparse)
plot(dense, interpolate(dense), label = ''InterpolatedUnivariateSpline'', linewidth = 2)

smoothing = ip.UnivariateSpline(sparse, fsparse)
plot(dense, smoothing(dense), label = ''UnivariateSpline'', color = ''k'', linewidth = 2)

ip1d = ip.interp1d(sparse, fsparse, kind = ''cubic'')
plot(dense, ip1d(dense), label = ''interp1d'')

ylim(.9, 3.3)

legend(loc = ''upper right'', frameon = False)

ylabel(''f(x)'')

#Plot the fractional error
subplot(2, 1, 2, sharex = ax)

plot(dense, smoothing(dense) / fdense - 1, label = ''UnivariateSpline'')
plot(dense, interpolate(dense) / fdense - 1, label = ''InterpolatedUnivariateSpline'')
plot(dense, ip1d(dense) / fdense - 1, label = ''interp1d'')

ylabel(''Fractional error'')
xlabel(''x'')
ylim(-.1,.15)

legend(loc = ''upper left'', frameon = False)

tight_layout()