Add files via upload

Updated the way forgetting rates are implemented. The forgetting rate parameter (omega) now takes values between 0 and 1, where higher values lead to greater forgetting. The updated implementation also only applies forgetting to concentration parameter values added to the initial prior values. This prevents the concentration parameters from moving to values that are implausibly low, which can cause numerical problems when running the code.

Author: rssmith33

Simplified_simulation_script.m

@@ -6,7 +6,7 @@
 % Application to Empirical Data
 
 % By: Ryan Smith, Karl J. Friston, Christopher J. Whyte
-
+% UPDATED: 8/28/2024 (modified forgetting rate implementation)
 rng('shuffle')
 close all
 clear
@@ -36,11 +36,11 @@
                % function starting on line 810. This includes, among
                % others (similar to in the main tutorial script):
 
-% prior beliefs about context (d): alter line 866
+% prior beliefs about context (d): alter line 876
 
-% beliefs about hint accuracy in the likelihood (a): alter lines 986-988
+% beliefs about hint accuracy in the likelihood (a): alter lines 996-998
 
-% to adjust habits (e), alter line 1145
+% to adjust habits (e), alter line 1155
 
 %% Specify Generative Model
 
@@ -336,7 +336,7 @@
             for modality = 1:NumModalities
                 % prior preferences about outcomes
                 predictive_observations_posterior = cell_md_dot(a{modality},Expected_states(:)); %posterior over observations
-                Gintermediate(policy) = Gintermediate(policy) + predictive_observations_posterior'*(C{modality}(:,timestep));
+                Gintermediate(policy) = Gintermediate(policy) + predictive_observations_posterior'*(C{modality}(:,t));
 
                 % Bayesian surprise about parameters 
                 if isfield(MDP,'a')
@@ -445,7 +445,7 @@
                 a_learning = spm_cross(a_learning,BMA_states{factor}(:,t));
             end
             a_learning = a_learning.*(MDP.a{modality} > 0);
-            MDP.a{modality} = MDP.a{modality}*omega + a_learning*eta;
+            MDP.a{modality} = (MDP.a{modality}-MDP.a_0{modality})*(1-omega) + MDP.a_0{modality} + a_learning*eta;
         end
     end 
 end 
@@ -454,13 +454,13 @@
 if isfield(MDP,'d')
     for factor = 1:NumFactors
         i = MDP.d{factor} > 0;
-        MDP.d{factor}(i) = omega*MDP.d{factor}(i) + eta*BMA_states{factor}(i,1);
+        MDP.d{factor}(i) = (1-omega)*(MDP.d{factor}(i)-MDP.d_0{factor}(i)) + MDP.d_0{factor}(i) + eta*BMA_states{factor}(i,1);
     end
 end
 
 % policies e (habits)
 if isfield(MDP,'e')
-    MDP.e = omega*MDP.e + eta*policy_posterior(:,T);
+    MDP.e = (1-omega)*(MDP.e - MDP.e_0) + MDP.e_0 + eta*policy_posterior(:,T);
 end
 
 % Free energy of concentration parameters
@@ -1163,9 +1163,9 @@
      eta = 1; % Default (maximum) learning rate
 
 % Omega: forgetting rate (0-1) controlling the magnitude of reduction in concentration
-% parameter values after each trial (if learning is enabled).
+% parameter values after each trial (if learning is enabled). NOTE THE FORM OF FORGETTING IMPLEMENTED HERE IS MODIFIED FROM THE DESCRIPTION IN THE PUBLISHED TUTORIAL FOR IMPROVED PERFORMANCE.
 
-     omega = 1; % Default value indicating there is no forgetting (values < 1 indicate forgetting)
+     omega = 0; % Default value indicating there is no forgetting (values approaching 1 indicate forgetting)
 
 % Beta: Expected precision of expected free energy (G) over policies (a 
 % positive value, with higher values indicating lower expected precision).
@@ -1194,13 +1194,13 @@
 mdp.B = B;                    % transition probabilities
 mdp.C = C;                    % preferred states
 mdp.D = D;                    % priors over initial states
-mdp.d = d;                    % enable learning priors over initial states
-
+mdp.d = d; mdp.d_0 = d;       % enable learning priors over initial states
+                              % d_0 is floor value for forgetting
 if Gen_model == 1
     mdp.E = E;                % prior over policies
 elseif Gen_model == 2
-    mdp.a = a;                % enable learning state-outcome mappings
-    mdp.e = e;                % enable learning of prior over policies
+    mdp.a = a; mdp.a_0 = a;   % enable learning state-outcome mappings and set floor value for forgetting (a_0)
+    mdp.e = e; mdp.e_0 = e;   % enable learning of prior over policies and set floor value for forgetting (e_0)
 end 
 
 mdp.eta = eta;                % learning rate
@@ -1225,4 +1225,4 @@
 
 MDP = mdp;
 
-end
+end

