My MATLAB psgnifit code is calculating the wrong values from the plot and I am unable to resolve it. Please help.

I am plotting a psychometric curve using psignifit toolbox in MATLAB, and need to calculate the values for DL25 (value of X when Y=0.25) andDL75(value of X when Y=0.75) from the plotted curve. I have defined a function getDL for the same.However, it is returning negative values (X ranges from 0-700 only) for DL25 and incorrect values for DL75 as evident from the plot. I have also defined a function PSE for calculating X when Y=0.5, with a similar function, and is called similarly, but that is functioning smoothly. I am attaching both codes as well as calls here for reference.

Any help in resolving this is greatly appreciated.

Thanks in advance

-Adithi Anil

function DL = getDL(result, targetY, verbose)
    % getDL - Find the X value corresponding to a specified Y value from the fitted curve

    % Syntax:DL = getPSE(result, targetY, verbose)
  
    % Inputs:
    %   result - struct containing the results from psignifit
    %   targetY - the Y value for which to find the corresponding X 
    %   verbose - optional flag for additional output (default is false)
  
    % Outputs:
    %   DL - the X value corresponding to the targetY value

    
    if nargin < 3
        verbose = false; % default verbose flag
    end
    
    % Extract the fit parameters
    params = result.Fit; % Fit parameters
    threshold = 0; % Starting at the lowest point
    width = params(2); % Width
    lambda = params(3); % Upper asymptote
    gamma = params(4); % Lower asymptote (guess rate)

    % Define the logistic sigmoid function
    sigmoidFunc = @(x) lambda + (1 - lambda - gamma) ./ (1 + exp(-(x - threshold) / width));

    % Use fzero to find the X value for which sigmoidFunc(X) = targetY
    DL = fzero(@(x) sigmoidFunc(x) - targetY, threshold); % Starting guess at threshold

    if verbose
        fprintf('DL (X value for Y=%.2f) = %.4f\n', targetY, DL);
    end
end

It is called as:

[DL75] = getDL(result, 0.7500, false); % DL75 = Difference Limen at .75 value of p(Long)

[DL25] = getDL(result, 0.2500, false); % DL25 = Difference Limen at .25 value of p(Long)

For reference, this is the similar function for finding X when Y=0.5, that is working smoothly

function PSE = getPSE(result, targetY, verbose)
    % getPSE - Find the X value corresponding to a specified Y value from the fitted curve
 
    % Syntax: PSE = getPSE(result, targetY, verbose)
   
    % Inputs:
    %   result - struct containing the results from psignifit
    %   targetY - the Y value for which to find the corresponding X (default is 0.5)
    %   verbose - optional flag for additional output (default is false)
  
    % Outputs:
    %   PSE - the X value corresponding to the targetY value

    if nargin < 2
        targetY = 0.5; % default target Y value
    end
    
    if nargin < 3
        verbose = false; % default verbose flag
    end
    
    % Extract the fit parameters
    params = result.Fit; % Fit parameters
    threshold = params(1); % PSE (threshold)
    width = params(2); % Width
    lambda = params(3); % Upper asymptote
    gamma = params(4); % Lower asymptote (guess rate)

    % Define the logistic sigmoid function
    sigmoidFunc = @(x) lambda + (1 - lambda - gamma) ./ (1 + exp(-(x - threshold) / width));

    % Use fzero to find the X value for which sigmoidFunc(X) = targetY
    PSE = fzero(@(x) sigmoidFunc(x) - targetY, threshold); % Starting guess at threshold

    if verbose
        fprintf('PSE (X value for Y=%.2f) = %.4f\n', targetY, PSE);
    end
end

Which is called as:

[PSE] = getPSE(result, 0.5000, false); % PSE= Point of Subjective Equality

Answer 1

[William Rose]

@Adithi Anil, I have not tried running your code. I notice that

gamma=0; % lower symptote, supposedly
lambda=10; % upper symptote, supposedly
threshold=0;  width=1;
%Define the logistic sigmoid function
sigmoidFunc = @(x) lambda + (1 - lambda - gamma) ./ (1 + exp(-(x - threshold) / width));

You say in the code comments that gamma and lambda are the lower and upper asymptotes, respectively. But your sigmoidFunc does not behave the way you want. See plot below.

A logistic function that DOES have gamma and lambda as the lower and upper symptotes, respectively, is defined as follows:

%Define another logistic sigmoid function
sigmoidFunc2 = @(x) gamma + (lambda - gamma) ./ (1 + exp(-(x - threshold) / width));

Plot both sigmoid functions:

x=-5:.2:5;
y1=sigmoidFunc(x);
y2=sigmoidFunc2(x);
figure; plot(x,y1,'-r.',x,y2,'-b.'); 
grid on; xlabel('X'); ylabel('Y')
legend('sigFunc','sigFunc2')

The plot shows that sigmoidFunction2() has the desired asymptotes, and your original sigmoidFunction() does not.

Does that fix your problem?

[Adithi Anil]

