Averaging scattered data

Question

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.

Answer 1

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()

Answer 2

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.

Answer 3

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.