Step_by_Step_AI_Guide.m

@@ -4,6 +4,7 @@
 % Application to Empirical Data
 
 % By: Ryan Smith, Karl J. Friston, Christopher J. Whyte
+% UPDATED: 8/28/2024 (modified forgetting rate implementation)
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
@@ -48,7 +49,7 @@
          % To reproduce fig. 11, use values of 3 or 4 (with Sim = 3)
          % This will have no effect on Sim = 4 or Sim = 5
 
-Sim = 1;
+Sim = 2;
 
 % When Sim = 5, if PEB = 1 the script will run simulated group-level
 % (Parametric Empirical Bayes) analyses.
@@ -358,7 +359,7 @@
 
 % Note that, expanded out, this means that the other C-matrices will be:
 
-% C{1} =      [0 0 0;     % No Hint
+% C{1} =      [0 0 0;    % No Hint
 %              0 0 0;    % Machine-Left Hint
 %              0 0 0];   % Machine-Right Hint
 % 
@@ -376,7 +377,18 @@
 % This will not be simulated here. However, this works by increasing the
 % preference magnitude for an outcome each time that outcome is observed.
 % The assumption here is that preferences naturally increase for entering
-% situations that are more familiar.
+% situations that are more familiar. To do so, you can specify starting
+% concentration parameters. For example:
+
+% c{1}      = zeros(No(1),T); % Hints
+% c{2}      = zeros(No(2),T); % Wins/Losses
+% c{3}      = zeros(No(3),T); % Observed Behaviors
+% 
+% c{2}(:,:) =    [1  1  1  ;  % Null
+%                 1  0  0.5;  % Loss
+%                 1  2  1.5]; % win
+
+% NOTE: These values must be non-negative; higher values = more preferred
 
 % Allowable policies: U or V. 
 %==========================================================================
@@ -461,12 +473,17 @@
 % degree to which newer experience can 'over-write' what has been learned
 % from older experiences. It is adaptive in environments where the true
 % parameters in the generative process (priors, likelihoods, etc.) can
-% change over time. A low value for omega can be seen as a prior that the
+% change over time. A high value for omega can be seen as a prior that the
 % world is volatile and that contingencies change over time.
 
-    omega = 1; % By default we here set this to 1 (indicating no forgetting, 
+  omega = 0.0; % By default we here set this to 0 (indicating no forgetting, 
                % but try changing its value to see how it affects model behavior. 
-               % Values below 1 indicate greater rates of forgetting.
+               % Values approaching 1 indicate greater rates of forgetting.
+               % NOTE: Trial 1 concentration parameter values are set as
+               % floor values (forgetting cannot reduce counts below those
+               % values - THIS IS MODIFIED FROM THE PUBLISHED TUTORIAL VERSION 
+               % SO THAT CONCENTRATION PARAMETERS ABOVE THE FLOOR VALUE 
+               % ARE MULTIPLIED BY 1-OMEGA)
 
 % Beta: Expected precision of expected free energy (G) over policies (a 
 % positive value, with higher values indicating lower expected precision).
@@ -577,8 +594,8 @@
 mdp.C = C;                    % preferred states
 mdp.D = D;                    % priors over initial states
 
-mdp.d = d;                    % enable learning priors over initial states
-    
+mdp.d = d; mdp.d_0 = mdp.d;   % enable learning priors over initial states
+                              %     and set lower bound on concentration paramaters (d_0)
 mdp.eta = eta;                % learning rate
 mdp.omega = omega;            % forgetting rate
 mdp.alpha = alpha;            % action precision
@@ -591,10 +608,10 @@
     % mdp.E = E;
 
 % or learning other parameters:
-    % mdp.a = a;                    
-    % mdp.b = b;
-    % mdp.c = c;
-    % mdp.e = e;         
+    % mdp.a = a;  mdp.a_0 = mdp.a;                  
+    % mdp.b = b;  mdp.b_0 = mdp.b;
+    % mdp.c = c;  mdp.c_0 = mdp.c; clear mdp.C = C;
+    % mdp.e = e;  mdp.e_0 = mdp.e;        
 
 % or specifying true states or outcomes:
 
@@ -1425,4 +1442,4 @@
 % now build a generative model of a task, run simulations, assess parameter
 % recoverability, do bayesian model comparison, and do hierarchical
 % bayesian group analyses. See the main text for further explanation of
-% other aspects of these steps.
+% other aspects of these steps.

spm_MDP_VB_X_tutorial.m

@@ -1,4 +1,5 @@
 function [MDP] = spm_MDP_VB_X_tutorial(MDP,OPTIONS)
+% UPDATED: 8/28/2024 (modified forgetting rate implementation)
 
 % active inference and learning using variational message passing
 % FORMAT [MDP] = spm_MDP_VB_X_tutorial(MDP,OPTIONS)
@@ -149,7 +150,7 @@
 
 % check MDP specification
 %--------------------------------------------------------------------------
-MDP = spm_MDP_check(MDP);
+ MDP = spm_MDP_check(MDP);
 
 % handle multiple trials, ensuring parameters (and posteriors) are updated
 %==========================================================================
@@ -225,7 +226,7 @@
 try, omega = MDP(1).omega; catch, omega = 1;    end % forgetting rate
 try, tau   = MDP(1).tau;   catch, tau   = 4;    end % update time constant
 try, chi   = MDP(1).chi;   catch, chi   = 1/64; end % Occam window updates
-try, erp   = MDP(1).erp;   catch, erp   = 4;    end % update reset
+try, erp   = MDP(1).erp;   catch, erp   = 4;    end % update reset 
 
 % preclude precision updates for moving policies
 %--------------------------------------------------------------------------
@@ -258,8 +259,8 @@
     end
     for g = 1:Ng(m)
         No(m,g) = size(MDP(m).A{g},1);     % number of outcomes
-    end
-    
+    end    
+
     % parameters of generative model and policies
     %======================================================================
 
@@ -306,7 +307,6 @@
                 sB{m,f}(:,:,j) = spm_norm(MDP(m).B{f}(:,:,j) );
                 rB{m,f}(:,:,j) = spm_norm(MDP(m).B{f}(:,:,j)');
             end
-            
         end
 
         % prior concentration paramters for complexity
@@ -341,7 +341,7 @@
     % priors over policies - concentration parameters
     %----------------------------------------------------------------------
     if isfield(MDP,'e')
-        E{m} = spm_norm(MDP(m).e);
+        E{m} = spm_norm(MDP(m).e); 
     elseif isfield(MDP,'E')
         E{m} = spm_norm(MDP(m).E);
     else
@@ -377,6 +377,7 @@
             end
         end
         C{m,g} = spm_log(spm_softmax(C{m,g}));
+
     end
 
     % initialise  posterior expectations of hidden states
@@ -1123,7 +1124,7 @@
                     da = spm_cross(da,X{m,f}(:,t));
                 end
                 da     = da.*(MDP(m).a{g} > 0);
-                MDP(m).a{g} = MDP(m).a{g}*omega + da*eta;
+                MDP(m).a{g} = (MDP(m).a{g}-MDP(m).a_0{g})*(1-omega) + MDP(m).a_0{g} + da*eta;
             end
         end
 
@@ -1135,7 +1136,7 @@
                     v   = V{m}(t - 1,k,f);
                     db  = u{m}(k,t)*x{m,f}(:,t,k)*x{m,f}(:,t - 1,k)';
                     db  = db.*(MDP(m).b{f}(:,:,v) > 0);
-                    MDP(m).b{f}(:,:,v) = MDP(m).b{f}(:,:,v)*omega + db*eta;
+                    MDP(m).b{f}(:,:,v) = (MDP(m).b{f}(:,:,v)-MDP(m).b_0{f}(:,:,v))*(1-omega) + MDP(m).b_0{f}(:,:,v) + db*eta;
                 end
             end
         end
@@ -1147,10 +1148,10 @@
                 dc = O{m}{g,t};
                 if size(MDP(m).c{g},2) > 1
                     dc = dc.*(MDP(m).c{g}(:,t) > 0);
-                    MDP(m).c{g}(:,t) = MDP(m).c{g}(:,t)*omega + dc*eta;
+                    MDP(m).c{g}(:,t) = (MDP(m).c{g}(:,t)-MDP(m).c_0{g}(:,t))*(1-omega) + MDP(m).c_0{g}(:,t) + dc*eta;
                 else
                     dc = dc.*(MDP(m).c{g}>0);
-                    MDP(m).c{g} = MDP(m).c{g}*omega + dc*eta;
+                    MDP(m).c{g} = (MDP(m).c{g}-c_0{g})*(1-omega) + c_0{g} + dc*eta;
                 end
             end
         end
@@ -1161,14 +1162,14 @@
     if isfield(MDP,'d')
         for f = 1:Nf(m)
             i = MDP(m).d{f} > 0;
-            MDP(m).d{f}(i) = MDP(m).d{f}(i)*omega + X{m,f}(i,1)*eta;
+            MDP(m).d{f}(i) = (MDP(m).d{f}(i)-MDP(m).d_0{f}(i))*(1-omega) + MDP(m).d_0{f}(i) + X{m,f}(i,1)*eta;
         end
     end
 
     % policies
     %----------------------------------------------------------------------
     if isfield(MDP,'e')
-        MDP(m).e = MDP(m).e*omega + eta*u{m}(:,T);
+        MDP(m).e = (MDP(m).e-MDP(m).e_0)*(1-omega) + MDP(m).e_0 + eta*u{m}(:,T);
     end
 
     % (negative) free energy of parameters (complexity): outcome specific
@@ -1682,4 +1683,3 @@
 
 end
 
-