Thank you for clarifying that part of the code and I have implemented that change. However, the problem with my DL25 and DL75 values still persist in that they are all returning NaN as the value.

[William Rose]

@Adithi Anil, please attach a script, inlcuding all functions called, that runs, and which returns NaNs for DL25 and DL75.

[Adithi Anil]

Sure. So this is the script for one of the datasets for which these functions are called. Hope this is what you are asking for.

% Data in the format (x | nCorrect | total)

data =[...
     100,   1,    15;...
     200,   0,    15;...
     300,   3,    15;...
    400,   10,    15;...
    500,   7,   15;...
    600,   13,   15;...
    700,   12,   15];

% Creating an options struct which specifies the kind of experiment and  other parameters of the fit

options             = struct;   % initialize as an empty struct

% Setting the different options with lines of the form [name].[field]

options.sigmoidName = 'logistic'; 
options.expType     = 'YesNo';   


% Now run psignifit which fits the function to the data. 
% You obtain a struct, which contains all the information about the
% fitted function and can be passed to the many other functions in this
% toolbox, to further process the results.

result = psignifit(data,options);


% Result is a struct which contains all information obtained from fitting your data. 

% Plot the psychometric function with the data:

plotPsych(result);


[PSE] = getPSE(result, 0.5000, false); % PSE= Point of Subjective Equality

[DL75] = getDL(result, 0.7500, false);  % DL75 = Difference Limen at .75 value 

[DL25] = getDL(result, 0.2500, false); % DL25 = Difference Limen at .25 value 

DL_Self = DL75-DL25;  % DL = Difference Limen (.75-.25) value 

WF_Self = DL_Self/PSE; % WF = Weber Fraction (.75-.25)/PSE 

title('Temporal Bisection Task');
xlabel('Stimulus duration (ms)');
ylabel('Proportion of long response p(Long)');
legend({'S-S Condition'});

% In case of Memory issues, drop the Posterior from the result

resultSmall = rmfield(result,{'Posterior','weight'});

Thank you.

[William Rose]

@Adithi Anil,

Thank you. Running your scrtipt:

% Data in the format (x | nCorrect | total)

data =[...
    100,   1,    15;...
    200,   0,    15;...
    300,   3,    15;...
    400,   10,    15;...
    500,   7,   15;...
    600,   13,   15;...
    700,   12,   15];
% Creating an options struct which specifies the kind of experiment and  other parameters of the fit
options             = struct;   % initialize as an empty struct
% Setting the different options with lines of the form [name].[field]
options.sigmoidName = 'logistic'; 
options.expType     = 'YesNo';   
% Now run psignifit which fits the function to the data. 
% You obtain a struct, which contains all the information about the
% fitted function and can be passed to the many other functions in this
% toolbox, to further process the results.
result = psignifit(data,options);

Unrecognized function or variable 'psignifit'.

% Result is a struct which contains all information obtained from fitting your data. 
% Plot the psychometric function with the data:
plotPsych(result);
[PSE] = getPSE(result, 0.5000, false); % PSE= Point of Subjective Equality
[DL75] = getDL(result, 0.7500, false);  % DL75 = Difference Limen at .75 value 
[DL25] = getDL(result, 0.2500, false); % DL25 = Difference Limen at .25 value 
DL_Self = DL75-DL25;  % DL = Difference Limen (.75-.25) value 
WF_Self = DL_Self/PSE; % WF = Weber Fraction (.75-.25)/PSE 
title('Temporal Bisection Task');
xlabel('Stimulus duration (ms)');
ylabel('Proportion of long response p(Long)');
legend({'S-S Condition'});
% In case of Memory issues, drop the Posterior from the result
resultSmall = rmfield(result,{'Posterior','weight'});

The code does not run because function psignifit() is not supplied. You will also need to supply funciton getDL().

You can format your code as code in this window, and run it, and then you will see what functions, if any, are missing. This will allow others to help you more effectively.

[Adithi Anil]

Sorry for the oversight.

The psignifit () function which comes as a default in the psignifit toolbox is attached below

function result=psignifit(data,options)
% main function for fitting psychometric functions
%function result=psignifit(data,options)
% This function is the user interface for fitting psychometric functions to
% data.
%
% pass your data in the n x 3 matrix of the form:
%       [x-value, number correct, number of trials]
%
% options should be a 1x1 struct in which you set the options for your fit.
% You can find a full overview over the options in demo002
%
% The result of this function is a struct, which contains all information
% the program produced for your fit. You can pass this as whole to all
% further processing function provided with psignifit. Especially to the
% plot functions.
% You can find an explanation for all fields of the result in demo003
%
%
% To get an introduction to basic useage start with demo001




%% input parsing
%% data
if all(data(:,2) <=1 & data(:,2) >= 0) && any(data(:,2) > 0 & data(:,2) < 1) % percent correct data
    data(:,2) = round(data(:,3).*data(:,2)); % we try to convert it to our notation
    % Note: No sanity checks here!
end


