[From Bill Powers (950729.1150 MDT)]
To Modelers:
Appended below is the source code for CYCLIC2.PAS, the program for fitting a
model with explicit collection time to the data take by Staddon and reported
by Bruce Abbott. The program creates a plot of the original Staddon data (
food delivery rate vs. mean behavior rate, red circles), the within-ratio
pressing rate calculated from the assumed collection time (food delivery rate
vs. pressing rate, green circles), the apparatus response to a given pressing
rate alternating with a collection time (white curves), and a model of the
rat with an assumed reference level, output gain, and error limit (yellow
curve).
This is a graphical method for solving the control system equations. The
apparatus equation is
r = 1/(m/p + c) where
r = rate of delivery of food pellets per hour
m = ratio, presses per pellet
p = pressing rate per hour, within ratio (excludes collection time)
c = collection time, seconds.
In this equation, p is an independent variable and r is a dependent variable.
p, however, is plotted on the vertical axis and r on the horizontal axis.
The model of the rat (the "rat function") is
error := ref - r where
ref = reference level for food delivery rate
r = actual food delivery rate
error = output of comparator
A nonlinear error-output function is used, of the form
p = ErrLimit * (1 - exp(-gain*error/ErrLimit)) where
ErrLimit = constant determining maximum possible output
gain = ratio of output to error near zero error
p = output pressing rate between collection pauses
This function is plotted over the plot for the apparatus. Since we consider
only the case for r < ref, the curve begins at zero when the food delivery
rate is at the reference level, and rises to the left as the actual food
delivery rate falls below the reference level.
Now the independent variable is the rate of food delivery (horizontal axis as
before) and the dependent variable is the rate of pressing (vertical axis).
Where the rat-function curve crosses the apparatus curves for each value of m
are the solutions of the whole system equation. These intersection points,
when the four adjustable parameters are optimized, should lie on the green
circles representing the corrected observations.
This form of the error function is arbitrary, chosen for a good fit with the
data when all parameters are adjusted for best fit. When data are available
for the actual collection time and the actual mean rate of pressing during
consecutive bar presses, the collection time and assumed error limit can be
eliminated as parameters, leaving only the reference level and the output
gain to be determined from the data.
···
---------------------------------------
Each rat's data and rat-function are plotted in a separate frame. At the
right side of the screen is a box containing the values of each of the four
adjustable parameters for each rat: collection time, reference level, gain,
and error limit. The up/down arrows keys position a pointer and the + and -
keys increase and decrease the indicated parameter. The Enter key causes the
data to be replotted. Because of a limitation in the "setparams" unit, all
four rat plots are replotted, even for the plotting frames that have not been
affected.
The collection time is included as a rat parameter, but it is an unmodeled
parameter. It affects primarily the form of the apparatus function and the
placement of the corrected data points, the green circles. The reason is that
the independent variable for the rat function is the rat's actual rate of
pressing, which does not take the collection time into account. The
correction of the observed data and the apparent form of the apparatus
function depend on the assumed collection time. In a simulation of the rat,
as opposed to this graphical solution, the collection time would be
introduced as an explicit pause in pressing, as in OPCOND5.PAS, and mean
behavior rates and food delivery rates could be compared directly against
observations.
---------------------------------------
The present model does not take the "value" of the food pellets into account.
What is actually controlled by the rat, presumably, is the product of food
delivery rate and the value of each delivery. If food value is varied, this
is equivalent to using value*rate (or some function of that product) for the
horizontal axis; the remainder of the model would be unchanged.
---------------------------------------
The model is initialized to the best values of collection time, gain,
reference level, and limit that I could achieve for each plot by eye. It
would be very useful to automate the matching of the model to the data by
adjustments of the parameters. Simply adding a display of the RMS error
between model and data would aid in manually adjusting the parameters for
best fit.
-----------------------------------------------------------------------
Best,
Bill P.
program cyclic2;
uses dos,crt,graph,grutils,setparam,frameplt;
const ratio: array[0..5] of real = (2.0,4.0,8.0,16.0,32.0,64.0);
type dataarray = array[1..4] of array[0..5] of real;
csystype = record
ref,gain,error,limit,output: real;
end;
var r,p,b,g,m,ref: real;
c: array[1..4] of real;
csys: array[1..4] of csystype;
i: integer;
rat: integer;
param: paramlisttype;
frame: array[1..4] of frametype;
maxx,maxy: integer;
ch: char;
b_per_hr,p_per_hr: dataarray;
procedure loadparams;
begin
with param[1] do
begin
legend := 'coll time 1';
kind := 'r';
rvinit := 4.4;
rvmin := 0.0;
rvmax := 100.0;
rvstep := 0.1;
rv := @c[1] ;
end;
with param[2] do
begin
legend := 'Ref level 1';
kind := 'r';
rvinit := 576.0;
rvmin := 0.0;
rvmax := 1500.0;
rvstep := 2.0;
rv := @csys[1].ref;
end;
with param[3] do
begin
legend := 'gain 1';
kind := 'r';
rvinit := 50.0;
rvmin := 0.0;
rvmax := 10000.0;
rvstep := 1.0;
rv := @csys[1].gain;
end;
with param[4] do
begin
legend := 'Err Limit 1';
kind := 'r';
rvinit := 5200.0;
rvmin := 0.0;
rvmax := 12000.0;
rvstep := 50.0;
rv := @csys[1].limit
end;
with param[5] do
begin
legend := 'coll time 2';
kind := 'r';
rvinit := 4.4;
rvmin := 0.0;
rvmax := 100.0;
rvstep := 0.1;
rv := @c[2] ;
end;
with param[6] do
begin
legend := 'Ref level 2';
kind := 'r';
rvinit := 584;
rvmin := 0.0;
rvmax := 1500.0;
rvstep := 2.0;
rv := @csys[2].ref;
end;
with param[7] do
begin
legend := 'gain 2';
kind := 'r';
rvinit := 59.0;
rvmin := 0.0;
rvmax := 10000.0;
rvstep := 1.0;
rv := @csys[2].gain;
end;
with param[8] do
begin
legend := 'Err Limit 2';
kind := 'r';
rvinit := 6250.0;
rvmin := 0.0;
rvmax := 12000.0;
rvstep := 50.0;
rv := @csys[2].limit
end;
with param[9] do
begin
legend := 'coll time 3';
kind := 'r';
rvinit := 3.9;
rvmin := 0.0;
rvmax := 100.0;
rvstep := 0.1;
rv := @c[3] ;
end;
with param[10] do
begin
legend := 'Ref level 3';
kind := 'r';
rvinit := 596.0;
rvmin := 0.0;
rvmax := 1500.0;
rvstep := 2.0;
rv := @csys[3].ref;
end;
with param[11] do
begin
legend := 'gain 3';
kind := 'r';
rvinit := 51.0;
rvmin := 0.0;
rvmax := 10000.0;
rvstep := 1.0;
rv := @csys[3].gain;
end;
with param[12] do
begin
legend := 'Err Limit 3';
kind := 'r';
rvinit := 7850.0;
rvmin := 0.0;
rvmax := 12000.0;
rvstep := 50.0;
rv := @csys[3].limit
end;
with param[13] do
begin
legend := 'coll time 4';
kind := 'r';
rvinit := 4.3;
rvmin := 0.0;
rvmax := 100.0;
rvstep := 0.1;
rv := @c[4] ;
end;
with param[14] do
begin
legend := 'Ref level 4';
kind := 'r';
rvinit := 600.0;
rvmin := 0.0;
rvmax := 1500.0;
rvstep := 2.0;
rv := @csys[4].ref;
end;
with param[15] do
begin
legend := 'gain 4';
kind := 'r';
rvinit := 103.0;
rvmin := 0.0;
rvmax := 10000.0;
rvstep := 1.0;
rv := @csys[4].gain;
end;
with param[16] do
begin
legend := 'Err Limit 4';
kind := 'r';
rvinit := 8700.0;
rvmin := 0.0;
rvmax := 12000.0;
rvstep := 50.0;
rv := @csys[4].limit
end;
end;
procedure loadframes;
begin
with frame[1] do
begin
numyvars := 1; mx := maxx; my := maxy;
xbase := 20; ybase := 270;
xsize := 180; ysize := 180;
numxgrid := 20; numygrid := 60;
xzero := 0; yzero := 0;
xmax := 800.0;
ymax[1] := 12000;
ylegend[1] := 'Press Rate' ;
xlegend := 'Deliv Rate' ;
color[1] := white;
yvar[1] := @b ;
xvar := @r ;
end;
with frame[2] do
begin
numyvars := 1; mx := maxx; my := maxy;
xbase := 250; ybase := 270;
xsize := 180; ysize := 180;
numxgrid := 20; numygrid := 60;
xzero := 0; yzero := 0;
xmax := 800.0;
ymax[1] := 12000;
ylegend[1] := 'Press Rate' ;
xlegend := 'Deliv Rate' ;
color[1] := white;
yvar[1] := @b ;
xvar := @r ;
end;
with frame[3] do
begin
numyvars := 1; mx := maxx; my := maxy;
xbase := 20; ybase := 30;
xsize := 180; ysize := 180;
numxgrid := 20; numygrid := 60;
xzero := 0; yzero := 0;
xmax := 800.0;
ymax[1] := 12000;
ylegend[1] := 'Press Rate' ;
xlegend := 'Deliv Rate' ;
color[1] := white;
yvar[1] := @b ;
xvar := @r ;
end;
with frame[4] do
begin
numyvars := 1; mx := maxx; my := maxy;
xbase := 250; ybase := 30;
xsize := 180; ysize := 180;
numxgrid := 20; numygrid := 60;
xzero := 0; yzero := 0;
xmax := 800.0;
ymax[1] := 12000;
ylegend[1] := 'Press Rate' ;
xlegend := 'Deliv Rate' ;
color[1] := white;
yvar[1] := @b ;
xvar := @r ;
end;
end;
procedure loadtable;
begin
b_per_hr[1,0] := 1007; b_per_hr[2,0] := 1050;
b_per_hr[3,0] := 1088; b_per_hr[4,0] := 1125;
b_per_hr[1,1] := 1767; b_per_hr[2,1] := 1825;
b_per_hr[3,1] := 1887; b_per_hr[4,1] := 2022;
b_per_hr[1,2] := 2711; b_per_hr[2,2] := 3026;
b_per_hr[3,2] := 3120; b_per_hr[4,2] := 3522;
b_per_hr[1,3] := 3660; b_per_hr[2,3] := 4146;
b_per_hr[3,3] := 4723; b_per_hr[4,3] := 5243;
b_per_hr[1,4] := 4255; b_per_hr[2,4] := 4820;
b_per_hr[3,4] := 5821; b_per_hr[4,4] := 6486;
b_per_hr[1,5] := 4711; b_per_hr[2,5] := 5893;
b_per_hr[3,5] := 6575; b_per_hr[4,5] := 7481;
end;
procedure control(rat,v: integer);
var i: integer;
u,p: real;
x,y: integer;
begin
for i := 0 to round(frame[rat].xmax) do
with csys[rat] do
begin
r := i;
error := ref - r;
if error < 0.0 then error := 0.0;
output := limit*(1 - exp(-gain*error/limit));
b := output;
with frame[rat] do
begin
u := output * yscale[v];
y := my - ybase - yzero - round(u);
x := xbase + xzero + round(r*xscale);
if (y > my - ybase - ysize) and (y < my - ybase) then
putpixel(x,y,yellow);
end
end;
end;
procedure placepoint(rat,v: integer);
var u: real;
x,y: integer;
begin
with frame[rat] do
begin
u := yvar[v]^ * yscale[v];
y := my - ybase - yzero - round(u);
x := xbase + xzero + round(xvar^*xscale);
if (y > my - ybase - ysize) and (y < my - ybase) then
circle(x,y,3);
end
end;
procedure showratdata(rat: integer);
var m: real;
i: integer;
begin
for i := 0 to 5 do
begin
m := ratio[i];
r := b_per_hr[rat,i]/m;
b := b_per_hr[rat,i]; { apparent behavior rate }
setcolor(lightred);
placepoint(rat,1);
b := m/(m/b_per_hr[rat,i] - c[rat]/3600.0);
setcolor(lightgreen);
placepoint(rat,1);
end;
end;
procedure showformula(rat: integer);
var m: real;
j,i: integer;
begin
for i := 0 to 5 do
begin
m := ratio[i];
for j := 1 to 200 do
begin
p := 60*j;
b := p;
r := 1.0/(m/p + c[rat]/3600.0);
plotvar(frame[rat]);
end;
end;
end;
begin
clrscr;
initgraphics;
maxx := getmaxx; maxy := getmaxy;
loadparams;
setupparam(maxx - 200,150,16,param);
loadframes;
for rat := 1 to 4 do
initframe(frame[rat]);
outtextxy(0,225,' RAT 1 RAT 2');
outtextxy(0,255,' RAT 3 RAT 4');
outtextxy(maxx - 170,maxy - 90,'"q" or Esc to quit');
outtextxy(maxx - 170,60,'Up-Dn arrow to select');
outtextxy(maxx - 170,100,'Return key to replot');
loadtable;
repeat
for rat := 1 to 4 do
begin
showratdata(rat);
showformula(rat);
control(rat,1);
end;
repeat
ch := changeparam(param);
until (ch = chr(13)) or (ch = 'q') or (ch = chr(27));;
for i := 1 to 4 do clrplot(frame[i]);