home *** CD-ROM | disk | FTP | other *** search
/ Turbo Toolbox / Turbo_Toolbox.iso / 1989 / 09 / hitech / linfit.pas < prev   
Encoding:
Pascal/Delphi Source File  |  1989-06-27  |  11.7 KB  |  353 lines
  1. (* ------------------------------------------------------ *)
  2. (*                    LINFIT.PAS                          *)
  3. (*         Robust Fitting in Turbo Pascal 4/5             *)
  4. (*           (c) TOOLBOX & Peter Kurzweil                 *)
  5. (* ------------------------------------------------------ *)
  6. PROGRAM RobustFitting;
  7. USES Crt;
  8. CONST
  9.   npt    = 200;                       { Zahl der Meßpunkte }
  10.   spread = 1.0;                                 { Streuung }
  11. TYPE
  12.   {$IFDEF CPU87}
  13.     REAL = Extended;
  14.   {$ENDIF}
  15.   Vector = ARRAY[1..npt] OF REAL;
  16. VAR
  17.   a, b, da, db, r, xm, ym, sx,
  18.   sy, xo, dy, phi, xmin, xmax : REAL;
  19.   x, y                        : Vector;
  20.   i, k, i1, i2                : WORD;
  21.  
  22.   FUNCTION Log10(x : REAL) : REAL;      { dek. Logarithmus }
  23.   BEGIN
  24.     IF x <> 0 THEN
  25.       Log10 := Ln(Abs(x)) / Ln(10.0)
  26.     ELSE
  27.       Log10 := 0;
  28.   END;
  29.  
  30.   FUNCTION Exp10(x : REAL) : REAL;             { 10 hoch x }
  31.   BEGIN
  32.     Exp10 := Exp(x * Ln(10));
  33.   END;
  34.  
  35.   PROCEDURE ProduceOutliers;         { Ausreißer aufprägen }
  36.   BEGIN
  37.     IF spread <> 0 THEN BEGIN
  38.       FOR k := 1 TO (npt DIV 10) DO BEGIN
  39.         i := Random(npt);             { Punkt auswählen... }
  40.         x[i] := ((i - 1) * (xmax - xmin) / npt + xmin) *
  41.                  (Random - 0.5);
  42.                                      { ...und verunstalten }
  43.         y[i] := Sin(x[i]) * (xmax - xmin) * b + a;
  44.       END;
  45.     END;
  46.   END;
  47.  
  48.   FUNCTION StdAbw(z : Vector; k1, k2 : WORD) : REAL;
  49.                                       { Standardabweichung }
  50.   VAR
  51.     s, zs, zz : REAL;
  52.     i, n      : WORD;
  53.   BEGIN
  54.     IF k1 > k2 THEN BEGIN
  55.       i := k1;  k1 := k2;  k2 := i;
  56.     END;
  57.     n  := k2 - k1 + 1;
  58.     zz := 0;
  59.     zs := 0;
  60.     FOR i := k1 TO k2 DO BEGIN
  61.       zs := zs + z[i];
  62.       zz := zz + Sqr(z[i]);
  63.     END;
  64.     s := Sqrt((zz - zs * zs/n) / (n - 1));
  65.     StdAbw := s;
  66.   END;
  67.  
  68.   FUNCTION Mittel(z : Vector; k1, k2 : WORD) : REAL;
  69.                                               { Mittelwert }
  70.   VAR
  71.     i, n : WORD;
  72.     s    : REAL;
  73.   BEGIN
  74.     IF k1 > k2 THEN BEGIN
  75.       i := k1;  k1 := k2;  k2 := i;
  76.     END;
  77.     n := k2 - k1 + 1;
  78.     s := 0;
  79.     FOR i := k1 TO k2 DO s := s + z[i];
  80.     Mittel := s/n;
  81.   END;
  82.  
  83.   PROCEDURE LinFit(x, y : Vector;  k1, k2     : WORD;
  84.                    VAR a, b, r, da, db,
  85.                        dy, xo, xm, ym, sx, sy : REAL);
  86.   VAR
  87.     xs, ys, xx, yy, xy,
  88.     r1, r2, a1, a2, b1, b2 : REAL;
  89.     i, n                   : WORD;
  90.   BEGIN
  91.     IF k1 > k2 THEN BEGIN
  92.       n := k1;  k1 := k2;  k2 := n;
  93.     END;
  94.     n := k2 - k1 + 1;                    { Zahl der Punkte }
  95.     da := 0;  db := 0;  r := 1;
  96.     sx := 0;  sy := 0;
  97.     CASE n OF
  98.       0 : Exit;                                  { sinnlos }
  99.       1 : BEGIN                    { Gerade durch Ursprung }
  100.             a  := y[k1];  b  := y[k1];
  101.             xm := x[k1];  ym := y[k1];
  102.             IF x[k1] <> 0 THEN b := y[k1] / x[k1];
  103.           END;
  104.       2 : BEGIN                 { Gerade durch zwei Punkte }
  105.             b1 := y[k2] - y[k1];
  106.             b2 := x[k2] - x[k1];
  107.             IF b2 = 0 THEN BEGIN
  108.               b := 99999;
  109.               a := 0;     xo := x[k1];         { Vertikale }
  110.             END ELSE BEGIN
  111.               b  := b1 / b2;
  112.               xm := (x[k1] + x[k2]) / 2;
  113.               ym := (y[k1] + y[k2]) / 2;
  114.               a  := ym - b * xm;
  115.               IF b <> 0 THEN xo := -a / b ELSE xo := -0;
  116.               sx := StdAbw(x, k1, k2);
  117.               sy := StdAbw(y, k1, k2);
  118.             END;
  119.           END;
  120.       ELSE BEGIN                        { Ausgleichsgerade }
  121.         xs := 0;  ys := 0;
  122.         xx := 0;  yy := 0;  xy := 0;
  123.         FOR i := k1 TO k2 DO BEGIN           { Gauß-Summen }
  124.           xs := xs + x[i];       ys := ys + y[i];
  125.           xx := xx + x[i]*x[i];  yy := yy + y[i]*y[i];
  126.           xy := xy + x[i]*y[i];
  127.         END;
  128.         a1 := ys*xx - xy*xs; a2 := n*xx - xs*xs;
  129.         b1 := n*xy - ys*xs;
  130.         IF a2 = 0 THEN BEGIN
  131.           a  := a1;  b := 99999;  r := 1;
  132.           xo := x[k1];                         { Vertikale }
  133.         END ELSE BEGIN
  134.           a  := a1/a2;
  135.           b  := b1/a2;
  136.           r1 := xy - xs * ys / n;
  137.           r2 := Abs((xx - xs * xs / n)*(yy - ys * ys / n));
  138.           IF r2 <> 0 THEN r  := r1 / Sqrt(r2) ELSE r  := 1;
  139.           IF b <> 0  THEN xo := -a / b        ELSE xo := 0;
  140.         END;
  141.         dy := Sqrt(Abs((yy - a*ys - b*xy)/(n-2)));
  142.                                                   { Fehler }
  143.         da := dy * Sqrt(Abs(xx/(n*xx - xs*xs)));
  144.         db := dy * Sqrt(Abs(n/(n*xx - xs*xs)));
  145.         xm := Mittel(x, k1, k2);
  146.         ym := Mittel(y, k1, k2);
  147.                                            { Statistisches }
  148.         sx := StdAbw(x, k1, k2);
  149.         sy := StdAbw(y, k1, k2);
  150.       END;
  151.     END;
  152.   END;
  153.  
  154.   PROCEDURE HeapSort(VAR a : Vector; n : WORD);
  155.                             { Zahlen aufsteigend sortieren }
  156.   VAR
  157.     l, j, k, i : WORD;
  158.     h          : REAL;
  159.   BEGIN
  160.     l := (n DIV 2) + 1;  k := n;
  161.     WHILE TRUE DO BEGIN
  162.       IF l > 1 THEN BEGIN
  163.         l := l - 1;
  164.         h := a[l];
  165.       END ELSE BEGIN
  166.         h := a[k];  a[k] := a[1];  k := k - 1;
  167.         IF k = 1 THEN BEGIN
  168.           a[1] := h;  Exit;
  169.         END;
  170.       END;
  171.       i := l;  j := l + l;
  172.       WHILE j <= k DO BEGIN
  173.         IF j < k THEN IF a[j] < a[j+1] THEN j := j + 1;
  174.         IF h < a[j] THEN BEGIN
  175.           a[i] := a[j];  i := j;  j := j + j;
  176.         END ELSE j := k + 1;
  177.       END;
  178.       a[i] := h;
  179.     END;
  180.   END;
  181.  
  182.   PROCEDURE RobustLinFit(x, y : Vector;  k1, k2 : WORD;
  183.                          VAR a, b, da, db, dy   : REAL);
  184.   LABEL 9;
  185.   VAR
  186.     i, k, n                                       : WORD;
  187.     xs, ys, xy, xx, f, f1, f2, b1, b2, d, chi2, h : REAL;
  188.  
  189.     FUNCTION Theory(b : REAL; VAR a, dy : REAL) : REAL;
  190.                                 { Rechte Seite von Gl. (1) }
  191.     VAR
  192.       d, s            : REAL;
  193.       i, k, n, n1, n2 : WORD;
  194.       f               : Vector;
  195.     BEGIN
  196.       n  := k2 - k1 + 1;
  197.       n2 := (n+1) DIV 2;
  198.       n1 := (n+1) - n2;   i := 0;
  199.       FOR k := k1 TO k2 DO BEGIN
  200.         i := i + 1; f[i] := y[k] - b * x[k];
  201.       END;
  202.       HeapSort(f, n);
  203.       a := 0.5*(f[n2] + f[n1]);   { Bestimmung des Medians }
  204.       s := 0; i:=0; dy:=0;
  205.       FOR i := k1 to k2 DO BEGIN
  206.         d  := y[i] - (b * x[i] + a);
  207.         dy := dy + Abs(d);
  208.         IF d > 0 THEN s := s + x[i] ELSE s := s - x[i];
  209.       END;
  210.       Theory := s;
  211.     END;
  212.  
  213.   BEGIN
  214.     IF k1 > k2 THEN BEGIN
  215.       n := k1;  k1 := k2;  k2 := n;
  216.     END;
  217.     n  := k2 - k1 + 1;                  { Stützstellenzahl }
  218.     xs := 0;  ys := 0;
  219.     xy := 0;  xx := 0;  chi2:=0;
  220.     FOR i := k1 TO k2 DO BEGIN               { Gauß-Summen }
  221.       xs := xs + x[i];  xy := xy + x[i] * y[i];
  222.       ys := ys + y[i];  xx := xx + x[i] * x[i];
  223.     END;
  224.     d := n*xx - xs*xs;
  225.     a := (xx*ys - xs*xy)/d;
  226.     b := (n*xy - xs*ys)/d;
  227.     FOR i := k1 TO k2 DO
  228.       chi2 := chi2 + Sqr(y[i] - (a + b * x[i]));
  229.                                       { Fehlerquadratsumme }
  230.     h  := Sqrt(chi2/d);
  231.     b1 := b;
  232.     f1 := Theory(b1, a, dy);
  233.     IF f1 > 0 THEN
  234.       b2 := b + Abs(3 * h)
  235.     ELSE
  236.       b2 := b - Abs(3 * h);
  237.     f2 := Theory(b2, a, dy);
  238.     WHILE ((f1 * f2 > 0) AND (f1 <> f2)) DO BEGIN
  239.                               { Löse Gl. (2) für versch. b }
  240.       b  := 2 * b2 - b1;
  241.       b1 := b2;  f1 := f2;
  242.       b2 := b;   f2 := Theory(b2, a, dy);
  243.     END;
  244.     h := 0.01 * h;
  245.     WHILE Abs(b2 - b1) > h DO BEGIN     { Nullstellensuche }
  246.       b := 0.5 * (b1 + b2);
  247.       IF ((b = b1) OR (b = b2)) THEN GOTO 9;
  248.       f := Theory(b, a, dy);
  249.       IF f * f1 >= 0 THEN BEGIN
  250.         f1 := f;  b1 := b;
  251.       END ELSE BEGIN
  252.         f2 := f;  b2 := b;
  253.       END;
  254.     END;
  255. 9:
  256.     dy := dy/n;                          { Fehler des Fits }
  257.     da := dy * Sqrt(Abs(xx / (n * xx - xs * xs)));
  258.     db := dy * Sqrt(Abs(n / (n * xx - xs * xs)));
  259.   END;
  260.  
  261.   PROCEDURE LinearFit;
  262.                { Treiberroutine und Ausgabe der Ergebnisse }
  263.   BEGIN
  264.     ClrScr;
  265.     WriteLn('A u s g l e i c h s g e r a d e',
  266.             '    y = b x + a');
  267.     WriteLn;
  268.     LinFit(x,y, i1,i2, a,b,r,da,db,dy,xo,xm,ym,sx,sy);
  269.     phi := ArcTan(b) * 180 / Pi;
  270.     WriteLn('Routine LINFIT',#13#10);
  271.     WriteLn('Steigung:                             b  = ',
  272.              b:12:6,' ± ',db:12:6);
  273.     WriteLn('-                                          ',
  274.              phi:12:6,'°');
  275.     WriteLn('Achsenabschnitt:                      a  = ',
  276.              a:12:6,' ± ',da:12:6);
  277.     WriteLn('Regression:                           r  = ',
  278.              (r*100):12:6,' %');
  279.     WriteLn('Absolute Abweichung der y-Werte:      dy = ',
  280.              dy:12:6);
  281.     WriteLn('Schnittpunkt mit der x-Achse:         xo = ',
  282.              xo:12:6);
  283.     WriteLn('x-Werte: Mittel, Standardabweichung:  xm = ',
  284.              xm:12:6,' ± ',sx:12:6);
  285.     WriteLn('y-Werte: -       -                    ym = ',
  286.              ym:12:6,' ± ',sy:12:6);
  287.     WriteLn;
  288.     WriteLn('Routine ROBUSTLINFIT:                 Einen ',
  289.             'Moment Geduld...',#13#10);
  290.     RobustLinFit(x,y, i1,i2, a,b,da,db,dy);
  291.     phi := ArcTan(b) * 180 / Pi;
  292.     WriteLn('Steigung:                             b  = ',
  293.              b:12:6,' ± ',db:12:6);
  294.     WriteLn('-                                          ',
  295.              phi:12:6,'°');
  296.     WriteLn('Achsenabschnitt:                      a  = ',
  297.              a:12:6,' ± ',da:12:6);
  298.     WriteLn('Absolute Abweichung der y-Werte:      dy = ',
  299.              dy:12:6);
  300.     Writeln('Schnittpunkt mit der x-Achse:         xo = ',
  301.              (-a/b):12:6);
  302.     WriteLn;
  303.     WriteLn('Weiter mit RETURN');
  304.     ReadLn;
  305.   END;
  306.  
  307. {----------------------- Beispiele ------------------------}
  308. BEGIN
  309.   Randomize;                                  { Beispiel 1 }
  310.   xmin := -10;
  311.   xmax := +10;              { x-Intervall der Stützstellen }
  312.   b    := 0.75;
  313.   a    := -2;              { Steigung b, Achsenabschnitt a }
  314.   FOR i := 1 TO npt DO BEGIN   { Verrauschter Meßdatensatz }
  315.     x[i] := (i - 1) * (xmax - xmin) / npt + xmin;
  316.     y[i] := (b * x[i] + a) + spread * 2 * (Random-0.5);
  317.   END;
  318.   ProduceOutliers;                   { Ausreißer aufprägen }
  319.   i1 := 1;                                  { Anfangspunkt }
  320.   i2 := npt;                                    { Endpunkt }
  321.   LinearFit;
  322.                        { Beispiel 2 - Radioaktiver Zerfall }
  323.   xmin := 0;  xmax := 10;                      { Zeitachse }
  324.   b := 1.5;                                   { b = ln 2/τ }
  325.   a := 100;                                 { Anfangsmenge }
  326.   FOR i := 1 TO npt DO BEGIN   { Verrauschter Meßdatensatz }
  327.     x[i] := (i - 1) * (xmax - xmin) / npt + xmin;
  328.     y[i] := a * Exp(-b * x[i]);           { Zerfallsgesetz }
  329.     y[i] := y[i] + spread * 2 * (Random - 0.5);
  330.                                        { Rauschen aufgeben }
  331.   END;
  332.   ProduceOutliers;                   { Ausreißer aufprägen }
  333.   i1 := 1;                                  { Anfangspunkt }
  334.   i2 := npt;                                    { Endpunkt }
  335.   FOR i := i1 TO i2 DO y[i] := Log10(y[i]); { Vorskalieren }
  336.   LinearFit;
  337.   ClrScr;
  338.   WriteLn('Ausgleichsfunktion                     y = a` + ',
  339.           'EXP(b` x) ');
  340.   WriteLn;
  341.   WriteLn('Steigung:                  b` = 2.303 b  = ',
  342.           (Ln(10)*b):12:6);
  343.   WriteLn('Achsenabschnitt:           a` = 10^a     = ',
  344.           (Exp10(a)):12:6);
  345.   WriteLn('Halbwertszeit:             τ  = -ln 2/b` = ',
  346.           (-Ln(2)/Ln(10)/b):12:6);
  347.   WriteLn('Anfangsmenge:              No = a`       = ',
  348.           (Exp10(a)):12:6);
  349.   WriteLn;
  350. END.
  351. (* ------------------------------------------------------ *)
  352. (*                Ende von LINFIT.PAS                     *)
  353.