%% options
if ~exist('options','var'),                  options=struct;                          end
if ~isfield(options,'sigmoidName'),          options.sigmoidName    = 'norm';         end
if ~isfield(options,'expType'),              options.expType        = 'YesNo';        end
if ~isfield(options,'estimateType'),         options.estimateType   = 'MAP';          end
if ~isfield(options,'confP'),                options.confP          = [0.95,0.9,.68]; end
if ~isfield(options,'instantPlot'),          options.instantPlot    = 0;              end
if ~isfield(options,'maxBorderValue'),       options.maxBorderValue = .00001;         end
if ~isfield(options,'moveBorders'),          options.moveBorders    = 1;              end
if ~isfield(options,'dynamicGrid'),          options.dynamicGrid    = 0;              end
if ~isfield(options,'widthalpha'),           options.widthalpha     = .05;            end
if ~isfield(options,'threshPC'),             options.threshPC       = .5;             end
if ~isfield(options,'CImethod'),             options.CImethod       = 'percentiles';  end
if ~isfield(options,'gridSetType'),          options.gridSetType    = 'cumDist';      end
if ~isfield(options,'fixedPars'),            options.fixedPars      = nan(5,1);       end
if ~isfield(options,'nblocks'),              options.nblocks        = 25;             end
if ~isfield(options,'useGPU'),               options.useGPU         = 0;              end
if ~isfield(options,'poolMaxGap'),           options.poolMaxGap     = inf;            end
if ~isfield(options,'poolMaxLength'),        options.poolMaxLength  = inf;            end
if ~isfield(options,'poolxTol'),             options.poolxTol       = 0;              end
if ~isfield(options,'betaPrior'),            options.betaPrior      = 10;             end
if ~isfield(options,'verbose'),              options.verbose        = 0;              end
if ~isfield(options,'stimulusRange'),        options.stimulusRange  = 0;              end
if ~isfield(options,'fastOptim'),            options.fastOptim      = false;          end



if strcmp(options.expType,'2AFC'),
    options.expType        = 'nAFC';
    options.expN           = 2;
end
if strcmp(options.expType,'3AFC'),
    options.expType        = 'nAFC';
    options.expN           = 3;        
end
if strcmp(options.expType,'4AFC'),
    options.expType        = 'nAFC';
    options.expN           = 4;        
end

if all(~isnan(options.fixedPars(3:4)))
    if options.fixedPars(3)>(1-options.fixedPars(4))
        error('You fixed the lapse rate and the guessing rate to values corresponding to a decreasing sigmoid, which is not supported');
    end
elseif any(~isnan(options.fixedPars(3:4)))
    if ~isnan(options.fixedPars(3))
    end
end

if strcmp(options.expType,'nAFC') && ~isfield(options,'expN');
    error('For nAFC experiments please also pass the number of alternatives (options.expN)'); end

switch options.expType
    case 'YesNo'
        if ~isfield(options,'stepN'),   options.stepN   = [40,40,20,20,20];  end
        if ~isfield(options,'mbStepN'), options.mbStepN = [25,30,10,10,15];  end
    case 'nAFC'
        if ~isfield(options,'stepN'),   options.stepN   = [40,40,20,1,20];   end
        if ~isfield(options,'mbStepN'), options.mbStepN = [30,40,10,1,20];   end
        assert((options.mbStepN(4) == 1) && (options.stepN(4) == 1), 'For nAFC experiments gamma is fixed. The number of gridpoints for it must be 1!')
    case 'equalAsymptote'
        if ~isfield(options,'stepN'),   options.stepN   = [40,40,20,1,20];   end
        if ~isfield(options,'mbStepN'), options.mbStepN = [30,40,10,1,20];   end
        assert((options.mbStepN(4) == 1) && (options.stepN(4) == 1), 'For equal asymptote experiments gamma is fixed equal to lambda. The number of gridpoints for it must be 1!')
    otherwise
        error('You specified an illegal experiment type')
end

assert(max(data(:,1)) > min(data(:,1)) , 'Your data does not have variance on the x-axis! This makes fitting impossible')

% log space sigmoids
% we fit these functions with a log transformed physical axis
% This is because it makes the paramterization easier and also the priors
% fit our expectations better then.
% The flag is needed for the setting of the parameter bounds in setBorders

if any(strcmpi(options.sigmoidName,{'Weibull','logn','weibull'})) % This is NOT run if options.sigmoidName is a handle here
    options.logspace = 1;
    assert(min(data(:,1)) >= 0, 'The sigmoid you specified is not defined for negative data points!');
else
    options.logspace=0;
end


% if range was not given take from data
if numel(options.stimulusRange)<=1
    if options.logspace
        options.stimulusRange = log([min(data(data(:,1)>0,1)),max(data(:,1))]);
    else
        options.stimulusRange = [min(data(:,1)),max(data(:,1))];
    end
    stimRangeSet = false;
else
    stimRangeSet = true;
    if options.logspace
        options.stimulusRange = log(options.stimulusRange);
    end
