forum.alglib.net

ALGLIB forum
It is currently Sat Nov 23, 2024 3:11 pm

All times are UTC


Forum rules


1. This forum can be used for discussion of both ALGLIB-related and general numerical analysis questions
2. This forum is English-only - postings in other languages will be removed.



Post new topic Reply to topic  [ 2 posts ] 
Author Message
 Post subject: Double Exponential curve fitting
PostPosted: Wed Mar 20, 2013 1:45 pm 
Offline

Joined: Wed Mar 20, 2013 11:21 am
Posts: 4
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 554 times
Top
 Profile  
 
 Post subject: Re: Double Exponential curve fitting
PostPosted: Thu Mar 21, 2013 6:24 pm 
Offline

Joined: Wed Mar 20, 2013 11:21 am
Posts: 4
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!


Top
 Profile  
 
Display posts from previous:  Sort by  
Post new topic Reply to topic  [ 2 posts ] 

All times are UTC


Who is online

Users browsing this forum: No registered users and 2 guests


You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot post attachments in this forum

Search for:
Jump to:  
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group