From Bill Powers (950813.0925 MDT)]
To modellers:
Using the data Bruce Abbott supplied for the Collier et. al. experiments, I
have matched a simple control model to the data for each of three rats.
Below: first the data set, then the results for each rat.
···
======================================================================
RAT RATIO WEIGHT INTAKE MEALRATE MEALSIZE (intake/rate)
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.87
Output gain = 1.04
Ref = 24.92
Mealsize correlation, real vs model = 0.98
Intake correlation, real vs model = 0.98
RAT 2
Corr, meal size vs intake = -0.91
Output gain = 1.20
Ref = 29.31
Mealsize correlation, real vs model = 0.99
Intake correlation, real vs model = 0.99
RAT 3
Corr, meal size vs intake = -0.47
Output gain = 0.40
Ref = 11.36
Mealsize correlation, real vs model = 0.93
Intake correlation, real vs model = 0.93
The basic procedure was as follows for each rat.
First the meal size was calculated (last column) by dividing total intake per
day (INTAKE) by number of meals per day (MEALRATE).
Then the correlation, regression, and intercept were calculated for
{mealrate:intake}. The input-output gain of the rat was taken to be the
negative of the regression, and the reference level as the intercept.
Then, using the number of meals per day (mealrate) as a multiplicative
disturbance in the feedback path, the model was run for each value of the
ratio. The ratio itself was not a parameter. The intake and meal size
correlations are identical because the multiplying disturbance was always the
same for the model as it was for the real rat. However, note that the models
intake and meal size were not identical to the real ones; they are both
produced by the model.
Some improvement in the model:real correlations could be obtained by a
successive-approximation fit with output gain and reference level being
varied. The above results came from using the slope and intercept of the
input-to-output relationship directly measured from the data. This does not
generally produce the best possible fit.
Here is the model:
Ref signal
+|
-------------->[COMP] ----error--->---
> - |
[x1] [Gain]
> >
> v
intake <-----[multiplier]<-----------meal size
^
>
>
----<---- meal rate
(multiplicative
disturbance)
Note that the loop gain is the output gain times the meal rate. The meal rate
ranges from 11 per day down to 0.6 per day, so the quality of control is
reasonably high at the low ratios and becomes very poor at the high ratios.
Note also that the control-system model produces a much higher correlation
with the actual output than the direct measurement of input-output regression
and slope do. This is because the control system model varies both the input
and the output of the system.
Below is the program for producing the above results. It reads in a data file
(which is a text file that can be created by copying the table above to a
file named COLLDATA after removing the blank lines). It outputs the table and
the results to a text file called COLLDATA.OUT. The program is specific to
this data set; to generalize it for other data sets would require some
rewriting.
----------------------------------------------------------------------
Best,
Bill P.
program Coll1;
uses dos, crt, printer,stats;
type colltype = record
rat: integer;
weight,ratio,intake,mealrate,mealsize,pressrate: real;
end;
var xvar,yvar: 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;
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,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 := 31; 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[0],@yvar[0],n2 - n1 + 1);
writeln(outfile,'Corr, meal size vs intake = ',corr:5:2);
writeln(outfile,' gain = ',-regression:5:2);
writeln(outfile,' Ref = ',intercept:5:2);
r := intercept;
g := -regression;
for i := n1 to n2 do
with coll[i] do
csys(n1,i);
for i := 0 to n2 - n1 + 1 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:5:2);
for i := 0 to n2 - n1 + 1 do
begin
xvar[i] := p[i];
yvar[i] := coll[i+ n1].intake;
end;
writeln(outfile,'Intake correlation, real vs model = ',corr:5:2);
writeln(outfile);
end;
begin
clrscr;
assign(colldata,'colldata');
reset(colldata);
assign(outfile,'colldata.out');
rewrite(outfile);
numlines := 0;
i := 0;
while not eof(colldata) do
begin
with coll[i] do
begin
readln(colldata,rat,ratio,weight,intake,mealrate);
mealsize := intake/mealrate;
end;
inc(i);
end;
close(colldata);
numlines := i;
writeln(outfile,
'RAT RATIO WEIGHT INTAKE MEALRATE MEALSIZE');
for i := 0 to numlines-1 do
with coll[i] do
writeln(outfile,
rat:3,ratio:9:1,weight:9:0,intake:9:2,mealrate:9:2,mealsize:9:2);
writeln(outfile);
matchmodel(1);
matchmodel(2);
matchmodel(3);
close(outfile);
end.
Here is the STATS.PAS unit that this program uses:
{$N+}
unit stats;
interface
var sx,sy,sx2,sy2,sxy,xbar,ybar,sigx,sigy,corr,regression,intercept: real;
procedure correl(ty: char; x,y: pointer ; datasize: integer);
{note: requires first argument of 'r' or 'i' to identify data type}
implementation
type dataarraytype = array[0..2047] of integer;
rdataarraytype = array[0..2047] of real;
dataptrtype = ^dataarraytype;
rdataptrtype = ^rdataarraytype;
procedure correl;
var n,u,v,w,z: real;
i: integer;
dxptr,dyptr: dataptrtype;
rdxptr,rdyptr: rdataptrtype;
begin
if ty = 'i' then
begin
dxptr := x; dyptr := y;
end
else
begin
rdxptr := x; rdyptr := y;
end;
sx := 0.0; sy := 0.0;
n := datasize;
for i := 0 to datasize - 1 do
begin
if ty = 'i' then
begin
u := dxptr^[i]; v := dyptr^[i];
end
else
begin
u := rdxptr^[i]; v := rdyptr^[i];
end;
sx := sx + u; sy := sy + v;
end;
xbar := sx/n; ybar := sy/n;
sx2 := 0.0; sy2 := 0.0; sxy := 0.0;
for i := 0 to datasize - 1 do
begin
if ty = 'i' then
begin
u := dxptr^[i]; v := dyptr^[i];
end
else
begin
u := rdxptr^[i]; v := rdyptr^[i];
end;
sx2 := sx2 + (u - xbar)*(u - xbar);
sy2 := sy2 + (v - ybar)*(v - ybar);
sxy := sxy + (u - xbar)*(v - ybar);
end;
sigx :=sqrt(sx2/n);
sigy :=sqrt(sy2/n);
if (abs(sigx*sigy) > 0.0001) then
z := sxy/(n * sigx * sigy)
else z := 0.0;
corr := z;
if abs(sigx) > 0.001 then regression := z * sigy/sigx
else regression := 9999.99;
intercept := ybar - regression*xbar;
end;
end.