← Back to branch summary

~ubuntu-branches/ubuntu/trusty/mricron/trusty-proposed

« back to all changes in this revision

Viewing changes to npm/dmath/fitiexpo.pas

  • Committer: Bazaar Package Importer
  • Author(s): Michael Hanke
  • Date: 2010-07-29 22:07:43 UTC
  • Revision ID: james.westby@ubuntu.com-20100729220743-q621ts2zj806gu0n
Tags: upstream-0.20100725.1~dfsg.1
Import upstream version 0.20100725.1~dfsg.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
npm/dmath/fitiexpo.pas
 
1
{ **********************************************************************
 
2
  *                         Unit FITIEXPO.PAS                          *
 
3
  *                            Version 1.2                             *
 
4
  *                    (c) J. Debord, August 2000                      *
 
5
  **********************************************************************
 
6
  This unit fits the increasing exponential :
 
7
 
 
8
                          y = A.[1 - exp(-k.x)]
 
9
 
 
10
  ********************************************************************** }
 
11
 
 
12
unit FitIExpo;
 
13
 
 
14
{$F+}
 
15
 
 
16
interface
 
17
 
 
18
uses
 
19
  FMath, Matrices, Stat, Regress;
 
20
 
 
21
function FuncName : String;
 
22
 
 
23
function FirstParam : Integer;
 
24
 
 
25
function LastParam : Integer;
 
26
 
 
27
function ParamName(I : Integer) : String;
 
28
 
 
29
function RegFunc(X : Float; B : PVector) : Float;
 
30
 
 
31
procedure DerivProc(X : Float; B, D : PVector);
 
32
 
 
33
function FitModel(Method : Integer; X, Y, W : PVector;
 
34
                  N : Integer; B : PVector) : Integer;
 
35
 
 
36
 
 
37
implementation
 
38
 
 
39
  function FuncName : String;
 
40
  { --------------------------------------------------------------------
 
41
    Returns the name of the regression function
 
42
    -------------------------------------------------------------------- }
 
43
  begin
 
44
    FuncName := 'y = A[1 - exp(-k.x)]';
 
45
  end;
 
46
 
 
47
  function FirstParam : Integer;
 
48
  { --------------------------------------------------------------------
 
49
    Returns the index of the first parameter to be fitted
 
50
    -------------------------------------------------------------------- }
 
51
  begin
 
52
    FirstParam := 0;
 
53
  end;
 
54
 
 
55
  function LastParam : Integer;
 
56
  { --------------------------------------------------------------------
 
57
    Returns the index of the last parameter to be fitted
 
58
    -------------------------------------------------------------------- }
 
59
  begin
 
60
    LastParam := 1;
 
61
  end;
 
62
 
 
63
  function ParamName(I : Integer) : String;
 
64
  { --------------------------------------------------------------------
 
65
    Returns the name of the I-th parameter
 
66
    -------------------------------------------------------------------- }
 
67
  begin
 
68
    case I of
 
69
      0 : ParamName := 'A';
 
70
      1 : ParamName := 'k';
 
71
    end;
 
72
  end;
 
73
 
 
74
  function RegFunc(X : Float; B : PVector) : Float;
 
75
  { --------------------------------------------------------------------
 
76
    Computes the regression function at point X
 
77
    B is the vector of parameters, such that :
 
78
 
 
79
    B^[0] = A     B^[1] = k
 
80
    -------------------------------------------------------------------- }
 
81
  begin
 
82
    RegFunc := B^[0] * (1.0 - Expo(- B^[1] * X));
 
83
  end;
 
84
 
 
85
  procedure DerivProc(X : Float; B, D : PVector);
 
86
  { --------------------------------------------------------------------
 
87
    Computes the derivatives of the regression function at point X
 
88
    with respect to the parameters B. The results are returned in D.
 
89
    D^[I] contains the derivative with respect to the I-th parameter.
 
90
    -------------------------------------------------------------------- }
 
91
  var
 
92
    E : Float;
 
93
  begin
 
94
    E := Expo(- B^[1] * X);  { exp(-k.x) }
 
95
    D^[0] := 1.0 - E;        { dy/dA = 1 - exp(-k.x) }
 
96
    D^[1] := B^[0] * X * E;  { dy/dk = A.x.exp(-k.x) }
 
97
  end;
 
98
 
 
99
  function FitModel(Method : Integer; X, Y, W : PVector;
 
100
                    N : Integer; B : PVector) : Integer;
 
101
  { --------------------------------------------------------------------
 
102
    Approximate fit of the increasing exponential by linear regression:
 
103
    Ln(1 - y/A) = -k.x
 
104
    --------------------------------------------------------------------
 
105
    Input :  Method = 0 for unweighted regression, 1 for weighted
 
106
             X, Y   = point coordinates
 
107
             W      = weights
 
108
             N      = number of points
 
109
    Output : B      = estimated regression parameters
 
110
    -------------------------------------------------------------------- }
 
111
  var
 
112
    Y1 : PVector;       { Transformed ordinates }
 
113
    W1 : PVector;       { Weights }
 
114
    A : PVector;        { Linear regression parameters }
 
115
    V : PMatrix;        { Variance-covariance matrix }
 
116
    K : Integer;        { Loop variable }
 
117
    ErrCode : Integer;  { Error code }
 
118
  begin
 
119
    DimVector(Y1, N);
 
120
    DimVector(W1, N);
 
121
    DimVector(A, 1);
 
122
    DimMatrix(V, 1, 1);
 
123
 
 
124
    { Estimation of A }
 
125
    B^[0] := 1.1 * Max(Y, 1, N);
 
126
 
 
127
    for K := 1 to N do
 
128
      begin
 
129
        Y1^[K] := Log(1.0 - Y^[K] / B^[0]);
 
130
        W1^[K] := Sqr(Y^[K] - B^[0]);
 
131
        if Method = 1 then W1^[K] := W1^[K] * W^[K];
 
132
      end;
 
133
 
 
134
    ErrCode := WLinFit(X, Y1, W1, N, A, V);
 
135
 
 
136
    if ErrCode = MAT_OK then
 
137
      B^[1] := - A^[1];
 
138
 
 
139
    FitModel := ErrCode;
 
140
 
 
141
    DelVector(Y1, N);
 
142
    DelVector(W1, N);
 
143
    DelVector(A, 1);
 
144
    DelMatrix(V, 1, 1);
 
145
  end;
 
146
 
 
147
  end.