forum.alglib.net
http://forum.alglib.net/

Double Exponential curve fitting
http://forum.alglib.net/viewtopic.php?f=2&t=800
Page 1 of 1

Author:  DarkDudae [ Wed Mar 20, 2013 1:45 pm ]
Post subject:  Double Exponential curve fitting

Hi there:

I am trying to solve a complex problem.

I have a Matrix with M values inside. Let?s call it: Comptes.

Matrix:
Channel Comptes
700 6
701 20
702 5
. .
. .
. .
800 7

M goes from 700 to 820. (length 120)

Well, I need that Summation from c = 700 to 820 of (Comptes_c - Comptes_Gauss_c)^2 = min

Where Comptes_Gauss = A*exp(-(c-Xk)^2/(2s^2))+B*exp(-(c-XL)^2/(2s^2))

(c = Channel)
s = sigma is a constant: 6.3

So I have to find the best suitable values for A, B, Xk and XL in order to minimize the given Summation.
I have initial approximate values to start using Min L/M iteration:
A=250
Xk=750
B=200
XL=775

I have tried to code this:

Code:
   Memo1.Clear;
    //We have 4 var: A, Xk, B, XLa
    N := 4;

    SetLength(S, N);

    //S: Our initial guess
    S[0]:=250; S[1]:=750; S[2]:=200; S[3]:=775;

    M := 120;
    SetLength(X1, M);
    SetLength(Y1, M);
    I:=0;
    while I<=M-1 do
    begin
        X1[I] := 700+I;
        //We read the 120 values from a memo (mValores)
        Y1[I] := strtofloat(mValores.Lines[i]);
        Inc(I);
    end;

    //
    // Now S stores starting point, X and Y store points being fitted.
    MinLMCreateFJ(N, M, S, State);
    MinLMSetCond(State, 0.0, 0.0, 0.001, 0);
    while MinLMIteration(State) do
    begin

        if State.NeedF then
        begin
            State.F := 0;
        end;
        I:=0;
        while I<=M-1 do
        begin

            //
            // "A" is stored in State.X[0]
            // "Xk" - State.X[1]
            // "B" - State.X[2]
            // "XLa" - State.X[3]
            //

            FI := Y1[I]-State.X[0]*exp((-1)*((AP_Sqr(X1[I]-State.X[1])/(2*AP_Sqr(6.3)))))-State.X[1]*exp((-1)*((AP_Sqr(X1[I]-State.X[3])/(2*AP_Sqr(6.3)))));

            if State.NeedF then
            begin
               //Summation here?
                State.F := State.F+AP_Sqr(FI);
            end;
            if State.NeedFiJ then
            begin
                //
                // Fi
                //
                State.Fi[I] := FI;

                //
                // dFi/dA = -exp^(-((x-Xk)^2)/2s^2)
                //
                State.J[I,0] := (-1)*exp((-1)*((AP_Sqr(X1[I]-State.X[1])/(2*AP_Sqr(6.3)))));
                //
                // dFi/dXk = (A*(xK-x)*exp(-(((x-Xk)^2)/(2*s^2))))))/(6.3^2))

                //
                State.J[I,1] := (State.X[0]*(State.X[1]-X1[I])*exp((-1)*((AP_Sqr(X1[I]-State.X[1])/(2*AP_Sqr(6.3))))))/(AP_Sqr(6.3));

                //
                // dFi/dB = -exp^(-((x-XL)^2)/2s^2)
                //
                State.J[I,2] := (-1)*exp((-1)*((AP_Sqr(X1[I]-State.X[3])/(2*AP_Sqr(6.3)))));

                //
                // dFi/dXLa = (B*(xL-x)*exp(-(((x-xL)^2)/(2*s^2))))))/(6.3^2))
                //
                State.J[I,3] := (State.X[2]*(State.X[3]-X1[I])*exp((-1)*((AP_Sqr(X1[I]-State.X[3])/(2*AP_Sqr(6.3))))))/(AP_Sqr(6.3));


            end;
            Inc(I);
        end;
    end;
    MinLMResults(State, S, Rep);

    //
    // output results
    //
    Memo1.Lines.Add(Format('A = %4.2f'#13#10'',[
        S[0]]));
    Memo1.Lines.Add(Format('Xk = %4.2f'#13#10'',[
        S[1]]));
    Memo1.Lines.Add(Format('B = %4.2f'#13#10'',[
        S[2]]));
    Memo1.Lines.Add(Format('XLa = %4.2f'#13#10'',[
        S[3]]));
    Memo1.Lines.Add(Format('TerminationType = %0d (should be 2 - stopping when step is small enough)'#13#10'',[
        Rep.TerminationType]));


Excel does a pretty good job with Solver and find the next values:
A=144,58
Xk=755,95
B=278,40
XL=772,68

I have uploaded the excel file I am using where you can see the real values and the used ecuations in a more graphic way.

However, when I use the Delphi implementation, I get wrong values for A, Xk, B, and XL. I don?t know what is wrong (maybe derivatives?)


Thank you in advance.

Attachments:
File comment: Excel file with text values and graphic ecuations
Channel_Values.xlsx [15.23 KiB]
Downloaded 355 times

Author:  DarkDudae [ Thu Mar 21, 2013 6:24 pm ]
Post subject:  Re: Double Exponential curve fitting

A partner has already found the mistake:

Code:
...   ))-State.X[1]*exp((-1)*((AP_S ...


Instead of "State.X[1]", it should have been "State.X[2]"

Thanks anyway!

Page 1 of 1 All times are UTC
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group
http://www.phpbb.com/