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.

Regards,

Bruce

···

----------------------------------------------------------------
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
                   abbott@cvax.ipfw.indiana.edu
}

uses DOS, CRT, Graph, GrUtils;

const
     MAXDATA = 3600;

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

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

procedure InitHeapVars;
begin
  new(h);
  new(mh);
  new(state);
  new(d);
end;

procedure DisposeHeapVars;
begin
  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: real;
var
  i: 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];
  for i := 1 to NPoints do
  begin
    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;
  end;
  modelrun := newsumsq;
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;
    if abs(delta) < 0.001 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/NPoints);
  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 - 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);
    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;
  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);
  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;
  rmsModel := sqrt(sumsq/NPoints);

  GraphData;
end;

begin
  InitScreen;
  InitHeapVars;
  RestoreCRTMode;
  ReadData(Success);
  if Success then PCTmodel;
  repeat
    ch := readkey;
  until ch = #27;
  DisposeHeapVars;
  CloseGraph;
end.