end

if ~isfield(options,'widthmin')
    if length(unique(data(:,1)))>1 && ~stimRangeSet
        if options.logspace
            options.widthmin  = min(diff(sort(unique(log(data(:,1))))));
        else
            options.widthmin  = min(diff(sort(unique(data(:,1)))));
        end
    else
        options.widthmin = 100*eps(options.stimulusRange(2));
    end
end

% check gpuOptions
if options.useGPU && ~gpuDeviceCount
    warning('You wanted to use your GPU but MATLAB does not recognize any useable GPU. We thus disabled GPU useage')
    options.useGPU=0;
end
if options.useGPU
    gpuDevice(options.useGPU);
end


% add priors
if options.threshPC  ~= .5 && ~isfield(options,'priors')
    warning('psignifit:ThresholdPCchanged','You changed the percent correct corresponding to the threshold\n please check that the prior is still sensible!')
end

if ~isfield(options,'priors')
    options.priors         = getStandardPriors(options);
else 
    if iscell(options.priors)
        priors = getStandardPriors(options);
        % if provided priors are to short
        if length(options.priors)<5 
            options.priors = [options.priors,cell(1,5-length(options.priors))];
        end
        for ipar = 1:5
            if isa(options.priors{ipar},'function_handle')
                % use the provided prior
            else
                options.priors{ipar} = priors{ipar};
            end
        end
    else
        error('if you provide your own priors it should be a cell array of function handles')
    end
    %check priors 
    checkPriors(data,options);
end

% for dynamic grid setting
if options.dynamicGrid && ~isfield(options,'GridSetEval'),   options.GridSetEval   = 10000; end
if options.dynamicGrid && ~isfield(options,'UniformWeight'), options.UniformWeight = 1;    end



%% initialize

% Warning if many blocks were measured -> adaptive sampling?
if length(unique(data(:,1))) >= 25 && numel(options.stimulusRange)==1
    warning('psignifit:probablyAdaptive', 'The data you supplied contained >= 25 stimulus levels.\n Did you sample adaptively?\n If so please specify a range which contains the whole psychometric function in options.stimulusRange.\n This will allow psignifit to choose an appropriate prior.\n For now we use the standard heuristic, assuming that the psychometric function is covered by the stimulus levels\n, which is frequently invalid for adaptive procedures!')
end

% pool data if necessary
% -> more than options.nblocks blocks or only 1 trial per block
if max(data(:, 3)) == 1 || size(data, 1) > options.nblocks
    warning('psignifit:pooling','We pooled your data, to avoid problems with n=1 blocks or to save time fitting because you have a lot of blocks\n You can force acceptence of your blocks by increasing options.nblocks');
    data = poolData(data,options);
    %options.nblocks = size(data, 1);
else
    %options.nblocks = size(data, 1);
end

% Warning if few trials per block -> adaptive sampling
if all(data(:,3)<=5) && numel(options.stimulusRange)==1
    warning('psignifit:probablyAdaptive', 'All provided data blocks contain <= 5 trials \n Did you sample adaptively?\n If so please specify a range which contains the whole psychometric function in options.stimulusRange.\n This will allow psignifit to choose an appropriate prior.\n For now we use the standard heuristic, assuming that the psychometric function is covered by the stimulus levels\n, which is frequently invalid for adaptive procedures!')
end

% create function handle to the sigmoid
if ~isfield(options,'sigmoidHandle')
    options.sigmoidHandle = getSigmoidHandle(options);
    if isa(options.sigmoidName, 'function_handle')
        options.sigmoidName = 'Custom Handle Provided';
    end
else 
    if strcmp(options.sigmoidName,'norm') % i.e. sigmoidName not set by user
        options.sigmoidName = 'Custom Handle Provided';
    end
end
% borders of integration
if isfield(options, 'borders')
    if any(options.borders(:,1)==options.borders(:,2))
        error('You specified equal upper and lower borders, please use options.fixedPars instead! This is caught here as only chaning the borders will mess with the MAP optimization.')
    elseif any(options.borders(:,2)<options.borders(:,1))
        error('The borders for one or more parameters are in the wrong order!')
    end
    borders = setBorders(options);
    options.borders(isnan(options.borders)) = borders(isnan(options.borders));
    options.borders(~isnan(options.fixedPars),1) = options.fixedPars(~isnan(options.fixedPars)); %fix parameter values
    options.borders(~isnan(options.fixedPars),2) = options.fixedPars(~isnan(options.fixedPars)); %fix parameter values
else
    options.borders = setBorders(options);
    options.borders(~isnan(options.fixedPars),1) = options.fixedPars(~isnan(options.fixedPars)); %fix parameter values
    options.borders(~isnan(options.fixedPars),2) = options.fixedPars(~isnan(options.fixedPars)); %fix parameter values
end

% normalize priors to first choice of borders

options.priors = normalizePriors(options);

if options.moveBorders
    options.borders = moveBorders(data,options);
