a 5273j.1,fmdekking

begin	comment fmd,corrcoef uit een tabel;
	integer nn,i,j,k,m,prob,n; real xmean,ymean,sxx,syy,sx,sy,sxsh,sysh,sxy,hx,hy;
	procedure ass(i,k,s,t); integer i,k; real s,t;
	  for i:=1 step 1 until k do s:=t;
	procedure p(s,a,b,x); string s; integer a,b; real x;
	  begin printtext(s); fixt(a,b,x); nlcr end;
	integer procedure ss;
	  begin integer x; ss:=x:=resym; prsym(x) end;

	nn:=119; for i:=1,2,3,4 do read1(ss,nn);
next:	prob:=read1(ss,nn); if prob=-1 then exit; k:=read1(ss,nn); m:=read1(ss,nn);
begin	integer array c,fx[1:k],d,fy[1:m],f[1:k,1:m]; real array xred[1:k],yred[1:m];
	ass(j,m,d[j],read1(ss,nn));
	for i:=1 step 1 until k do begin c[i]:=read1(ss,nn); ass(j,m,f[i,j],read1(ss,nn)) end;
	n:=sum(i,1,k,sum(j,1,m,f[i,j])); hx:=c[2]-c[1]; hy:=d[2]-d[1];
	ass(i,k,fx[i],sum(j,1,m,f[i,j])); ass(j,m,fy[j],sum(i,1,k,f[i,j]));
--------------------------------
 nlcr;
	printtext("x randtotalen"); for i:=1 step 1 until k do absfixt(4,0,fx[i]); nlcr;
	printtext("y randtotalen"); for j:=1 step 1 until m do absfixt(4,0,fy[j]); nlcr;
	xmean:=sum(i,1,k,c[i]*fx[i])/n; ymean:=sum(j,1,m,d[j]*fy[j])/n;
	ass(i,k,xred[i],c[i]-xmean); ass(j,m,yred[j],d[j]-ymean);
	sxx:=sum(i,1,k,fx[i]*xred[i]*xred[i])/n; syy:=sum(j,1,m,fy[j]*yred[j]*yred[j])/n;
	sx:=sqrt(sxx); sy:=sqrt(syy); sxsh:=sqrt(sxx-hx*hx/12); sysh:=sqrt(syy-hy*hy/12);
	sxy:=sum(i,1,k,sum(j,1,m,f[i,j]*xred[i]*yred[j]))/n;
	p("n",8,0,n); p("xmean",4,6,xmean); p("ymean",4,6,ymean);
	p("sxx",6,6,sxx); p("syy",6,6,syy); p("sxy",6,6,sxy);
--------------------------------
	p("sx",7,6,sx); p("sy",7,6,sy); p("sxshep",3,6,sxsh); p("syshep",3,6,sysh);
	p("corrco",3,6,sxy/(sx*sy)); p("corrsh",3,6,sxy/(sxsh*sysh))
end; newpage; goto next
end
1]7