[From Bill Powers (950310.1130 MST)]

Rick Marken, Bruce Abbott, Tom Bourbon, other modelers:

Rick, you're the champ right now. Thanks for the confirming data.

Appended below is a version of the SD3an analysis called SD3anp1 for Powers

modification 1. This version has two main new features.

First, the parameters of the tracking loop are calculated using data

starting some time after the condition switch, so the data during the delay

and switchover transient are excluded from the calculation. The model is

run as usual, but the squared error is accumulated only for the periods not

excluded. The k parameter is adjusted to minimize this squared error

instead of the squared error over the whole run.

This brought the RMS error for my second run down from 24 to 6 pixels, and

the correlation of model and real mouse position up from 0.978 to 0.999

(for the nonexcluded periods). The excluded period was set at 80 iterations

after the condition switch. The lag can still be changed as usual, because

the exclusion period is calculated independently and affects only the times

during which squared error is accumulated.

Temporary arrays temp1 and temp2 are used to store model and real hndle

positions so the correlation can be computed using only those values.

The value of k obtained from the above process is applicable to the

tracking process between switching transients. Since in the model the

handle position is

h = k*integral(error),

it follows that

error = (dh/dt)/k.

The perceptual signal is

p = h + disturbance, so

error = reference - (h + disturbance), and

reference = error + (h + disturbance), or

reference = (dh/dt)/k + (h + disturbance).

Once we have a measure of k from the initial fit of the model to the data

from the periods between switching transients, we can use the original

experimental handle positions to compute a table of reference-signal values

for the whole run (in the array calcref^).

This is the same method I used in the article "Quantitative measurement of

volition", Ch. 13 in Wayne Hershberger's book _Volitional Action: conation

and control_. (Amsterdam, North-Holland, 1989). Part of the data is

obtained under normal tracking conditions, with a period during the

experiment where a different condition holds. The model found from the

normal part of the data is used to find a value of the integration factor

k, and then the raw data plus this derived integration factor is used to

calculate the behavior of the reference signal during the whole run, under

the hypothesis that all tracking errors not accounted for by the basic

control model are due to variations in the reference signal.

This method was used in the article to recover a pattern of reference-

signal changes which the subject was instructed to carry out during the

middle third of a run. The pattern intended was a stairstep of reference

positions going upward and then downward. The stairstep was not visible in

the resulting data, either in the cursor positions or the handle positions,

because there was a large and fast enough disturbance to mask the pattern.

However, the calculation of the running values of the reference signal,

when smoothed, produced a very clear stairstep.

In the present application, this method is used to deduce the behavior of

the reference signal entering the control system that is doing the

tracking.

The variable "lag" is now located in the main begin..end at the bottom of

the program. Just below it, the excluded period "exclude" is calculated as

80 - lag. Thus no matter what value is used for "lag", the tracking data

counted for computing squared error starts 80 iterations after the change

of the "state^" variable. I picked 80 after a few quick trials; other

values can be explored to see their effects.

The ModelRun function is entered with the exclusion period as an argument

"ex". If the exclusion period is zero, the reference levels are computed as

before, depending on the changes in the state^ variable. Otherwise the

reference signal is taken from the calcref^ table.

All this, needless to say, needs to be checked by someone else to see if

the logic is what I think it is.

You run the program in the usual way. First you get the numerical results,

and then a graphic plot of the results as before. You can adjust "lag" to

measure the delay in the outer control loop. The value of "lag" does not

have any effect on computing the reference signal.

When you hit the escape key, the reference-signal table is calculated, and

then is plotted over the existing display. You will see that the computed

reference signal jumps very rapidly, although not instantly, to its new

value, with a pronounced overshoot on most transitions. It then comes to a

nearly constant value centered on the yellow line representing the ideal

setting of the reference signal. You will also see that this transition

begins after the ideal reference signal jumps to the other position, which

jump happens exactly at the moment the conditions change. So the delay can

be measured here, too.

This method asserts that all errors not accounted for by the basic tracking

model come from variations in the reference signal. When you look at the

model run using the array calcref as the reference signal, which happens

when you call the ModelRun function with an argument of 0, you will find a

fantastically good fit. This is because all errors, wherever they really

originate, are put into the reference signal. If all the calculations were

done with infinite accuracy, you would get an RMS error of exactly 0. The

present program does no calculations of model fit when the running

reference signal is used.

So to some extent the reference signal array probably includes some

variations that belong inside the tracking control system. However, when

you look at the plot of the derived reference signal, you will see that the

initial transient is probably real, because neither handle position nor

cursor position shows any overshoot of that size. The small noisy

fluctuations may or may not reflect noise in the reference signal; some of

it may be noise in the tracking control system. I don't see just now how to

determine where the variations really come from.

However, it's not a bad guess to say that the big overshoot in the

reference signal is a logical way of designing a system so its output will

produce a quick change when used to drive a system that is somewhat

sluggish. It's not a bad guess, either, to say that this plot of the

reference signal gives us a better picture of what the next level up is

doing than does the assumption of an instantaneous switch between exactly

