⟦1a34873a8⟧

TextFile

funcfit=algol list.yes index.no xref.no   
\f




funcfit

begin
comment
      Fitting of functionparameters to a set of points
      by least squares method.
      GOS:  december 1979.

   In the function   f(x1,x2, . . . xp)  certain parameters
pm(1), pm(2), . . pm(n) may be fitted. The function is defined
by the user in an external real procedure with formal parameters
(x,pm,dy,b) :
      x    (call value, array) The independent variables of a
                 point.
      pm   (call value, array) The values of the parameters.
      f    (return value, real) The value of the function calcu-
                 lated from x and pm.
      dy   (return value, array) The values of the partial deri-
                 vatives of y with respect to the parameters,
                 calculated from x and pm.
      b    (call value, boolean) If b = false, the dy-values
                 should not be evaluated.
The name of the translated external must be <:funclinda:>.

The input requirements of the main program are:
   1)  A textstring in < >, max. 72 characters.
   2)  The number of points, N.
   3)  The number of independent variables, p.
   4)  The number of parameters to be fitted, n.
   5)  The number of fixed parameters.
   6)  Start values of the parameters.
   7)  Indices of the fixed parameters according to their
       order of specification (point 6).
   8)  The N points, each specified by
           x1, x2, . . . xp,  y , sigma(y)
       with sigma(y) as the uncertainty of y.
   9)  The function used in fitting, or any other function deri-
       ved from the parameters, may be plotted (with 3 sigma con-
       fidence limits). If (at input 2) an N < 0 is specified
       funcfit stores the necessary data on disc for use by
       the drawing program, funcdraw. Funcdraw requieres the
       same external procedure funclinda during compilation
       as used with funcfit;
integer N, n, nt, nf, p, i, j, k, l, m, lines;
boolean B, continue, yzero, draw;
array head(1:12);
connectlso;
readhead(in,head,1); read(in,N,p,nt,nf);  n:= nt-nf;
draw:= N<0; N:= abs N; i:= 1;
write(out,<:<12>:>,nl,3,string head(increase(i)),nl,3);

begin
real S, s, h, y, Sold, Smn, fx, fy;
integer array no,cif,length,exp(1:nt), space(1:p+2), td(1:8);
array point(1:N,1:p+2),x(1:p),pm,dy,lo(1:nt),a(1:N,1:n+1),
      dr,dd(1:nf),A(1:if n>0 then n else 1,1:n+1),layout(1:p+2);
read(in,pm);
for i:=n+1 step 1 until nt do read(in,no(i)); j:= 0;
for i:=1 step 1 until nt do begin
   cif(i):= 0; continue:= true; k:= n;
   for k:=k+1 while k<=nt and continue do
         continue:= i<>no(k);
   if continue then begin j:=j+1; no(j):= i end
end;
read(in,point); continue:= n>0; S:= 1; Sold:= 0;
m:= 0;
for i:=1 step 1 until N do
  if point(i,p+2)<>0 then m:= m+1;
B:= m>n and n>0;
if m<n then goto slut;
if -,B then begin
  for i:=1 step 1 until N do
    if point(i,p+2)<>0 then point(i,p+2):= 1
end;

comment Evaluation of point layout starts here;

begin
real man; integer b,h,d,fn,fe,s,ex;
yzero:= true;
write(out,<:     b     h     d    fn     s    fe   space:>,nl,2);
for j:=1 step 1 until p+2 do begin
h:= b:= 1; d:= fn:= 0;
for i:=1 step 1 until N do begin
   man:= point(i,j); if man<>0 then begin
   if j=p+1 then yzero:= false;
   ex:= entier(ln(abs man)/ln10); man:= man*10**(-ex);
   if man<0 then fn:= 1; 
   man:= abs man; y:= man - round man; k:= 0;
   for k:=k+1 while abs y>10**(k-10) do begin
      man:= (man - entier man)*10; y:= man - round man end;
   if k>b then b:= k;
   s:= if ex<0 then 1 else ex+1;
   if s>h then h:= s;
   s:= if ex<0 then k-ex-1 else if s>k then 0 else k-s;
   if s>d then d:= s
end end i;
if j=p+1 then begin
  fn:= 1; Smn:= 10**(if d>0 then -d else h-b);
end;
comment b,h,d determine a layout of the form
        <<ddd.ddd0000> for column j;
if h+d>b+3 then begin
   comment exponent layout;
   s:= b-d-1; fe:= if s<0 then 1 else 0;
   s:= if abs s>9 or h-1>9 then 2 else 1;
   h:= 1; d:= b-1; space(j):= b+fn+s+fe+2;
end else begin
   s:= fe:= 0; space(j):= h+fn+(if d>0 then d+1 else 0)
