Get the formula of a interpolation function created by scipy

Question

I have done some work in Python, but I'm new to scipy. I'm trying to use the methods from the interpolate library to come up with a function that will approximate a set of data.

I've looked up some examples to get started, and could get the sample code below working in Python(x,y):

import numpy as np
from scipy.interpolate import interp1d, Rbf
import pylab as P

# show the plot (empty for now)
P.clf()
P.show()

# generate random input data
original_data = np.linspace(0, 1, 10)

# random noise to be added to the data
noise = (np.random.random(10)*2 - 1) * 1e-1

# calculate f(x)=sin(2*PI*x)+noise
f_original_data = np.sin(2 * np.pi * original_data) + noise

# create interpolator
rbf_interp = Rbf(original_data, f_original_data, function='gaussian')

# Create new sample data (for input), calculate f(x) 
#using different interpolation methods
new_sample_data = np.linspace(0, 1, 50)
rbf_new_sample_data    = rbf_interp(new_sample_data)

# draw all results to compare
P.plot(original_data, f_original_data, 'o', ms=6, label='f_original_data')
P.plot(new_sample_data, rbf_new_sample_data, label='Rbf interp')
P.legend()

The plot is displayed as follows:

Now, is there any way to get a polynomial expression representing the interpolated function created by Rbf (i.e. the method created as rbf_interp)?

Or, if this is not possible with Rbf, any suggestions using a different interpolation method, another library, or even a different tool are also welcome.

Answer 1

The answer is no, there is no "nice" way to write down the formula, or at least not in a short way. Some types of interpolations, like RBF and Loess, do not directly search for a parametric mathematical function to fit to the data and instead they calculate the value of each new data point separately as a function of the other points.

These interpolations are guaranteed to always give a good fit for your data (such as in your case), and the reason for this is that to describe them you need a very large number of parameters (basically all your data points). Think of it this way: you could interpolate linearly by connecting consecutive data points with straight lines. You could fit any data this way and then describe the function in a mathematical form, but it would take a large number of parameters (at least as many as the number of points). Actually what you are doing right now is pretty much a smoothed version of that.

If you want the formula to be short, this means you want to describe the data with a mathematical function that does not have many parameters (specifically the number of parameters should be much lower than the number of data points). Such examples are logistic functions, polynomial functions and even the sine function (that you used to generate the data). Obviously, if you know which function generated the data that will be the function you want to fit.

Answer 2

The RBF uses whatever functions you ask, it is of course a global model, so yes there is a function result, but of course its true that you will probably not like it since it is a sum over many gaussians. You got:

 rbf.nodes   # the factors for each of the RBF (probably gaussians)
 rbf.xi      # the centers.
 rbf.epsilon # the width of the gaussian, but remember that the Norm plays a role too

So with these things you can calculate the distances (with rbf.xi then pluggin the distances with the factors in rbf.nodes and rbf.epsilon into the gaussian (or whatever function you asked it to use). (You can check the python code of __call__ and _call_norm)

So you get something like sum(rbf.nodes[i] * gaussian(rbf.epsilon, sqrt((rbf.xi - center)**2)) for i, center in enumerate(rbf.nodes)) to give some funny half code/formula, the RBFs function is written in the documentation, but you can also check the python code.

Answer 3

RBF likely stands for Radial Basis Function. I wouldn't be surprised if scipy.interpolate.Rbf was the function you're looking for.

However, I doubt you'll be able to find a polynomial expression to represent your result.

If you want to try different interpolation methods, check the corresponding Scipy documentation, that gives link to RBF, splines...

Answer 4

I don’t think SciPy’s RBF will give you the actual function. But one thing that you could do is sample the function that SciPy’s RBF gave you (ie 100 points). Then use Lagrange interpretation with those points. This will generate a polynomial function for you. Here is an example on how this would look. If you do not want to use Lagrange interpolation, You can also use “Newton’s dividend difference method” to generate a polynomial function.

Answer 5

My answer is based on numpy only :

import matplotlib.pyplot as plt
import numpy as np
x_data = [324, 531, 806, 1152, 1576, 2081, 2672, 3285, 3979, 4736]
y_data = [20, 25, 30, 35, 40, 45, 50, 55, 60, 65]
x =  np.array(x_data)
y = np.array(y_data)
model = np.poly1d(np.polyfit(x, y, 2))
ynew = model(x)
plt.plot(x, y, 'o', x, ynew, '-' , )
plt.ylabel( str(model).strip())
plt.show()