SD3AN.PAS: Analysis Program for SDTEST3

[From Bruce Abbott (950307.19:40 EST)]

Below is SD3AN.PAS, an analysis program for data from the SDTEST3.EXE
discriminated operant experiment.

The program fits a one-level PCT model to the data, using the actual target
displacements as the reference levels of the low-level system controlling
cursor position. In reality, these reference values would come from
another, higher-level system controlling the perception of cursor over left
target when cursor is green, cursor over right target when cursor is red. A
yet higher-level system would be controlling the rate of point increase,
with reference at some high level. Or at least that seems plausible to me.

I'm not sure how to go about fitting multi-level models, so this one just
provides the lowest level and assumes that higher-level systems are setting
the reference positions to the correct values. I'll rely on our PCT experts
for advice on how to proceed further.

Meanwhile, here are the results of three runs of SDTEST3 by yours truly:

Run k rms r best
005 0.0283 45.68 0.919 45.38
006 0.0205 47.23 0.908 46.40
007 0.0186 54.88 0.903 53.79

k = integration constant for the model
rms = root mean square error, model vs. actual
r = Pearson r, model vs. actual
best = root mean square error, regression of model on actual

The correlations are low by PCT standards. This may reflect the simplified
nature of the model as well as the early stage of learning represented in
the data.




program SD3an;

{ 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

uses DOS, CRT, Graph, GrUtils;

     MAXDATA = 3600;

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

  h, mh, d, state: dataptr;
  MaxX, MaxY, NPoints: integer;
  rms, k: real;
  Success: boolean;
  ch: char;

procedure InitHeapVars;

procedure DisposeHeapVars;

procedure InitScreen;
  MaxX := GetMaxX; MaxY := GetMaxY;

function FileExists(Filename: PathStr): Boolean;
  TextFile: Text;
  Assign(TextFile, Filename);
  FileExists := (IOResult = 0);

procedure ReadData(var FileRead: boolean);
  Filename: Pathstr;
  DataFile: Text;
  i, j, result: integer;
  Filename := ParamStr(1);
  if Filename = '' then
      gotoXY(1, 5); Write('Enter Data Filename or type QUIT to exit: ');
    until Filename <> '';
  for i := 1 to length(Filename) do
    Filename[i] := upcase(filename[i]);
  if Filename = 'QUIT' then
      FileRead := False;
      writeln('QUIT: No files read or written');
  gotoXY(1, 6); write('Filename = ', Filename);
  if FileExists(Filename) then
    i := 0;
    Assign(DataFile, Filename);
    Result := IOResult;
    if Result = 0 then
        gotoXY(1, 7); write('Reading ',Filename);
        while (NOT EOF(DataFile)) do
            Readln(DataFile, h^[i], d^[i], state^[i]);
            gotoXY(1, 8); write(i:5, ' data points read');
        FileRead := true;
        if i > MAXDATA then NPoints := MAXDATA else NPoints := i;
        FileRead := false;
        gotoXY(1, 8); write('Unable to read data file...');
      FileRead := false;
      gotoXY(1, 8); write('Unable to find ', Filename);

procedure correlation(x, y: dataptr; n: word;
             var a, b, rms, r, Xbar, Ybar: real);
  i: integer;
  sx, sy, ssx, ssy, sxy, sigx, sigy: real;
sx := 0.0; sy := 0.0; ssx := 0.0; ssy := 0.0; sxy := 0.0;
for i :=1 to n do
   sx := sx + x^[i];
   sy := sy + y^[i];
  xbar := sx/n;
  ybar := sy/n;
for i :=1 to n do
   ssx := ssx + (x^[i] - xbar) * (x^[i] - xbar);
   ssy := ssy + (y^[i] - ybar) * (y^[i] - ybar);
   sxy := sxy + (x^[i] - xbar) * (y^[i] - ybar);
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));

function ModelRun: real;
  i: integer;
  handle, newsumsq, error, p, temp,
  ref, ref1, ref2: real;
  ref1 := 1.0*(MaxX div 3 - MaxX div 2);
  ref2 := -ref1;
  newsumsq := 0.0;
  handle := h^[1];
  for i := 1 to NPoints do
    p := handle + d^[i];
    if state^[i] = 1 then ref := ref1 else ref := ref2;
    error := p - ref;
    handle := handle - k*error;
    mh^[i] := round(handle);
    temp := 1.0*h^[i] - handle;
    newsumsq := newsumsq + temp*temp;
  modelrun := newsumsq;

function fitmodel: real;
var i, iter: integer;
    delta, sumsq, savesum: real;
    done: boolean;
  gotoXY(30, 1); write('FITTING MODEL TO DATA');
  k := 0.5;
  delta := -0.25;
  iter := 0;
  sumsq := 1e30;
  done := false;
    iter := iter + 1;
    gotoxy(23, 3); write('k = ', k:10:6, ' delta = ', delta:10:6);
    savesum := sumsq;
    sumsq := modelrun;
    if abs(delta) < 0.001 then done := true
        if sumsq > savesum then delta := -delta/2.0;
        k := k + delta;
        if k <= 0 then
            delta := - delta;
            k := k + delta;
  until done or (iter >= 50);
  fitmodel := sqrt(sumsq/NPoints);

procedure GraphData;
var i, ref1, ref2, temp: integer;
  outtextxy(0,40,'PARTICIPANT''S HANDLE');
  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
      putpixel(i div 6,maxy div 2 + d^[i] + h^[i],lightred);
      putpixel(i div 6,maxy div 2 - mh^[i], lightgreen);
      putpixel(i div 6,maxy div 2 - h^[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);

procedure PCTmodel;
  sumsq, rmsModel, a, b, r, rms, Xbar, Ybar: real;
  rms := fitmodel;
  writeln('Value of k: ', k:7:4,' RMS prediction error, pixels: ',rms:3:2);
  Correlation(mh, h, NPoints, a, b, rms, r, Xbar, Ybar);
  writeln('Model vs Actual Handle Correlation = ', r:8:3);
  writeln(' RMS error = ', rms:8:3);
  'Fitting model handle. Enter value of k from a previous run? Y/N: ');
  ch := readkey; writeln(ch);
  if (ch = 'Y') or (ch = 'y') then
      write(' Enter value: '); readln(k);
  write('Fitting model to data, please wait...');
  sumsq := modelrun;
  rmsModel := sqrt(sumsq/NPoints);


  if Success then PCTmodel;
    ch := readkey;
  until ch = #27;