Recent

Author Topic: Lmath ulogifit samples  (Read 1753 times)

erik2727

  • New Member
  • *
  • Posts: 14
Lmath ulogifit samples
« on: September 03, 2020, 02:58:21 am »
Hi,

I was looking for some samples on how to use ulogifit for 4 parameters logistic curve fit,
Could someone share some tips or guide please? Below are my code snippet, not sure if the TVector B parameter result output is correct or not ?


var

 X:array [0..9] of Float=(0,1,2,3,4,5,6,7,8,9,10);
 Y:array [0..9] of Float=(10,11,25,36,82,90,106,201,301,403);
 TX,TY:TVector;
  vecB: array of Float;
  B: TVector;
  MatV: TMatrix;
  beta,alpha:Float;
  A1,B1,N:Float;

begin

 TX:=X;
TY:=Y;
DimVector(TX,10);
DimVector(TY,10);
  SetLength(vecB,5);
  B := vecB;
 
  DimMatrix(MatV, 0, 4);
  LogiFit(TX, TY, 1, 9,true,true, 1000, 0.001, B, MatV); 

A1:=B[0];
B1:=B[1];
alpha:=B[2];
beta:=B[3];
N:=B[4];
« Last Edit: September 05, 2020, 04:12:06 am by erik2727 »

Paolo

  • Jr. Member
  • **
  • Posts: 68
Re: Lmath ulogifit samples
« Reply #1 on: September 13, 2020, 01:06:07 pm »
I don't know about LogiFit, but  I see something wrong in your code:

