[From Bill Powers (950814.1115 MDT)]
To Modellers:
I've added a successive-approximation routine to the COLLIER program
(COLLIER1) to iteratively determine both gain and reference level, minimizing
the RMS error between predicted and actual meal size. As Bruce guessed, the
reference level for the third rat was in the 20s (the statistical intercept
did not give a good estimate). Also now listed is the RMS prediction error,
and all the predicted and actual values of meal size.
The program takes longer to run now, since the successive approximations for
gain and reference level are run alternatively a total of 30 times to make
sure of reaching the actual final value (with fewer repeats it was sensitive
to the starting values). The output goes to COLLDATA.PAS so you can
view it easily with TP-7.0. Some encouragement is shown on the screen as the
program runs.
I have truncated the number of data points used for rat 3 to 10. The last two
points are anomalous. Using only the first 10 points (to a ratio of 1280),
the correlation is over 0.98.
Here is the output of the program:
···
=======================================================================
RAT RATIO WEIGHT INTAKE MEALRATE MEALSIZE
1 1.0 378 20.00 8.90 2.25
1 5.0 423 22.00 7.90 2.78
1 10.0 405 19.90 6.80 2.93
1 20.0 416 20.90 7.50 2.79
1 40.0 421 18.80 5.20 3.62
1 80.0 425 19.80 5.80 3.41
1 160.0 419 22.90 4.60 4.98
1 320.0 430 17.70 3.40 5.21
1 640.0 430 14.80 2.00 7.40
1 1280.0 416 15.80 1.20 13.17
1 2560.0 397 12.80 1.20 10.67
1 5120.0 402 9.80 0.60 16.33
2 1.0 404 23.20 9.10 2.55
2 5.0 423 20.90 8.30 2.52
2 10.0 432 20.90 7.00 2.99
2 20.0 442 20.90 6.40 3.27
2 40.0 455 20.90 5.40 3.87
2 80.0 452 21.00 4.20 5.00
2 160.0 427 17.90 1.60 11.19
2 320.0 401 12.90 1.00 12.90
3 1.0 334 17.70 11.40 1.55
3 5.0 330 22.90 9.70 2.36
3 10.0 347 16.70 9.30 1.80
3 20.0 356 15.80 5.80 2.72
3 40.0 366 17.60 6.00 2.93
3 80.0 377 17.70 4.80 3.69
3 160.0 377 22.10 6.80 3.25
3 320.0 382 15.90 3.80 4.18
3 640.0 387 13.80 2.20 6.27
3 1280.0 371 14.80 2.20 6.73
3 2560.0 382 16.90 1.40 12.07
3 5120.0 351 8.80 1.10 8.00
RAT 1
Corr, meal size vs intake = -0.8741 Sigma = 4.46
Mealsize correlation, real vs model = 0.9899
Ref = 23.00 Gain = 1.26
RMS prediction error = 0.66 N = 12
MEAL SIZE:
PREDICTED OBSERVED
2.37 2.25
2.64 2.78
3.03 2.93
2.77 2.79
3.84 3.62
3.49 3.41
4.26 4.98
5.48 5.21
8.23 7.40
11.53 13.17
11.53 10.67
16.48 16.33
RAT 2
Corr, meal size vs intake = -0.9137 Sigma = 3.86
Mealsize correlation, real vs model = 0.9943
Ref = 25.68 Gain = 1.09
RMS prediction error = 0.46 N = 8
MEAL SIZE:
PREDICTED OBSERVED
2.56 2.55
2.79 2.52
3.24 2.99
3.51 3.27
4.07 3.87
5.02 5.00
10.20 11.19
13.40 12.90
RAT 3
Corr, meal size vs intake = -0.5430 Sigma = 1.66
Mealsize correlation, real vs model = 0.9849
Ref = 21.73 Gain = 0.85
RMS prediction error = 0.30 N = 10
MEAL SIZE:
PREDICTED OBSERVED
1.73 1.55
2.00 2.36
2.07 1.80
3.11 2.72
3.03 2.93
3.64 3.69
2.72 3.25
4.37 4.18
6.43 6.27
6.43 6.73
The revised program follows at the end.
Best,
Bill P.
program Coll1a;
uses dos, crt, printer,stats;
type colltype = record
rat: integer;
weight,ratio,intake,mealrate,mealsize,pressrate: real;
end;
var xvar,yvar,zvar: array[0..11] of real;
coll: array[0..100] of colltype;
g,r: real;
p,o: array[0..11] of real;
i,j,numlines: integer;
dummy: string;
colldata: text;
outfile: text;
ch: char;
delta, bestg, bestr,sumsqr,lastsum: real;
procedure csys(n,i: integer);
var k: real;
j,h: integer;
begin
h := i - n;
o[h] := 0;
for j := 1 to 2 do
begin
k := coll[i].mealrate;
p[h] := k*o[h];
o[h] := o[h] + (g*(r - p[h]) - o[h])/(1 + g*k);
end;
end;
procedure matchmodel(ratnum: integer);
var i,k,n1,n2: integer;
begin
case ratnum of
1: begin n1 := 0; n2 := 11; end;
2: begin n1 := 12; n2 := 19; end;
3: begin n1 := 20; n2 := 29; end;
end;
for i := n1 to n2 do
with coll[i] do
begin
yvar[i - n1] := mealsize;
xvar[i - n1] := intake;
end;
writeln(outfile,'RAT ',ratnum);
correl('r',@xvar,@yvar,n2 - n1 + 1);
writeln(outfile,'Corr, meal size vs intake = ',corr:7:4, ' Sigma =
',sigy:5:2);
g := -regression;
r := 30;
for k := 1 to 30 do
begin
lastsum := 1e38;
delta := -0.1;
repeat
for i := n1 to n2 do
with coll[i] do csys(n1,i);
for i := 0 to n2 - n1 do
begin
xvar[i] := o[i];
yvar[i] := coll[i + n1].mealsize;
end;
correl('r',@xvar,@yvar,n2 - n1 + 1);
sumsqr := 0.0;
for i := 0 to n2 - n1 do
sumsqr := sumsqr + sqr(yvar[i] - xvar[i]);
if sumsqr <= lastsum then
begin
bestg := g;
end
else
delta := - delta/2.0;
lastsum := sumsqr;
g := g + delta;
until abs(delta) < 0.0005;
g := bestg;
delta := -5;
lastsum := 1e38;
repeat
for i := n1 to n2 do
with coll[i] do csys(n1,i);
for i := 0 to n2 - n1 do
begin
xvar[i] := o[i];
yvar[i] := coll[i + n1].mealsize;
end;
correl('r',@xvar,@yvar,n2 - n1 + 1);
sumsqr := 0.0;
for i := 0 to n2 - n1 do
sumsqr := sumsqr + sqr(yvar[i] - xvar[i]);
if sumsqr <= lastsum then
begin
bestr := r;
end
else
delta := - delta/2.0;
lastsum := sumsqr;
r := r + delta;
until abs(delta) < 0.005;
r := bestr;
end;
for i := n1 to n2 do
with coll[i] do csys(n1,i);
for i := 0 to n2 - n1 do
begin
xvar[i] := o[i];
yvar[i] := coll[i + n1].mealsize;
end;
correl('r',@xvar,@yvar,n2 - n1 + 1);
writeln(outfile,'Mealsize correlation, real vs model = ',corr:7:4);
writeln(outfile,'Ref = ',r:5:2,' Gain = ',g:5:2);
sumsqr := sqrt(sumsqr/(n2 - n1));
writeln(outfile,' RMS prediction error = ',sumsqr:4:2, ' N = ',(n2 - n1 +
1):3);
writeln(outfile);
writeln(outfile,' MEAL SIZE:');
writeln(outfile,'PREDICTED OBSERVED');
for i := 0 to n2 - n1 do
writeln(outfile,xvar[i]:6:2,' ',yvar[i]:6:2);
writeln(outfile);
end;