end;
write(out,<<dddddd>,b,h,d,fn,s,fe,space(j),nl,2);
k:= b shift 4 add h shift 4 add d shift 4 add fn
      shift 2 add s shift 4 add fe;
layout(j):= 0.0 add (if j=p+1 then 127 else 124) shift 41 add k
end j;
end layout;

repeter:
if -, continue then begin
   j:= 1;
   write(out,<:<12>:>,nl,2,string head(increase(j)),nl,2)
end;
for i:=1 step 1 until nt do length(i):= 9;
l:= m:= 0; j:= k:= 1;
rep:  write(out,nl,2,<:Parameters::>);
repwr:  m:= m+1;
if m<>k or m=1 then begin
  y:= pm(m); i:= cif(m);
  exp(m):= if i>9 then 9 else i;
  h:= format(length(m),exp(m),y,lo(m));
  l:= length(m)+l
end else l:= length(m);
if l<=60 then begin
  write(out,string h,y);
  if m<nt then goto repwr
end else m:= m-1;
if Sold<>0 and B then begin
  write(out,nl,1,<:delta(pm) ::>);
  for i:=k step 1 until m do
  if i=no(j) and j<=n then begin
    write(out,string lo(i),sqrt(S*A(j,j))*10**(-exp(i)));
    if exp(i)<>0 then write(out,<:':>,<<d>,exp(i));
    j:= j+1
  end else write(out,false add 32,length(i)-5,<:fixed:>)
end;
k:= m+1; if k<=nt then goto rep;
if -,continue and B then begin array C(1:n,1:n);
   for i:=1 step 1 until n do
   for j:=1 step 1 until i do C(i,j):= A(i,j)*S;
   write(out,nl,3,<:Error-matrix::>);
   writemat(out,<<__-d.dddddd'-dd>,C,n,4);
   for i:=1 step 1 until n do begin
      dy(i):= sqrt(A(i,i)); C(i,i):= cif(no(i));
      for j:=i-1 step -1 until 1 do C(i,j):=
        if dy(i)=0 or dy(j)=0 then 0 else A(i,j)/dy(i)/dy(j)
   end;
   write(out,nl,1,<:Correlation-matrix::>);
   writemat(out,<<-dd.ddd>,C,n,10);
   for j:=1 step 1 until nf do dd(j):= dr(j):= 0;
end;
if -,continue or -,B then begin
   write(out,nl,if B then 1 else 3);
   l:= 4;
   for j:=1 step 1 until p do begin
      k:= (space(j)-1)//2;
      write(out,sp,k+l,<:x:>,<<d>,j);
      l:= space(j)-k+2 end;
   if yzero then write(out,sp,l,
            <:   deviation         delta      t:>,nl,3)
   else begin
     k:= (space(p+1)-2)//2;
     write(out,sp,k+l+2,<:obs.:>);
     l:= space(p+1)-k; k:= (space(p+2)-1)//2-2;
     write(out,sp,k+l,<:delta:>,sp,
           space(p+1)+space(p+2)-k-6,<:obs-calc.    t:>,nl,3);
end end;
S:= 0; m:= 0;
for j:=1 step 1 until 8 do td(j):= 0;
for i:=1 step 1 until N do begin
   h:= point(i,p+2);
   for j:=1 step 1 until p do x(j):= point(i,j);
   if h<>0 or -,continue or -,B then
   s:= point(i,p+1) - funclinda(x,pm,dy,true);
   if -,continue or -,B then begin real t;
      t:= 0;
      if B and Sold>0 then begin
      for j:=1 step 1 until n do begin y:= dy(no(j));
      for k:=j step 1 until n do
         t:= t+y*A(j,k)*dy(no(k))*(if j=k then 1 else 2)
      end;
      if h<>0 then begin t:= h*h-t;
        t:= if t<=0 then 9.9 else s/sqrt(Sold*t);
        j:= entier(t+t);
        j:= (if j<-4 then -4 else
             if j=-4 then -3 else
             if j= 3 then  2 else
             if j> 3 then  3 else j)+5;
        td(j):= td(j)+1;
        for j:=1 step 1 until nf do begin y:= dy(no(n+j))/h;
          dr(j):= dr(j)+y*s/h; dd(j):= dd(j)+y*y end;
      end else t:= sqrt(t*Sold)
      end;
      for j:=1 step 1 until p do
        write(out,string layout(j),x(j));
      if yzero then begin
        write(out,sp,5,<<-d.dddd'-dd>,s);
        if h=0 then write(out,sp,6,<<-d.d'-dd>,t)
               else write(out,sp,10-space(p+2),
                          string layout(p+2),h)
      end else begin
        if point(i,p+1)=0 and h=0 then begin
          k:= entier(ln(abs(s))/ln10)-entier(ln(t)/ln10)+1;
          j:= 0; fy:=format(j,k,-s,fx);
          write(out,string fy,
                sp,space(p+1)-j+6,-s,sp,space(p+2)-j+4,t)
        end else begin
          write(out,string layout(p+1),point(i,p+1));
          write(out,string layout(p+2),if h=0 then t else h)
        end;
        if point(i,p+1)<>0 then write(out,string layout(p+1),s);
      end;
      if h<>0 then write(out,<<    -d.d>,t);
      write(out,nl,1); forceout(out);
   end;
   if h<>0 then begin
   s:= s/h; m:= m+1;
   if continue then begin
      for j:=1 step 1 until n do a(m,j):= dy(no(j))/h;
      a(m,n+1):= s
   end;
   S:= S+ s*s
end end i;
k:= if B then (m-n) else 1; S:= S/k;
s:= sqrt(S); i:= 4; j:= 3; fx:= format(i,j,s,fy);
write(out,nl,2,<:Standard deviation:  :>,string fx,s,
      <:  +-:>,sqrt(S/(k+k)),nl,if continue then 4 else 2);
if continue then begin
  continue:= Sold=0 and S<>0;
  for j:=1 step 1 until n do begin
    y:= 0;
    for i:=1 step 1 until m do begin
      h:= abs a(i,j); if h>y then y:= h end;
    lo(j):= y
  end;
  ortho(a,A,m,n);
  y:= n/(if B then s else Smn);
  for j:=1 step 1 until n do begin
    k:= no(j);  h:= lo(j)*y;
    pm(k):= pm(k)+A(j,n+1);
    continue:= continue or h*abs A(j,n+1) > 1;
    i:= if abs pm(k)=0 or h=0 then 0 else
        entier(ln(abs pm(k))/ln10)
       -entier(-ln(h*2)/ln10)+1;
    cif(k):= if i<1 then 1 else i;
  end;
  continue:= (continue or (B and Sold>S+S)) and (Sold>S or Sold=0);
  Sold:= S;
  goto repeter
end;
if B then begin
  write(out,
    sp,17,<:-2      -1     -0.5      0      0.5      1       2:>,
    nl,1,<:t-dist::>);
  for j:=1 step 1 until 8 do write(out,<<dddddddd>,td(j));
  write(out,nl,1,<:theor.:  :>); fx:= 2; fy:= 0; h:= 0.5;
  for j:=8 step -1 until 6 do begin
    y:= student(fx,fy,k,'-4); h:= h-y;
    fy:= fx; fx:= fx-(if j=8 then 1 else 0.5);
    td(j):= td(9-j):= y*m*10 end;
  td(4):= td(5):= h*m*10;
  for j:=1 step 1 until 8 do write(out,<<dddddd.d>,td(j)/10);
  write(out,nl,1,<:
Significance estimate for fixed constants::>,nl,1);
  for j:=1 step 1 until nf do begin
    if j mod 4 = 1 then write(out,nl,1);
    write(out,<<ddddd>,no(n+j),<:::>,
          << -d.dddd'-d>,dr(j)/sqrt(dd(j)*S*(m-n)));
end end;
write(out,nl,1);

if draw then begin
integer array tail(1:10);
zone res(128,1,stderror);
l:= 14+nt+(n+3)//4+n*n+(p+1)*N;
tail(1):= (l+127)//128;
for i:=2 step 1 until 10 do tail(i):= 0;
open(res,4,<:fdraw:>,0);
if monitor(42,res,0,tail)=0 then begin
  write(out,<:<10>***fdraw reserved<10>:>);
  goto slut
end;
monitor(40,res,0,tail); monitor(50,res,15,tail);
outrec(res,if l>128 then 128 else l);
for i:=1 step 1 until 12 do res(i):= head(i);
res(13):= 0.0 add nt shift  8 add n shift 8
              add p  shift 24 add N;
res(14):= S;
for i:=1 step 1 until n do begin
  k:= (i+3)//4+14;
  res(k):= res(k) shift 12 add no(i)
end;
for i:=1 step 1 until nt do res(k+i):= pm(i);
k:= k+nt; p:= p+1;
for i:=1 step 1 until n do
for j:=1 step 1 until n do begin
  if k=128 then begin
    l:= l-128; k:= 0;
    outrec(res,if l>128 then 128 else l)
  end;
  k:= k+1; res(k):= A(i,j)
end;
for i:=1 step 1 until N do
for j:=1 step 1 until p do begin
  if k=128 then begin
    l:= l-128; k:= 0;
    outrec(res,if l>128 then 128 else l)
  end;
  k:= k+1; res(k):= point(i,j)
end;
close(res,true);
end disc-lagring;
end arrays;

slut: end
▶EOF◀
DataMuseum.dk

RC4000/8000/9000

⟦1a34873a8⟧ TextFile

Derivation

TextFile