end

% Warning for high confidence values
if any(options.confP > 0.95)
    warning('psingifit:confPLarge','You requested a confidence with higher confidence than 95%%. This was not thoroughly tested.');
end


%% core

result = psignifitCore(data, options);


%% after processing

%check that the marginals go to nearly 0 at the borders of the grid
if options.verbose > -5
    if result.marginals{1}(1).* result.marginalsW{1}(1) > .001
        warning('psignifit:borderWarning',...
           ['The marginal for the threshold is not near 0 at the lower border\n',...
            'This indicates that smaller Thresholds would be possible'])
    end
    if result.marginals{1}(end).* result.marginalsW{1}(end) > .001
        warning('psignifit:borderWarning',...
           ['The marginal for the threshold is not near 0 at the upper border\n',...
            'This indicates that your data is not sufficient to exclude much higher thresholds.\n',...
            'Refer to the paper or the manual for more info on this topic'])
    end
    if result.marginals{2}(1).* result.marginalsW{2}(1) > .001
        warning('psignifit:borderWarning',...
           ['The marginal for the width is not near 0 at the lower border\n',...
            'This indicates that your data is not sufficient to exclude much lower widths.\n',...
            'Refer to the paper or the manual for more info on this topic'])
    end
    if result.marginals{2}(end).* result.marginalsW{2}(end) > .001
        warning('psignifit:borderWarning',...
           ['The marginal for the width is not near 0 at the lower border\n',...
            'This indicates that your data is not sufficient to exclude much higher widths.\n',...
            'Refer to the paper or the manual for more info on this topic'])
    end
end

result.psiHandle = @(x) result.Fit(4)+ (1-result.Fit(3)-result.Fit(4))*result.options.sigmoidHandle(x,result.Fit(1),result.Fit(2));

[result.devianceResiduals,result.deviance] = getDeviance(result);

result.Cov = getCov(result);
result.Cor = getCor(result);

result.timestamp = datestr(now);

if options.instantPlot
    plotPsych(result);
    plotBayes(result);
end

The function getPSE() is defined as:

function PSE = getPSE(result, targetY, verbose)
    % getPSE - Find the X value corresponding to a specified Y value from the fitted curve
    %
    % Syntax: PSE = getPSE(result, targetY, verbose)
    %
    % Inputs:
    %   result - struct containing the results from psignifit
    %   targetY - the Y value for which to find the corresponding X (default is 0.5)
    %   verbose - optional flag for additional output (default is false)
    %
    % Outputs:
    %   PSE - the X value corresponding to the targetY value

    if nargin < 2
        targetY = 0.5; % default target Y value
    end
    
    if nargin < 3
        verbose = false; % default verbose flag
    end
    
    % Extract the fit parameters
    params = result.Fit; % Fit parameters
    threshold = params(1); % PSE (threshold)
    width = params(2); % Width
    lambda = params(3); % Upper asymptote
    gamma = params(4); % Lower asymptote (guess rate)

    % Define the logistic sigmoid function
     sigmoidFunc = @(x) gamma + (lambda - gamma) ./ (1 + exp(-(x - threshold) / width));


    % Use fzero to find the X value for which sigmoidFunc(X) = targetY
    PSE = fzero(@(x) sigmoidFunc(x) - targetY, threshold); % Starting guess at threshold

    if verbose
        fprintf('PSE (X value for Y=%.2f) = %.4f\n', targetY, PSE);
    end
end

And the function getDL(), the one which I am encountering problems with is as follows:

function DL = getDL(result, targetY, verbose)
    % getDL - Find the X value corresponding to a specified Y value from the fitted curve

    % Syntax:DL = getPSE(result, targetY, verbose)
  
    % Inputs:
    %   result - struct containing the results from psignifit
    %   targetY - the Y value for which to find the corresponding X 
    %   verbose - optional flag for additional output (default is false)
  
    % Outputs:
    %   DL - the X value corresponding to the targetY value

    
    if nargin < 3
        verbose = false; % default verbose flag
    end
    
    % Extract the fit parameters
    params = result.Fit; % Fit parameters
    threshold = 0; % Starting at the lowest point
    width = params(2); % Width
    lambda = params(3); % Upper asymptote
    gamma = params(4); % Lower asymptote (guess rate)

    %Define another logistic sigmoid function
    sigmoidFunc = @(x) gamma + (lambda - gamma) ./ (1 + exp(-(x - threshold) / width));

    % Use fzero to find the X value for which sigmoidFunc(X) = targetY
    DL = fzero(@(x) sigmoidFunc(x) - targetY, threshold); % Starting guess at threshold

    if verbose
        fprintf('DL (X value for Y=%.2f) = %.4f\n', targetY, DL);
    end
end

The script for one of the datasets is reattached to avoid any confusion.

% Data in the format (x | nCorrect | total)

data =[...
     100,   1,    15;...
     200,   0,    15;...
     300,   3,    15;...
    400,   10,    15;...
    500,   7,   15;...
    600,   13,   15;...
    700,   12,   15];

