Python: generating a "curve fit score"
Question: I am working on a project in which I am trying to model the movement of an
object in a [kymograph](http://en.wikipedia.org/wiki/Kymograph). In order to
do so, I fit a curve to each line of pixels in an image, and append the
location of the vertex to approximately model the location of the object in
the image. Below is a sample image.

![Image to be analyzed](http://i.stack.imgur.com/D1CyB.png)

As you can see, early in the time series (at the top of the image) the
position of the object is nicely focused and easily modeled with a Gaussian
curve. However, closer to the end of the time series (at the bottom of the
image), the peak is much more diffuse. I suspect that the data at the bottom
of the image will be fit much more closely by a curve modeling a Poisson
distribution (image below, right) while the data at the top/middle of the
image will be fit much more closely by a Gaussian or polynomial curve (image
below, left).

![good and bad examples](http://i.stack.imgur.com/SKu1a.png)

Is there any way to, for each line of pixels, fit more than one curve to the
same data and then score each for a least-squares fit? This way, I could
(hopefully) switch models midway through an image to accommodate changing
behaviors of the object that I am trying to track. My current code is below:

    
    
    from PIL import Image
    
    def populateData(picture) :
        """Open an image and populate a list of lists with the grayscale value"""
        im = Image.open(picture)           
        size = im.size                    
        width = size[0]                   
        height = size[1]                   
        allPixels = list(im.getdata())      
        pixelList = [allPixels[width*i :    
                           width * (i+1)] for i in range(height)] 
        return(pixelList)
    
    rawData = populateData("testTop.tif")   
    
    import numpy as np
    from scipy.optimize import curve_fit
    
    def findVertex(listOfRows) :
    
        xList = []
        for row in listOfRows :
            x = np.arange(len(row))
            ffunc = lambda x, a, x0, s: a*np.exp(-0.5*(x-x0)**2/s**2) 
            p, _ = curve_fit(ffunc, x, row, p0=[100,5,2])
            x0 = p[1] 
            xList.append(x0)
        xArray = np.array(xList)    
        return(xArray)
    
    xValues = findVertex(rawData)
    
    def buildRows(listOfRows) :
        yArray = np.arange(len(listOfRows))
        return(yArray)
    
    yValues = buildRows(rawData)
    
    from matplotlib import pyplot as plt
    from scipy import ndimage
    
    image = ndimage.imread("testTop.tif",flatten=True)
    fig = plt.figure()
    axes = fig.add_subplot(111)
    axes.imshow(image)
    axes.plot(xValues, yValues, 'k-')
    
    axes.set_title('testLine')
    axes.grid()
    axes.set_xlabel('x')
    axes.set_ylabel('time')
    plt.show()
    

![program output](http://i.stack.imgur.com/La2sc.png)

EDIT: This is the file I used as an input (testTop.tif) ![input
image](http://i.stack.imgur.com/Gh5dB.png)


Answer: You will need to work out some form of [goodness of
fit](http://en.wikipedia.org/wiki/Goodness_of_fit) between the fit and your
data. Taking the sum of the squared differences between your current fit (a
Gaussian) and your data divided by the variance.

    
    
    sumerrsq = 0.
    for i in range(yValues.shape[0]):
        sumerrsq += np.power(yValues[i] - xValues[i],2)
    
    goodfit = np.sqrt(sumerrsq/var)
    

I think you can use use the second output from curve fit (the covariance) to
get the variance,

    
    
    p, pcov = curve_fit(ffunc, x, row, p0=[100,5,2])
    var = np.diag(pcov)
    

You can then check the value of goodfit and if it is not sufficient, switch to
a different distribution. In using a different distribution, you may need to
use a different estimation of error (this assumes the errors are normally
distributed).

Note, without the data (and not being sure what array was which) I couldn't
test any of this code...

