home *** CD-ROM | disk | FTP | other *** search
/ Turbo Toolbox / Turbo_Toolbox.iso / spezial / 15 / graphen / reg / roblinfi.pas < prev    next >
Encoding:
Pascal/Delphi Source File  |  1989-05-18  |  3.1 KB  |  98 lines
  1. { --------------------------------------------------------------- }
  2. {                   File : ROBLINFI.PAS                           }
  3. {          Copyright (c) : 1989  Peter Kurzweil, H.H  & TOOLBOX   }
  4. {                Sprache : TURBO PASCAL 4.0  (oder 3.0)           }
  5. {              Ein Programm zur Auswertung von Meßwerten          }
  6. {                   mit Hilfe von GRAPH.PAS                       }
  7. { --------------------------------------------------------------- }
  8.  
  9. FUNCTION RobustLinFit ( VAR x,y : Vector; k1,k2 : word;
  10.                         VAR a,b,da,db,dy : real)    : Boolean ;
  11.                         
  12. LABEL 9;
  13. VAR i,k,n:word; xs,ys,xy,xx,f,f1,f2,b1,b2,d,chi2,h:real;
  14.  
  15. PROCEDURE HeapSort(VAR a:Vector; n:word);     { Zahlen aufsteigend sortieren }
  16. VAR l,j,k,i:word; h:real;
  17. BEGIN
  18.   l:=(n DIV 2)+1; k:=n;
  19.   WHILE true DO BEGIN
  20.     IF l>1 THEN BEGIN l:=l-1; h:=a[l]; END
  21.     ELSE BEGIN
  22.       h:=a[k]; a[k]:=a[1]; k:=k-1;
  23.       IF k=1 THEN BEGIN a[1]:=h; Exit; END
  24.     END;
  25.     i:=l; j:=l+l;
  26.     WHILE j<=k DO BEGIN
  27.       IF j<k THEN IF a[j]<a[j+1] THEN j:=j+1;
  28.       IF h<a[j] THEN BEGIN a[i]:=a[j]; i:=j; j:=j+j; END ELSE j:=k+1;
  29.     END;
  30.     a[i]:=h;
  31.   END;
  32. END;
  33.  
  34.   FUNCTION Theory(b:real; VAR a,dy:real):real;    { Rechte Seite von Gl. (1) }
  35.   VAR d,s:real; i,k,n,n1,n2:word; f:Vector;
  36.   BEGIN
  37.     n:=k2-k1+1; n2:=(n+1) DIV 2; n1:=(n+1)-n2; i:=0;
  38.     FOR k:=k1 to k2 DO BEGIN i:=i+1; f[i]:=y[k]-b*x[k]; end;
  39.     HeapSort(f,n);
  40.     a:=0.5*(f[n2]+f[n1]);                           { Bestimmung des Medians }
  41.     s:=0; i:=0; dy:=0;
  42.     FOR i:=k1 to k2 DO BEGIN
  43.       d:=y[i]-(b*x[i]+a);
  44.       dy:=dy+abs(d);
  45.       IF d>0 THEN s:=s+x[i] ELSE s:=s-x[i];
  46.     END;
  47.     Theory:=s;
  48.   END;
  49.  
  50. BEGIN
  51.   IF k1>k2 THEN BEGIN n:=k1; k1:=k2; k2:=n; END;
  52.   n:=k2-k1+1;                                            { Stuetzstellenzahl }
  53.   xs:=0; ys:=0; xy:=0; xx:=0; chi2:=0;
  54.  
  55.   FOR i:=k1 to k2 DO                                     { Gaußsummen }
  56.   BEGIN
  57.        xs := xs + x[i]     ; ys := ys + y[i];
  58.        xy := xy + x[i]*y[i]; xx := xx + sqr (x[i]);
  59.   END;
  60.  
  61.   d:=n*xx-xs*xs;
  62.  
  63.   IF d = 0 THEN                                   { nur ein Punkt eingegeben }
  64.   BEGIN
  65.        RobustLinFit := FALSE;
  66.        EXIT;
  67.   END;
  68.  
  69.   RobustLinFit := TRUE;
  70.   a:=(xx*ys-xs*xy)/d;
  71.   b:=(n*xy-xs*ys)/d;
  72.   FOR i:=k1 to k2 DO chi2:=chi2+sqr(y[i]-(a+b*x[i]));   { Fehlerquadratsumme }
  73.   h:=sqrt(chi2/d);
  74.   b1:=b;
  75.   f1:=Theory(b1,a,dy);
  76.  
  77.   IF f1>0 THEN b2:=b+abs(3*h) ELSE b2:=b-abs(3*h);
  78.   f2:=Theory(b2,a,dy);
  79.   WHILE ((f1*f2>0) and (f1<>f2)) DO BEGIN     { Loese Gl. (2) fuer versch. b }
  80.     b :=2*b2-b1;
  81.     b1:=b2; f1:=f2;
  82.     b2:=b;  f2:=Theory(b2,a,dy);
  83.   END;
  84.   h:=0.01*h;
  85.   WHILE abs(b2-b1)>h DO BEGIN                             { Nullstellensuche }
  86.     b:=0.5*(b1+b2);
  87.     IF ((b=b1) OR (b=b2)) THEN GOTO 9;
  88.     f:=Theory(b,a,dy);
  89.     IF f*f1>=0 THEN BEGIN f1:=f; b1:=b; END
  90.     ELSE BEGIN f2:=f; b2:=b; END;
  91.   END;
  92. 9:
  93.   dy:=dy/n;                                                { Fehler des Fits }
  94.   da:=dy*sqrt(abs(xx/(n*xx-xs*xs)));
  95.   db:=dy*sqrt(abs(n/(n*xx-xs*xs)));
  96. END;
  97.  
  98.