-
Notifications
You must be signed in to change notification settings - Fork 6
Home
- The COIN class
- Running a COIN model simulation
- Performing multiple runs of a COIN model simulation
- Plotting COIN model variables
- Fitting the parameters of the COIN model to data
- Storing COIN model variables for offline analysis
- Averaging COIN model variables over particles and runs
The COIN model is implemented as a class in MATLAB called COIN. An object of the class can be created by calling the class name:
obj = COIN
This object has properties that define the model (e.g., model parameters) and the paradigm (e.g., perturbations, sensory cues). Some properties have default values (set at the top of COIN.m), which can be overwritten. Data stored in properties can be operated on by methods (functions) of the class.
For a full list of properties of the COIN class, call help COIN (see 'PROPERTIES' section).
To run a COIN model simulation, first define a trial-by-trial sequence of perturbations via the 'perturbations' property (use NaN to indicate a channel trial). For example, for a spontaneous recovery paradigm:
obj.perturbations = [zeros(1,50) ones(1,125) -ones(1,15) NaN(1,150)];
If an experiment involves sensory cues, these can also be provided via the 'cues' property. Next call the simulate_COIN method on the class object 'obj':
OUTPUT = obj.simulate_COIN;
The simulation results are contained in the structure OUTPUT. The state feedback and motor output of the COIN model can be plotted as follows:
figure
plot(OUTPUT.runs{1}.state_feedback,'.')
hold on
plot(OUTPUT.runs{1}.motor_output)
legend('state feedback','motor output')
Note that the state feedback observations are generated randomly within the simulate_COIN method by adding observation noise (sampled from the prior) to the user-defined perturbations.
Multiple calls to the simulate_COIN method will produce different results, even for the same perturbation sequence and parameters, as inference in the COIN model depends on state feedback that is generated randomly within the simulate_COIN method. Such variability in model behaviour is often undesirable. To produce consistent model behaviour for a given perturbation sequence and set of parameters, multiple runs of a simulation can be performed, each conditioned on a different state feedback sequence, and the results of the runs averaged. In general, averaging over more runs will produce cleaner, less variable results.
To specify the number of runs to perform, use the 'runs' property. For example, to perform 10 runs:
obj = COIN
obj.perturbations = [zeros(1,50) ones(1,125) -ones(1,15) NaN(1,150)];
obj.runs = 10;
OUTPUT = obj.simulate_COIN;
OUTPUT.runs is a cell array with one cell per run and OUTPUT.weights is a vector specifying the relative weight of each run (weights are non-negative and sum to 1). These weights should be used to compute a weighted average of COIN model variables over runs:
for run = 1:obj.runs
motor_output(:,run) = OUTPUT.runs{run}.motor_output;
end
plot(sum(OUTPUT.weights.*motor_output,2))
In this example, each run is assigned an equal weight, as no adaptation data was passed to the class object. However, if the parameters of the COIN model have been first fit to adaptation data, and this data is passed to the class object before calling the simulate_COIN method, each run will be assigned a weight based on the likelihood of the adaptation data (see Performing multiple runs of a COIN model simulation conditioned on the adaptation data of a participant).
The time it takes to execute multiple runs of a COIN model simulation in series (e.g., using a for loop) can be prohibitively long if the number of runs is large. To execute multiple runs in parallel across different CPU cores, and thus reduce the model runtime, specify the maximum number of CPU cores available using the 'max_cores' property (the default setting of max_cores is 0, which executes multiple runs in series).
The distributions or expected values of variables inferred by the COIN model can be plotted by activating the relevant properties of the class object. For a full list of variables that can be plotted and how to plot them, call help COIN, and see the 'plot flags' and 'plot inputs' subsections of the 'PROPERTIES' section. As an example, to plot the predicted state distribution for each context and the predicted probabilities, call:
obj = COIN
obj.perturbations = [zeros(1,50) ones(1,125) -ones(1,15) NaN(1,150)];
obj.runs = 10;
obj.plot_state_given_context = true;
obj.state_values = linspace(-1.5,1.5,500); % points to evaluate the ‘state | context’ distribution at
obj.plot_predicted_probabilities = true;
OUTPUT = obj.simulate_COIN;
Note that for distributions of continuous variables, the points at which to evaluate the distribution can be specified (as shown in the example above). Next call the plot_COIN method:
PLOTS = obj.plot_COIN(OUTPUT);
This will generate the requested plots—a state | context plot and a predicted probabilities plot in this example. The structure PLOTS contains the data that is plotted (to see how the data is plotted, view the generate_figures method in COIN.m). Note that these plots may take some time to generate, as they require contexts to be relabelled in multiple particles and multiple runs on each trial. Once the contexts have been relabelled, the relevant distributions or expected values of variables are automatically averaged over particles and runs (as described in Averaging COIN model variables over particles and runs). In general, using more runs will result in cleaner, less variable plots (see Performing multiple runs of a COIN model simulation).
The parameters of the COIN model are fit to data by minimising the negative log of the likelihood function. For an optimiser to perform this minimisation, it must have access to a function that can evaluate the negative log-likelihood. In the COIN model, this function is provided by the objective_COIN method:
negative_log_likelihood = obj.objective_COIN;
The objective_COIN method returns a stochastic estimate of the negative log-likelihood. The estimate is stochastic because it is calculated by simulating one or more runs of the COIN model conditioned on randomly-generated state feedback observations. To aid parameter optimisation, the variance of the estimate of the negative log-likelihood can be reduced by increasing the number of runs via the 'runs' property; this is best done in conjunction with parallel computing to avoid excessive runtimes. Because the estimate of the negative log-likelihood is stochastic, an optimiser that can handle a stochastic objective function should be used to fit the parameters of the COIN model to data (e.g., BADS).
To view a working example of how to fit the parameters of the COIN model to the data of an individual participant, see fit_model_to_data_individual.m (BADS installation required).
Performing multiple runs of a COIN model simulation conditioned on the adaptation data of a participant
Once the parameters of the COIN model have been fit to the adaptation data of a participant, each run of a COIN model simulation can be assigned a weight based on the likelihood of the adaptation data. In general, these weights will not be equal, although they can be. To generate a set of weighted runs conditioned on the adaptation data of a participant, follow the instructions described in Performing multiple runs of a COIN model simulation but in addition specify the adaptation data of the participant via the 'adaptation' property (the data should be in vector form with one element per trial—use NaN on trials where adaptation was not measured) and set the parameters of the COIN model to their fitted values via the corresponding properties.
To view a working example of how to fit the parameters of the COIN model to the average data of a group of participants, see fit_model_to_data_group.m (BADS installation required). Note that this method assumes that the average adaptation data over participants can be calculated without first reordering the adaptation data of each participant (i.e., that the nth average adaptation data point of the group is the average nth adaptation data point of each participant).
To fit the parameters of the COIN model to data from multiple experiments simultaneously (e.g., spontaneous and evoked recovery), the negative log-likelihood should be calculated for each experiment individually and the sum of the log-likelihoods passed to an optimiser.
Inference in the COIN model is performed in a sequential manner by recursively propagating the joint posterior distribution from one trial to the next after each new set of observations is made. These recursive updates depend only on the most recent posterior distribution, due to the Markov property, rather than on all past posterior distributions. Therefore, to reduce memory requirements, the COIN model does not store all past inferences in memory by default but only on request. To make such a request, add the names of any variables to be stored in memory to the 'store' property. For example, to store the context responsibilities and the context-specific Kalman gains:
obj = COIN
obj.perturbations = [zeros(1,50) ones(1,125) -ones(1,15) NaN(1,150)];
obj.store = {'responsibilities','Kalman_gains'};
OUTPUT = obj.simulate_COIN;
The stored variables are contained in the cell array OUTPUT. For a full list of the names of variables that can be stored, call help COIN and see the 'VARIABLES' section. Note that some variables (the predicted probabilities and the means and variances of the predicted state and state feedback distributions for each context) are stored before the resampling step of the COIN model inference algorithm (so that they don’t depend on the latest state feedback observation), whereas others are stored after the resampling step. This has implications for equating particles across stored COIN model variables.
Each particle of each run of a COIN model simulation generates its own set of inferences. For plotting purposes, it is useful to compute the expected value of these inferences by averaging over both particles and runs. This is achieved by computing a simple average of each COIN model variable over particles within each run (as all particles of the same run have the same weight) and a weighted average over runs (as different runs can have different weights, see Performing multiple runs of a COIN model simulation). For example, to compute the expected value of the Kalman gain of the context with the highest responsibility (c*):
obj = COIN
obj.perturbations = [zeros(1,50) ones(1,125) -ones(1,15) NaN(1,150)];
obj.runs = 10;
obj.store = {'Kalman_gains','responsibilities'};
OUTPUT = obj.simulate_COIN;
for trial = 1:numel(obj.perturbations)
for run = 1:obj.runs
for particle = 1:obj.particles
[~,j] = max(OUTPUT.runs{run}.responsibilities(:,particle,trial));
Kalman_gain(trial,run,particle) = OUTPUT.runs{run}.Kalman_gains(j,particle,trial);
end
end
end
plot(sum(OUTPUT.weights.*mean(Kalman_gain,3),2)) % weighted average over runs (dimension 2), simple average over particles (dimension 3)
ylabel('Kalman gain | c^*')
Note that this code is not optimised for performance; nested for loops have been used to enhance the readability of the code.
The resample step of the COIN model inference algorithm changes the indices of particles for each COIN model variable. Therefore, if two COIN model variables are stored at times separated by one or more resampling steps (e.g. before vs. after the resampling step on the same trial), the same set of indices may map to different sets of particles in each of the two stored variables. To find equivalent particle indices for any two stored variables, use the map_particles method:
i_particles = obj.map_particles(OUTPUT,variable_names,trials);
This method takes as input the 'OUTPUT' of a COIN model simulation, a cell array 'variable_names' that contains the names of the two stored variables whose particle indices are to be equated and a two-element vector 'trials' that specifies the trial on which each of the two variables was stored. The output of the method contains, for each run, a set of indices for each of the two variables that map to the same set of particles. As an example, to compute the expected value of the mean of the predicted state distribution (stored before the resample step) for the context with the highest responsibility (stored after the resample step):
obj = COIN;
obj.perturbations = [zeros(1,50) ones(1,125) -ones(1,15) NaN(1,150)];
obj.runs = 10;
obj.store = {'state_mean','responsibilities'};
OUTPUT = obj.simulate_COIN;
for trial = 1:numel(obj.perturbations)
variable_names = {'state_mean','responsibilities'};
trial_state_mean = trial;
trial_responsibilities = trial;
trials = [state_mean_trial responsibilities_trial];
i_particles = obj.map_particles(OUTPUT,variable_names,trials);
for run = 1:obj.runs
for particle = 1:obj.particles
[~,j] = max(OUTPUT.runs{run}.responsibilities(:,i_particles.runs{run}.responsibilities(particle),trial_responsibilities));
state(trial,run,particle) = OUTPUT.runs{run}.state_mean(j,i_particles.runs{run}.state_mean(particle),trial_state_mean);
end
end
end
figure
plot(sum(OUTPUT.weights.*mean(state,3),2))
ylabel('E[state | c^*]')
Note that, in the previous example, the map_particles method was not used to compute the Kalman gain of the context with the highest responsibility, as the Kalman_gains and responsibilities variables are both stored after the resampling step, and hence they share the same particle indices when they are analysed on the same trial.
Contexts that have the same index in different particles cannot, in general, be equated due to the 'label-switching problem' (see Section 2.9.2 of [link to SM]). This problem can be addressed in one of two ways. The first way, demonstrated in this section, is to analyse a single context in each particle regardless of its index (e.g., the context with the largest predicted probability or responsibility). The second way, employed by the plot_COIN method in the Plotting COIN model variables section, is to optimally permute the indices of contexts such that they can be treated as equivalent across particles.