commiting code

PBups_psychometrics/bootstrap_conditioned_psychometrics.m

@@ -0,0 +1,30 @@
+% wrapper for computing which parameters are significantly modulated by
+% trial history:
+
+nboots = 100;
+
+path_pkgdata = '~/ondrive/analysisDG/PBups_trialhistory/data/packaged_reaction_time_data/';
+files = dir([path_pkgdata, 'sess_rawdata_*.mat']);
+filenames = arrayfun(@(x) x.name,files,'UniformOutput',false);
+filenames = filenames(~contains(filenames, 'wholeStim'));
+
+
+for i = 1:length(filenames)
+    load([path_pkgdata, filenames{i}]);
+
+    data = [];
+    data.pokedR = avgdata.pokedR;
+    data.gamma = avgdata.gamma;
+    data.hits = avgdata.hits;
+
+    tic()
+    boot_rc = bootstrap_psych_fit(data, 'nboots', nboots, 'h_ind', 1, 'pR_ind', 1, 'trial_back', 1);
+    boot_lc = bootstrap_psych_fit(data, 'nboots', nboots, 'h_ind', 1, 'pR_ind', 0, 'trial_back', 1);
+    toc()
+
+    for pn = 1:4
+        [p,h] = ranksum(boot_rc(:,pn), boot_lc(:,pn))
+    end
+
+end
+

PBups_psychometrics/bootstrap_psych_fit.m

@@ -0,0 +1,76 @@
+function[boot_prm] = bootstrap_psych_fit(data_master, varargin)
+
+p = inputParser;
+addParameter(p,'h_ind', []);
+addParameter(p,'pR_ind', []);
+addParameter(p,'trial_back', +1);
+% +1 means 1 trial after the conditioned trial
+addParameter(p, 'xreg','gamma'); % gamma/ Delta
+addParameter(p, 'nboots', 1000);
+addParameter(p, 'ax', []);
+addParameter(p, 'plot_bootstraps', false)
+addParameter(p, 'seed', 1)
+parse(p,varargin{:});
+
+xreg = p.Results.xreg;
+rng(p.Results.seed);
+
+if ~isempty(p.Results.h_ind) & ~isempty(p.Results.pR_ind)
+    ids = find(data_master.hits(1:end) == p.Results.h_ind & data_master.pokedR(1:end) == p.Results.pR_ind) + p.Results.trial_back;
+    ids(ids<1) = [];
+    ids(ids>length(data_master.(xreg))) = [];
+elseif ~isempty(p.Results.h_ind)
+    ids = find(data_master.hits(1:end) == p.Results.h_ind) + p.Results.trial_back;
+    ids(ids<1) = [];
+    ids(ids>length(data_master.(xreg))) = [];
+elseif isempty(p.Results.h_ind) & isempty(p.Results.pR_ind)
+    ids = 1:length(data_master.(xreg));
+end
+
+pokedR = data_master.pokedR;
+x = data_master.(xreg);
+ax = p.Results.ax;
+
+if p.Results.plot_bootstraps
+    plotfit = true;
+    plotdata = true;
+    ploterrorbar = true;
+    if isempty(ax)
+        figure();
+        ax = gca();
+        hold on;
+    end
+else
+    plotfit = false;
+    plotdata = false;
+    ploterrorbar = false;
+    ax = [];
+
+end
+
+% data_master
+boot_prm = zeros(p.Results.nboots, 4);
+parfor b = 1:p.Results.nboots  % resample data
+
+    id_sample = datasample(ids, length(ids));
+    data = [];
+    data.pokedR = pokedR(id_sample);
+    data.(xreg) = x(id_sample); 
+    prm = pbups_psych_gamma(data, ...
+        'xreg', p.Results.xreg, ...
+        'plotfit', plotfit, ...
+        'plotdata', plotdata,...
+        'ploterrorbar', ploterrorbar,...
+        'fitLineColor', [0.8, 0.8, 0.8],...
+        'axHandle', ax);
+    drawnow;
+
+%     sens = prm.beta(3);
+%     bias = prm.beta(4);
+%     left_lapse = prm.beta(1);
+%     right_lapse = 1 - prm.beta(1) - prm.beta(2);
+
+    boot_prm(b,:) = prm.beta;
+end
+
+

