// Display mode
mode(0);

// Display warning for floating point exception
ieee(1);

// Exercise 2.1
// 
// This Scilab script accompanies the chapter "Distinguishing between models
// of perceptual decision making" by J. Ditterich in the book "An
// Introduction to Model-Based Cognitive Neuroscience". It demonstrates that
// an integration-to-threshold mechanism can make accurate decisions based
// on a relatively small number of required samples. Please note that this
// script has been optimized for readability, not execution speed.
// 
// J. Ditterich, 12/2012

clear;
xdel(winsid());

// Parameters

// Adjust the decision threshold (between 5 and 20) until an accuracy of
// approximately 90% is obtained!
thr         = 15; // decision threshold (A)

dist_mean   = 1; // mean of normal distribution from which random samples are drawn
dist_sd     = 3; // standard deviation of normal distribution from which random samples are drawn
trials      = 20000; // number of trials to be simulated

// Let's go...
dt_vec = []; // vector of decision times
num_correct = 0; // number of correct decisions
disp("Starting simulation...");

for i = 1:trials // simulate each trial
    cur_sum = 0; // initialize current sum
    num_samples = 0; // initialize number of drawn samples

    while 1 // loop for sequential sampling
        cur_samp = grand(1,1,"nor",dist_mean,dist_sd); // draw a new random number
        num_samples = num_samples + 1; // update number of random samples (reflecting decision time)
        cur_sum = mtlb_a(cur_sum,cur_samp); // update the sum

        if mtlb_logic(cur_sum,">",thr) then // Did we exceed the upper threshold?
            dt_vec = [dt_vec,num_samples]; // store the decision time
            num_correct = num_correct + 1; // keep track of accuracy
            break; // terminate the loop for this trial
        end;

        if mtlb_logic(cur_sum,"<",-thr) then // Did we cross the lower threshold?
            dt_vec = [dt_vec,num_samples]; // store the decision time
            break; // terminate the loop for this trial
        end;
    end;

    if pmodulo(i,1000)==0 then
        disp(string(i)+" trials completed.");
    end;
end;

disp(" ");
disp("Accuracy: "+msprintf("%.1f",(num_correct/trials)*100)+"%");
disp("Average number of evaluated samples: "+msprintf("%.1f",mean(dt_vec,"m")));
