This file is indexed.

/usr/share/psychtoolbox-3/PsychDemos/MinExpEntStairDemo.m is in psychtoolbox-3-common 3.0.11.20131230.dfsg1-1build1.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
% The demo and MinExpEntStair use nested functions internally, something
% not supported by Octave, so this is a no-go unless somebody rewrites this
% stuff:
if IsOctave
    fprintf('Sorry, this demo does not yet work on GNU/Octave, neither does the MinExpEntStair function used.\n');
    return;
end

close all; clear all

% stair input
probeset    = -15:0.5:15;       % set of possible probe values
meanset     = -10:0.2:10;       % sampling of pses, doesn't have to be the same as probe set
slopeset    = [.5:.1:5].^2;     % set of slopes, quad scale
lapse       = 0.05;             % lapse/mistake rate
guess       = 0.50;             % guess rate / minimum correct response rate (for detection expt: 1/num_alternative, 0 for discrimination expt)

% general settings
ntrial  = 40;
qpause  = false;    % pause after every iteration? (press any key to continue)
qplot   = false;    % plot information about each trial? (this pauses as well, regardless of whether you specified qpause as true)

% model observer parameters 
qusemodel = true;   % use model observer to get responses? Or, if false, input responses by hand (0/1)
truepse = 0;        % inflection point (50% if guess rate is 0)
truedl  = 4;        % (75%-25% point)/2 if guess rate is 0. In general, take the position of the halfway points between the inflection point and the upper and lower asymptotes, then its the distance between them

% model observer definition. uses a cum normal for psychometric function,
% the formula for which is equivalent to what is used by the staircase
% internally (if it is set up to use a cumulative Gaussian)
if guess==0
    resp    = @(probe)  lapse/2 + (1-lapse)      *normcdf(probe,truepse,truedl/sqrt(2)/erfinv(0.5)) > rand;
else
    resp    = @(probe)   guess  + (1-lapse-guess)*normcdf(probe,truepse,truedl/sqrt(2)/erfinv(0.5)) > rand;
end


% Create staircase instance. If you want to interleave staircases, or
% otherwise run multiple conditions, create one staircase per condition.
% You can store these in a cell-array and of course use different settings
% for each as needed.
stair = MinExpEntStair('v2');
% use lookup table to cache pvalues and avoid unnecessary evaluations of
% psychometric function? Can require lots of memory, especially when
% stepsize of probeset and meanset is not equal. Call before calling
% stair.init.
stair.set_use_lookup_table(true);
% option: use logistic instead of default cumulative normal. best to call
% before stair.init
% stair('set_psychometric_func','logistic');
% init stair
stair.init(probeset,meanset,slopeset,lapse,guess);
% option: use a subset of all data for choosing the next probe, use
% proportion of available data (good idea for robustness - see docs)
stair.toggle_use_resp_subset_prop(10,.9);
% option: instead of choosing a random value for the first probe,
% you can set which value is to be tested first.
stair.set_first_value(3);

for ktrial = 1:ntrial
    % trial
    [p,entexp,ind]  = stair.get_next_probe();      % get next probe to test
    fprintf('%d, new sample point: %f\nexpect ent: %f\n', ...
        ktrial,p,entexp(ind));
    
    if qusemodel % set whether model creates response or you do by manual input
        % get observer response from model observer
        r = resp(p);
        fprintf('response: %d\n',r);
    else
        % make the response yourself, provide either 0 or 1 (actually,
        % anything below and including 0 or anything above 0)
        r = input(sprintf('r(%d): ',ktrial));
        qpause = false;
    end
    stair.process_resp(r);                        % store response in staircase
    % end trial
    
    if ktrial == ntrial || qplot
        
        [m,s,loglik]    = stair.get_fit();
        [ps,rs]         = stair.get_history();
        
        figure(1);
        subplot(1,3,1);
        imagesc(exp(.5*loglik));
        set(gca,'YTick',1:4:length(slopeset));
        set(gca,'YTickLabel',slopeset(1:4:end));
        set(gca,'XTick',1:5:length(meanset));
        set(gca,'XTickLabel',meanset(1:5:end));
        title('estimated likelihood function');
        xlabel('mean (a)')
        ylabel('slope (b)')
        
        subplot(1,3,2);
        hold off;
        if ~isempty(entexp)
            plot(probeset,entexp,'k-o');
            hold on;
            plot(ps(ktrial)*[1,1],[min(entexp),max(entexp)],'r-');
        else
            plot(ps(ktrial)*[1,1],[0,1],'r-');
        end
        title('expected entropy vs probe');
        xlabel('possible probe values')
        xlim([min(probeset) max(probeset)]);
        
        subplot(1,3,3);
        pind = find(rs>0);
        nind = setdiff(1:length(ps),pind);
        plot(1:length(ps),ps,'k-',pind,ps(pind),'bo',nind,ps(nind),'ro');
        ylim([min(probeset) max(probeset)]);
        title('probe-resp history');
    end
    
    % pause if needed
    if (ktrial ~= ntrial) && (qpause || qplot)
        pause;
    end
    
end % loop over trials

%%% show final results
%  [mu,sigma] = stair('get_fit');    % get fitted parameters of cumulative Gaussian
%  DL = sigma*erfinv(.5)*sqrt(2)     % convert sigma (std) to DL (75%-25% point)/2
% get DL from staircase directly, NB: the space of the outputted
% loglikelihood is the mean/slope space as defined atop this script, its
% not a PSE/DL space
[PSEfinal,DLfinal,loglikfinal]  = stair.get_PSE_DL();
finalent                        = sum(-exp(loglikfinal(:)).*loglikfinal(:));
fprintf('final estimates:\nPSE: %f\nDL: %f\nent: %f\n',PSEfinal,DLfinal,finalent);
% for actual offline fitting of your data, you would probably want to use a
% dedicated toolbox such as Prins, N & Kingdom, F. A. A. (2009) Palamedes:
% Matlab routines for analyzing psychophysical data.
% http://www.palamedestoolbox.org.
% Also note that while the staircase runs far more rebust when a small
% lapse rate is assumed, it is common to either fit the psychometric
% function without a lapse rate, or otherwise with the lapse rate as a free
% parameter (possibily varying only over subjects, but not over conditions
% within each subject).

figure(2);
imagesc(exp(.5*loglikfinal));
set(gca,'YTick',1:4:length(slopeset));
set(gca,'YTickLabel',slopeset(1:4:end));
set(gca,'XTick',1:5:length(meanset));
set(gca,'XTickLabel',meanset(1:5:end));
xlabel('$\mu$','interpreter','latex')
switch stair.get_psychometric_func()
    case 'cumGauss'
        title('estimated likelihood function - cumulative Gaussian')
        ylabel('$\sigma$','interpreter','latex')
    case 'logistic'
        title('estimated likelihood function - logistic')
        ylabel('$s$','interpreter','latex')
end