PBups_psychometrics/extract_cond_ids.m

@@ -0,0 +1,37 @@
+function[ids] = extract_cond_ids(d, cond)
+    % this function extracts the ids of the trial for which
+    % condition cond is satisfied. session boundaries are respected
+    % d is a structure with following fields:
+    % d.sessid: vector of session ids corr to each trial 
+    % d.hits: 1 hit 0 err 
+    % d.pokedR: 1 poked right 0 poked L
+
+
+    last_trial_ind = find(diff([d.sessid, 1]));
+
+    switch cond
+        case 'last_trial'
+            ids = last_trial_ind;
+        case 'not_last_trial'
+            ids = setdiff(1:length(d.sessid), last_trial_ind);
+        case 'corr'
+            ids = find(d.hits == 1);
+        case 'err'
+            ids = find(d.hits == 0);
+        case 'post_corr'
+            ids = setdiff(find(d.hits == 1), last_trial_ind)+1;
+        case 'post_err'
+            ids = setdiff(find(d.hits == 0), last_trial_ind)+1;
+        case 'post_right_corr'
+            ids = setdiff(find(d.hits == 1 & d.pokedR == 1), last_trial_ind)+1;
+        case 'post_left_corr'
+            ids = setdiff(find(d.hits == 1 & d.pokedR == 0), last_trial_ind)+1;
+        case 'post_right_err'
+            ids = setdiff(find(d.hits == 0 & d.pokedR == 1), last_trial_ind)+1;
+        case 'post_left_err'
+            ids = setdiff(find(d.hits == 0 & d.pokedR == 0), last_trial_ind)+1;
+
+
+    end
+
+end

PBups_psychometrics/fit_logistic2.m

@@ -0,0 +1,60 @@
+function [P varargout] = fit_logistic2(pokedR, xregressor)
+
+
+if size(pokedR) ~= size(xregressor)
+    error("one choice for each regressor, something is wrong!")
+end
+x_uniq = unique(xregressor);
+
+nPokedR = [];
+nPokes = [];
+
+for c = 1:numel(x_uniq)
+
+    idx = x_uniq(c) == xregressor;
+
+    nPokedR(c,1) = sum(pokedR(idx)); 
+    nPokes(c,1) = sum(idx);
+end
+
+% *** bounds ***
+lowerbound = [0, -range(xregressor)]; % lapse [0,1], sensitivity [0 5], bias [-40, 40]
+upperbound = [20,  range(xregressor)];
+% lowerbound = [];
+% upperbound = [];
+
+
+% *** initial parameters ***
+param0 = [0.5, 0]; % no lapse, mediocre sensitivity, no bias
+
+warning('off', 'MATLAB:nchoosek:LargeCoefficient')
+
+% ***************** %
+%       Fmincon     %
+% ***************** %
+nLL0 = negLogLike_logistic2(param0, nPokes, nPokedR, x_uniq);
+if isinf(nLL0) || isnan(nLL0)
+    error('Failure to evaluate initial values')
+end
+
+[x_fmin, f_fmin, exitflag, output, ~, grad, hessian] = ...
+    fmincon(@(param) negLogLike_logistic2(param, nPokes, nPokedR, x_uniq), ...
+            param0, ...
+            [], [], [], [], ...
+            lowerbound, upperbound, ...
+            [], ...
+            optimset('DiffMinChange', 0.00001, ...
+                     'MaxIter', 1000, ...
+                     'Algorithm', 'interior-point', ...
+                     'Display', 'off')     );
+P = struct;
+P.sens   = x_fmin(1);
+P.bias   = x_fmin(2);
+
+% *** vargout ***
+varargout{1} = f_fmin;
+varargout{2} = exitflag;
+varargout{3} = output;
+varargout{4} = grad;
+varargout{5} = hessian;
+end

PBups_psychometrics/fit_logistic3.m