1) you declare X : array [0..9], this means that you have 10 elements in the array, but you initialize 11 elements ! from 0 to 10, I think this should not compile !
2) why declare vecB : arrray of float and not as TVector !
3) by using DimMatrix(MatV, 0, 4) are you saying that the matrix has one dimension = 0 ! I think also this should not complie, it should be DimMAtrix(MatV, 5, 5) according to matrix 5x5 = [0..4,0..4]. (Maybe the martix is just an output so it is the LogiFit that domesnions it properly, I don't know.

erik2727

  • New Member
  • *
  • Posts: 14
Re: Lmath ulogifit samples
« Reply #2 on: September 17, 2020, 06:07:42 am »
Hi,

If you don't know lmath especially logistics mathematical functions in detail, i suggest you don't simply give your
 comments, pascal blueprint are based on Fortran!


First X element should be 1 not zero, rest of the
code compile successfully, readup Lmath what is underlying principle and basis on model construction, I has compared the result with famous C++ non linear solver just like to know why this algorithm assume highest k carrying capacity is the actual high end number and deduce  rest of parameters bases on fixed preset as I have no access to Original Fortran 90 solver
« Last Edit: September 17, 2020, 11:59:06 am by erik2727 »

Paolo

  • Jr. Member
  • **
  • Posts: 68
Re: Lmath ulogifit samples
« Reply #3 on: September 21, 2020, 01:51:50 pm »
Hi, I admit to have replaied too hasty (in the attempt to help in someway).
Bye.

alanphys

  • New Member
  • *
  • Posts: 46
Re: Lmath ulogifit samples
« Reply #4 on: September 21, 2020, 03:02:13 pm »
Hi

Have a look at the regexlin.pas demo in the demo/bgi/regmodel directory of LMath.

Regards
Alanphys
Fedora 32 + KDE/QT5, Tinycore 8, Windows XP-10
https://github.com/alanphys

wp

  • Hero Member
  • *****
  • Posts: 7739
Re: Lmath ulogifit samples
« Reply #5 on: September 21, 2020, 04:55:59 pm »
I was looking for some samples on how to use ulogifit for 4 parameters logistic curve fit,
Could someone share some tips or guide please?
Enclosed is a demo project performing a logistic fit to synthesized data. In the "Input parameters" box the coefficients of the logistic function can be specified for creation of data. In the "Start parameters" you must provide estimates for these coefficients with which the fit should start, they should not be off too much, otherwise numerical overflows can occur depending on the choice of the fitting algorithm (I found that the "Simplex algorithm" seems to be more tolerant than the default "Marquardt algorithm"). Click "Create data" to create the data to be fitted, then "Calc fit" to execute the fit. The fitted curve is drawn over the input data, and the fitted parameters are listed.

var
 X:array [0..9] of Float=(0,1,2,3,4,5,6,7,8,9,10);
 Y:array [0..9] of Float=(10,11,25,36,82,90,106,201,301,403);
...
DimVector(TX,10);
DimVector(TY,10);
... 
  DimMatrix(MatV, 0, 4);
  LogiFit(TX, TY, 1, 9,true,true, 1000, 0.001, B, MatV); 
I agree with Paolo that there is a problem with the indices. The 3rd and 4th parameters of the LogiFit procedure are the start and end index of the data array provided as first and second parameter. Since DimVector() internally simply calls Setlength() which creates a 0-based array the 3rd LogiFit parameter must be 0 rather than 1; the 4th parameter (last array index) is correct. And the matrix MatV is an m x m matrix where m is the count of fit parameters (5, here); so, the DimMatrix() call must be DimMatrix(MatV, 5, 5); each index i,j in MatV[i,j] then runs from 0 to 4
« Last Edit: September 21, 2020, 05:07:58 pm by wp »
Mainly Lazarus trunk / fpc 3.2.0 / all 32-bit on Win-10, but many more...

erik2727

  • New Member
  • *
  • Posts: 14
Re: Lmath ulogifit samples
« Reply #6 on: September 22, 2020, 03:17:53 pm »
Hi,

Thank you so much  to those who have responded, I  pursuit another non linear fit  solution to my test data with the following code snippets  using modify logistic growth model equation, the output obtained was comparable and almost identical to that of the well known C++ nonlinear solver yield  , the most challenging task of using this method is the proper coding of  partial derivative apply to function parameters required by  Jacobian matrix as well as the best guess input estimate for the parameters   



function RegFunc(X : Float; B1 : TVector) : Float;
{B1[1]=k,B1[2]=C,B1[3]=b }
{y=k/1+(X/C)^b}

begin
  RegFunc := B1[1] / (1+Math.Power(X/B1[2],B1[3]));
end;

procedure DerivProc(X, Y : Float; B1, D : TVector); {Jacobian Matrix}
{B1[1]=k,B1[2]=C,B1[3]=b }

begin
  D[1]:=1/ (1+Math.Power(X/B1[2],B1[3])); {pdk}
  D[2]:=B1[1]*Math.Power(X,B1[3])*B1[3]*Math.Power(B1[2],(-B1[3]-1))/Math.Power(1+Math.Power(X/B1[2],B1[3]),2);  {pdc}
  D[3]:=-B1[1]*Math.Power((X/B1[2]),B1[3])*Ln(X/B1[2])/Math.Power(1+Math.Power(X/B1[2],B1[3]),2); {pdb}

end;


procedure Multi_Fit();


  FirstPar = 1;       { Index of first fitted parameter }
  LastPar  = 3;       { Index of last fitted parameter }
  MaxIter  = 1000;    { Max. number of iterations }
  Tol      = 1.0E-3;  { Required precision }
  Alpha    = 0.05;    { Significance level }


  var
    B1:TVector;
    V:TMatrix;
    N1:integer ;      { Number of points }

begin
  DimVector(B1,LastPar);
  DimMatrix(V, LastPar, LastPar);
  N1:= ??    {set the number of data point to be fitted}
  // set algorithm
  SetOptAlgo(NL_SIMP);

  //initial estimate input to solver
  B1[1]:= ???
  B1[2]:= ???
  B1[3]:= ???


  NLFit(@RegFunc, @DerivProc, TX, TY, 1, N1, MaxIter, Tol,
          B1, FirstPar, LastPar, V);

   ShowMessage('B1[0..N] = ' + FloatToStr(B1[1]) + '/' + FloatToStr(B1[2]) +
    '/' + FloatToStr(B1[3]));


end;
« Last Edit: September 24, 2020, 10:32:03 am by erik2727 »

wp

  • Hero Member
  • *****
  • Posts: 7739
Re: Lmath ulogifit samples
« Reply #7 on: September 22, 2020, 04:35:44 pm »
Is this a question?

I doubt that your procedure is correct. I must say that it's been a long time that I had been able to calculate partial derivatives by myself. Today I call Wolfram Alpha and enter the expression to be differentiated: "derivative  of <expression>". Your expression is
Code: [Select]
y := k / (1+(x/c)^b)Since Wolfram Alpha only calculates a derivative with respect to a variable named x I do a replacement: When I want the partial derivative with respect to b I replace your b by x and your x by a new variable A:
Code: [Select]
derivative of k / (1 + (A/c)^x)The result is
Code: [Select]
dy/dx = -k * ((A/C)^x * ln(A/C)) / ( (A/C)^x + 1)^2or, with the original variables
Code: [Select]
dy/db = -k * (x/C)^b * ln(x/C)) / ( (x/C)^b + 1)^2Part of it is looking like your pdb =  D[3], but yours is lacking the denominator...

