added non-decision time

fit_bandit.m

@@ -21,7 +21,7 @@
     param(1).ub = 20;   % upper bound
 
     % learning rate
-    param(2).name = 'lr_pos';
+    param(2).name = 'lr';
     param(2).hp = [1.2 1.2];    % hyperparameters of beta prior
     param(2).logpdf = @(x) sum(log(betapdf(x,param(2).hp(1),param(2).hp(2))));
     param(2).lb = 0;
@@ -33,11 +33,12 @@
     param(3).lb = 1e-3;
     param(3).ub = 20;
 
+    % non-decision time
+    param(4).name = 'T';
+    param(4).logpdf = @(x) 0;
+    param(4).lb = 0;
+    param(4).ub = 1;
+
     % fit model
     f = @(x,data) likfun_bandit(x,data);    % log-likelihood function
-    results = mfit_optimize(f,param,data);
-
-    % get latent variables
-    for i = 1:length(data)
-        [~,results.latents(i)] = likfun_bandit(results.x(i,:),data(i));
-    end
+    results = mfit_optimize(f,param,data);

fit_discount.m

@@ -33,6 +33,12 @@
     param(3).lb = 1e-3;
     param(3).ub = 20;
 
+    % non-decision time
+    param(4).name = 'T';
+    param(4).logpdf = @(x) 0;
+    param(4).lb = 0;
+    param(4).ub = 1;
+
     % fit model
     f = @(x,data) likfun_discount(x,data);    % log-likelihood function
     results = mfit_optimize(f,param,data);

fit_gonogo.m

@@ -41,11 +41,12 @@
     param(5).lb = 1e-3;
     param(5).ub = 20;
 
+    % non-decision time
+    param(6).name = 'T';
+    param(6).logpdf = @(x) 0;
+    param(6).lb = 0;
+    param(6).ub = 1;
+
     % fit model
     f = @(x,data) likfun_gonogo(x,data);    % log-likelihood function
-    results = mfit_optimize(f,param,data);
-
-    % get latent variables
-    for i = 1:length(data)
-        [~,results.latents(i)] = likfun_gonogo(results.x(i,:),data(i));
-    end
+    results = mfit_optimize(f,param,data);

fit_multi.m

@@ -26,11 +26,12 @@
     param(2).lb = 1e-3;
     param(2).ub = 40;
 
+    % non-decision time
+    param(3).name = 'T';
+    param(3).logpdf = @(x) 0;
+    param(3).lb = 0;
+    param(3).ub = 1;
+
     % fit model
     f = @(x,data) likfun_multi(x,data);    % log-likelihood function
-    results = mfit_optimize(f,param,data);
-
-    % get latent variables
-    for i = 1:length(data)
-        [~,results.latents(i)] = likfun_bandit(results.x(i,:),data(i));
-    end
+    results = mfit_optimize(f,param,data);

likfun_bandit.m

@@ -9,6 +9,7 @@
     %       x(1) - drift rate differential action value weight (b)
     %       x(2) - learning rate for state-action values (alpha)
     %       x(3) - decision threshold (a)
+    %       x(4) - non-decision time (T)
     %   data - structure with the following fields
     %           .c - [N x 1] choices
     %           .r - [N x 1] rewards
@@ -20,20 +21,23 @@
     % OUTPUTS:
     %   lik - log-likelihood
     %   latents - structure with the following fields:
-    %           .P - [N x 1] action probability
     %           .v - [N x 1] drift rate
+    %           .P - [N x 1] probability of chosen option
+    %           .RT_mean - [N x 1] mean response time for chosen option
     %
-    % Sam Gershman, Nov 2015
+    % Sam Gershman, Aug 2016
 
     % set parameters
     b = x(1);           % drift rate differential action value weight
     lr = x(2);          % learning rate
     a = x(3);           % decision threshold
+    T = x(4);           % non-decision time
 
     % initialization
     lik = 0; C = data.C;
     S = length(unique(data.s)); % number of states
     Q = zeros(S,C);    % initial state-action values
+    data.rt = data.rt - T;
 
     for n = 1:data.N
 
@@ -46,11 +50,8 @@
         v = b*(Q(s,2)-Q(s,1));
 
         % accumulate log-likelihod
-        if data.c(n) == 2
-            P = wfpt(data.rt(n),-v,a);  % Wiener first passage time distribution
-        else
-            P = wfpt(data.rt(n),v,a);
-        end
+        if data.c(n) == 1; v = -v; end
+        P = wfpt(data.rt(n),-v,a);  % Wiener first passage time distribution
         if isnan(P) || P==0; P = realmin; end % avoid NaNs and zeros in the logarithm
         lik = lik + log(P);
 