% Creating an options struct which specifies the kind of experiment and  other parameters of the fit

options             = struct;   % initialize as an empty struct

% Setting the different options with lines of the form [name].[field]

options.sigmoidName = 'logistic'; 
options.expType     = 'YesNo';   


% Now run psignifit which fits the function to the data. 
% You obtain a struct, which contains all the information about the
% fitted function and can be passed to the many other functions in this
% toolbox, to further process the results.

result = psignifit(data,options);


% Result is a struct which contains all information obtained from fitting your data. 

% Plot the psychometric function with the data:

plotPsych(result);


[PSE] = getPSE(result, 0.5000, false); % PSE= Point of Subjective Equality

[DL75] = getDL(result, 0.7500, false);  % DL75 = Difference Limen at .75 value 

[DL25] = getDL(result, 0.2500, false); % DL25 = Difference Limen at .25 value 

DL_Self = DL75-DL25;  % DL = Difference Limen (.75-.25) value 

WF_Self = DL_Self/PSE; % WF = Weber Fraction (.75-.25)/PSE 

title('Temporal Bisection Task');
xlabel('Stimulus duration (ms)');
ylabel('Proportion of long response p(Long)');
legend({'S-S Condition'});

% In case of Memory issues, drop the Posterior from the result

resultSmall = rmfield(result,{'Posterior','weight'});

Hope this conveys all the required information.

Thank you

[Adithi Anil]

Additionally, the plotPsych() is also a predefined function that is defined as:

function [hline,hdata] = plotPsych(result,plotOptions)
% plot your data with the fitted function
%function plotPsych(result,plotOptions)
% This function produces a plot of the fitted psychometric function with
% the data.
% In the struct plotOptions you may set a lot of options, which are listed 
% here with their default values:
% 
%plotOptions.dataColor      = [0,105/255,170/255];  % color of the data
%plotOptions.plotData       = 1;                    % plot the data?
%plotOptions.lineColor      = [0,0,0];              % color of the PF
%plotOptions.lineWidth      = 2;                    % lineWidth of the PF
%plotOptions.xLabel         = 'Stimulus Level';     % xLabel
%plotOptions.yLabel         = 'Percent Correct';    % yLabel
%plotOptions.labelSize      = 15;                   % font size labels
%plotOptions.fontSize       = 10;                   % font size numbers
%plotOptions.fontName       = 'Helvetica';          % font
%plotOptions.tufteAxis      = false;                % use special axis
%plotOptions.plotAsymptote  = true;                 % plot Asympotes 
%plotOptions.plotThresh     = true;                 % plot Threshold Mark
%plotOptions.aspectRatio    = false;                % set aspect ratio
%plotOptions.extrapolLength = .2;                   % extrapolation percentage
%plotOptions.CIthresh       = false;                % draw CI on threhold
%plotOptions.dataSize       = 10000./sum(result.data(:,3)) % size of the data-dots

if ~exist('plotOptions','var'),         plotOptions           = struct;               end
assert(isstruct(plotOptions),'If you pass an option file it must be a struct!');

if isfield(plotOptions,'h')
    h = plotOptions.h;
else 
    h = gca;
end
assert(ishandle(h),'Invalid axes handle provided to plot in.')
axes(h);

if ~isfield(plotOptions,'dataColor'),      plotOptions.dataColor      = [0,105/255,170/255]; end
if ~isfield(plotOptions,'plotData'),       plotOptions.plotData       = 1;                   end
if ~isfield(plotOptions,'lineColor'),      plotOptions.lineColor      = [0,0,0];             end
if ~isfield(plotOptions,'lineWidth'),      plotOptions.lineWidth      = 2;                   end
if ~isfield(plotOptions,'xLabel'),         plotOptions.xLabel         = 'Stimulus Level';    end
if ~isfield(plotOptions,'yLabel'),         plotOptions.yLabel         = 'Proportion Correct';end
if ~isfield(plotOptions,'labelSize'),      plotOptions.labelSize      = 15;                  end
if ~isfield(plotOptions,'fontSize'),       plotOptions.fontSize       = 10;                  end
if ~isfield(plotOptions,'fontName'),       plotOptions.fontName       = 'Helvetica';         end
if ~isfield(plotOptions,'tufteAxis'),      plotOptions.tufteAxis      = false;               end
if ~isfield(plotOptions,'plotAsymptote'),  plotOptions.plotAsymptote  = true;                end
if ~isfield(plotOptions,'plotThresh'),     plotOptions.plotThresh     = true;                end
if ~isfield(plotOptions,'aspectRatio'),    plotOptions.aspectRatio    = false;               end
if ~isfield(plotOptions,'extrapolLength'), plotOptions.extrapolLength = .2;                  end
if ~isfield(plotOptions,'CIthresh'),       plotOptions.CIthresh       = false;               end

if isnan(result.Fit(4)),                   result.Fit(4)=result.Fit(3);    end

