6

I'm having a bit of trouble with fitting a curve to some data, but can't work out where I am going wrong.

In the past I have done this with numpy.linalg.lstsq for exponential functions and scipy.optimize.curve_fit for sigmoid functions. This time I wished to create a script that would let me specify various functions, determine parameters and test their fit against the data. While doing this I noticed that Scipy leastsq and Numpy lstsq seem to provide different answers for the same set of data and the same function. The function is simply y = e^(l*x) and is constrained such that y=1 at x=0.

Excel trend line agrees with the Numpy lstsq result, but as Scipy leastsq is able to take any function, it would be good to work out what the problem is.

import scipy.optimize as optimize
import numpy as np
import matplotlib.pyplot as plt

## Sampled data
x = np.array([0, 14, 37, 975, 2013, 2095, 2147])
y = np.array([1.0, 0.764317544, 0.647136491, 0.070803763, 0.003630962,     0.001485394,     0.000495131])

# function
fp = lambda p, x: np.exp(p*x)

# error function
e = lambda p, x, y: (fp(p, x) - y)

# using scipy least squares
l1, s =  optimize.leastsq(e, -0.004, args=(x,y))
print l1
# [-0.0132281]


# using numpy least squares
l2 = np.linalg.lstsq(np.vstack([x, np.zeros(len(x))]).T,np.log(y))[0][0]
print l2
# -0.00313461628963 (same answer as Excel trend line)

# smooth x for plotting
x_ = np.arange(0, x[-1], 0.2)

plt.figure()
plt.plot(x, y, 'rx', x_, fp(l1, x_), 'b-', x_, fp(l2, x_), 'g-')
plt.show()

Edit - additional information

The MWE above includes a small sample of the dataset. When fitting the actual data the scipy.optimize.curve_fit curve presents an R^2 of 0.82, while the numpy.linalg.lstsq curve, which is the same as that calculated by Excel, has an R^2 of 0.41.

StacyR
  • 193
  • 1
  • 7

2 Answers2

4

You are minimizing different error functions.

When you use numpy.linalg.lstsq, the error function being minimized is

np.sum((np.log(y) - p * x)**2)

while scipy.optimize.leastsq minimizes the function

np.sum((y - np.exp(p * x))**2)

The first case requires a linear dependency between the dependent and independent variables, but the solution is known analitically, while the second can handle any dependency, but relies on an iterative method.

On a separate note, I cannot test it right now, but when using numpy.linalg.lstsq, I you don't need to vstack a row of zeros, the following works as well:

l2 = np.linalg.lstsq(x[:, None], np.log(y))[0][0]
Jaime
  • 65,696
  • 17
  • 124
  • 159
  • Thanks @Jaime - great answer! Unfortunately my maths knowledge is not that great; is one write or wrong [also see edit above], or are they just fundamentally different...? What are the implications for other functions, for example, if I wanted test the fit of a Sigmoid or Gompertz curve to the same data? – StacyR Jan 16 '13 at 21:14
  • @StacyR I don't have the knowledge to properly answer your question, but I am pretty sure that fitting an exponential as you did with `np.linalg.lstsq` is just a quick'n'dirty trick that doesn't compute errors properly. There's some discussion (hard for me to follow) here: http://mathworld.wolfram.com/LeastSquaresFittingExponential.html If you don't want to dive really deep into this stuff, I'd go with scipy's method for everything: it should give better fits, and your results will be consistent for all functions. – Jaime Jan 17 '13 at 21:26
  • thanks again! I have done some more research on this and, as you mentioned, have found that the `np.linalg.lstsq` method overly weights y-errors at low x values. The link you shared, and some other resources I found, allowed me to derive one other analytical method (the thing that makes it tricky is the constraint --- all the books describe the method for y=a*e^b*x rather than y=e^b*x), however, this also produces a worse fitting curve than the iterative `scipy.optimize.leastsq`. – StacyR Jan 18 '13 at 02:35
1

To expound a bit on Jaime's point, any non-linear transformation of the data will lead to a different error function and hence to different solutions. These will lead to different confidence intervals for the fitting parameters. So you have three possible criteria to use to make a decision: which error you want to minimize, which parameters you want more confidence in, and finally, if you are using the fitting to predict some value, which method yields less error in the interesting predicted value. Playing around a bit analytically and in Excel suggests that different kinds of noise in the data (e.g. if the noise function scales the amplitude, affects the time-constant or is additive) leads to different choices of solution.

I'll also add that while this trick "works" for exponential decay to 0, it can't be used in the more general (and common) case of damped exponentials (rising or falling) to values that cannot be assumed to be 0.

user3117404
  • 11
  • 2