Take the 2-minute tour ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

I have multiple sets of measured data that can easily be visualized using a scatter plot (red and black points in the figure). If my measurements were perfect, the red and black points should lie on a single curve. However they are not, thus I would like to "average"/fit the data points in such a way that I obtain the black line in the figure. The black line could be a line or a set of discrete points with minimum distance to all the measurement data.

Can someone help me get started? Ideally there is already a python library that can handle this for me, but I'm also happy with algorithm descriptions.

Scattered data and averaged/fitted curve

share|improve this question
1  
do I understand it correctly that you want to have all red dots on the "same side" of the fit and the black dots on the other side? –  GertVdE Feb 26 '13 at 18:57
    
@GertVdE In case of two measurement sets, yes. For three measurement sets, imagine an extra set of green points with a similar topology. Then some of the points could "jump" sides. Is that clear? –  akid Feb 26 '13 at 19:52
1  
@GertVdE I removed my answer because it was only valid if one could assume the same number of points on both sides, which is not the case. –  Dr_Sam Feb 27 '13 at 13:44
1  
@Dr_Sam: ok, but removing it was a bit drastic. Your approach was valid for the case where you would have a certain number of measurements at the same $x$ values. You wouldn't necessarily need the same number of points above and below... –  GertVdE Feb 27 '13 at 13:47
1  
@GertVdE Ok, it was maybe rough to delete the post, but I did not want to mislead future readers of this question. I revived my answer just for the sake of helping this beta :-) –  Dr_Sam Feb 27 '13 at 14:03

3 Answers 3

up vote 4 down vote accepted

Based on your last comment (the fact that you might have multiple measurements and you don't care whether a certain set of measurements is above or below the fit), I think what you are looking for is a spline fit. You can do this using the scipy.interpolate B-spline routines. The script below generates three sets of data based on a function (that you can consider to model your system). For the first set, it adds some random error to the function "above" (red dots), for the second below (blue dots) and for the third just above and below (green dots). The black line is a plot of the exact function (without errors).

And then you use the scipy.interpolate.splrep function to generate a B-spline representation for this data. A B-spline is a piecewise combination of polynomial functions, typically cubic polynomials (cubic B-splines) which have some nice properties. If you want to know more about them, I would strongly suggest to read works by C. de Boor and P. Dierckx. The scipy.interpolate.splrep function returns a tuple of the knots, the coefficients and the order. Using the scipy.interpolate.splev function, you can then evaluate the B-spline representation in a number of points in the interval of the fit.

The script below plots an evaluation of the B-spline in magenta.

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

def f(x):
    return np.sin(np.pi*2*x)*np.exp(-2*x)

N=20
err=0.15

data = np.zeros([3*N,2])

#red dots above curve
data[:N,0] = np.sort(np.random.rand(N))
data[:N,1] = f(data[:N,0])+err*np.random.rand(N)

#blue dots below curve
data[N:2*N,0] = np.sort(np.random.rand(N))
data[N:2*N,1] = f(data[N:2*N,0])-err*np.random.rand(N)

#green dots above and below curve
data[2*N:,0] = np.sort(np.random.rand(N))
data[2*N:,1] = f(data[2*N:,0])-2*err*(np.random.rand(N)-0.5)

plt.plot(data[:N,0],data[:N,1],'ro')
plt.plot(data[N:2*N,0],data[N:2*N,1],'bo')
plt.plot(data[2*N:,0],data[2*N:,1],'go')

data = data[data[:,0].argsort()]
x = data[:,0]
y = data[:,1]
y_exact = f(x)


plt.plot(x,y_exact,'k')

w = np.ones([len(x),1])
spl = ip.splrep(x,y,w)

xn = np.linspace(0,1.0,100)
sple = ip.splev(xn,spl)

plt.plot(xn,sple,'m')

plt.show()
share|improve this answer
    
I'm missing SciPy on my workstation (which surprised me, because I do have NumPy), so I need to wait for IT to install it before I can test your solution. –  akid Feb 27 '13 at 19:38
    
@akid: ok. let me know if this is what you expected... –  GertVdE Feb 27 '13 at 20:16
    
Thanks, this seems to be what I need! –  akid Mar 9 '13 at 15:22

To reproduce your particular figure given the red and black labeling of the points, you could use a piece of a voronoi diagram. In other words, your black line will be the set of points that is equidistant between the nearest red and the nearest black point.

But assuming that you want to do something less contrived, you could use some kind of nonlinear svm instead. For example as shown in figures here.

share|improve this answer
    
SVM looks interesting, but @GertVDE's answer is simpler and also does the job I think. –  akid Mar 9 '13 at 15:20

A simple idea to obtain a continuous piecewise linear approximation would be to make the average of each pair of points (a red and the corresponding black) and make a piecewise linear interpolation of these averages. That would produce something close to what you show. (However, it does not make use a regression)

Would that fit your needs?

For the implementation, I think that you can do that manually, the method is very simple.

share|improve this answer
    
It is just by chance that I have an identical number of points for each set in the example. I might have 30 values in one set and 50 in the other, and 25 in a third, so the pairing is a bit of a problem. –  akid Feb 26 '13 at 17:11
    
Also, if you look at the central bump in the example sketch, I think you would end up with a different curve depending on how you pair... –  akid Feb 26 '13 at 17:12
    
Ok, this wasn't clear from the question. I'll try to update. –  Dr_Sam Feb 27 '13 at 7:08

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.