if isempty(result.data)
    return
end



if ~isfield(plotOptions,'dataSize')
    plotOptions.dataSize  = 10000./sum(result.data(:,3));
end


switch result.options.expType
    case {'nAFC'}
        ymin = 1./result.options.expN;
        ymin = min(ymin,min(result.data(:,2)./result.data(:,3)));
    otherwise
        ymin = 0;
end

% backcompatibility
if isfield(plotOptions,'plotPar') && ~plotOptions.plotPar
    plotOptions.plotAsymptote  = false;
    plotOptions.plotThresh     = false;
end

%% plot data
holdState = ishold(h);
if ~holdState
    cla(h);
end
hold on
if exist ('OCTAVE_VERSION', 'builtin')
    hdata = zeros(size(result.data,1),1);
else
if verLessThan('matlab', '8.1')
    hdata = zeros(size(result.data,1),1);
else
    hdata = gobjects(size(result.data,1),1);
end
end
if plotOptions.plotData
    for i=1:size(result.data,1)
        hdata(i) = plot(result.data(i,1),result.data(i,2)./result.data(i,3),'.','MarkerSize',sqrt(plotOptions.dataSize*result.data(i,3)),'Color',plotOptions.dataColor);
    end
end

%% plot fitted function
%for dashed ends:
if result.options.logspace
    xOK = result.data(result.data(:,1)>0,1);
    xlength   = log(max(xOK))-log(min(xOK));
    x         = exp(linspace(log(min(xOK)),log(max(xOK)),1000));
    xLow      = exp(linspace(log(min(xOK))-plotOptions.extrapolLength*xlength,log(min(xOK)),100));
    xHigh     = exp(linspace(log(max(xOK)),log(max(xOK))+plotOptions.extrapolLength*xlength,100));
else
    xlength   = max(result.data(:,1))-min(result.data(:,1));
    x         = linspace(min(result.data(:,1)),max(result.data(:,1)),1000);
    xLow      = linspace(min(result.data(:,1))-plotOptions.extrapolLength*xlength,min(result.data(:,1)),100);
    xHigh     = linspace(max(result.data(:,1)),max(result.data(:,1))+plotOptions.extrapolLength*xlength,100);
end
fitValuesLow    = (1-result.Fit(3)-result.Fit(4))*arrayfun(@(x) result.options.sigmoidHandle(x,result.Fit(1),result.Fit(2)),xLow)+result.Fit(4);
fitValuesHigh   = (1-result.Fit(3)-result.Fit(4))*arrayfun(@(x) result.options.sigmoidHandle(x,result.Fit(1),result.Fit(2)),xHigh)+result.Fit(4);

fitValues = (1-result.Fit(3)-result.Fit(4))*arrayfun(@(x) result.options.sigmoidHandle(x,result.Fit(1),result.Fit(2)),x)+result.Fit(4);
hline = plot(x,     fitValues,          'Color', plotOptions.lineColor,'LineWidth',plotOptions.lineWidth);
plot(xLow,  fitValuesLow,'--',  'Color', plotOptions.lineColor,'LineWidth',plotOptions.lineWidth)
plot(xHigh, fitValuesHigh,'--', 'Color', plotOptions.lineColor,'LineWidth',plotOptions.lineWidth)

if result.options.logspace
    set(gca,'XScale','log')
end

%% plot parameter Illustrations
% threshold
if plotOptions.plotThresh
    if result.options.logspace
        plot([exp(result.Fit(1)),exp(result.Fit(1))],[ymin,result.Fit(4)+(1-result.Fit(3)-result.Fit(4)).*result.options.threshPC],'-','Color',plotOptions.lineColor);
    else
        plot([result.Fit(1),result.Fit(1)],[ymin,result.Fit(4)+(1-result.Fit(3)-result.Fit(4)).*result.options.threshPC],'-','Color',plotOptions.lineColor);
    end
end
if plotOptions.plotAsymptote
    % asymptotes
    plot([min(xLow),max(xHigh)],[1-result.Fit(3),1-result.Fit(3)],':','Color',plotOptions.lineColor);
    plot([min(xLow),max(xHigh)],[result.Fit(4),result.Fit(4)],':','Color',plotOptions.lineColor);
end

if plotOptions.CIthresh
    if result.options.logspace
        result.conf_Intervals(1,:,1) = exp(result.conf_Intervals(1,:,1));
    end
    plot(result.conf_Intervals(1,:,1),repmat(result.Fit(4)+result.options.threshPC*(1-result.Fit(3)-result.Fit(4)),1,2),'Color',plotOptions.lineColor)
    plot(repmat(result.conf_Intervals(1,1,1),1,2),repmat(result.Fit(4)+result.options.threshPC*(1-result.Fit(3)-result.Fit(4)),1,2)+[-.01,+.01],'Color',plotOptions.lineColor)
    plot(repmat(result.conf_Intervals(1,2,1),1,2),repmat(result.Fit(4)+result.options.threshPC*(1-result.Fit(3)-result.Fit(4)),1,2)+[-.01,+.01],'Color',plotOptions.lineColor)