@@ -60,7 +61,8 @@
         % store latent variables
         if nargout > 1
             latents.v(n,1) = v;
-            latents.P(n,1) = P;
+            latents.P(n,1) = 1/(1+exp(-a*v));
+            latents.RT_mean(n,1) = (0.5*a/v)*tanh(0.5*a*v)+T;
         end
 
     end

likfun_discount.m

@@ -9,6 +9,7 @@
     %       x(1) - drift rate differential action value weight (b)
     %       x(2) - discount parameter (k)
     %       x(3) - decision threshold (a)
+    %       x(4) - non-decision time (T)
     %   data - structure with the following fields
     %           .c - [N x 1] choices
     %           .r - [N x 2] reward for each option
@@ -20,19 +21,22 @@
     % OUTPUTS:
     %   lik - log-likelihood
     %   latents - structure with the following fields:
-    %           .lik - [N x 1] likelihood
-    %           .P - [N x 1] probability of choosing option 1
-    %           .RT_mean - [N x 1] mean response time
     %           .v - [N x 1] drift rate
+    %           .P - [N x 1] probability of chosen option
+    %           .RT_mean - [N x 1] mean response time for chosen option
     %
     % Sam Gershman, Nov 2015
 
     % set parameters
     b = x(1);           % drift rate differential action value weight
     k = x(2);           % discount parameter
     a = x(3);           % decision threshold
+    T = x(4);           % non-decision time
 
+    % initialization
     lik = 0;
+    data.rt = data.rt - T;
+
     for n = 1:data.N
 
         % data for current trial
@@ -44,20 +48,16 @@
         v = b*(diff(V));
 
         % accumulate log-likelihod
-        if data.c(n) == 2
-            P = wfpt(data.rt(n),-v,a);  % Wiener first passage time distribution
-        else
-            P = wfpt(data.rt(n),v,a);
-        end
+        if data.c(n) == 1; v = -v; end
+        P = wfpt(data.rt(n),-v,a);  % Wiener first passage time distribution
         if isnan(P) || P==0; P = realmin; end % avoid NaNs and zeros in the logarithm
         lik = lik + log(P);
 
         % store latent variables
         if nargout > 1
             latents.v(n,1) = v;
-            latents.lik(n,1) = P;
             latents.P(n,1) = 1/(1+exp(-a*v));
-            latents.RT_mean(n,1) = (0.5*a/v)*tanh(0.5*a*v);
+            latents.RT_mean(n,1) = (0.5*a/v)*tanh(0.5*a*v)+T;
         end
 
     end

likfun_gonogo.m

@@ -33,12 +33,14 @@
     b3 = x(3);          % drift rate Pavlovian bias weight
     lr = x(4);          % learning rate
     a = x(5);           % decision threshold
+    T = x(6);           % non-decision time
 
     % initialization
     lik = 0; C = data.C;
     S = length(unique(data.s)); % number of states
     Q = zeros(S,C);    % initial state-action values
     V = zeros(S,1);    % initial state values
+    data.rt = data.rt - T;
     mx = max(data.rt)+0.1;  % max RT
 
     for n = 1:data.N
@@ -68,6 +70,8 @@
         % store latent variables
         if nargout > 1
             latents.v(n,1) = v;
+            latents.P(n,1) = 1/(1+exp(-a*v));
+            latents.RT_mean(n,1) = (0.5*a/v)*tanh(0.5*a*v)+T;
         end
 
     end

likfun_multi.m

@@ -9,6 +9,7 @@
     %   x - parameters:
     %       x(1) - drift rate value weight (b)
     %       x(2) - decision threshold (a)
+    %       x(3) - non-decision time (T)
     %   data - structure with the following fields
     %           .c - [N x 1] choices
     %           .V - [N x C] values
@@ -19,17 +20,18 @@
     % OUTPUTS:
     %   lik - log-likelihood
     %   latents - structure with the following fields:
-    %           .logP - [N x 1] action log probability
     %           .v - [N x C] drift rates
     %
     % Sam Gershman, Nov 2015
 
     % set parameters
     b = x(1);           % drift rate action value weight
     a = x(2);           % decision threshold
+    T = x(3);           % non-decision time
 
     % initialization
     lik = 0; C = data.C;
+    data.rt = data.rt - T;
 
     for n = 1:data.N
 
@@ -51,7 +53,6 @@
         % store latent variables
         if nargout > 1
             latents.v(n,:) = v;
-            latents.logP(n,1) = logP;
         end
 
     end