[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;