Exponential curve b*exp(-a*x) fit problem

french · **Joined:** Sat May 29, 2010 9:27 am **Posts:** 2

Hi,

I'm trying to fit exponential curve b*exp(-a*x) using LMfit. I adapted source code from example. First, I tried to fit it to set of noised data, with no success. After, I generated ideal exp function, but the result was the same. I was using lsfitnonlinearfgh and also lsfitnonlinearfg. I've changed fit to b*exp(a*x). I tried different combinations, but fitting doesn't work. Fit result is identical as starting values. When I use matlab curve fitting tool with LM algorithm (alglib use LM too), I've got pretty fine results (a = 0.003231, b = 610.1 - fit to data from attached signal.txt in range [350,1000]).
What am i doing wrong? I would be grateful for any help.

This is my source code:

Code:

using System;
using System.IO;
using alglib;

public class ExpFit
{
    public static int Main(string[] args)
    {
        lsfit.lsfitreport rep = new lsfit.lsfitreport();
        lsfit.lsfitstate state = new lsfit.lsfitstate();

        int begin = 350;
        int end = 1000;
        int n = end - begin;
        int m = 1;
        int k = 2;

        double a = 0.003231;
        double b = 610.1;

        //
        // Prepare task matrix
        //
        double[] y = new double[n];
        double[,] x = new double[n, m];
        double[] c = new double[k];

        double[] data = LoadFile("signal.txt");
        for(int i = 0; i < n; i++)
        {
            x[i,0] = i + begin;
            y[i] = data[i + begin];
            //y[i] = b * Math.Exp(-a * x[i,0]);
        }
        double epsf = 0.0;
        double epsx = 0.0;
        int maxits = 0;
        int info = 0;
        int j = 0;

        c[0] = 500;
        c[1] = 0.001;
        
        //
        // Solve
        //
        lsfit.lsfitnonlinearfgh(ref x, ref y, ref c, n, m, k, ref state);
        //lsfit.lsfitnonlinearfg(ref x, ref y, ref c, n, m, k, false, ref state);
        lsfit.lsfitnonlinearsetcond(ref state, epsf, epsx, maxits);
        while( lsfit.lsfitnonlineariteration(ref state) )
        {
            //
            // F(x) = b*exp(-a*x)
            //
            state.f = state.c[0] * Math.Exp(-(state.c[1] * state.x[0]));

            //
            // F(x)   = b*exp(-a*x)
            // dF/da  = -b*x*exp(-a*x)
            // dF/db  = exp(-a*x)
            //
            if (state.needfg | state.needfgh)
            {
                state.g[0] = -state.c[0] * state.x[0] * Math.Exp(-(state.c[1] * state.x[0]));
                state.g[1] = Math.Exp(-(state.c[1] * state.x[0]));
            }

            //
            // F(x)      = b*exp(-a*x)
            // d2F/da2   = b*x^2*exp(-a*x)
            // d2F/dadb  = -x*exp(a*x)
            // d2F/db2   = 0
            //
            if (state.needfgh)
            {
                state.h[0, 0] = state.c[0] * AP.Math.Sqr(state.x[0]) * Math.Exp(-(state.c[1] * state.x[0]));
                state.h[0, 1] = -state.x[0] * Math.Exp(-(state.c[1] * state.x[0]));
                state.h[1, 0] = -state.x[0] * Math.Exp(-(state.c[1] * state.x[0]));
                state.h[1, 1] = 0.0;
            }

            j++;
        }
        lsfit.lsfitnonlinearresults(ref state, ref info, ref c, ref rep);
        System.Console.Write("a:       ");
        System.Console.Write("{0,0:F8}",c[1]);
        System.Console.WriteLine();
        System.Console.Write("b:       ");
        System.Console.Write("{0,0:F8}",c[0]);
        System.Console.WriteLine();
        System.Console.Write("rms.err: ");
        System.Console.Write("{0,0:F3}",rep.rmserror);
        System.Console.WriteLine();
        System.Console.Write("max.err: ");
        System.Console.Write("{0,0:F3}",rep.maxerror);
        System.Console.WriteLine();
        System.Console.Write("iterations: ");
        System.Console.Write("{0}", j);
        System.Console.WriteLine();
        System.Console.Write("Termination type: ");
        System.Console.Write("{0,0:d}",info);
        System.Console.WriteLine();
        System.Console.WriteLine();
        System.Console.WriteLine();
        return 0;
    }