the right levels, as assumed in the initial phase of the analysis. The

reference signal changes rapidly, yes. But it does not seem to change in

the way we would expect if the higher system were merely pointing to a

stationary target and saying "track that one." That would not involve an

overshoot.

Comments?

## ···

-----------------------------------------------------------------------

Best to all,

# Bill P.

program SD3anp1;

{ Fits a simple one-level PCT model to data from the SDTEST3.EXE

discriminated operant procedure and displays the results. The

program may be called from the command line by entering

SD3an filename.ext

or by entering the filename at the prompt after the program is

running. SDTEST3.EXE data files are named SDDATA.XXX, where XXX is

a three-digit number.

Written by Bruce Abbott

Psychological Sciences

Indiana University - Purdue University

Fort Wayne, IN 46805-1499

(219) 481-6399

abbott@cvax.ipfw.indiana.edu

Modified by WTP 950310 to find best tracking model using data between

switches, then deduce reference signal from that.

}

uses DOS, CRT, Graph, GrUtils;

const

MAXDATA = 3600;

type datalist = array[1..MAXDATA] of integer;

dataptr = ^datalist;

var

h, mh, d, state,temp1,temp2,calcref,dhdt: dataptr;

MaxX, MaxY, NPoints,exclude,lag: integer;

rms, k, trakdata: real;

Success: boolean;

ch: char;

procedure InitHeapVars;

begin

new(h);

new(mh);

new(state);

new(d);

new(temp1);

new(temp2);

new(calcref);

new(dhdt);

end;

procedure DisposeHeapVars;

begin

dispose(dhdt);

dispose(calcref);

dispose(temp2);

dispose(temp1);

dispose(d);

dispose(state);

dispose(mh);

dispose(h);

end;

procedure InitScreen;

begin

ClrScr;

InitGraphics;

MaxX := GetMaxX; MaxY := GetMaxY;

ClearViewPort;

end;

function FileExists(Filename: PathStr): Boolean;

var

TextFile: Text;

begin

{$I-}

Assign(TextFile, Filename);

Reset(TextFile);

Close(TextFile);

{$I+}

FileExists := (IOResult = 0);

end;

procedure ReadData(var FileRead: boolean);

var

Filename: Pathstr;

DataFile: Text;

i, j, result: integer;

begin

Filename := ParamStr(1);

if Filename = '' then

repeat

gotoXY(1, 5); Write('Enter Data Filename or type QUIT to exit: ');

Readln(Filename);

until Filename <> '';

for i := 1 to length(Filename) do

Filename[i] := upcase(filename[i]);

if Filename = 'QUIT' then

begin

FileRead := False;

writeln('QUIT: No files read or written');

Exit;

end;

gotoXY(1, 6); write('Filename = ', Filename);

if FileExists(Filename) then

begin

i := 0;

{$I-}

Assign(DataFile, Filename);

Reset(DataFile);

Result := IOResult;

{$I+}

if Result = 0 then

begin

gotoXY(1, 7); write('Reading ',Filename);

while (NOT EOF(DataFile)) do

begin

inc(i);

Readln(DataFile, h^[i], d^[i], state^[i]);

gotoXY(1, 8); write(i:5, ' data points read');

end;

Close(DataFile);

FileRead := true;

if i > MAXDATA then NPoints := MAXDATA else NPoints := i;

end

else

begin

FileRead := false;

gotoXY(1, 8); write('Unable to read data file...');

end;

end

else

begin

FileRead := false;

gotoXY(1, 8); write('Unable to find ', Filename);

end;

end;

procedure correlation(x, y: dataptr; n: word;

var a, b, rms, r, Xbar, Ybar: real);

var

i: integer;

sx, sy, ssx, ssy, sxy, sigx, sigy: real;

begin

sx := 0.0; sy := 0.0; ssx := 0.0; ssy := 0.0; sxy := 0.0;

for i :=1 to n do

begin

sx := sx + x^[i];

sy := sy + y^[i];

end;

xbar := sx/n;

ybar := sy/n;

for i :=1 to n do

begin

ssx := ssx + (x^[i] - xbar) * (x^[i] - xbar);

ssy := ssy + (y^[i] - ybar) * (y^[i] - ybar);

sxy := sxy + (x^[i] - xbar) * (y^[i] - ybar);

end;

sigx := sqrt(ssx/n);

sigy := sqrt(ssy/n);

r := sxy/(n*sigx*sigy);

b := sxy/ssx;

a := ybar - b*xbar;

rms := sqrt((1-r*r)*ssy/(n-2));

end;

function ModelRun(ex: integer): real;

var

i, j, oldstate, count, change: integer;

handle, newsumsq, error, p, temp,

ref, ref1, ref2: real;

begin

ref1 := 1.0*(MaxX div 3 - MaxX div 2);

ref2 := -ref1;

newsumsq := 0.0;

handle := h^[1];

if ex > 0 then

begin if state^[1] = 1 then ref := ref1 else ref := ref2; end

else

ref := calcref^[1];

oldstate := 0;

