Building Your Own model¶

Foreword¶

This section is a guide for writing a Dynamic Causal Model for your own application. Be warned that the code is currently under development, and that these instructions are likely to change significantly as time goes on.

Overview¶

Creating your model¶

The dcEmb Dynamic Causal Model library provides support for creating arbitary DCM models through inheritance. The core idea is that the only things you should need to write for your own application are:

Your forward (generative) model, and
Estimates of parameter/hyperparameter distributions.

To write your own DCM you simply need to create a new class that inherits from the dynamic_model class. The Dynamic Model class contains two virtual methods that you will be forced to overwrite. This are:

virtual Eigen::VectorXd get_observed_outcomes();

A function that returns a vector of all observations you are training the model based on, subject to any suitable transformations.

virtual std::function<Eigen::VectorXd (Eigen::VectorXd)> get_forward_model_function();

A function that return a wrapped link to your generative function. The wrapped function must take a vector of parameters, and return a vector of predicted observations.

Once you have overwritten these methods in your class, the only thing you need to do is create an instance of your class, populate the parameters specifying parameter/hyperparameter distributions (plus a few others) and call new_model.invert(). The current parameters that are mandatory are:

1. Eigen::VectorXd prior_parameter_expectations;, a vector of the prior parameter expectations

2. Eigen::MatrixXd prior_parameter_covariances;, the covariance Matrix containing the prior parameter covariances

3. Eigen::MatrixXd prior_hyper_covariances;, a vector of the prior “hyperparameter” expectations

4. Eigen::MatrixXd prior_hyper_covariances;, the covariance Matrix containing the prior “hyperparameter” covariances

5. Eigen::MatrixXd response_vars;, the matrix of observations over (potentially serveral variables) to train the model on

6. int num_samples;, the number of timesteps in the observation matrix

7. int num_response_vars;, the number of different variables in the observation matrix

8. Eigen::VectorXi select_response_vars;, a vector enumerating which columns of the outcome matrix generated by our generative model correspond to the columns of the observation matrix we are training on.

Advanced functions¶

There are several functions, for example Bayesian Model Reduction (BMR), or Parametric Empirical Bayes (PEB) which are called on dynamic causal models themselves. These functions have been made generic and should be able to be called over any DCM-like class through the use of templating. For example, to create a PEB class in the COVID example, we do:

peb_model<dynamic_COVID_model> PEB;

Some of these classes may require additional parameters to be specified before they can be run. PEB, for example will require a matrix specifying random effects, as well as a between-subject design matrix.

A worked example for these functions is still a work in progress. We suggest anyone interested in pursuing this further references the documentation for these functions.

Worked Example: dcm_3body¶

In this section, we work through the steps we’ve just described in an example creating our three body problem DCM example.

The goal of this example is to use the DCM method to recover a known set of results from a physical system. The system in question is the three body problem, the problem of finding the trajectories of three planetary bodies in space. There exist several known stable orbits for the three body problem, which are periodic. In this example, we will take one of these stable solutions (a “figure of 8” orbit), and the positions, masses and velocities of one of the planets in this orbit, and use the DCM method to (hopefully accurately) estimate the positions, masses and velocities of the other two.

Creating the base model¶

In this section, we walk through creating a new DCM, and placing an appropriate three body problem forward model inside it.

Chronologically, the first thing we need to do following the instructions above, is to create a new class that inherits from the dynamic_model class. To do this we create a new folder in the include subdirectory of the project, called “3body”, and inside this create a new file called “dynamic_3body_model.hh”, with the following content:

include/3body/dynamic_3body_model.hh

#include "dynamic_model.hh"
#pragma once

class dynamic_3body_model : public dynamic_model {
  dynamic_3body_model(); // Class Constructor
};

Continuing to follow the above instructions, we also need to populate this class with at least the following functions: get_observed_outcomes() and get_forward_model_function().

include/3body/dynamic_3body_model.hh

#include "dynamic_model.hh"
#pragma once

class dynamic_3body_model : public dynamic_model {
 public:
  Eigen::VectorXd get_observed_outcomes();
  std::function<Eigen::VectorXd(Eigen::VectorXd)>get_forward_model_function();
  dynamic_3body_model();  // Class Constructor
};

The above code realizes the class definition for dynamic_3body_model. We now need to go away and create the implementations of these functions we have just defined.

In this case, our forward model will take the initial state of our three planets (positions, masses, and velocities) and produce a time series of the state over a given period. We can do this by using newtons law of gravitation to create equations of motion dictating the planets movements over time (see here for a simple walkthrough). These equations of motion lack an analytical solution, so we solve them using numerical methods. We suggest using the runge-kutta method, of which an implementation is provided in dCEmbs utility module.

Putting this together speaks to the creation of two functions, an “equation of motion” function evaluating the rate of change of the state, given the current state, and a main “generative model” function that iteratively applies the runge-kutta method to this to produce a time series. In the current version of the code, we also implement a forward model() function that explicitly wraps the entire forward model. We should implement these functions as member functions of our dynamic_3body_model class. To do this, we will first need to change our dynamic_3body_model.hh header file to include definitions for three new functions eval_generative(), forward model(), and differential_eq():

include/3body/dynamic_3body_model.hh

class dynamic_3body_model : public dynamic_model {
public:
  parameter_location_3body parameter_locations;
  int num_bodies;
  Eigen::VectorXd get_observed_outcomes();
  std::function<Eigen::VectorXd(Eigen::VectorXd)> get_forward_model_function();
  Eigen::VectorXd forward_model(
      const Eigen::VectorXd& parameters,
      const int& timeseries_length,
      const Eigen::VectorXi& select_response_vars);
  Eigen::MatrixXd eval_generative(
      const Eigen::VectorXd& parameters,
      const int& timeseries_length);
  Eigen::MatrixXd eval_generative(
      const Eigen::VectorXd& parameters,
      const int& timeseries_length,
      const Eigen::VectorXi& select_response_vars);
  Eigen::VectorXd differential_eq(const Eigen::VectorXd& state,
                                  const int& num_bodies);
  dynamic_3body_model();
};

As well as the parameters() and timeseries_length() parameters, which we use in our function definitions for the parameters our generative model will use, and the length of the timeseries it produces, we also pass it select_response_vars(). This is a variable that we will use later in order to match up our generative models outputs to our training data.

We also need to create implementations for each of these functions. To do this, we need to create a new folder in the src subdirectory of our project, called 3body, and create a new file in it called dynamic_3body_model.cc. This file should have each of the following blocks of content in it: