home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / l / l042 / 2.ddi / CHAP9.ARC / LEAST1.INC < prev    next >
Encoding:
Text File  |  1987-12-30  |  40.3 KB  |  862 lines
  1. {----------------------------------------------------------------------------}
  2. {-                                                                          -}
  3. {-     Turbo Pascal Numerical Methods Toolbox                               -}
  4. {-     Copyright (c) 1986, 87 by Borland International, Inc.                -}
  5. {-                                                                          -}
  6. {----------------------------------------------------------------------------}
  7.  
  8. function ModuleName{(Fit : FitType) : TNString40};
  9. begin
  10.   case Fit of
  11.     Poly    : ModuleName := '       Polynomial Least Squares Fit';
  12.     Fourier : ModuleName := ' Finite Fourier Series Least Squares Fit';
  13.     Power   : ModuleName := '       Power Law Least Squares Fit';
  14.     Expo    : ModuleName := '      Exponential Least Squares Fit';
  15.     Log     : ModuleName := '      Logarithmic Least Squares Fit';
  16.     User    : ModuleName := '  User''s Fit - currently powers of X';
  17.   end;
  18. end; { function ModuleName }
  19.  
  20. {-------------------------------------------------------------------------}
  21. {-                                                                       -}
  22. {-            Chebyshev Polynomials                                      -}
  23. {-                                                                       -}
  24. {- This module creates basis vectors to find a least squares solution    -}
  25. {- of the form    f(X) = SUM from i=1 to n of (a[i] * Ti[X]), where Ti   -}
  26. {- is the ith Chebyshev polynomial.  The coefficients of the Ti[X] are   -}
  27. {- converted to coefficients of X^(i-1).                                 -}
  28. {-                                                                       -}
  29. {- The function  ModuleName identifies this as the polynomial fit.       -}
  30. {- The procedure TransformPoly translates and scales the XData to the    -}
  31. {-               interval [-1, 1]. The YData is unchanged.               -}
  32. {- The procedure InverseTransform doesn't do anything in this module.    -}
  33. {- The procedure CreateBasisFunctions creates above basis vectors.       -}
  34. {- The procedure TransformSolution changes the solution vector from      -}
  35. {- coefficients of the Chebyshev polynomials to coefficients of a power  -}
  36. {- series, including adjusting for the shifted data done in TransformPoly-}
  37. {-                                                                       -}
  38. {-------------------------------------------------------------------------}
  39. procedure PolyTransform(NumPoints  : integer;
  40.                     var XData      : TNColumnVector;
  41.                     var YData      : TNColumnVector;
  42.                     var Multiplier : Float;
  43.                     var Constant   : Float;
  44.                     var WData      : TNColumnVector;
  45.                     var ZData      : TNColumnVector;
  46.                     var Error      : byte);
  47.  
  48. {---------------------------------------------------------------}
  49. {- Input : NumPoints, XData, YData,                            -}
  50. {- Output: WData, ZData, Error                                 -}
  51. {-                                                             -}
  52. {- This procedure maps the XData linearly into the interval    -}
  53. {- [-1, 1] returning the transformed data in WData.  The YData -}
  54. {- is passed to ZData unchanged.                               -}
  55. {---------------------------------------------------------------}
  56.  
  57. var
  58.   XDataMin, XDataMax : Float;
  59.   Row : integer;
  60.  
  61. begin
  62.   XDataMin := XData[1];
  63.   XDataMax := XData[1];
  64.   for Row := 1 to NumPoints do
  65.   begin
  66.     if XDataMin > XData[Row] then
  67.       XDataMin := XData[Row];
  68.     if XDataMax < XData[Row] then
  69.       XDataMax := XData[Row];
  70.   end;
  71.   Multiplier := 2.0 / (XDataMax - XDataMin);
  72.   Constant := - Multiplier * (XDataMax + XDataMin) / 2.0;
  73.   for Row := 1 to NumPoints do
  74.     WData[ Row ] := Multiplier * XData[ Row ] + Constant;
  75.   ZDAta := YData;
  76. end; { procedure PolyTransform }
  77.  
  78. procedure PolyInverseTransform(Multiplier : Float;
  79.                                Constant   : Float;
  80.                            var YFit       : Float);
  81.  
  82. {---------------------------------------------------}
  83. {- Input: Multiplier, Constant, YFit               -}
  84. {- Output: YFit                                    -}
  85. {-                                                 -}
  86. {- This procedure undoes the transformation of     -}
  87. {- the YFit values.  Here, no inverse transform    -}
  88. {- is performed because there was no               -}
  89. {- transformation of Y values in procedure         -}
  90. {- PolyTransform.                                  -}
  91. {---------------------------------------------------}
  92.  
  93. begin
  94. end; { procedure PolyInverseTransform }
  95.  
  96. procedure PolyCreateBasisFunctions(NumPoints : integer;
  97.                                    NumTerms  : integer;
  98.                                var WData     : TNColumnVector;
  99.                                var Basis     : TNmatrix);
  100.  
  101. {---------------------------------------------------------}
  102. {- Input: NumPoints, NumTerms, WData                     -}
  103. {- Output: Basis                                         -}
  104. {-                                                       -}
  105. {- This procedure creates a matrix of basis vectors.     -}
  106. {- The basis vectors are the CHEBYSHEV POLYNOMIALS.      -}
  107. {-                                                       -}
  108. {- The elements of the matrix are:                       -}
  109. {-     Basis[i, j] = T[j](WData[i])                      -}
  110. {- where T[j](WData[i]) is the jth basis vector          -}
  111. {- evaluated at the value WData[i].                      -}
  112. {-                                                       -}
  113. {- The vectors are:                                      -}
  114. {-   T[1] = 1                                            -}
  115. {-   T[2] = X                                            -}
  116. {-   T[3] = 2x*X - 1                                     -}
  117. {-   T[4] = (4x*X - 3)*X                                 -}
  118. {-   T[5] = (8x*X - 8)*X*X + 1                           -}
  119. {-   etc.                                                -}
  120. {-                                                       -}
  121. {- The Chebyshev Polynomials can be defined recursively: -}
  122. {-   T[1] = 1, T[2] = X                                  -}
  123. {-   T[j] = 2x * T[j - 1] - T[j - 2]                     -}
  124. {-                                                       -}
  125. {---------------------------------------------------------}
  126.  
  127. var
  128.   Row, Column : integer;
  129.  
  130. begin
  131.   for Row := 1 to NumPoints do
  132.   begin
  133.     Basis[Row, 1] := 1;
  134.     Basis[Row, 2] := WData[Row];
  135.     for Column := 3 to NumTerms do
  136.       Basis[Row, Column] := 2 * WData[Row] * Basis[Row, Column - 1]
  137.                             - Basis[Row, Column - 2];
  138.   end;
  139. end; { procedure PolyCreateBasisFunctions }
  140.  
  141. procedure PolyTransformSolution(NumTerms    : integer;
  142.                             var OldSolution : TNRowVector;
  143.                                 Multiplier  : Float;
  144.                                 Constant    : Float;
  145.                             var NewSolution : TNRowVector);
  146.  
  147. {--------------------------------------------------------------}
  148. {- Input: NumTerms, OldSolution, Multiplier, Constant         -}
  149. {- Output: NewSolution                                        -}
  150. {-                                                            -}
  151. {- The least squares solution will be more useful if it is    -}
  152. {- expressed as a linear expansion of powers of X, rather     -}
  153. {- than as a linear expansion of Chebyshev polynomials.       -}
  154. {-                                                            -}
  155. {- This procedure converts the coefficients of the Chebyshev  -}
  156. {- polynomials to coefficients of powers of X.                -}
  157. {- The vectors ConversionVec and OldConversionVec store       -}
  158. {- information about the relationship between these two sets  -}
  159. {- of coefficients.  This relationship is defined recursively -}
  160. {- above in the procedure PolyCreateBasisFunctions.           -}
  161. {- The parameters Multiplier and Constant define the linear   -}
  162. {- transformation of the XData, so this is accounted for in   -}
  163. {- finding the polynomial coefficients.                       -}
  164. {--------------------------------------------------------------}
  165.  
  166. var
  167.   Index, Term : integer;
  168.   Sum : Float;
  169.   OldConversionVec, ConversionVec : TNRowVector;
  170.  
  171. begin
  172.   FillChar(OldConversionVec, SizeOf(OldConversionVec), 0);
  173.   for Index := 1 to NumTerms do
  174.   begin
  175.     Sum := 0;
  176.     if Index > 1 then
  177.       ConversionVec[Index - 1] := 0;
  178.     for Term := Index to NumTerms do
  179.     begin
  180.       if Term = 1 then
  181.         ConversionVec[Term] := 1.0
  182.       else
  183.         if Term = 2 then
  184.           begin
  185.             if Index = 1 then
  186.               ConversionVec[Term] := Constant
  187.             else
  188.               ConversionVec[Term] := Multiplier
  189.           end
  190.         else
  191.           ConversionVec[Term] := 2 * Multiplier * OldConversionVec[Term - 1]
  192.                                  + 2 * Constant * ConversionVec[Term - 1]
  193.                                  - ConversionVec[Term - 2];
  194.       Sum := Sum + ConversionVec[Term] * OldSolution[Term];
  195.     end;
  196.     NewSolution[Index] := Sum;
  197.     OldConversionVec := ConversionVec;
  198.   end;
  199. end; { procedure PolyTransformSolution }
  200.  
  201. {-------------------------------------------------------------------------}
  202. {-                                                                       -}
  203. {-            Fourier series                                             -}
  204. {-                                                                       -}
  205. {- This module creates basis vectors to find a least squares solution    -}
  206. {- of the form:  f(x) = a[0] + SUM from i=1 to n/2 of (a[i] - Cos(ix) +  -}
  207. {- a[i+1] - Sin(ix)).  A least squares fit with basis vectors 1, Cos(x), -}
  208. {- Sin(x), Cos(2x), Sin(2x), etc. is made to the data (x, y).            -}
  209. {-                                                                       -}
  210. {- The function ModuleName identifies this as the finite Fourier         -}
  211. {-                              series fit.                              -}
  212. {- The procedure Transform doesn't do anything in this module.           -}
  213. {- The procedure InverseTransform doesn't do anything in this module.    -}
  214. {- The procedure CreateBasisFunctions creates the above basis vectors.   -}
  215. {- The procedure TransformSolution doesn't do anything in this module.   -}
  216. {- The first element of the solution vector will be the constant,        -}
  217. {- the second element will be the coefficient of Cos(x), the third       -}
  218. {- element will be the coefficient of Sin(x), the fourth element will    -}
  219. {- be the coefficient of Cos(2x), etc.                                   -}
  220. {-                                                                       -}
  221. {-------------------------------------------------------------------------}
  222.  
  223. procedure FourierTransform(NumPoints       : integer;
  224.                        var XData           : TNColumnVector;
  225.                        var YData           : TNColumnVector;
  226.                            DummyMultiplier : Float;
  227.                            DummyConstant   : Float;
  228.                        var WData           : TNColumnVector;
  229.                        var ZData           : TNColumnVector;
  230.                        var Error           : byte);
  231.  
  232. {---------------------------------------------------------------}
  233. {- Input : NumPoints, XData, YData, DummyMultiplier,           -}
  234. {-         DummyConstant                                       -}
  235. {- Output: WData, ZData, Error                                 -}
  236. {-                                                             -}
  237. {- No transformations are needed for Fourier Series            -}
  238. {---------------------------------------------------------------}
  239.  
  240. begin
  241.   WData := XData;
  242.   ZData := YData;
  243. end; { procedure FourierTransform }
  244.  
  245. procedure FourierInverseTransform(DummyMultiplier : Float;
  246.                                   DummyConstant   : Float;
  247.                               var YFit            : Float);
  248.  
  249. {---------------------------------------------------}
  250. {- Input: DummyMultiplier, DummyConstant, YFit     -}
  251. {- Output: YFit                                    -}
  252. {-                                                 -}
  253. {- This procedure undoes the transformation of     -}
  254. {- the YFit values.  Here, no inverse transform    -}
  255. {- is performed because there was no               -}
  256. {- transformation in procedure Transform.          -}
  257. {---------------------------------------------------}
  258.  
  259. begin
  260. end; { procedure FourierInverseTransform }
  261.  
  262. procedure FourierCreateBasisFunctions(NumPoints : integer;
  263.                                       NumTerms  : integer;
  264.                                   var WData     : TNColumnVector;
  265.                                   var Basis     : TNmatrix);
  266.  
  267. {-------------------------------------------------------------}
  268. {- Input: NumPoints, NumTerms, WData                         -}
  269. {- Output: Basis                                             -}
  270. {-                                                           -}
  271. {- This procedure creates a matrix of basis vectors.         -}
  272. {- The basis vectors are the FOURIER SERIES.                 -}
  273. {-                                                           -}
  274. {- The elements of the matrix are:                           -}
  275. {-     Basis[i, j] = F[j](WData[i])                          -}
  276. {- where F[j](WData[i]) is the jth basis vector              -}
  277. {- evaluated at the value WData[i].                          -}
  278. {-                                                           -}
  279. {- The vectors are:                                          -}
  280. {-      F[1] = 1                                             -}
  281. {-      F[2] = Cos(x);                                       -}
  282. {-      F[3] = Sin(x);                                       -}
  283. {-      F[4] = Cos(2x);                                      -}
  284. {-      F[5] = Sin(2x);                                      -}
  285. {-      F[6] = Cos(3x);                                      -}
  286. {-      F[7] = Sin(3x);                                      -}
  287. {-      etc.                                                 -}
  288. {-                                                           -}
  289. {- This series is defined recursively by:                    -}
  290. {-      F[1] = 1, F[2] = Cos(x),  F[3] = Sin(x)              -}
  291. {-      F[j] = F[2] - F[j - 2] - F[3] - F[j - 1]  for even j -}
  292. {-      F[j] = F[3] - F[j - 3] + F[2] - F[j - 2]  for odd j  -}
  293. {-------------------------------------------------------------}
  294.  
  295. var
  296.   Row, Column : integer;
  297.  
  298. begin
  299.   for Row := 1 to NumPoints do
  300.   begin
  301.     Basis[Row, 1] := 1;
  302.     Basis[Row, 2] := Cos(WData[Row]);
  303.     Basis[Row, 3] := Sin(WData[Row]);
  304.     for Column := 4 to NumTerms do
  305.       if Odd(Column) then
  306.         Basis[Row, Column] := Basis[Row, 3] * Basis[Row, Column-3]
  307.                             + Basis[Row, 2] * Basis[Row, Column-2]
  308.       else
  309.         Basis[Row, Column] := Basis[Row, 2] * Basis[Row, Column-2]
  310.                             - Basis[Row, 3] * Basis[Row, Column-1];
  311.   end;
  312. end; { procedure FourierCreateBasisFunctions }
  313.  
  314.  
  315. procedure FourierTransformSolution(NumTerms        : integer;
  316.                                var OldSolution     : TNRowVector;
  317.                                    DummyMultiplier : Float;
  318.                                    DummyConstant   : Float;
  319.                                var NewSolution     : TNRowVector);
  320.  
  321. {----------------------------------------------------}
  322. {- Input: NumTerms, OldSolution, DummyMultiplier,   -}
  323. {-        DummyConstant                             -}
  324. {- Output: NewSolution                              -}
  325. {-                                                  -}
  326. {- No need to change the coefficients of the        -}
  327. {- Fourier series.                                  -}
  328. {----------------------------------------------------}
  329.  
  330. begin
  331.   NewSolution := OldSolution    { no transformation  }
  332. end; { procedure FourierTransformSolution }
  333.  
  334. {----------------------------------------------------------------------------}
  335. {-                                                                          -}
  336. {-            Y = AX^B                                                      -}
  337. {-                                                                          -}
  338. {- This module creates basis vectors to find a least squares solution       -}
  339. {- of the form    f(X) = A * X^B.  Taking the logarithm of both sides:      -}
  340. {- Ln(f(X)) = Ln(A) + B * Ln(X).  A linear least squares fit                -}
  341. {- (i.e. with basis vectors Ln(X) and 1) is then made to the data           -}
  342. {- (Ln(X), Ln(Y)). The slope will be B, and the intercept will be Ln(A).    -}
  343. {-                                                                          -}
  344. {- The function ModuleName identifies this as the power law fit.            -}
  345. {- The procedure Transform converts the data from (X, Y) to (Ln(X), Ln(Y)). -}
  346. {- The procedure InverseTransform converts from YFit to Exp(YFit)           -}
  347. {- The procedure CreateBasisFunctions creates the vectors Ln(X) and 1.      -}
  348. {- The procedure TransformSolution changes the solution vector from         -}
  349. {- (Ln(A), B) to (A, B).  Therefore, the first coefficient is A,            -}
  350. {- and the second coefficient is B.                                         -}
  351. {-                                                                          -}
  352. {----------------------------------------------------------------------------}
  353.  
  354. procedure PowerTransform(NumPoints     : integer;
  355.                      var XData         : TNColumnVector;
  356.                      var YData         : TNColumnVector;
  357.                      var Multiplier    : Float;
  358.                          DummyConstant : Float;
  359.                      var WData         : TNColumnVector;
  360.                      var ZData         : TNColumnVector;
  361.                      var Error         : byte);
  362.  
  363. {---------------------------------------------------------------}
  364. {- Input : NumPoints, XData, YData, Multiplier, DummyConstant  -}
  365. {- Output: WData, ZData, Error                                 -}
  366. {-                                                             -}
  367. {- This procedure transforms the X and Y values to their       -}
  368. {- logarithms.  A linear least squares fit will then be made   -}
  369. {- to to Ln(Y) = B * Ln(X) + Ln(A). If the Y values are of     -}
  370. {- differing sign, Error = 3 is returned.                      -}
  371. {---------------------------------------------------------------}
  372.  
  373. var
  374.   Index : integer;
  375.   YPoint : Float;
  376.  
  377. begin
  378.   Index := 0;
  379.   if YData[1] < 0 then
  380.     Multiplier := -1
  381.   else
  382.     Multiplier := 1;
  383.   while (Index < NumPoints) and (Error = 0) do
  384.   begin
  385.     Index := Succ(Index);
  386.     if XData[Index] <= 0 then
  387.       Error := 3
  388.     else
  389.       begin
  390.         YPoint := Multiplier * YData[Index];
  391.         if YPoint <= 0 then { The data must all have the same sign  }
  392.           Error := 3
  393.         else
  394.           begin
  395.             WData[Index] := Ln(XData[Index]);
  396.             ZData[Index] := Ln(YPoint);
  397.           end;
  398.       end;
  399.   end; { while }
  400. end; { procedure PowerTransform }
  401.  
  402. procedure PowerInverseTransform(Multiplier    : Float;
  403.                                 DummyConstant : Float;
  404.                             var YFit          : Float);
  405.  
  406. {-------------------------------------------------}
  407. {- Input: Multiplier, DummyConstant, YFit        -}
  408. {- Output: YFit                                  -}
  409. {-                                               -}
  410. {- This procedure undoes the transformation of   -}
  411. {- the YFit values.  Here, the function          -}
  412. {- YFit := Exp(YFit) is performed to undo the    -}
  413. {- Ln transformation in procedure Transform.     -}
  414. {-------------------------------------------------}
  415.  
  416. begin
  417.   YFit := Multiplier * Exp(YFit);
  418. end; { procedure PowerInverseTransform }
  419.  
  420. procedure PowerCreateBasisFunctions(NumPoints : integer;
  421.                                 var NumTerms  : integer;
  422.                                 var WData     : TNColumnVector;
  423.                                 var Basis     : TNmatrix);
  424.  
  425. {---------------------------------------------------------}
  426. {- Input: NumPoints, NumTerms, WData                     -}
  427. {- Output: Basis                                         -}
  428. {-                                                       -}
  429. {- This procedure creates a matrix of basis vectors.     -}
  430. {- The elements of the matrix are:                       -}
  431. {-     Basis[i, j] = C[j](WData[i])                      -}
  432. {- where C[j](WData[i]) is the jth basis vector          -}
  433. {- evaluated at the value WData[i].                      -}
  434. {-                                                       -}
  435. {-                                                       -}
  436. {- The vectors are:                                      -}
  437. {-   C[1] = 1                                            -}
  438. {-   C[2] = X                                            -}
  439. {---------------------------------------------------------}
  440.  
  441. var
  442.   Row : integer;
  443.  
  444. begin
  445.   NumTerms := 2;  { This is only a straight line least squares  }
  446.   for Row := 1 to NumPoints do
  447.   begin
  448.     Basis[Row, 1] := 1;
  449.     Basis[Row, 2] := WData[Row];
  450.   end;
  451. end; { procedure PowerCreateBasisFunctions }
  452.  
  453. procedure PowerTransformSolution(NumTerms      : integer;
  454.                              var OldSolution   : TNRowVector;
  455.                                  Multiplier    : Float;
  456.                                  DummyConstant : Float;
  457.                              var NewSolution   : TNRowVector);
  458.  
  459. {--------------------------------------------------------------}
  460. {- Input: NumTerms, OldSolution, Multiplier, DummyConstant    -}
  461. {- Output: NewSolution                                        -}
  462. {-                                                            -}
  463. {- The least squares solution will be more useful if it is    -}
  464. {- expressed in terms of Y = AX^B, rather than in terms       -}
  465. {- of Ln(Y) = B * Ln(X) + Ln(A).                              -}
  466. {--------------------------------------------------------------}
  467.  
  468. begin
  469.   NewSolution[1] := Multiplier * Exp(OldSolution[1]);
  470. end; { procedure PowerTransformSolution }
  471.  
  472. {-------------------------------------------------------------------------}
  473. {-                                                                       -}
  474. {-            Y = A * Exp(bx)                                            -}
  475. {-                                                                       -}
  476. {- This module creates basis vectors to find a least squares solution    -}
  477. {- of the form    f(X) = A * Exp(bx).  Taking the logarithm of both      -}
  478. {- sides:   Ln(f(X)) = Ln(A) + B * X. A linear least squares fit         -}
  479. {- (i.e. with basis vectors X and 1) is then made to the data            -}
  480. {- (X, Ln(Y)). The slope will be B, and the intercept will be Ln(A).     -}
  481. {-                                                                       -}
  482. {- The function ModuleName identifies this as the exponential fit        -}
  483. {- The procedure Transform converts the data from (X, Y) to (X, Ln(Y)).  -}
  484. {- The procedure InverseTransform converts from YFit to Exp(YFit)        -}
  485. {- The procedure CreateBasisFunctions creates the vectors 1 and X.       -}
  486. {- The procedure TransformSolution changes the solution vector from      -}
  487. {- (Ln(A), B) to (A, B).  Therefore, the first coefficient is a,         -}
  488. {- and the second coefficient is B.                                      -}
  489. {-                                                                       -}
  490. {-------------------------------------------------------------------------}
  491.  
  492. procedure ExpoTransform(NumPoints     : integer;
  493.                     var XData         : TNColumnVector;
  494.                     var YData         : TNColumnVector;
  495.                     var Multiplier    : Float;
  496.                         DummyConstant : Float;
  497.                     var WData         : TNColumnVector;
  498.                     var ZData         : TNColumnVector;
  499.                     var Error         : byte);
  500.  
  501. {---------------------------------------------------------------}
  502. {- Input : NumPoints, XData, YData, Multiplier, DummyConstant  -}
  503. {- Output: WData, ZData, Error                                 -}
  504. {-                                                             -}
  505. {- This procedure transforms the Y values to their             -}
  506. {- logarithms.  A linear least squares fit will then be made   -}
  507. {- to Ln(Y) = bx + Ln(A). If the Y values are of different     -}
  508. {- sign, then Error = 3 is returned.                           -}
  509. {---------------------------------------------------------------}
  510.  
  511. var
  512.   Index : integer;
  513.   YPoint : Float;
  514.  
  515. begin
  516.   WData := XData;
  517.   if  YData[1] < 0 then
  518.     Multiplier := -1
  519.   else
  520.     Multiplier := 1;
  521.   Index := 0;
  522.   while (Index < NumPoints) and (Error = 0) do
  523.   begin
  524.     Index := Succ(Index);
  525.     YPoint := Multiplier * YData[Index];
  526.     if YPoint <= 0 then
  527.       Error := 3   { The Y values must all have the same sign }
  528.     else
  529.       ZData[Index] := Ln(YPoint);
  530.   end;
  531. end; { procedure ExpoTransform }
  532.  
  533. procedure ExpoInverseTransform(Multiplier    : Float;
  534.                                DummyConstant : Float;
  535.                            var YFit          : Float);
  536.  
  537. {-------------------------------------------------}
  538. {- Input: Multiplier, DummyConstant, YFit        -}
  539. {- Output: YFit                                  -}
  540. {-                                               -}
  541. {- This procedure undoes the transformation of   -}
  542. {- the YFit values.  Here, the function          -}
  543. {- YFit := Exp(YFit) is performed to undo the    -}
  544. {- Ln transformation in procedure Transform.     -}
  545. {-------------------------------------------------}
  546.  
  547. begin
  548.   YFit := Multiplier * Exp(YFit);
  549. end; { procedure ExpoInverseTransform }
  550.  
  551. procedure ExpoCreateBasisFunctions(NumPoints : integer;
  552.                                var NumTerms  : integer;
  553.                                var WData     : TNColumnVector;
  554.                                var Basis     : TNmatrix);
  555.  
  556. {---------------------------------------------------------}
  557. {- Input: NumPoints, NumTerms, WData                     -}
  558. {- Output: Basis                                         -}
  559. {-                                                       -}
  560. {- This procedure creates a matrix of basis vectors.     -}
  561. {-                                                       -}
  562. {- The elements of the matrix are:                       -}
  563. {-     Basis[i, j] = C[j](WData[i])                      -}
  564. {- where C[j](WData[i]) is the jth basis vector          -}
  565. {- evaluated at the value WData[i].                      -}
  566. {-                                                       -}
  567. {- The vectors are:                                      -}
  568. {-   C[1] = 1                                            -}
  569. {-   C[2] = X                                            -}
  570. {---------------------------------------------------------}
  571.  
  572. var
  573.   Row : integer;
  574.  
  575. begin
  576.   NumTerms := 2;  { This is only a straight line least squares }
  577.   for Row := 1 to NumPoints do
  578.   begin
  579.     Basis[Row, 1] := 1;
  580.     Basis[Row, 2] := WData[Row];
  581.   end;
  582. end; { procedure ExpoCreateBasisFunctions }
  583.  
  584. procedure ExpoTransformSolution(NumTerms      : integer;
  585.                             var OldSolution   : TNRowVector;
  586.                                 Multiplier    : Float;
  587.                                 DummyConstant : Float;
  588.                             var NewSolution   : TNRowVector);
  589.  
  590. {--------------------------------------------------------------}
  591. {- Input: NumTerms, OldSolution, Multiplier, DummyConstant    -}
  592. {- Output: NewSolution                                        -}
  593. {-                                                            -}
  594. {- The least squares solution will be more useful if it is    -}
  595. {- expressed in terms of Y = A - Exp(bx), rather than in      -}
  596. {- terms of Ln(Y) = B - X + Ln(A).                            -}
  597. {--------------------------------------------------------------}
  598.  
  599. begin
  600.   NewSolution[1] := Multiplier * Exp(OldSolution[1]);
  601. end; { procedure ExpoTransformSolution }
  602.  
  603. {-------------------------------------------------------------------------}
  604. {-                                                                       -}
  605. {-            Y = A * Ln(bx)                                             -}
  606. {-                                                                       -}
  607. {- This module creates basis vectors to find a least squares solution    -}
  608. {- of the form    f(X) = A * Ln(bx).  Rewriting the right side of the    -}
  609. {- equation:  f(X) = A * Ln(X) + A * Ln(B).  A linear least squares fit  -}
  610. {- (i.e. with basis vectors Ln(X) and 1) is then made to the data        -}
  611. {- (Ln(X), Y). The slope will be A, and the intercept will be A * Ln(B). -}
  612. {-                                                                       -}
  613. {- The function ModuleName identifies this as the logarithmic fit        -}
  614. {- The procedure Transform converts the data from (X, Y) to (Ln(X), Y).  -}
  615. {- The procedure InverseTransform doesn't do anything in this module.    -}
  616. {- The procedure CreateBasisFunctions creates the vectors Ln(X) and 1.   -}
  617. {- The procedure TransformSolution changes the solution vector from      -}
  618. {- (A, A * Ln(B)) to (A, B).  Therefore, the first coefficient is A,     -}
  619. {- and the second coefficient is B.                                      -}
  620. {-                                                                       -}
  621. {-------------------------------------------------------------------------}
  622.  
  623. procedure LogTransform(NumPoints     : integer;
  624.                    var XData         : TNColumnVector;
  625.                    var YData         : TNColumnVector;
  626.                    var Multiplier    : Float;
  627.                        DummyConstant : Float;
  628.                    var WData         : TNColumnVector;
  629.                    var ZData         : TNColumnVector;
  630.                    var Error         : byte);
  631.  
  632. {---------------------------------------------------------------}
  633. {- Input : NumPoints, XData, YData, Multiplier, DummyConstant  -}
  634. {- Output: WData, ZData, Error                                 -}
  635. {-                                                             -}
  636. {- This procedure transforms the X values to their             -}
  637. {- logarithms.  A linear least squares fit will then be made   -}
  638. {- to Y = ALn(X) + ALn(B).  If the X values are of different   -}
  639. {- sign, then Error = 3 is returned.                           -}
  640. {---------------------------------------------------------------}
  641.  
  642. var
  643.   Index : integer;
  644.   XPoint : Float;
  645.  
  646. begin
  647.   ZData := YData;
  648.   if  XData[1] < 0 then
  649.     Multiplier := -1
  650.   else
  651.     Multiplier := 1;
  652.   Index := 0;
  653.   while (Index < NumPoints) and (Error = 0) do
  654.   begin
  655.     Index := Succ(Index);
  656.     XPoint := Multiplier * XData[Index];
  657.     if XPoint <= 0 then
  658.       Error := 3   { The X values must all have the same sign  }
  659.     else
  660.       WData[Index] := Ln(XPoint);
  661.   end;
  662. end; { procedure LogTransform }
  663.  
  664. procedure LogInverseTransform(Multiplier    : Float;
  665.                               DummyConstant : Float;
  666.                           var YFit          : Float);
  667.  
  668. {---------------------------------------------------}
  669. {- Input: Multiplier, DummyConstant, YFit          -}
  670. {- Output: YFit                                    -}
  671. {-                                                 -}
  672. {- This procedure undoes the transformation of     -}
  673. {- the YFit values.  Here, no inverse transform    -}
  674. {- is performed because the was no transformation  -}
  675. {- in procedure Transform.                         -}
  676. {---------------------------------------------------}
  677.  
  678. begin
  679. end;{ procedure LogInverseTransform }
  680.  
  681. procedure LogCreateBasisFunctions(NumPoints : integer;
  682.                               var NumTerms  : integer;
  683.                               var WData     : TNColumnVector;
  684.                               var Basis     : TNmatrix);
  685.  
  686. {---------------------------------------------------------}
  687. {- Input: NumPoints, NumTerms, WData                     -}
  688. {- Output: Basis                                         -}
  689. {-                                                       -}
  690. {- This procedure creates a matrix of basis vectors.     -}
  691. {-                                                       -}
  692. {- The elements of the matrix are:                       -}
  693. {-     Basis[i, j] = C[j](WData[i])                      -}
  694. {- where C[j](WData[i]) is the jth basis vector          -}
  695. {- evaluated at the value WData[i].                      -}
  696. {-                                                       -}
  697. {- The vectors are:                                      -}
  698. {-   C[1] = X                                            -}
  699. {-   C[2] = 1                                            -}
  700. {---------------------------------------------------------}
  701.  
  702. var
  703.   Row : integer;
  704.  
  705. begin
  706.   NumTerms := 2; { This is only a straight line least squares  }
  707.   for Row := 1 to NumPoints do
  708.   begin
  709.     Basis[Row, 1] := WData[Row];
  710.     Basis[Row, 2] := 1;
  711.   end;
  712. end; { procedure LogCreateBasisFunctions }
  713.  
  714. procedure LogTransformSolution(NumTerms      : integer;
  715.                            var OldSolution   : TNRowVector;
  716.                                Multiplier    : Float;
  717.                                DummyConstant : Float;
  718.                            var NewSolution   : TNRowVector);
  719.  
  720. {--------------------------------------------------------------}
  721. {- Input: NumTerms, OldSolution, Multiplier, DummyConstant    -}
  722. {- Output: NewSolution                                        -}
  723. {-                                                            -}
  724. {- The least squares solution will be more useful if it is    -}
  725. {- expressed in terms of Y = A - Ln(bx), rather than in       -}
  726. {- terms of Y = ALn(X) + ALn(B).                              -}
  727. {--------------------------------------------------------------}
  728.  
  729. begin
  730.   NewSolution[2] := Multiplier * Exp(OldSolution[2]/OldSolution[1]);
  731. end; { procedure LogTransformSolution }
  732.  
  733. {-------------------------------------------------------------------------}
  734. {-                                                                       -}
  735. {-            User Defined function                                      -}
  736. {-                                                                       -}
  737. {- This module provides the format for the user to create her own basis  -}
  738. {- vectors.                                                              -}
  739. {-                                                                       -}
  740. {- The function ModuleName identifies this as the user's Module.         -}
  741. {-   This function should be changed to identify the user's basis.       -}
  742. {- The procedure Transform converts the data from (X, Y) to an           -}
  743. {-   appropriate format (e.g. (Ln(X), Ln(Y)) ). If no transformation     -}
  744. {-   is needed, this procedure should not be changed.                    -}
  745. {- The procedure InverseTransform undoes the transformation of the       -}
  746. {-   Y-coordinate.  In the above example, the procedure would perform    -}
  747. {-   the function YFit := Exp(YFit).  This allows comparison between     -}
  748. {-   least squares approximation and the actual Y-values.                -}
  749. {- The procedure CreateBasisFunctions creates the  basis vectors. The    -}
  750. {-   least squares solution will be coefficients of these basis vectors. -}
  751. {-   Currently the basis vectors are powers of X.                        -}
  752. {- The procedure TransformSolution transforms the solution vector to     -}
  753. {-   an appropriate format.  This usually undoes the transformation made -}
  754. {-   in procedure Transform. If no transformation is needed, this        -}
  755. {-   procedure should not be changed.                                    -}
  756. {-                                                                       -}
  757. {-------------------------------------------------------------------------}
  758.  
  759. procedure UserTransform(NumPoints       : integer;
  760.                     var XData           : TNColumnVector;
  761.                     var YData           : TNColumnVector;
  762.                     var DummyMultiplier : Float;
  763.                     var DummyConstant   : Float;
  764.                     var WData           : TNColumnVector;
  765.                     var ZData           : TNColumnVector;
  766.                     var Error           : byte);
  767.  
  768. {---------------------------------------------------------------}
  769. {- Input : NumPoints, XData, YData, DummyMultiplier,           -}
  770. {-         DummyConstant                                       -}
  771. {- Output: WData, ZData, Error                                 -}
  772. {-                                                             -}
  773. {- This procedure transforms the input data to an appropriate  -}
  774. {- format.  The transformed (or possibly unchanged) data is    -}
  775. {- returned in WData, ZData.                                   -}
  776. {---------------------------------------------------------------}
  777.  
  778. var
  779.   Index : integer;
  780.  
  781. begin
  782.   WData := XData;     { No transformation  }
  783.   ZData := YData;     { No transformation  }
  784. end; { procedure UserTransform }
  785.  
  786. procedure UserInverseTransform(DummyMultiplier : Float;
  787.                                DummyConstant   : Float;
  788.                            var YFit            : Float);
  789.  
  790. {---------------------------------------------------}
  791. {- Input: DummyMultiplier, DummyConstant, YFit     -}
  792. {- Output: YFit                                    -}
  793. {-                                                 -}
  794. {- This procedure undoes the transformation of     -}
  795. {- the YFit values.  No inverse transformation     -}
  796. {- may be necessary.                               -}
  797. {---------------------------------------------------}
  798.  
  799. begin
  800. end; { procedure UserInverseTransform }
  801.  
  802. procedure UserCreateBasisFunctions(NumPoints : integer;
  803.                                var NumTerms  : integer;
  804.                                var WData     : TNColumnVector;
  805.                                var Basis     : TNmatrix);
  806.  
  807. {---------------------------------------------------------}
  808. {- Input: NumPoints, NumTerms, WData                     -}
  809. {- Output: Basis                                         -}
  810. {-                                                       -}
  811. {- This procedure creates a matrix of basis vectors.     -}
  812. {- The user must modify this procedure.                  -}
  813. {-                                                       -}
  814. {- The elements of the matrix must be:                   -}
  815. {-     Basis[i, j] = Bj(WData[i])                        -}
  816. {- where Bj(WData[i]) is the jth basis vector evaluated  -}
  817. {- at the value WData[i].                                -}
  818. {-                                                       -}
  819. {- For example, the basis vector might be powers of X:   -}
  820. {-   B1 = 1                                              -}
  821. {-   B2 = X                                              -}
  822. {-   B3 = X^2                                            -}
  823. {-   B4 = X^3                                            -}
  824. {-   etc.                                                -}
  825. {-                                                       -}
  826. {- These vectors can be defined recursively:             -}
  827. {-   B1 = 1                                              -}
  828. {-   Bj = X * B[j - 1]                                   -}
  829. {---------------------------------------------------------}
  830.  
  831. var
  832.   Row, Column : integer;
  833.  
  834. begin
  835.   for Row := 1 to NumPoints do
  836.   begin
  837.     Basis[Row, 1] := 1;
  838.     for Column := 2 to NumTerms do
  839.     Basis[Row, Column] := WData[Row] * Basis[Row, Column - 1];
  840.   end;
  841. end; { procedure UserCreateBasisFunctions }
  842.  
  843. procedure UserTransformSolution(NumTerms        : integer;
  844.                             var OldSolution     : TNRowVector;
  845.                                 DummyMultiplier : Float;
  846.                                 DummyConstant   : Float;
  847.                             var NewSolution     : TNRowVector);
  848.  
  849. {--------------------------------------------------------------}
  850. {- Input: NumTerms, OldSolution, DummyMultiplier,             -}
  851. {-        DummyConstant                                       -}
  852. {- Output: NewSolution                                        -}
  853. {-                                                            -}
  854. {- This procedure transforms the solution into an appropriate -}
  855. {- form.  The transformed (or possibly unchanged) solution    -}
  856. {- is returned in NewSoluition.                               -}
  857. {--------------------------------------------------------------}
  858.  
  859. begin
  860.   NewSolution := OldSolution; { No transformation  }
  861. end; { procedure UserTransformSolution }
  862.