    public static double[] LoadFile(string path)
    {
        StreamReader sr = new StreamReader(path);
        uint n = uint.Parse(sr.ReadLine());
        uint dummy = uint.Parse(sr.ReadLine());
        double[] data = new double[n];
        for (uint i = 0; i < n; i++)
        {
            data[i] = double.Parse(sr.ReadLine());
        }
        sr.Close();
        return data;
    }
}

And console output:

Code:

a:       0.00100000
b:       500.00000000
rms.err: 176.992
max.err: 229.619
iterations: 34450
Termination type: 2

Sergey.Bochkanov · **Joined:** Fri May 07, 2010 7:06 am **Posts:** 927

Code:

//
// F(x) = b*exp(-a*x)
//
state.f = state.c[0] * Math.Exp(-(state.c[1] * state.x[0]));

//
// F(x)   = b*exp(-a*x)
// dF/da  = -b*x*exp(-a*x)
// dF/db  = exp(-a*x)
//
if (state.needfg | state.needfgh)
{
    state.g[0] = -state.c[0] * state.x[0] * Math.Exp(-(state.c[1] * state.x[0]));
    state.g[1] = Math.Exp(-(state.c[1] * state.x[0]));
}

//
// F(x)  = b*exp(-a*x)
// d2F/da2   = b*x^2*exp(-a*x)
// d2F/dadb  = -x*exp(a*x)
// d2F/db2   = 0
//
if (state.needfgh)
{
    state.h[0, 0] = state.c[0] * AP.Math.Sqr(state.x[0]) * Math.Exp(-(state.c[1] * state.x[0]));
    state.h[0, 1] = -state.x[0] * Math.Exp(-(state.c[1] * state.x[0]));
    state.h[1, 0] = -state.x[0] * Math.Exp(-(state.c[1] * state.x[0]));
    state.h[1, 1] = 0.0;
}

There is mistake in your derivative calculation code.

When you calculate F(x,a,b), you assume that c[0]=b, c[1]=a. But when you calculate grad(F) and hess(F), you treat c[0] as a, c[1] as b. So algorithm tries to search in a wrong direction and stops (because it is NOT descent direction).

french · **Joined:** Sat May 29, 2010 9:27 am **Posts:** 2

Thank you very much for your reply and pointing mistake.

serrrgio · **Joined:** Tue Nov 27, 2012 7:27 pm **Posts:** 1

I was very excited when I found libraries of computational algorithms for programs written in Delphi. It seemed to offer exactly what I needed for my research project.
I would like to fit experimental data by two-exponential function y=a*exp(-x/b)+c*exp(-x/d)+k
I attempted using function LSFitNonlinearFG (i.e. approximation without Gaussian), but encountered the following problems. Approximation fails even in case of rather smooth data without noise. When I ran my program without limiting the number of iterations, it continued for hours without returning any result. If I limit the number of iterations, then it stops after executing indicated number of steps, but found parameters of 2-exponential approximation do not fit the results.
I used your library version 2.6.0 for Delphi, since I don’t have experience with programming languages corresponding to the later version 3.4.0. I’m slightly familiar with C++, and therefore it might be another opportunity, if my Delphi program could not be made working with any earlier libraries.
All methods on your website refer to the library v. 3.4.0, and I couldn’t find any description for v.2.6.0 Does it still exist?
I kindly ask for your advice in this regard.
Could any other function utilizing limitations of the number of iterations improve multi-exponential approximation?
Do you believe that Gaussian methodology might be of advantage for this purpose?
Could your libraries v.3.40 for C++ be used in a program written in Delphi?

Here’s my code:

Code:

unit Unit1;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls, LsFit ,AP;

