2011-11-21 20 views
9

Estoy convirtiendo un programa MatLab en Python, y tengo problemas para entender por qué scipy.interpolate.interp1d está dando resultados diferentes que MatLab interp1.Los resultados de SciPy interp1d son diferentes que MatLab interp1

En MatLab el uso es ligeramente diferente:

yi = interp1(x,Y,xi,'cubic') 

SciPy:

f = interp1d(x,Y,kind='cubic') 
yi = f(xi) 

Para un ejemplo trivial de los resultados son los mismos: MatLab:

interp1([0 1 2 3 4], [0 1 2 3 4],[1.5 2.5 3.5],'cubic') 
    1.5000 2.5000 3.5000 

Python:

interp1d([1,2,3,4],[1,2,3,4],kind='cubic')([1.5,2.5,3.5]) 
    array([ 1.5, 2.5, 3.5]) 

embargo, para un ejemplo del mundo real no son los mismos:

x = 0.0000e+000 2.1333e+001 3.2000e+001 1.6000e+004 2.1333e+004 2.3994e+004 
Y = -6 -6 20 20 -6 -6 
xi = 0.00000 11.72161 23.44322 35.16484... (2048 data points) 

Matlab:

-6.0000e+000 
-1.2330e+001 
-3.7384e+000 
    ... 
7.0235e+000 
7.0028e+000 
6.9821e+000 

SciPy:

array([[ -6.00000000e+00], 
     [ -1.56304101e+01], 
     [ -2.04908267e+00], 
     ..., 
     [ 1.64475576e+05], 
     [ 8.28360759e+04], 
     [ -5.99999999e+00]]) 

Cualquier pensamiento en cuanto a cómo puedo obtener resultados que sean consistentes con MatLab?

Editar: Entiendo que haya una cierta flexibilidad en la implementación de los algoritmos de interpolación cúbica que probablemente explica las diferencias que estoy viendo. También parece que el programa MatLab original que estoy convirtiendo debería haber utilizado la interpolación lineal, por lo que la pregunta probablemente sea discutible.

Respuesta

11

El método de interpolación subyacente que scipy.interpolate.interp1d y interp1 son diferentes. Scipy usa las rutinas netlib fitpack, que producen splines cúbicos continuos estándar C2. El argumento "cúbico" en interp1 utiliza polinomios de interpolación de hermite cúbicos por partes, que no son C2 continuos. Vea here para una explicación de lo que hace Matlab.

Sospecho que es la fuente de la diferencia que está viendo.

0

En el actual scipy use http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.PchipInterpolator.html Esto creará la interpolación cúbica monotónica de la y = f (x) pasada, y utiliza el algoritmo pchip para determinar las pendientes en los puntos.

Por lo tanto, para cada sección, (x, y) es pasado por usted, el algoritmo pchip calculará (x, dy/dx), y solo hay cúbico pasando por 2 puntos con derivada conocida en estos puntos. Por construcción, continuará con la primera derivada continua.

Cuestiones relacionadas