I have a very large data set (roughly 11,000 points) that follow a roughly exponential curve with random variation. Here's a small sample of my data set:
The underlying exponential function follows a curve a*b^x, not a*e^(b*x). In the case of the graph above, the equation is something like 0.05 * 1.195^x (I don't know the exact values)
My aim is to find the two parameters (a and b) that best fit the data. If the equation were of the form a*e^(b*x) then I could use standard linear regression techniques, but that's not the case.
So far I've taken the approach of computing the sum of squared errors (for all 11,000 data points) and attempting to minimize this error function. I've had varying degrees of success with this:
- Using Excel I can set up two "input" cells to be my parameters
aandb, add a column for "computed value" beside my data, add a column for "squared error", then add an "output" cell that sums this squared error column. I then use Solver to minimize the output cell by changing the input cells and it performs very well. On my raw data set I get a sum squared error around 48,000 - Scipy's curve_fit utterly failed when passed my data set, giving a nonsensical answer. Scipy's minimize equally failed when passed an error function that computed the sum of squared errors but gave a meaningful message that overflow was encountered and so the desired error was not achieved due to precision loss. The error was in the range of 1e147
- OpenOffice also has a Solver, like Excel, but it utterly failed and generated nonsense answers. The error was in the range of 1e9
- Using both the genetic optimization and Powell optimization from optimization-js returned nonsense answers. The error was in the range of 1e9. I wasn't sure if I could take advantage of the L-BFGS or gradient descent options since I don't know how to compute the derivative for my error function
- I tried to use liboptim but couldn't figure out how to get all of the dependencies working (Armadillo, LAPACK, etc)
- I wrote my own very naive nonlinear solver which uses a pseudo-binary search for the first parameter and steps through all possible digits for the second parameter, stopping when it hits an inflection point. This seemed to work pretty well. It returned an error of 55,000 (not the 48,000 Excel got, but far better than I could do by hand)
In trying to research the mathematics involved in non-linear optimization so I can improve my naive optimizer, I keep stumbling upon recurring terms that I don't fully understand. Like whether a function is "Lipschitz", or whether the function is "convex".
My first question is: Given the definition of my error function (the sum of squared differences between an exponential function and a data set that is roughly exponential), what properties would my function have? Is it convex? Is it "Lipschitz"?
My second question is: Am I overdoing this? Is there an easier solution?
Best Answer
Your model is a pure exponential $$y=a\, b^x=a\, e^{x\log(b)}=a\, e^{cx}$$ but it is nonlinear with respect to its parameters; so you need some reasonable guesses to start.
Keeping your formulation, in a first step, linearize the model $$y=a\, b^x \implies \log(y)=\log(a)+x \log(b)=\alpha + \beta x$$ A first linear regression gives $\alpha$ and $\beta$ and then $a=e^{\alpha}$ and $b=e^{\beta}$. Now, start the nonlinear regression.
Edit
You could even reduce the problem to one equation in $b$ $$a=\frac{\sum_{i=1}^n y_i b^{x_i} } {\sum_{i=1}^n b^{2x_i} }$$ and then $$f(b)=\frac{\sum_{i=1}^n y_i b^{x_i} } {\sum_{i=1}^n b^{2x_i} }-\frac{\sum_{i=1}^n x_iy_i b^{x_i} } {\sum_{i=1}^n x_ib^{2x_i} }=0$$ SInce you have the estimate from $\beta$, even plotting will give you the result