@@ -0,0 +1,61 @@
+function [P varargout] = fit_logistic3(pokedR, xregressor)
+
+
+if size(pokedR) ~= size(xregressor)
+    error("one choice for each regressor, something is wrong!")
+end
+x_uniq = unique(xregressor);
+
+nPokedR = [];
+nPokes = [];
+
+for c = 1:numel(x_uniq)
+
+    idx = x_uniq(c) == xregressor;
+
+    nPokedR(c,1) = sum(pokedR(idx)); 
+    nPokes(c,1) = sum(idx);
+end
+
+% *** bounds ***
+lowerbound = [0, 0, -range(xregressor)]; % lapse [0,1], sensitivity [0 5], bias [-40, 40]
+upperbound = [1, 20,  range(xregressor)];
+% lowerbound = [];
+% upperbound = [];
+
+
+% *** initial parameters ***
+param0 = [0.1, 0.5, 0]; % no lapse, mediocre sensitivity, no bias
+
+warning('off', 'MATLAB:nchoosek:LargeCoefficient')
+
+% ***************** %
+%       Fmincon     %
+% ***************** %
+nLL0 = negLogLike_logistic3(param0, nPokes, nPokedR, x_uniq);
+if isinf(nLL0) || isnan(nLL0)
+    error('Failure to evaluate initial values')
+end
+
+[x_fmin, f_fmin, exitflag, output, ~, grad, hessian] = ...
+    fmincon(@(param) negLogLike_logistic3(param, nPokes, nPokedR, x_uniq), ...
+            param0, ...
+            [], [], [], [], ...
+            lowerbound, upperbound, ...
+            [], ...
+            optimset('DiffMinChange', 0.00001, ...
+                     'MaxIter', 1000, ...
+                     'Algorithm', 'interior-point', ...
+                     'Display', 'off')     );
+P = struct;
+P.lapse  = x_fmin(1);
+P.sens   = x_fmin(2);
+P.bias   = x_fmin(3);
+
+% *** vargout ***
+varargout{1} = f_fmin;
+varargout{2} = exitflag;
+varargout{3} = output;
+varargout{4} = grad;
+varargout{5} = hessian;
+end

PBups_psychometrics/fit_logistic4.m

@@ -0,0 +1,66 @@
+function [P varargout] = fit_logistic4(pokedR, xregressor)
+
+rng shuffle
+if size(pokedR) ~= size(xregressor)
+    error("one choice for each regressor, something is wrong!")
+end
+x_uniq = unique(xregressor);
+
+nPokedR = [];
+nPokes = [];
+
+for c = 1:numel(x_uniq)
+
+    idx = x_uniq(c) == xregressor;
+
+    nPokedR(c,1) = sum(pokedR(idx)); 
+    nPokes(c,1) = sum(idx);
+end
+
+% *** bounds ***
+lowerbound = [0, 0, 0.0001*range(xregressor), -range(xregressor)]; % lapse [0,1], sensitivity [0 5], bias [-40, 40]
+upperbound = [1.0, 1.0, 100*range(xregressor), range(xregressor)];
+
+
+function [c,ceq] = boundcon(x)
+    c(1) = x(1) + x(2) - 1; % gamm0 + gamm1  <= 1
+    c(2) = -(x(1) + x(2));  % gamm0 + gamm1 >= 0
+    ceq = [];
+end
+
+% *** initial parameters ***
+param0 = [0.01, 0.9, 0.01, 0]; % low lapse, mediocre sensitivity, no bias
+
+warning('off', 'MATLAB:nchoosek:LargeCoefficient')
+
+% ***************** %
+%       Fmincon     %
+% ***************** %
+nLL0 = negLogLike_logistic4(param0, nPokes, nPokedR, x_uniq);
+if isinf(nLL0) || isnan(nLL0)
+    error('Failure to evaluate initial values')
+end
+
+[x_fmin, f_fmin, exitflag, output, ~, grad, hessian] = ...
+    fmincon(@(param) negLogLike_logistic4(param, nPokes, nPokedR, x_uniq), ...
+            param0, ...
+            [], [], [], [], ...
+            lowerbound, upperbound, ...
+            @boundcon, ...
+            optimset('DiffMinChange', 1e-12, ...
+                     'MaxIter', 5000, ...
+                     'Algorithm', 'interior-point', ...
+                     'Display', 'off')     );
+P = struct;
+P.gamma0 = x_fmin(1);
+P.gamma1 = x_fmin(2);
+P.sens   = x_fmin(3);
+P.bias   = x_fmin(4);
+
+% *** vargout ***
+varargout{1} = f_fmin;
+varargout{2} = exitflag;
+varargout{3} = output;
+varargout{4} = grad;
+varargout{5} = hessian;
+end