Skip to content

Commit

Permalink
first release/upload to github
Browse files Browse the repository at this point in the history
  • Loading branch information
Heiko Schütt committed May 21, 2014
1 parent 88b01f9 commit fd65945
Show file tree
Hide file tree
Showing 30 changed files with 2,042 additions and 3 deletions.
6 changes: 3 additions & 3 deletions README.md
View file
@@ -1,9 +1,9 @@
psignifit
=========

Toolbox for Bayesian inference for psychometric functions
Copyright (C) 2014 Heiko Schuett
Toolbox for Bayesian psychometric function estimation

Copyright (C) 2014 Heiko Schütt

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
Expand Down
106 changes: 106 additions & 0 deletions demo_001.m
View file
@@ -0,0 +1,106 @@
%% DEMO_001 basic useage
% The Psifnifit 101

%% save data in right format

% First we need the data in the format (x | nCorrect | total)
% see demo_XXX to find ways how to obtain this
% As an example we use the following dataset from a 2AFC experiment with 90
% trials at each stimulus level. This dataset comes from a simple signal
% detection experiment.

data = [...
0.0010, 45.0000, 90.0000;...
0.0015, 50.0000, 90.0000;...
0.0020, 44.0000, 90.0000;...
0.0025, 44.0000, 90.0000;...
0.0030, 52.0000, 90.0000;...
0.0035, 53.0000, 90.0000;...
0.0040, 62.0000, 90.0000;...
0.0045, 64.0000, 90.0000;...
0.0050, 76.0000, 90.0000;...
0.0060, 79.0000, 90.0000;...
0.0070, 88.0000, 90.0000;...
0.0080, 90.0000, 90.0000;...
0.0100, 90.0000, 90.0000];

% remark: This format differs slightly from the format used in older
% psignifit versions. You can convert your data by using the reformat
% comand. If you are a user of the older psignifits, see demo_XXX for an
% overview over the differences in use.

%% construct an options struct
% To start psignifit you need to pass a struct, which specifies, what kind
% of experiment you did and any other parameters of the fit you might want
% to set:

% You can create a struct by simply calling [name]=struct

options = struct; % initialize as an empty struct

%Now you can set the different options with lines of the form
%[name].[field] as in the following lines:

options.sigmoidName = 'norm'; % choose a cumulative Gauss as the sigmoid
options.expType = '2AFC'; % choose 2-AFC as the paradigm of the experiment
% this sets the guessing rate to .5 and
% fits the rest of the parameters

% There are 3 other types of experiments supported out of the box:
% n alternative forces choice. The guessing rate is known.
% options.expType = "nAFC"
% options.expN = [number of alternatives]
% Yes/No experiments a free guessing and lapse rate is estimated
% options.expType = "YesNo"
% equal asymptote, as Yes/No, but enforces that guessing and lapsing occure
% equally often
% options.expType = "equalAsymptote"

% Out of the box psignifit supports the following sigmoid functions,
% choosen by:
% options.sigmoidName = ...
%
% 'norm' a cummulative gauss distribution
% 'logistic' a logistic function
% 'gumbel' a cummulative gumbel distribution
% 'rgumbel' a reversed gumbel distribution
% 'tdist' a t-distribution with df=1 as a heavytail distribution
%
% for positive stimulus levels which make sence on a log-scale:
% 'logn' a cumulative lognormal distribution
% 'Weibull' a Weibull function

% There are many other options you can set in the options-file. You find
% them in demo_005


%% Now run psignifit
% Now we are ready to run the main function, 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.

res = psignifit(data,options);

%% visualize the results
% For example you can use the result struct res to plot your psychometric
% function with the data:

plotPsych(res);

% You find different visualizations in demo_XXX and an introduction to
% the statistical tests on psychometric functions in demo_XXX. All of these
% are started with the result struct as an argument.



%% remark for insufficient memory issues
% especially for YesNo experiments the result structs can become rather
% large. If you run into Memory issues you can drop the Posterior from the
% result with the following command.

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

% without these fields you will not be able to use the 2D Bayesian plots
% anymore and the equalityTest will fail. All other functions work without
% it.
18 changes: 18 additions & 0 deletions getColormap.m
View file
@@ -0,0 +1,18 @@
function MAP = getColormap()
% returns our standard Tuebingen Colormap
%function getColormap()

%midBlue = [0 105 170]./255;
midBlue = [165 30 55]./255;
%lightBlue = [80 170 200]./255;
%lightBlue = [125 165 75]./255;
lightBlue = [210 150 0]./255;

steps = linspace(0,1,100)';

%MAP = bsxfun(@times,steps,midBlue);
MAP = [];
MAP = [MAP; bsxfun(@times,steps,lightBlue)+bsxfun(@times,1-steps,midBlue)];
MAP = [MAP; bsxfun(@times,steps,[1,1,1])+bsxfun(@times,1-steps,lightBlue)];

end
112 changes: 112 additions & 0 deletions getConfRegion.m
View file
@@ -0,0 +1,112 @@
function [conf_Intervals, confRegion]=getConfRegion(result)
% get confidence intervals and region for parameters
% function [conf_Intervals, confRegion]=getConfRegion(result)
% This function returns the conf_intervals for all parameters and a
% confidence region on the whole parameter space.
%
% Useage
% pass the result obtained from psignifit
% additionally in confP the confidence/ the p-value for the interval is required
% finally you can specify in CImethod how to compute the intervals
% 'project' -> project the confidence region down each axis
% 'stripes' -> find a threshold with (1-alpha) above it
% 'percentiles' -> find alpha/2 and 1-alpha/2 percentiles
% (alpha = 1-confP)


mode = result.options.CImethod;
d = length(result.X1D);

if nargout > 1 || strcmp(mode,'project')
[~,order] = sort(result.Posterior(:),'descend');

Mass = result.Posterior.*result.weight;
Mass = cumsum(Mass(order));

confRegion = true(size(result.Posterior));
confRegion(order(Mass>=result.options.confP)) = false;
end
%% get confIntervals for each paramerter-> marginalize

conf_Intervals=zeros(d,2);

switch mode
case 'project'
for id=1:d
confRegionM = confRegion;
for id2=1:d
if id~=id2
confRegionM = any(confRegionM,id2);
end
end
start = result.X1D{id}(find(confRegionM,1,'first'));
stop = result.X1D{id}(find(confRegionM,1,'last'));
conf_Intervals(id,:) = [start,stop];
end
case 'stripes'
for id=1:d
[margin,x,weight1D] = marginalize(result,id);
[~,order] = sort(margin,'descend');

Mass = margin.*weight1D;
MassSort = cumsum(Mass(order));
% find smallest possible percentage above confP
confP1 = min(MassSort(MassSort>result.options.confP));
confRegionM = true(size(margin));
confRegionM(order(MassSort>confP1))=false;

% Now we have the confidence regions
% put the borders between the nearest contained and the first
% not contained point

% we move in from the outer points to collect the half of the
% leftover confidence from each side

startIndex=find(confRegionM,1,'first');
pleft = confP1-result.options.confP;
if startIndex>1,
start = (x(startIndex)+x(startIndex-1))/2;
start = start + pleft/2/margin(startIndex);
else start = x(startIndex); end
stopIndex=find(confRegionM,1,'last');
if stopIndex < length(x)
stop = (x(stopIndex)+x(stopIndex+1))/2;
stop = stop - pleft/2/margin(stopIndex);
else stop = x(stopIndex); end

conf_Intervals(id,:)=[start,stop];
end
case 'percentiles'
for id=1:d
[margin,x,weight1D]=marginalize(result,id); % marginalize
if length(x)==1 % catch when a dimension was fixed
start=x;
stop=x;
else
Mass = margin.*weight1D; % probability Mass for each point
cumMass = cumsum(Mass); % cumulated p-Mass

% in the intervall are all points between the alpha and the
% 1-alpha percentile.-> where alpha<cumMass<(1-alpha)
confRegionM = cumMass > ((1-result.options.confP)/2) & cumMass < (1-(1-result.options.confP)/2);

% Now we have the confidence regions
% put the borders between the nearest contained and the first
% not contained point
alpha = (1-result.options.confP)/2;
startIndex = find(confRegionM,1,'first');
if startIndex > 1
start = x(startIndex-1) +(alpha-cumMass(startIndex-1))/(cumMass(startIndex) - cumMass(startIndex-1))*(x(startIndex)-x(startIndex-1));
else start = x(startIndex); end

stopIndex=find(confRegionM,1,'last');
if stopIndex < length(x)
stop = x(stopIndex) + (1-alpha-cumMass(stopIndex))/(cumMass(stopIndex+1) - cumMass(stopIndex))*(x(stopIndex+1)-x(stopIndex));
else stop = x(stopIndex); end
end

conf_Intervals(id,:)=[start,stop];
end
otherwise
error('You specified an invalid mode')
end
11 changes: 11 additions & 0 deletions getCor.m
View file
@@ -0,0 +1,11 @@
function Cor = getCor(result)
% calculates the correlation matrix
%function Cor = getCor(result)
% it uses getCov to get the covariance matrix
% result should be a result struct which contains the full posterior
% Cor then is the correlation matrix

Cor = getCov(result);
std = sqrt(diag(Cor));
Cor = bsxfun(@rdivide,Cor,std);
Cor = bsxfun(@rdivide,Cor,reshape(std,1,[]));
48 changes: 48 additions & 0 deletions getCov.m
View file
@@ -0,0 +1,48 @@
function Cov = getCov(result)
% this function calculates the covariance matrix
%function Cov = getCov(result)
% result should be a standard result struct which contains the whole
% posterior
% the function then returns the covariance matrix of the posterior believe
% about the parameters
% If result has a precalculated covariance matrix it is returned instead of
% a newly calculated one

if isfield(result,'Cov')
Cov = result.Cov;
else
Cov = nan(5,5);

for i = 1:5
for j = 1:5
if j < i
Cov(i,j) = Cov(j,i);
else
switch i
case 1, x = reshape(result.X1D{i} ,[],1);
case 2, x = reshape(result.X1D{i} ,1,[]);
case 3, x = reshape(result.X1D{i} ,1,1,[]);
case 4, x = reshape(result.X1D{i} ,1,1,1,[]);
case 5, x = reshape(result.X1D{i} ,1,1,1,1,[]);
end
switch j
case 1, y = reshape(result.X1D{j} ,[],1);
case 2, y = reshape(result.X1D{j} ,1,[]);
case 3, y = reshape(result.X1D{j} ,1,1,[]);
case 4, y = reshape(result.X1D{j} ,1,1,1,[]);
case 5, y = reshape(result.X1D{j} ,1,1,1,1,[]);
end
if length(x) == 1 || length(y)==1
Cov(i,j) = 0;
else
[marginal,~,weight] = marginalize(result,[i,j]);
Mass = marginal .*weight;
Exy = sum(sum(bsxfun(@times,bsxfun(@times,Mass,x),y),i),j);
Ex = sum(sum(bsxfun(@times,Mass,x),i),j);
Ey = sum(sum(bsxfun(@times,Mass,y),j),i);
Cov(i,j) = Exy-Ex*Ey;
end
end
end
end
end
60 changes: 60 additions & 0 deletions getSeed.m
View file
@@ -0,0 +1,60 @@
function Seed = getSeed(data,options)
% calculation of a first guess for the parameters
%function Seed=getSeed(data,options)
% This function finds a seed for initial optimization by logistic
% regression of the data
% Additional parameter is the options.widthalpha which specifies the
% scaling of the width by
% width= psi^(-1)(1-alpha) - psi^(-1)(alpha)
% where psi^(-1) is the inverse of the sigmoid function.

%input parsing
alpha0 = options.widthalpha;
if options.logspace, data(:,1) = log(data(:,1)); end

x = data(:,1);
y = data(:,2)./data(:,3);

% lower and upper asymptote taken simply from min/max of the data
lambda = 1-max(data(:,2)./data(:,3));
gamma = min(data(:,2)./data(:,3));

% rescale y
y = (y - gamma) ./ (1 - lambda - gamma);


% prevent 0 and 1 as bad input for the logit
% this moves the data in from the borders by .25 of a trial from up and
% .25 of a trial from the bottom
factor = .25 ./ data(:,3);
y = factor + (1-2 .* factor) .* y;

% logit transform
y = log(y ./ (1-y));


% fit robust if possible
if exist('robustfit')
fit = robustfit(x, y);
else % Without statistics toolbox do standard linear regression
fit = polyfit(x, y, 1);
fit = [fit(2);fit(1)]; % change to format as it is returned from robustfit
end


% threshold at the zero of the linear fit
% fit(2)*alpha+fit(1)=0
alpha = -fit(1) / fit(2);

% width of the function difference between x'es where the logistic is alpha
% and where it is 1-alpha
beta = (log((1-alpha0) ./ alpha0) - log(alpha0 ./ (1-alpha0))) ./ fit(2);

% varscale is set to almost 0-> we start with the binomial model
varscale = exp(-20);


Seed = [alpha; beta; lambda; gamma; varscale];


end

0 comments on commit fd65945

Please sign in to comment.