-
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
- More advanced information
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 using the dot notation:
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')
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 set of parameters and perturbation sequence, as inference in the COIN model depends on state feedback, which is generated randomly within the simulate_COIN method. Such variability in model behaviour is often undesirable. To produce consistent model behaviour for a given set of parameters and perturbation sequence, 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 example, for the motor output:
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). You may want to use a computer cluster.
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 example, to plot the predicted state distribution for each context, the predicted probabilities and the overall predicted state distribution:
obj = COIN
obj.perturbations = [zeros(1,50) ones(1,125) -ones(1,15) NaN(1,150)];
obj.runs = 10;
obj.max_cores = feature('numcores');
obj.plot_state_given_context = true;
obj.plot_predicted_probabilities = true;
obj.plot_state = true;
obj.state_values = linspace(-1.5,1.5,500); % points to evaluate the ‘state | context’ distribution at
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. In general, using more runs will result in cleaner, less variable plots. 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.
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 as a function of the model parameters. In the COIN model, this function is performed 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, each conditioned on a different random sequence of 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 used to calculate the negative log-likelihood via the 'runs' property; to avoid excessive runtimes, multiple runs should be performed in parallel. 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 code showing how to fit the parameters of the COIN model to the data of an individual participant, see fit_parameters_to_data.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 pass the adaptation data of the participant to the class object 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 maximum-likelihood (fitted) values via the appropriate properties.
To view view code showing how to fit the parameters of the COIN model to the average data of a group of participants, see fit_parameters_to_data.m (specify the number of participants in the group at the top of the file; BADS installation required). Note that this method assumes that the average adaptation data across participants can be calculated without reordering the adaptation data of each participant (the nth average adaptation data point of the group is calculated as 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 of each experiment should be calculated and the sum of the log-likelihoods passed to an optimiser.
The above information covers the main features of the COIN model and is sufficient to reproduce the simulation results in the COIN paper. For users who would like to explore the COIN model more themselves and compute additional quantities not shown in the COIN paper, the following technical details are provided.
Inference in the COIN model is performed in a sequential manner by recursively propagating the joint posterior distribution of all latent quantities from one trial to the next after each new set of observations is made. Due to the Markov property, these recursive updates depend only on the most recent posterior distribution and not on all past posterior distributions. Therefore, to reduce memory requirements, the COIN model does not store all past inferences in memory unless a request is made to do so. To make such a request, add the names of variables to be stored 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 structure OUTPUT. For a full list of the names of variables that can be stored in memory, call help COIN and see the 'VARIABLES' section. Note that variables related to predictive distributions (specifically, the predicted probabilities and the means and variances of the predicted state and state feedback distributions for each context) are stored before particles are resampled in the COIN model (so that they don’t depend on the current state feedback), whereas other variables are stored after particles are resampled; 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 of each COIN model variable over runs (as different runs can have different weights). 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.max_cores = feature('numcores');
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 merely 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 resample steps (e.g. before vs. after the resample step on the same trial), the same set of indices in each of the two stored variables may map to different sets of particles. The map_particles method can be used to find two sets of indices (one per variable) that map to the same set of particles:
i_particles = obj.map_particles(OUTPUT,variable_names,trials);
This method takes as input the 'OUTPUT' of a COIN model simulation, a cell array called 'variable_names' that contains the names of the two stored variables and a vector called 'trials' that specifies, for each variable, the trial on which it was stored. The output contains, for each run, a set of indices for each of the two variables that map to the same set of particles. For 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.max_cores = feature('numcores');
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 = [trial_state_mean trial_responsibilities];
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 resample step, and hence they share the same particle indices (when analysed on the same trial).
Contexts that have the same label in different particles cannot, in general, be assumed to represent the same entity due to the 'label-switching problem' (see Section 2.9.2 of [link to SM]). The label-switching 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 the context labels (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 labels of contexts in each particle such that these labels can be treated as equivalent across particles.