Initial public commit

.gitignore

@@ -0,0 +1,8 @@
+/ModelRecoveryTests/*
+/schedule/*
+/logs/*
+/Other people's code/*
+/Plots/*
+Old/
+*Dirs.mat
+RepDirs.mat

README.md

@@ -0,0 +1,22 @@
+# distractor-stats-analysis
+
+This repository contains code for the analysis in the paper [_Explaining the effects of distractor statistics in visual search_](https://doi.org/10.1101/2020.01.03.893057). Most code is written for Matlab.
+
+Contact: Joshua Calder-Travis, [email protected]
+
+_Those sections of the code used in associated paper have been carefully checked. Nevertheless, no guarantee or warranty is provided for any part of the code._
+
+## Main functions
+- `runMultipleLogisicRegressions.m`: Code for running the logistic regressions reported in the paper
+- `fitModels.m`: Code for running fits, or for packing all relevant information so that fits can be run on a computer cluster
+- `mT_runOneJob.sh`: Code used for submitting jobs on a computer cluster
+- `makePlotsForPaper.m`: Code used for making the plots in the paper 
+- `attemptModelRecovery.m`: Code for simulating data with different models and then fitting the simulated datasets
+- `runAllTests.m`: Code for running various code checks
+- `runFullCollationPipeline.m`: Code used for collecting and managing data from the experiment
+
+These functions rely on code in the various subfolders. Therefore, the subfolders need to be on the Matlab path.
+
+## The DSet structure
+Where a dataset is requested (`DSet`) it should be in the format described in the
+README for the [modellingTools repositotry](https://github.com/jCalderTravis/mat-comp-model-tools).

analysisFuns/addNoiseToStim.m

@@ -0,0 +1,86 @@
+function percepts = addNoiseToStim(Model, ParamStruct, Data, stim, sampleShortcut)
+% Add noise to observed stimulus. Different noise can be added for different
+% trials. The code treats all orieations with the same first index, i.e.
+% stim(i, :, :) as coming from the same trial.
+
+% INPUT
+% Model: See 'findDefaultModelSettings' for a description of all model 
+% variants avaliable.
+% ParamStruct: Structure containing a field for each fitted parameter.
+% Data: Strcuture containing field for each feature of the
+% stimulus/behaviour. Fields contain arrays which are num trials long,
+% along the first dimention.
+% stim: [num trials X num locations X nDraws] array describing the orientation 
+% of the Gabor patches
+% sampleShortcut: Instead of sampling from the von Mises. Draw a limited
+% number of samples for each level of observer precision and
+% then just sample from these. Currently only works when SetSizePrec is
+% 'variable'.
+
+
+% Produce array of simulated von Mises noise. Set up an array where rows are
+% trials, and the third dimention indexes different simulations of the trial 
+% (the different nDraws).
+noise = NaN(size(stim));
+
+% Only need to simulate noise for locations in which a stimulus was presented.
+presentedLocs = ~isnan(stim);
+
+% Is the precision of the noise fixed for all trials or does it depend on set
+% size?
+if strcmp(Model.SetSizePrec, 'fixed')
+    assert(isequal(shape(ParamStruct.Kappa_x), [1, 1]))
+    relKappaX = ParamStruct.Kappa_x;
+
+    noise(presentedLocs) ...
+        = qrandvm(0, relKappaX, sum(presentedLocs(:)));
+
+elseif strcmp(Model.SetSizePrec, 'variable')
+    if ~sampleShortcut
+
+        % Find the relevant Kappa_x value. This depends on the set size.
+        assert(isequal(length(unique(Data.SetSizeCond)), ...
+            length(ParamStruct.Kappa_x)) || ...
+            isequal(length(unique(Data.SetSizeCond)), 1))
+
+        relKappaX = ParamStruct.Kappa_x(Data.SetSizeCond);
+        assert(isequal(size(relKappaX), size(Data.SetSizeCond)))
+
+        relKappaXShaped = repmat(relKappaX, 1, size(stim, 2));
+
+        if length(size(stim)) == 3
+            relKappaXShaped = repmat(relKappaXShaped, 1, 1, size(stim, 3));
+        end
+        assert(isequal(size(relKappaXShaped), size(stim)))
+
+        relKappaXShaped(~presentedLocs) = NaN;
+
+        noise ...
+            = qrandvm(0, relKappaXShaped, size(noise));
+
+    elseif sampleShortcut
+        % Loop through the observer noise levels
+        for iNoise = 1 : length(ParamStruct.Kappa_x)
+
+            numSamples = 5000;
+            vonMisesSamples = qrandvm(0, ParamStruct.Kappa_x(iNoise), numSamples);
+
+            relevantLocs = presentedLocs & (Data.SetSizeCond == iNoise);
+
+            noise(relevantLocs) ...
+                = vonMisesSamples( ...
+                randi(numSamples, [sum(relevantLocs(:)), 1, 1]));
+        end
+    end
+end
+
+
+% Add to the stimuli and then map them back in range if they are outside
+% [-pi pi]
+assert(isequal(size(stim), size(noise)))
+percepts = stim + noise;
+percepts = vS_mapBackInRange(percepts, -pi, pi);
+
+% Check we have NaNs in the same place as we started
+assert(isequal(isnan(stim), isnan(percepts)))
+

analysisFuns/assignDefaultParamValues.m

@@ -0,0 +1,67 @@
+function Params = assignDefaultParamValues(ModelName)
+% Assign some default values to all the parameters in the model ready for a
+% simulation. Parameters are produced in the 'unpacked' form (see 
+% modellingTools repository README).
+
+% INPUT
+% ModelName: One of the "model name" structures produced by enumerateModels
+
+
+% All models use kappa_x, but how many values of kappa_x depends on the model
+if strcmp(ModelName.SetSizePrec, 'variable')
+    logKappaX = log([ 2 3 4 5 ]');
+    logKappaX = logKappaX + (0.2 * randn); % Add randomness to parameter selection
+    Params.Kappa_x = exp(logKappaX);
+
+elseif strcmp(ModelName.SetSizePrec, 'fixed')
+    logKappaX = log(3.5) + (0.2 * randn);
+    Params.Kappa_x = exp(logKappaX);
+
+end
+
+
+% The 'min' inference model also uses thresholds. The number of these depends on
+% whether the observer uses different thresholds for different set sizes and
+% block types.
+if strcmp(ModelName.Inference, 'min')
+    if strcmp(ModelName.SetSizeThresh, 'variable')
+        Params.Thresh = [4*pi/16 : -pi/16 : pi/16]';
+        Params.Thresh = Params.Thresh - (pi/32) + ((pi/16)*rand(1));
+
+    elseif strcmp(ModelName.SetSizeThresh, 'fixed')
+        Params.Thresh = 2*pi/16;
+        Params.Thresh = Params.Thresh - (2*pi/32) + ((2*pi/16)*rand(1));
+    end
+
+    if strcmp(ModelName.BlockTypes, 'use')
+        Params.Thresh = ...
+            [Params.Thresh, (Params.Thresh - (pi/32))];
+    end
+
+    assert(all(Params.Thresh(:)>0))
+    assert(all(Params.Thresh(:)<pi))
+end
+
+
+% Are lapses being modelled?
+if strcmp(ModelName.Lapses, 'yes')    
+    Params.LapseRate = 0.15;
+    Params.LapseRate = Params.LapseRate - 0.1 + (0.2 * rand(1));
+end
+
+
+% The Bayesian model has some parameters unique to it
+if strcmp(ModelName.Inference, 'bayes')
+
+    % Does the observer use the true prior?
+    if strcmp(ModelName.Prior, 'biased')
+        Params.ObserverPrior = 0.65;
+        Params.ObserverPrior = Params.ObserverPrior - 0.3 + (0.6 * rand(1));
+    end
+
+    % Does the observer use the true kappa_s?
+    if strcmp(ModelName.BlockTypes, 'ignore')
+        Params.ObserverKappaS = 0.5;
+        Params.ObserverKappaS = Params.ObserverKappaS - 0.25 + (2 * rand(1));
+    end
+end

analysisFuns/circ_vmrnd_fixed.m

@@ -0,0 +1,16 @@
+function samples = circ_vmrnd_fixed(mu, kappa, shape)
+% For kappa_s ~ 0, circ_vmrnd such does not reshape the output. Do the
+% reshaping.
+
+assert(isequal(size(mu), [1, 1]))
+assert(isequal(size(kappa), [1, 1]))
+
+samples = circ_vmrnd(mu, kappa, shape);
+
+if kappa < 1e-6
+    samples = reshape(samples, shape);    
+end
+
+if any(samples(:) > pi) || any(samples(:) < -pi)    
+    error('Von Misses function is returning unexpected values')
+end

analysisFuns/combineSingleParticipantData.m

@@ -0,0 +1,64 @@
+function ConvertedDataFinal = combineSingleParticipantData(Data)
+% The function simSingleCond simulates data for a signle participant and single 
+% condition. This function combines all the simulations for a participant.
+
+% INPUT
+% Data: [num block types x num set sizes] struct array. Each structure is
+% produced by simSingleCond
+
+% Make the orientation data arrays the same size regardless of how many stimuli
+% were in each condition
+for iBlockType = 1 : size(Data, 1)
+    for iSetSize = 1 : size(Data, 2)
+        orientationDataSize = size(Data(iBlockType, iSetSize).Orientation);
+        assert(orientationDataSize(2) <= 8) 
+        % Code currently only deals with this case
+
+        % Pad the orientation data with NaNs until it is 8 columns wide
+        Data(iBlockType, iSetSize).Orientation = ...
+            [Data(iBlockType, iSetSize).Orientation, ...
+            NaN(orientationDataSize(1), 8-orientationDataSize(2))];
+    end
+end
+
+assert(~isfield(Data, 'Title'))
+assert(~isfield(Data, 'NumItems'))
+
+
+% Collect all orientation data, as this has to be treated slightly differently
+% in the combination of data because it is not a vector
+allOrientations = [];
+
+for iStruct = 1 : length(Data(:))
+    allOrientations = [allOrientations; Data(iStruct).Orientation];
+end
+
+% Now we just need to combine all the other data together
+Data = rmfield(Data, 'Orientation');
+dataTable = table();
+
+for iStruct = 1 : length(Data(:))    
+    strcutData = struct2table(Data(iStruct));
+    dataTable = [dataTable; strcutData];
+end
+
+
+% Combine everything 
+ConvertedData = table2struct(dataTable, 'ToScalar', true);
+ConvertedData.Orientation = allOrientations;
+
+for iPtpnt = 1 : length(ConvertedData)
+    vars = fieldnames(ConvertedData);
+
+    for iVar = 1 : length(vars)
+        if isfield(ConvertedData(iPtpnt), vars{iVar})
+
+            ConvertedDataFinal(iPtpnt).Raw.(vars{iVar}) = ...
+                ConvertedData(iPtpnt).(vars{iVar});
+        end
+    end
+end
+
+
+
+

analysisFuns/computeLogBesseliForDuplicatedValues.m

@@ -0,0 +1,12 @@
+function logBesseliVals = computeLogBesseliForDuplicatedValues(vals)
+% Compute log of besseli of order zero. Useful when vals contains many duplicated
+% values. In this case compuation may be sped up by this function. 
+
+uniqueVals = unique(vals);
+uniqueResults = log(besseli(0, uniqueVals));
+
+logBesseliVals = NaN(size(vals));
+
+for iVal = 1 : length(uniqueVals)    
+    logBesseliVals(vals == uniqueVals(iVal)) = uniqueResults(iVal);
+end

analysisFuns/computeStimStats.m

@@ -0,0 +1,113 @@
+function Data = computeStimStats(Data, statType, incTarget)
+% Compute some statistics of the simuli.
+
+% INPUT
+% Data              Struct array or structure with fields, 'Orientation',
+%                   'Target', and 'SetSize'. Additional fields will be added.
+% incTarget         If true then the target is treated as a distractor for the
+%                   computation of the distractor statistics.
+
+% Joshua Calder-Travis 
+% [email protected]
+% GitHub: jCalderTravis
+
+
+for iElement = 1 : length(Data(:))
+    if ~incTarget
+        % First need to find the indcies in Data.Orientation which correspond 
+        % to targets
+        targetIndex = sub2ind([size(Data(iElement).Orientation)], ...
+            find(Data(iElement).Target), ...
+            Data(iElement).TargetLoc(logical(Data(iElement).Target)));
+
+        % Retrieve the stimuli orientations and remove the target orientations 
+        % from these
+        distractorAngle = Data(iElement).Orientation;
+
+        % Code checks
+        targetVals = distractorAngle(targetIndex);
+        if isempty(targetVals)
+            warning('No target trials found.')
+        elseif length(unique(targetVals)) ~= 1
+            error('Bug')
+        end
+        assert(isequal(unique(targetVals), 0))
+
+        distractorAngle(targetIndex) = NaN;
+
+        distractorAngle = mat2cell(distractorAngle, ...
+            ones(size(distractorAngle, 1), 1));
+        assert(isequal(size(distractorAngle), [length(Data(iElement).Target), 1]))
+
+        distractorAngle = cellfun(@(vector) removeNaNs(vector), distractorAngle, ...
+            'UniformOutput', false);
+        assert(isequal(size(distractorAngle), [length(Data(iElement).Target), 1]))
+
+        % Check the resulting vectors are all the expected length
+        for iRow = 1 : length(distractorAngle)
+            expectedLength = Data(iElement).SetSize(iRow) ...
+                - double(Data(iElement).Target(iRow));
+
+            if length(distractorAngle{iRow}) ~= expectedLength; error('bug'); end
+        end
+
+        % If any cells are now empty replace them with NaNs
+        emptyCells = cellfun(@(vector) isempty(vector), distractorAngle);
+        distractorAngle(emptyCells) = {NaN};
+
+    elseif incTarget
+        distractorAngle = Data(iElement).Orientation;
+
+        distractorAngle = mat2cell(distractorAngle, ...
+            ones(size(distractorAngle, 1), 1));
+        assert(isequal(size(distractorAngle), [length(Data(iElement).Target), 1]))
+    end
+
+
+    % Compute some statistics of the stimuli
+    if strcmp(statType, 'circ')
+        [~, Data(iElement).DistractorStd] = cellfun( @(angles) circ_std(angles'), ...
+            distractorAngle);
+        assert(isequal(size(Data(iElement).DistractorStd), ...
+            size(Data(iElement).Target)))
+
+        Data(iElement).DistractorVar = cellfun( @(angles) circ_var(angles'), ...
+            distractorAngle);
+        assert(isequal(size(Data(iElement).DistractorVar), ...
+            size(Data(iElement).Target)))
+
+        Data(iElement).DistractorMean = cellfun( @(angles) abs(circ_mean(angles')), ...
+            distractorAngle);
+        assert(isequal(size(Data(iElement).DistractorMean), ...
+            size(Data(iElement).Target)))
+
+    elseif strcmp(statType, 'real')
+        Data(iElement).DistractorStd = cellfun( @(angles) std(angles'), ...
+            distractorAngle);
+
+        Data(iElement).DistractorVar = cellfun( @(angles) var(angles'), ...
+            distractorAngle);
+
+        Data(iElement).DistractorMean = cellfun( @(angles) abs(mean(angles')), ...
+            distractorAngle);
+    else
+        error('Bug')
+    end
+
+    Data(iElement).MostSimilarDistr = cellfun( @(angles) min(abs(angles)), ...
+        distractorAngle);
+    assert(isequal(size(Data(iElement).MostSimilarDistr), ...
+            size(Data(iElement).Target)))
+end
+
+end
+
+
+function vector = removeNaNs(vector)
+% Remove the NaNs from a vector
+
+vector(isnan(vector)) = [];
+
+end
+
+