-
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 inferred variables
- Integrating out particles and runs
- Plotting inferred variables
- Fitting the 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). These default values 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 simulate learning on a spontaneous recovery paradigm, define a series of perturbations (use NaN to indicate a channel trial):
obj.perturbations = [zeros(1,50) ones(1,125) -ones(1,15) NaN(1,150)];
Then run the model by calling the simulate_COIN method on object obj:
S = obj.simulate_COIN;
The adaptation and state feedback (perturbation plus random observation noise) are contained in S.runs{1}. They can be plotted as follows:
figure
plot(S.runs{1}.y,'.')
hold on
plot(S.runs{1}.yHat)
legend('state feedback','adaptation')
The actual observation noise and hence state feedback that a participant perceives is hidden from the perspective of the experimenter. Therefore, rather than perform one run of a simulation, conditioned on one particular state feedback sequence, multiple runs of a simulation can be performed, each conditioned on a different state feedback sequence (sampled from the prior). The state feedback can then be integrated out by computing a weighted average of the runs, with each run assigned a weight based on the likelihood of the adaptation data of the participant (if available). Formally, this corresponds to approximating the expected value of a function of the state feedback (e.g. the inferences made by the COIN model) using importance sampling. Use the 'runs' property to specify the number of runs to perform. For example, to perform 2 runs:
obj.runs = 2;
S = obj.simulate_COIN;
S.runs is a cell array with one cell per run and S.weights is a vector specifying the relative weight of each run. In this example, each run is assigned an equal weight as no adaptation data was passed to the object. However, if adaptation data is passed to the object before running the model, 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, as weights are reset whenever runs are resampled during sequential importance sampling). To generate a set of weighted runs, set the model parameters to their maximum likelihood estimates (the model should be fit to the adaptation data first, see Fitting the model to data), pass the adaptation data to the class object via the adaptation property and call the simulate_COIN method. The resultant weights should be used to compute a weighted average of the inferred variables across runs (as in Analysing stored variables). In general, using more runs will result in cleaner, less variable results.
Note that the model runtime scales linearly with the number of runs, as each run is executed in series (by default). To reduce the runtime by executing multiple runs in parallel see Parallel computing.
Inference in the COIN model is performed in an online (i.e. recursive) manner, which does not require all of the past values of all inferred variables to be stored in memory. Therefore, to reduce memory requirements, the past values of variables are only stored in memory (e.g. for later analysis) if an explicit request is made; the state feedback (model input) and adaptation (model output), which are always stored, are the only exceptions to this. To store the values of certain variables on all trials, add the names of these variables to the store property; for a full list of the names of variables that can be stored, call help COIN (see 'VARIABLES' section). For example, to store the Kalman gains ('k') and responsibilities ('cFilt'), set store to
obj.store = {'k','cFilt'};
S = obj.simulate_COIN;
The store property must be set before calling the simulate_COIN method (as shown above).
The stored variables are contained in the cell array S and can be analysed after the model has run. For example, to compute the Kalman gain of the context with the highest responsibility:
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(S.runs{run}.cFilt(:,particle,trial));
Kalman_gain(trial,run,particle) = S.runs{run}.k(j,particle,trial);
end
end
end
Inference in the COIN model is performed using a sequential Monte Carlo method known as particle learning that uses an ensemble of equally-weighted particles. Each run is associated with an ensemble of particle of particles that performs inference. Within each run, a simple average across particles can be computed (as all particles within the same run have the same weight), whereas across runs, a weighted average should be computed:
figure
plot(sum(S.weights.*mean(Kalman_gain,3),2)) % weighted average over runs (dim 2) and simple average over particles (dim 3)
ylabel('Kalman gain | c^*')
Some variables or their distributions can be plotted by activating the relevant flags in the properties of the class object. For a full list of variables or distributions 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. To plot a continuous distribution, also specify the points at which to evaluate the distribution. For example, to plot the distribution of the state of each context and the predicted probabilities:
obj.plot_state_given_context = true;
obj.state_values = linspace(-1.5,1.5,500); % values to evaluate state | context at
obj.plot_predicted_probabilities = true;
S = obj.simulate_COIN;
These properties should be set before calling the simulate_COIN method (as shown above), so that the variables needed to generate the requested plots can be stored in memory (this is done automatically).
After the model has run, call the plot_COIN method:
P = obj.plot_COIN(S);
This will generate the requested plots—a state | context plot and a predicted probabilities plot in this example. The structure P contains the data that is plotted (to see how the data in P 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 variables or distributions are averaged across 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.