Levenberg-Marquadt Algorithm (GSL implementation)

[pisuke]

Hello there,

I have to fit a one term exponential curve: y = a*exp(-b*t), and I have to implement the fitting procedure in a C program. Some days ago I discovered the GSL (GNU Scientific Library), and I told myself better later than never ;-)

When reading the GSL documentation and the implementation examples at http://www.gnu.org/software/gsl/man...-for-Nonlinear-Least_002dSquares-Fitting.html I gladly found that one of the examples is an exponential curve of the type y = a* exp(-b*t) + c, it means that with just a slight modification of the code I will get what I am looking for.

However, in the code there's something that I don't understand, the fitting procedure adds a term called sigma with a value of 0.1 in every calculation of the function, here is the explanation on the documentation:

The following example program fits a weighted exponential model with background to experimental data, Y = A \exp(-\lambda t) + b. The first part of the program sets up the functions expb_f and expb_df to calculate the model and its Jacobian. The appropriate fitting function is given by,

f_i = ((A \exp(-\lambda t_i) + b) - y_i)/\sigma_i

Here is the code doing the stuff:

...
sigma[i] = 0.1;
...

for (i = 0; i < n; i++)
         {
           /* Model Yi = A * exp(-lambda * i) + b */
           double t = i;
           double Yi = A * exp (-lambda * t) + b;
           gsl_vector_set (f, i, (Yi - y[i])/sigma[i]);
         }
...

I don't know the reason and what is more important I don't know if I have to use such term in my own program.

I am not a mathematician, I come from the biologic field and I need to fit data which comes from light absorption measurements.
Hopefully can someone here give me a hand on that issue.

Thanks,
Pisuke

 

[dodo]

While waiting for somebody more experienced for a proper answer, here's a snack:

You probably don't need the L-M algorithm at all; your equation can be fitted using logarithmic regression.
That is, given y = a* exp(-b*t), applying a logarithm to both sides, you get: ln y = ln a - b * t, which is a simple linear equation, which you can fit by the least squares method, for example.
Just use, instead of your (t,y) data pairs, the pairs (x', y') = (t, ln y), fit the linear equation y' = p + q * x', and then use b = -q and a = exp(p) in your original equation.

The "small" difference of a similar equation with an added constant +c in the L-M case example you cite, is due precisely because, after adding a constant, log regression is no longer applicable.

 

[pisuke]

Hello Dodo,

Thanks for your answer, but applying the natural logarithm to both sides of the equation does not solve my problem, what I have really is a double exponential decay of the form:

y = a*exp(-b*t) + c*exp(-d*t)

I approximate the double exponential to one exponential and I accept the deviation I get, because what I need is the factor, which multiplies the t in a single exponential decay.

y = a*exp(-b*t)

I have already tried your ln approach and I have compared the results with a direct LM fitting of the "deviated" one term exponential curve (I have used MATLAB fitting tool), the LM result is 10% better than the regression line result, that's why I need to fit directly the exponential decay.

Actually what would be ideal, would be a fitting procedure using the double exponential, and then somehow converge the parameters b and d in one, but that as long as I know is not possible.

Pisuke

 

[sk98]

As stated in the text that accompanies the example

The main part of the program sets up a Levenberg-Marquardt solver and some simulated random data. The data uses the known parameters (1.0,5.0,0.1) combined with gaussian noise (standard deviation = 0.1) over a range of 40 timesteps

So I think the extra sigma term in the example is the noise with a standard deviation of 0.1.

Hope it helps!

 

[lucbel]

hy,
i've read your question about the two exponential fitting of the function:

y = a*exp(-b*t) + c*exp(-d*t)


i've the same problem but i must implement the solving algorithm in excel

can you help me finding code example or other material?

you've write:

"I approximate the double exponential to one exponential and I accept the deviation I get, because what I need is the factor, which multiplies the t in a single exponential decay.

y = a*exp(-b*t)

I have already tried your ln approach and I have compared the results with a direct LM fitting of the "deviated" one term exponential curve (I have used MATLAB fitting tool), the LM result is 10% better than the regression line result, that's why I need to fit directly the exponential decay"

is it really a good approximation?

my email address is :
Lucabellesi@yahoo.it

best regards

 

[pisuke]

Hi!

hy,
i've read your question about the two exponential fitting of the function:

y = a*exp(-b*t) + c*exp(-d*t)


i've the same problem but i must implement the solving algorithm in excel

can you help me finding code example or other material?

In your case I would implement the algorithm using visual basic for applications. you can find the source code of a C implementation on the GSL (GNU Scientific Library): http://linux.duke.edu/~mstenner/free-docs/gsl-ref-1.1/gsl-ref_36.html

Unfortunately I have not experience in visual basic for applications, so I can't give you hints on translating the code to basic.

is it really a good approximation?

Buff... that depends on the function parameters, in my case it could be better but it could also be worse. However you should realize that excel, even offering a exponential fitting does not do an actual non linear fitting procedure. Excel does a linearization applying the natural logarithm to your data and afterwards a linear fitting procedure. And the results are radically different if your function is a complex exponential and not a single one.

I recommend you to play a little bit with Matlab's curve fitting tool (it has a fantastic help). Matlab also offers an Excel link tool, which allow you to link both programs, you can generate a function in matlab, which is a high level language very very flexible and then call that function from Excel. The MATLAB function can get as input the values which are in some cells of an excel sheet, and write the output in another cell of the same excel sheet. The shortcoming is that you need both programs, the advantage is that is easier than using visual basic, specially if you don't have previous experience.
Again the MATLAB help for the excel link tool is outstanding.

Regards,
Pisuke

 

[hotvette]

The 2nd attachment (Nonlinear Least Squares Example) in the below link might help.

https://www.physicsforums.com/showthread.php?t=97391

There are two ways this can be solved with Excel. The easiest is to set up the equations within the spreadsheet and solve using built-in Excel matrix functions. The second is to write VBA routine, which would still use the built-in Excel matrix functions.

 

[hotvette]

Attached is an example of how one might solve such a problem in Excel. The data I just made up. In this example, the matrix calculations are done with Excel built-in spreadsheet functions (i.e. MMULT, Transpose, Minverse) and the iterative portion is handled by some simple VBA code that does a copy/paste operation on the solution vector and adjusts the L-M parameter (Lambda) based on how the sum of sqaures is behaving. Of course, one can get much fancier by using matrix decomposition techniques such as QR instead of calculating a matrix inverse and using a line search to determine updated parameters instead of a fixed multiplication factor (mu), but additional VBA code would be needed. This particular example is tricky because it is ill conditioned and actually is unstable.