type
  TForm1 = class(TForm)
    Button1: TButton;
    Label1: TLabel;
    procedure Button1Click(Sender: TObject);
    procedure FormCreate(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
  end;

var
  Form1: TForm1;

implementation

{$R *.dfm}

procedure TForm1.Button1Click(Sender: TObject);
var
    M : AlglibInteger;
    N : AlglibInteger;
    K : AlglibInteger;
    Y : TReal1DArray;
    X : TReal2DArray;
    C : TReal1DArray;
    Rep : LSFitReport;
    State : LSFitState;
    Info : AlglibInteger;
    EpsF : Double;
    EpsX : Double;
    MaxIts : AlglibInteger;
    I : AlglibInteger;
    J : AlglibInteger;
    A : Double;
    B : Double;
begin
    Label1.Caption:='Ждите';
    Application.ProcessMessages;
    //
    // апроксимация 100exp(-x/10)+ 20exp(-x/200)+1 двумя exp:
    // * without Hessian (gradient only)
    // * using alpha=1 as initial value
    // * using 30 uniformly distributed points to fit to
    //
    // Notes:
    // * N - number of points
    // * M - dimension of space where points reside
    // * K - number of parameters being fitted
    //
    N := 30;
    M := 1;
    K := 5;
    A := 0.1;
    B := 500;
    
    //
    // Подготовка матрицы задачи
    //
    SetLength(Y, N);
    SetLength(X, N, M);
    SetLength(C, K);
    I:=0;
    while I<=N-1 do
    begin
        X[I,0] := A+(B-A)*I/(N-1);
        Y[I] := 100*exp(-X[I,0]/30)+ 2*exp(-X[I,0]/200)+1;
        Inc(I);
    end;
    C[0] := 2.0;
    C[1] := 1.0;
    C[2] := 1.0;
    C[3] := 2.0;
    C[4] := 1.0;
    EpsF := 0;
    EpsX := 0;
    MaxIts := 21474;//83647;
    
    //
    // Решение
    //
    LSFitNonlinearFG(X, Y, C, N, M, K, True, State);
    LSFitNonlinearSetCond(State, EpsF, EpsX, MaxIts);
    LSFitNonlinearSetStpMax(State, 0.3);
    while LSFitNonlinearIteration(State) do
    begin
        if State.NeedF then
        begin
            
            //
            // F(x) = a*exp(-x/b)+ cexp(-x/d)+k
            //
            State.F := C[0]*Exp(-State.X[0]/State.C[1])+C[2]*Exp(-State.X[0]/State.C[3])+C[4];
        end;
        if State.NeedFG then
        begin
            
            //
            // F(x)  = a*exp(-x/b)+ cexp(-x/d)+k
            // dF/da = exp(-x/b)
            // dF/db = a*x*exp(-x/b)/b^2
            // dF/dc = exp(-x/d)
            // dF/dd = c*x*exp(-x/d)/d^2
            // dF/dk = 1
            State.F := C[0]*Exp(-State.X[0]/State.C[1])+C[2]*Exp(-State.X[0]/State.C[3])+C[4];
            State.G[0] := Exp(-State.X[0]/State.C[1]);
            State.G[1] := C[0]*State.X[0]*Exp(-State.X[0]/State.C[1])/(State.C[1]*State.C[1]);
            State.G[2] := Exp(-State.X[0]/State.C[2]);
            State.G[3] := C[2]*State.X[0]*Exp(-State.X[0]/State.C[3])/(State.C[3]*State.C[3]);
            State.G[4] := 1;
        end;
    end;
  //Вывод результатов
    LSFitNonlinearResults(State, Info, C, Rep);
    Label1.Caption:='a:'+FloatToStr(C[0])+#13+#10
    +'b:'+FloatToStr(C[1])+#13+#10
    +'c:'+FloatToStr(C[2])+#13+#10
    +'d:'+FloatToStr(C[3])+#13+#10
    +'k:'+FloatToStr(C[4])+#13+#10
    +'rms.err:'+FloatToStr(Rep.RMSError)+#13+#10
    +'max.err:'+FloatToStr(Rep.MaxError)+#13+#10
    +'Termination type:'+FloatToStr(Info);
end;                              


procedure TForm1.FormCreate(Sender: TObject);
begin
 Label1.Caption:='Апроксимация 100exp(-x/30)+ 2exp(-x/200)+1 на отрезке [0.1,500] двумя exp';
end;

end.

and sample of my fitted data by other programs

And output:

Sergey.Bochkanov · **Joined:** Fri May 07, 2010 7:06 am **Posts:** 927

You have an error in your function/gradient calculation code - sometimes you use State.C, and sometimes just C (local array, not fields of State structure). This local array is not updated with recent changes in parameters, so algorithm is confused and tries to make wrong steps.

Sergey.Bochkanov · **Joined:** Fri May 07, 2010 7:06 am **Posts:** 927

BTW, I recommend you to increase StpMax to 10-20.... Your algorithm has to perform too many small steps in order to reach optimum parameters.

forum.alglib.net

Forum rules

Exponential curve bexp(-ax) fit problem

Who is online