Interpolation in SciPy: Finding X that produces Y

Question

Is there a better way to find which X gives me the Y I am looking for in SciPy? I just began using SciPy and I am not too familiar with each function.

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

x = [70, 80, 90, 100, 110]
y = [49.7, 80.6, 122.5, 153.8, 163.0]
tck = interpolate.splrep(x,y,s=0)
xnew = np.arange(70,111,1)
ynew = interpolate.splev(xnew,tck,der=0)
plt.plot(x,y,'x',xnew,ynew)
plt.show()
t,c,k=tck
yToFind = 140
print interpolate.sproot((t,c-yToFind,k)) #Lowers the spline at the abscissa

Furbish

Can you elaborate on what you want to be better? Performance, accuracy, conciseness? – Shattuck

Answer 1

The UnivariateSpline class in scipy makes doing splines much more pythonic.

x = [70, 80, 90, 100, 110]
y = [49.7, 80.6, 122.5, 153.8, 163.0]
f = interpolate.UnivariateSpline(x, y, s=0)
xnew = np.arange(70,111,1)

plt.plot(x,y,'x',xnew,f(xnew))

To find x at y then do:

yToFind = 140
yreduced = np.array(y) - yToFind
freduced = interpolate.UnivariateSpline(x, yreduced, s=0)
freduced.roots()

I thought interpolating x in terms of y might work but it takes a somewhat different route. It might be closer with more points.

Alisonalissa

Wouldn't this require twice the amount of CPU calculations since are interpolating practically the same data set two times? – Furbish

@JcMaco, the first use of UnivariateSpline is just to make a pretty plot. The second usage is what actually gives the values. – Harrietharriett

Craig is right, can you correct it on your example as it's otherwise great! – Grevera

Fixed the typo. Thanks Craig. – Alisonalissa

Answer 2

If all you need is linear interpolation, you could use the interp function in numpy.

Cootch

I'd prefer spline interpolation. How would the interp function help me solve my problem more easily? – Furbish

Your question didn't specify what type of interpolation you needed -- if linear isn't good enough for your problem, then I don't think interp will help. – Cootch

Answer 3

I may have misunderstood your question, if so I'm sorry. I don't think you need to use SciPy. NumPy has a least squares function.

#!/usr/bin/env python

from numpy.linalg.linalg import lstsq



def find_coefficients(data, exponents):
    X = tuple((tuple((pow(x,p) for p in exponents)) for (x,y) in data))
    y = tuple(((y) for (x,y) in data))
    x, resids, rank, s = lstsq(X,y)
    return x

if __name__ == "__main__":
    data = tuple((
        (1.47, 52.21),
        (1.50, 53.12),
        (1.52, 54.48),
        (1.55, 55.84),
        (1.57, 57.20),
        (1.60, 58.57),
        (1.63, 59.93),
        (1.65, 61.29),
        (1.68, 63.11),
        (1.70, 64.47),
        (1.73, 66.28),
        (1.75, 68.10),
        (1.78, 69.92),
        (1.80, 72.19),
        (1.83, 74.46)
    ))
    print find_coefficients(data, range(3))

This will return [ 128.81280358 -143.16202286 61.96032544].

>>> x=1.47 # the first of the input data
>>> 128.81280358 + -143.16202286*x + 61.96032544*(x**2)
52.254697219095988

0.04 out, not bad

Mundy

The problem is finding which number will give me 52.21. Of course, there can be many solutions if the interpolation is quadratic (or higher power). – Furbish