-
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
- Storing variables
- Integrating over particles and runs
- Plotting variables
- Fitting the parameters of the COIN model to data
- Parallel computing
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 (use NaN to indicate a channel trial). For example, to simulate 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 object obj:
OUTPUT = obj.simulate_COIN;
The results of the simulation are contained in the structure OUTPUT. The state feedback and adaptation can be plotted as follows:
figure
plot(OUTPUT.runs{1}.y,'.')
hold on
plot(OUTPUT.runs{1}.yHat)
legend('state feedback','adaptation')
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 parameters and perturbation sequence), as inference in the COIN model depends on state feedback that is generated randomly within the simulate_COIN method. To produce consistent results that are independent of the state feedback, multiple runs of a simulation can be performed, each conditioned on a different state feedback sequence, with the results averaged across runs. In general, averaging across more runs will produce cleaner, less variable results.
To specify the number of runs to perform, use the 'runs' property. For example, to perform 2 runs:
obj.runs = 2;
obj.perturbations = [zeros(1,50) ones(1,125) -ones(1,15) NaN(1,150)];
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). The weights should be used to compute a weighted average of the simulation results across runs (see Integrating over particles and runs). 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 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. 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:
- set the parameters of the COIN model to their maximum likelihood values via the corresponding properties (the parameters of the COIN model should first be fit to the adaptation data of the participant, see Fitting the parameters of the COIN model to data),
- pass the adaptation data of the participant to the class object via the 'adaptation' property (the data should be a vector with one element per trial—use NaN on trials where adaptation was not measured),
- specify the number of runs to perform via the 'runs' property,
- and call the simulate_COIN method.
By default, the model runtime scales linearly with the number of runs, as each run is executed in series. To reduce the model runtime by executing multiple runs in parallel, see Parallel computing.
Inference in the COIN model is performed in an online manner by recursively updating the most recent value of each inferred variable. Hence, to run a simulation, the full trial-by-trial sequence of each inferred variable does not need to be stored in memory. Therefore, to reduce memory requirements, the full sequence of each variable is only stored in memory if a request is made to do so. The state feedback and adaptation variables are the only exceptions to this, as they are stored in memory as standard and are provided as part of the basic model output, as demonstrated in Running a COIN model simulation. To store the full trial-by-trial sequence of one or more variables, add the names of these variables to the 'store' property. For example, to store the Kalman gains (named 'k') and responsibilities (named ‘cFilt’), call:
obj.store = {'k','cFilt'};
OUTPUT = obj.simulate_COIN;
The 'store' property must be set before the simulate_COIN method is called (as shown here). Stored variables are contained in the cell array OUTPUT and can be analysed after the model has run (see Integrating over particles and runs for one example analysis). For a full list of the names of variables that can be stored, call help COIN and see the 'VARIABLES' section.
The COIN performs inference using an ensemble of particles, each of which is associated with its own set of inferences. For analysis and plotting purposes, it is useful to integrate over both particles and runs to compute an expected value of the set of inferences. All particles within a run have equal weights, whereas different runs can have different weights (see Performing multiple runs of a COIN model simulation). Therefore, to integrate over particles and runs, a simple average of each inferred variable should be computed across particles within each run and a weighted average across runs.
STATE WHERE YOU CAN FIND PARTICLE DIMENSION IN GENERAL
For example, to compute the expected value of the Kalman gain of the context with the highest responsibility (c*) by integrating over both particles and runs:
obj.runs = 2;
obj.store = {'k','cFilt'};
OUTPUT = obj.simulate_COIN;
for run = 1:obj.runs % loop over runs
for trial = 1:numel(obj.perturbations) % loop over trials
for particle = 1:obj.particles % loop over particles
[~,j] = max(OUTPUT.runs{run}.cFilt(:,particle,trial));
Kalman_gain(trial,run,particle) = OUTPUT.runs{run}.k(j,particle,trial);
end
end
end
figure
plot(sum(OUTPUT.weights.*mean(Kalman_gain,3),2)) % weighted average over runs (dimension 2) and simple average over particles (dimension 3)
ylabel('Kalman gain | c^*')
The distributions or expected values of the main variables inferred by the COIN model can be plotted by activating the relevant properties of the class object. For a full list of inferred 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.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;
The relevant properties must be set before the simulate_COIN method is called (as shown here), so that the variables needed to generate the requested plots can be stored in memory (this is done automatically). Note that for distributions of continuous variables, the points at which to evaluate the distribution can be specified (as shown in the example above).
After the model has run, call the plot_COIN method:
PLOTS = obj.plot_COIN(S);
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 in PLOTS 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 are automatically averaged across particles and runs, as described in Integrating over particles and runs. In general, using more runs will result in cleaner, less variable results (see Performing multiple runs of a COIN model simulation).
The COIN model can be fit either to the data of an individual participant or to the average data of a group of participants. Here the latter case is considered, as an individual participant is a special case of a group with 1 member.
The COIN model is fit to data by finding the parameters that minimise the negative log of the likelihood function. To calculate the negative log-likelihood, create an object array with one object per participant:
S = 8; % number of participants in group
obj(1,S) = COIN; % array of S objects (1 object per member of the group)
For each object, define the model parameters, the paradigm (perturbations and, if applicable, cues) and the adaptation data of the participant using the properties of the object. The 'adaptation' property should be a vector with one element per trial (use NaN on trials where adaptation was not measured). Note that the number of adaptation measurements should be the same for all participants (the nth average adaptation measurement is calculated as the average nth adaptation measurement across participants).
After the properties of the object array have been set, call the objective_COIN method:
negative_log_likelihood = obj.objective_COIN;
This returns a stochastic estimate of the negative log-likelihood, which can be passed to an optimiser. The estimate is stochastic because it is calculated from simulations that are conditioned on random observation noise (state feedback). To aid parameter optimisation, the variance of the estimate can be reduced by increasing the number of runs using the 'runs' property (see Performing multiple runs of a COIN model simulation). This is best done in conjunction with parallel computing to avoid excessive runtimes. An optimiser that can handle a stochastic objective function should be used (e.g., BADS).
Note that to fit the model to multiple paradigms simultaneously (e.g., spontaneous and evoked recovery), the negative log-likelihood should be calculated for each paradigm separately and their sum passed to the optimiser.
To see how to fit the COIN model to synthetic force-field adaptation data with sensory cues, see fit_model_to_data.m (note that this example requires installation of BADS).
The time it takes to execute multiple runs in series (e.g. using a for loop) can be prohibitively long if the number of runs is large, as the computational complexity of the model scales linearly with the number of runs. Therefore, in practice, it is advisable to use the smallest number of runs that achieves desired performance. In addition, runtime can be reduced by executing multiple runs in parallel (e.g., across different CPU cores on a computer cluster). To execute multiple runs in parallel, specify the maximum number of CPU cores available using the 'maxCores' property. The default setting of maxCores is 0, which executes multiple runs in series.
Not that the parallel computing functionality can only be used when simulating the COIN model or when calculating the negative log-likelihood, that is with the simulate_COIN and objective_COIN methods, respectively. It cannot be used when generating plots of variables, that is with the plot_COIN method.