For calculating the derivative with respect to C I enter
Code: [Select]
derivative of k / (1 + (A/x)^b)and get the result
Code: [Select]
dy/dx = b * k *  (A/x)^b / ( x * ((A/x)^b - 1)^2 )or, with our variables
Code: [Select]
dy/dC = b * k * (x/C)^b / (C * ((x/C)^b - 1)^2)Again, this is different from your equation for pdc = D{2].

I don't insist on that my results are correct either, there may be typos everywhere... I just want to say: please check your result carefully.

Another (possible) issue refers to the indices, again. This is the declaration of LMath for DimVector (taken from unit utypes):
Code: Pascal  [Select][+][-]
  1. { Creates floating point vector V[0..Ub] }
  2. procedure DimVector(var V : TVector; Ub : Integer); overload;
As the comment says, the procedure creates a 0-based vector, the first element has index 0, the last element has index Ub. So, your function requires 3 arguments, k, C, and b. Since you call DimVector(B, lastPar), with lastPar 3, you are creating a vector with 4 elements, not three, i.e. you ignore the 1st element with index 0. As long as you create more elements then needed there is no problem with this (except for the wasted memory). But you really must keep in mind that this practice is against the Lazarus/FPC  convention. Almost all arrays (except for strings) begin with index 0 and expect the array element with this index to contain valid information. Therefore, I'd expect that you will sooner or later run into problems when stuffing 1-based arrays into 0-based memory allocation.

Suppose you want to display the fit parameters in a listbox. For this purpose you iterate through your pseudo-1-based array B and add the numbers as strings to the listbox:

Code: Pascal  [Select][+][-]
  1. for i := 1 to lastPar do
  2.   Listbox.Items.Add(FormatFloat('0.00000', B[i]);

Now you want to display the second parameter in a label; since the parameter is already converted to string you take the copy from the list box:

Code: Pascal  [Select][+][-]
  1.   Label1.Caption := Listbox.Items[2];

No! The items in the listbox are 0-based, but you think 1-based, you must use ListBox.Items[1]! Whenever I mixed 0-based and 1-based storage I ran into lots of trivial, but hard-to-find problems.

Therefore, I would like you to reconsider this practice and forget your obvious FORTRAN roots. The Pascal-compatible way would be to really use the 0-based array:
Code: Pascal  [Select][+][-]
  1. const
  2.   lastPar = 2;
  3. ...
  4.   DimVector(B, lastPar);
  5. ...
  6.   B[0] := k;
  7.   B[1] := C;
  8.   B[2] := b;  // of course, this will not compile...
  9. ...
  10.   D[0] := ... {pdk}
  11.   D[1] := ... {pdc}
  12.   D[2] := ... {pdb}
« Last Edit: September 22, 2020, 04:48:20 pm by wp »
Mainly Lazarus trunk / fpc 3.2.0 / all 32-bit on Win-10, but many more...

erik2727

  • New Member
  • *
  • Posts: 14
Re: Lmath ulogifit samples
« Reply #8 on: September 22, 2020, 06:03:08 pm »
Hi ,

I have vetted my result with few defacto  c++ nonlinear solver of which one is the GNU GSL library, the final result is what matter most, have I done  a wrong derivative somehow ,  I may not be able to obtain a default optimum solution , by the way , my main programming language root is C++!!! and I am not a rookie, the reason why I put up this is Lmath or Dmath manual are not very comprehensive and confusing in presenting the solution. 

As long as the result is vetted independently  , , I don't have time to   verify your so called computed derivative, I don't see why computed dy/dk still got k variable when
k*(Constant) is just  constant as far as simple Jacobain matrix of ith row is concern  .

ps.  missing denominator added to orginal

Matlab diff result as below:=

f=k/((x^b/c^b) + 1)

∂f/∂ k=
1/(x^b/c^b + 1)

∂f/∂c=
(b*k*x^b)/(c^(b + 1)*(x^b/c^b + 1)^2)


∂f/∂b=
-(k*((x^b*log(x))/c^b - (x^b*log(c))/c^b))/(x^b/c^b + 1)^2

      k
f=  ──────────
         -b   b
  1 + c  ⋅x

∂f/∂k=


         1
 ──────────
         -b   b
  1 + c  ⋅ x

∂f/∂c=


        -b      b
     b⋅c  ⋅k ⋅x
  ───────────────
                      2
     ⎛       -b  b  ⎞
  c⋅⎝1 + c  ⋅x    ⎠


∂f/∂b=

     ⎛ -b  b            -b  b          ⎞
  k⋅⎝c  ⋅x ⋅log(c) - c  ⋅x ⋅log(x) ⎠
  ─────────────────────────────────
                             2
            ⎛      -b   b  ⎞
            ⎝1 + c  ⋅x    ⎠




 I also could not understand why someone  just like keep on commenting on the index array LB.. UB , read what LMath manual was stating ,

 TMatrix  := total number of independent variable
0..2 is two elements  or 3 elements in total ???????

parameters; V = inverse matrix. Its dimentions must be [Lb..Ub, Lb..Ub],
and also look into the sample programme exhibit ,


Extract of the Lmath sample from regnlin.pas

function RegFunc(X : Float; B : TVector) : Float;
begin
  RegFunc := B[1] * Exp(- B[2] * X);
end;

procedure DerivProc(X, Y : Float; B, D : TVector);
begin
  D[1] := Exp(- B[2] * X);
  D[2] := - B[1] * X * D[1];
end;

 
« Last Edit: September 24, 2020, 10:45:48 am by erik2727 »

glorfin

  • Full Member
  • ***
  • Posts: 109
Re: Lmath ulogifit samples
« Reply #9 on: September 22, 2020, 06:29:23 pm »
Hi

Have a look at the regexlin.pas demo in the demo/bgi/regmodel directory of LMath.

Regards
By the way, I have committed BGI demo projects which compile and work under Windows 64. See Trunk.

glorfin

  • Full Member
  • ***
  • Posts: 109
Re: Lmath ulogifit samples
« Reply #10 on: September 22, 2020, 06:38:38 pm »
I agree with Paolo that there is a problem with the indices. The 3rd and 4th parameters of the LogiFit procedure are the start and end index of the data array provided as first and second parameter. Since DimVector() internally simply calls Setlength() which creates a 0-based array the 3rd LogiFit parameter must be 0 rather than 1; the 4th parameter (last array index) is correct. And the matrix MatV is an m x m matrix where m is the count of fit parameters (5, here); so, the DimMatrix() call must be DimMatrix(MatV, 5, 5); each index i,j in MatV[i,j] then runs from 0 to 4
DimVector() and DimMatrix(), unlike SetLength(), take highest index of an array as a parameter. That is, DimVector(M) calls SetLength(M+1). Same is true for DimMatrix. So, if you call DimMatrix(5,5), you can use indices from 1 to 5.

erik2727

  • New Member
  • *
  • Posts: 14
Re: Lmath ulogifit samples
« Reply #11 on: September 22, 2020, 06:48:47 pm »
hi,

glorfin (Viatcheslav V. Nesterov)

Is there a procedure in Lmath  to do away with Jacobian matrix for nonlinear regression fit similar to c++ gsl library where jacobian matrix can be optional or not?

Just realise and Look as if the jacobian matrix can be coded as empty procedure for max of 3 parameters  in Lmath afterall and yet still can get final output, great!!
« Last Edit: September 23, 2020, 10:58:21 am by erik2727 »

wp

  • Hero Member
  • *****
  • Posts: 7739
Re: Lmath ulogifit samples
« Reply #12 on: September 22, 2020, 08:00:53 pm »
So, if you call DimMatrix(5,5), you can use indices from 1 to 5.
Sorry to be insistent: While I know that this is caused by the math literature where all vectors and matrices begin with index 1 I think it is a bad concept to abuse standard dynamic arrays for 1-based arrays. As a consequence, the user is on his own, there is no help from the compiler or from run-time errors. As an example, the user of LMath must never mix LMath routines with the routines of the standard Math unit unless he is very aware what he is doing.

The following program calculates the average value of the components of a 5-dim vector. The vector is assumed to be 1-based and is created by LMath. The mean is calculated by the function Mean() in unit uMeanSD. Fine.

But now imagine that this a part of a large program, and somebody somewhere else has to calculate the mean again. He is not familiar with LMath and tries to avoid it. He checks the type TVector and sees: it is an array of float. Ah - there is a ready-made routine in the Math unit. But the result is wrong because Math.Mean runs over all elements of the array, including the 0-th element. There is nothing which tells this user that he must not access the index 0. No warning, no error messsage.

Another issue is that being an array of float the function Length() can be applied to a TVector - the result will be off by 1 for a 1-based vector. Same with the function Low().

Code: Pascal  [Select][+][-]
  1. program Project1;
  2.  
  3. uses
  4.   math,
  5.   uTypes, umeansd;
  6.  
  7. var
  8.   v: TVector;
  9.   i: Integer;
  10.  
  11. begin
  12.   RandSeed := 0;
  13.   DimVector(v, 5);
  14.  
  15.   for i := 1 to 5 do
  16.     v[i] := Random;
  17.  
  18.   WriteLn('Mean(LMath) = ', umeansd.Mean(v, 1, 5):0:3);  // ---> 0.661
  19.   WriteLn('Mean(math)  = ', math.Mean(v):0:3);           // ---> 0.551
  20.  
  21.   ReadLn;
  22.  
  23. end.  
« Last Edit: September 22, 2020, 10:11:54 pm by wp »
Mainly Lazarus trunk / fpc 3.2.0 / all 32-bit on Win-10, but many more...

erik2727

  • New Member
  • *
  • Posts: 14
Re: Lmath ulogifit samples
« Reply #13 on: September 22, 2020, 08:47:00 pm »
Hi

Please Don't simply change LMath coded routine and break all the compatibility with previous version just
because  1  and not zero being  set as  starting index option instead, it would be a big nightmares to re-programme existing coded block  from scratch !!!

In my view , TVector.TMatrix  is also  not advisable to  assume as a  regular float array block with 0 as initial, it is much better to  abide by classical mathematics principle i,j concept
« Last Edit: September 23, 2020, 12:38:28 pm by erik2727 »

glorfin

  • Full Member
  • ***
  • Posts: 109
Re: Lmath ulogifit samples
« Reply #14 on: September 23, 2020, 11:29:26 am »
As an example, the user of LMath must never mix LMath routines with the routines of the standard Math unit unless he is very aware what he is doing.

Yes, that is exactly why I try to make LMath self-contained without a necessety to use Math together with it.
In part I agree with you. But problem is that a lot of functions in DMath (which later became LMath) were translated from Pascal, and rewriting them for 0-based arrays can lead not only to broken compatibility with older versions, but also introduce a lot of bugs in process.

Therefore, internally LMath will remain largely 1-based.

However, in many cases there are Lb, Ub parameters, and in such cases one can freely set Lb = 0 and use effectively 0-based TVector or TMatrix. In uMeanSD you can freely call uMeanSD(V,0,5).


If there are no such params, look at the documentation for array base.

 

TinyPortal © 2005-2018