end

%% axis settings
axis tight
set(gca,'FontSize',plotOptions.fontSize)
ylabel(plotOptions.yLabel,'FontName',plotOptions.fontName,'FontSize', plotOptions.labelSize);
xlabel(plotOptions.xLabel,'FontName',plotOptions.fontName,'FontSize', plotOptions.labelSize);
if plotOptions.aspectRatio
    set(gca,'PlotBoxAspectRatio',[(1+sqrt(5))/2,1,1])
end
if plotOptions.tufteAxis
    if result.options.logspace
        tufteaxis([min(result.data(:,1)),exp(result.Fit(1)),max(result.data(:,1))],[ymin,1]);
    else
        tufteaxis([min(result.data(:,1)),result.Fit(1),max(result.data(:,1))],[ymin,1]);
    end
    set(get(gca,'xLabel'),'visible','on')
    set(get(gca,'yLabel'),'visible','on')
else 
    ylim([ymin,1]);
    set(gca,'TickDir','out')
    box off
end

%% toggle back hold state
if ~holdState
    hold off
end

[Adithi Anil]

And since I am noticing that even these functions are calling multiple other functions from the same toolbox, I am also attaching the link for psignifit toolbox: https://github.com/wichmann-lab/psignifit/archive/master.zip

[William Rose]

@Adithi Anil,

This is turning into a long set of functions that call other (long) functions. I will not install someone else's toolbox of functions. I can't run psignifit() because it calls psignifitCore(). And so on.

Is this a good summary:

You use psignifit() to fit a signmoid to response data. psignifit() returns a struture variable, response, which has the fitted parameters and a bunch of other info about the fit. Functions getDL() and getPSE() return the x value of the sigmoid that corresponds to a given y value. The only difference between getDL and getPSE is that getPSE defaults to a y value of 0.5, if no y-value is supplied as an argument. You say getPSE() works fine. You say getDL() returns incorect results when y=0.25 and when y=0.75.

Why write getPSE()? Why not just call getDL(), with y=0.5?

You have two definitions of PSE. Your function getPSE() includes the line

threshold = params(1); % PSE (threshold)

Exanmnation of the logistic function shows that "threshold", labelled PSE in the comment in the code, is the x-value at which y is halfway between the lower and upper asymptotes. However, you say in your explanation above that PSE is the x-value where the y-value=0.5. These are two different definitions. They are only the same if the upper and lower ansymptotes are symmetric about 0.5, which is not guaranteed. Which definition of PSE is right: "x value when y=0.5" or "x value when y is halfway between lower and upper asymptote"?

In the example data you provide,

data =[100, 1,15; 200, 0,15;...
       300, 3,15; 400,10,15;...
       500, 7,15; 600,13,15; 700,12,15];
frac=data(:,2)./data(:,3);
plot(data(:,1),frac,'-r.')
grid on; xlabel('x'); ylabel('Frac. Correct'); 

Have you checked the fitted parameters to make sure the upper fitted asymptote exceeds 0.75? If the upper fitted asymptote is less than 0.75, then of course it will be an error when you try to obtain the x vaue corresponding to y=0.75.

You are using 7 data points to fit a curve (a logistic function) with 4 parameters. (The parameters are upper and lower asymptotes, theshold (=x-value of PSE), and width. (the x-value corresponding to perceived subjective equality) I am confident that the uncertainties in the 4 fitted parameters of the logistic function will be very large, for the example data above. You need more data points in order to have reaosnable confidence in the fitted values. Talk to a statistician about how to obtain confidence intervals for fited parameters. Or read the Matlab help for model fitting functions. .

Finally: Use psignifit() to compute result. Save result as a .mat file. Upload the .mat file to this chat using the paper clip icon. Then others will be able to access result and use it as input to your functions, getDL() and getPSE(). PLease provide the data set used by psignifit(), so we can confirm that the fitted values are reasonable.

[William Rose]

@Adithi Anil, I just noticed another difference between getDL() and getPSE().

fragment of getDL()

% Extract the fit parameters
params = result.Fit; % Fit parameters
threshold = 0; % Starting at the lowest point
width = params(2); % Width
lambda = params(3); % Upper asymptote
gamma = params(4); % Lower asymptote (guess rate)

fragment of getPSE()

% Extract the fit parameters
params = result.Fit; % Fit parameters
threshold = params(1); % PSE (threshold)
width = params(2); % Width
lambda = params(3); % Upper asymptote
gamma = params(4); % Lower asymptote (guess rate)

The line "threshold=0" in getDL() is wrong. Change it to the line in getPSE(). Maybe that will solve your problems with getDL().

You could still get negative x-values back for DL25, depending on the data set and the fit, because there's nothing obvious that prevents the possibility of a fitted logistic funciton that has the value y=0.25 when x is <0.

You can still expect a weird value, or NaN, for DL75, if the upper asymptote is <0.75.