count := 0;

change := 0; trakdata := 0.0; j := 1;

for i := 1 to NPoints do

begin

p := handle + d^[i];

if oldstate <> state^[i] then

begin

count := ex;

oldstate := state^[i];

end;

if (i > lag) then

begin

if (ex > 0) then

begin if state^[i - lag] = 1 then ref := ref1 else ref := ref2; end

else

ref := calcref^[i];

end;

error := p - ref;

handle := handle - k*error;

mh^[i] := round(handle);

if count > 0 then dec(count)

else

begin

temp1^[j] := round(handle);

temp2^[j] := h^[i];

temp := 1.0*h^[i] - handle;

inc(j);

newsumsq := newsumsq + temp*temp;

trakdata := trakdata + 1.0;

end;

end;

modelrun := newsumsq/trakdata;

end;

function fitmodel: real;

var i, iter: integer;

delta, sumsq, savesum: real;

done: boolean;

begin

gotoXY(30, 1); write('FITTING MODEL TO DATA');

k := 0.5;

delta := -0.25;

iter := 0;

sumsq := 1e30;

done := false;

repeat

iter := iter + 1;

gotoxy(23, 3); write('k = ', k:10:6, ' delta = ', delta:10:6);

savesum := sumsq;

sumsq := modelrun(exclude + lag);

if abs(delta) < 0.0001 then done := true

else

begin

if sumsq > savesum then delta := -delta/2.0;

k := k + delta;

if k <= 0 then

begin

delta := - delta;

k := k + delta;

end;

end;

until done or (iter >= 50);

fitmodel := sqrt(sumsq);

writeln;

end;

procedure GraphData;

var i, ref1, ref2, temp: integer;

begin

setgraphmode(graphmode);

setcolor(lightred);

outtextxy(0,0,'CURSOR');

setcolor(yellow);

outtextxy(0,20,'REFERENCE');

setcolor(lightgreen);

outtextxy(0,40,'PARTICIPANT''S HANDLE');

setcolor(lightgray);

outtextxy(0,60,'MODEL''S HANDLE');

outtextxy(50, MaxY-50, 'Press ESCape key to quit...');

ref1 := MaxX div 2 - MaxX div 3;

ref2 := -ref1;

for i := 1 to NPoints do

begin

putpixel(i div 6,maxy div 2 + d^[i] + h^[i],lightred);

putpixel(i div 6,maxy div 2 - h^[i], lightgreen);

putpixel(i div 6,maxy div 2 - mh^[i], lightgray);

if state^[i] = 1 then temp := ref1 else temp := ref2;

putpixel(i div 6, maxy div 2 - temp, yellow);

putpixel(i div 6,maxy div 2,white);

end;

end;

procedure PCTmodel;

var

sumsq, rmsModel, a, b, r, rms, Xbar, Ybar: real;

begin

ClrScr;

rms := fitmodel;

writeln;

writeln('Value of k: ', k:7:4,' RMS prediction error, pixels:

',rms:3:2);

writeln('Number of evaluation points ',round(trakdata));

writeln;

Correlation(temp1, temp2, round(trakdata), a, b, rms, r, Xbar, Ybar);

writeln('Model vs Actual Handle Correlation = ', r:8:3);

writeln(' RMS error = ', rms:8:3);

writeln;

write(

'Fitting model handle. Enter value of k from a previous run? Y/N: ');

ch := readkey; writeln(ch);

if (ch = 'Y') or (ch = 'y') then

begin

write(' Enter value: '); readln(k);

end;

writeln;

write('Fitting model to data, please wait...');

sumsq := modelrun(exclude + lag);

rmsModel := sqrt(sumsq/trakdata);

GraphData;

end;

procedure showref;

var i: integer;

begin

for i := 1 to 600 do

putpixel(i,maxy div 2,white);

moveto(1,maxy div 2 + calcref^[1]);

setcolor(white);

for i := 2 to 600 do

begin

lineto(i,maxy div 2 + calcref^[i*6]);

end;

outtextxy(150,20,'CALCULATED REFERENCE SIGNAL (press esc to exit)');

end;

procedure getref;

var i: integer;

begin

fillchar(calcref^,sizeof(calcref^),0);

for i := 1 to NPoints do

calcref^[i] := round(dhdt^[i]/k) + (h^[i] + d^[i]);

end;

procedure makedhdt;

var i: integer;

begin

for i := 2 to NPoints - 1 do

dhdt^[i] := (h^[i+1] - h^[i-1]) div 2;

dhdt^[1] := dhdt^[2];

dhdt^[NPoints] := dhdt^[NPoints - 1];

end;

begin

InitScreen;

InitHeapVars;

RestoreCRTMode;

ReadData(Success);

Makedhdt;

lag := 40;

exclude := 80 - lag;

if Success then PCTmodel else halt;

repeat

ch := readkey;

until ch = #27;

getref;

ModelRun(0); { note: this run is not shown or analyzed}

showref;

repeat

ch := readkey;

until ch = #27;

DisposeHeapVars;

CloseGraph;

end.