-
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 over 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), 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 define a spontaneous recovery paradigm, call:
obj.perturbations = [zeros(1,50) ones(1,125) -ones(1,15) NaN(1,150)];
After specifying the sequence of perturbations, call the simulate_COIN method on object obj to run a simulation:
OUTPUT = obj.simulate_COIN;
The structure OUTPUT contains the simulation results. Note that the COIN model performs inference conditioned on a sequence of state feedback observations. These observations are generated by adding random observation noise to the user-defined perturbations. The state feedback and adaptation are contained in OUTPUT.runs{1} and can be plotted as follows:
figure
plot(OUTPUT.runs{1}.y,'.')
hold on
plot(OUTPUT.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 generated by adding random observation noise sampled from the prior to the user-defined perturbations . Use the 'runs' property to specify the number of runs to perform. For example, to perform 2 runs:
obj.runs = 2;
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. 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 Integrating over particles and runs). In general, using more runs will result in cleaner, less variable results.
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.
Note that the model runtime scales linearly with the number of runs, as each run is executed in series (by default). To execute multiple runs in parallel and thus reduce the runtime 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. 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.
Each particle in each run of the COIN model generates its own set of inferences. For analysis and plotting purposes, it is useful to integrate this set of inferences over both particles and runs to compute an expected value. Although all particles within a run have equal weights, all runs do not have equal weights in general. Therefore, integrating over particles and runs is performed by computing a simple average of each inferred variable 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 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.