Designing and Using Models

Introduction

The concept of a model is key to the use of QInfer. A model defines the probability distribution over experimental data given hypotheses about the system of interest, and given a description of the measurement performed. This distribution is called the likelihood function, and it encapsulates the definition of the model.

In QInfer, likelihood functions are represented as classes inheriting from either Model, when the likelihood function can be numerically evaluated, or Simulatable when only samples from the function can be efficiently generated.

Using Models and Simulations

Basic Functionality

Both Model and Simulatable offer basic functionality to describe how they are parameterized, what outcomes are possible, etc. For this example, we will use a premade model, SimplePrecessionModel. This model implements the likelihood function

\[\begin{split}\Pr(d | \omega; t) = \begin{cases} \cos^2 (\omega t / 2) & d = 0 \\ \sin^2 (\omega t / 2) & d = 1 \end{cases},\end{split}\]

as can be derived from Born’s Rule for a spin-½ particle prepared and measured in the \(\left|+\right\rangle\propto\left|0\right\rangle+\left|1\right\rangle\) state, and evolved under \(H = \omega \sigma_z / 2\) for some time \(t\).

In this way, we see that by defining the likelihood function in terms of the hypothetical outcome \(d\), the model parameter \(\omega\), and the experimental parameter \(t\), we can reason about the experimental data that we would extract from the system.

In order to use this likelihood function, we must instantiate the model that implements the likelihood. Since SimplePrecessionModel is provided with QInfer, we can simply import it and make an instance.

>>> from qinfer import SimplePrecessionModel
>>> m = SimplePrecessionModel()

Once a model or simulator has been created, you can query how many model parameters it admits and how many outcomes a given experiment can have.

>>> print(m.n_modelparams)
1
>>> print(m.modelparam_names)
['\\omega']
>>> print(m.is_n_outcomes_constant)
True
>>> print(m.n_outcomes(expparams=0))
2

Model and Experiment Parameters

The division between unknown parameters that we are trying to learn (\(\omega\) in the SimplePrecessionModel example) and the controls that we can use to design measurements (\(t\)) is generic, and is key to how QInfer handles the problem of parameter estimation. Roughly speaking, model parameters are real numbers that represent properties of the system that we would like to learn, whereas experiment parameters represent the choices we get to make in performing measurements.

Model parameters are represented by NumPy arrays of dtype float and that have two indices, one representing which model is being considered and one representing which parameter. That is, model parameters are defined by matrices such that the element \(X_{ij}\) is the \(j^{\text{th}}\) parameter of the model parameter vector \(\vec{x}_i\).

By contrast, since not all experiment parameters are best represented by the data type float, we cannot use an array of homogeneous dtype unless there is only one experimental parameter. The alternative is to use NumPy’s record array functionality to specify the heterogeneous type of the experiment parameters. To do so, instead of using a second index to refer to specific experiment parameters, we use fields. Each field then has its own dtype.

For instance, a dtype of [('t', 'float'), ('basis', 'int')] specifies that an array has two fields, named t and basis, having dtypes of float and int, respectively. Such arrays are initialized by passing lists of tuples, one for each field:

>>> eps = np.array([
...     (12.3, 2),
...     (14.1, 1)
... ], dtype=[('t', 'float'), ('basis', 'int')])
>>> print(eps)
[(12.3, 2) (14.1, 1)]
>>> eps.shape == (2,)
True

Once we have made a record array, we can then index by field names to get out each field as an array of that field’s value in each record, or we can index by record to get all fields.

>>> print(eps['t'])
[ 12.3  14.1]
>>> print(eps['basis'])
[2 1]
>>> print(eps[0])
(12.3, 2)

Model classes specify the dtypes of their experimental parameters with the property expparams_dtype. Thus, a common idiom is to pass this property to the dtype keyword of NumPy functions. For example, the model class BinomialModel adds an int field specifying how many times a two-outcome measurement is repeated, so to specify that we can use its expparams_dtype:

>>> from qinfer import BinomialModel
>>> bm = BinomialModel(m)
>>> print(bm.expparams_dtype)
[('x', 'float'), ('n_meas', 'uint')]
>>> eps = np.array([
...     (5.0, 10)
... ], dtype=bm.expparams_dtype)

Model Outcomes

Given a specific vector of model parameters \(\vec{x}\) and a specific experimental configuration \(\vec{c}\), the experiment will yield some outcome \(d\) according to the model distribution \(\Pr(d|\vec{x},\vec{c})\).

In many cases, such as SimplePrecessionModel discussed above, there will be a finite number of outcomes, which we can label by some finite set of integers. For example, we labeled the outcome \(\left|0\right\rangle\) by \(d=0\) and the outcome \(\left|1\right\rangle\) by \(d=1\). If this is the case for you, the rest of this section will likely not be very relevant, and you may assume your outcomes are zero-indexed integers ending at some value.

In other cases, there may be an infinite number of possible outcomes. For example, if the measurement returns the total number of photons measured in a time window, which can in principle be arbitrarily large, or if the measuremnt is of a voltage or current, which can be any real number. Or, we may have outcomes which require fancy data types. For instance, perhaps the output of a single experiment is a tuple of numbers rather than a single number.

To accomodate these possible situations, and to have a systematic way of testing whether or not all possibe outcomes can be enumerated, Simulatable (and subclasses like Model) has a method domain which for every given experimental parameter, returns a Domain object. One major benifit of explicitly storing these objects is that certain quantities (like bayes_risk) can be computed much more efficiently when all possible outcomes can be enumerated. Domain has attributes which specify whether or not it are finite, how many members it has and what they are, what data type they are, and so on.

For the BinomialModel defined above, there are n_meas+1 possible outcomes, with possible values the integers between 0 and n_meas inclusive.

>>> bdomain = bm.domain(eps)[0]
>>> bdomain.n_members
11
>>> bdomain.values
array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10])
>>> bdomain.dtype == np.int
True

We need to extract the \(0^\text{th}\) element of bm.domain(eps) above because eps is a vector of length \(1\) and domain always returns one domain for every member of eps. In the case where the domain is completely independent of eps, it should be possible to call m.domain(None) to return the unique domain of the model m.

The MultinomialModel requires a fancy datatype so that outcomes can be tuples of integers. In the following a single experiment of the model mm consists of throwing a four sided die n_meas times and recording how many times each side lands facing down.

>>> from qinfer import MultinomialModel, NDieModel
>>> mm = MultinomialModel(NDieModel(n=4))
>>> mm.expparams_dtype
[('exp_num', 'int'), ('n_meas', 'uint')]
>>> mmeps = np.array([(1, 3)], dtype=mm.expparams_dtype)
>>> mmdomain = mm.domain(mmeps)[0]
>>> mmdomain.dtype
dtype([('k', '<i8', (4,))])
>>> mmdomain.n_members
20
>>> print(mmdomain.values)
[([3, 0, 0, 0],) ([2, 1, 0, 0],) ([2, 0, 1, 0],) ([2, 0, 0, 1],)
 ([1, 2, 0, 0],) ([1, 1, 1, 0],) ([1, 1, 0, 1],) ([1, 0, 2, 0],)
 ([1, 0, 1, 1],) ([1, 0, 0, 2],) ([0, 3, 0, 0],) ([0, 2, 1, 0],)
 ([0, 2, 0, 1],) ([0, 1, 2, 0],) ([0, 1, 1, 1],) ([0, 1, 0, 2],)
 ([0, 0, 3, 0],) ([0, 0, 2, 1],) ([0, 0, 1, 2],) ([0, 0, 0, 3],)]

We see here all \(20\) possible ways to roll this die four times.

Note

Model inherits from Simulatable, and FiniteOutcomeModel inherits from Model. The subclass FiniteOutcomeModel is able to concretely define some methods (like simulate_experiment) because of the guarantee that all domains have a finite number of elements. Therefore, it is generally a bit less work to construct a FiniteOutcomeModel than it is to construct a Model.

Additionally, FiniteOutcomeModel automatically defines the domain corresponding to the experimental parameter ep by looking at n_outcomes, namely, if nep=n_outcomes(ep), then the corresponding domain has members 0,1,...,nep by default.

Finally, make note of the slightly subtle role of the method n_outcomes. In principle, n_outcomes is completely independent of domain. For FiniteOutcomeModel, it will almost always hold that m.n_outcomes(ep)==domain(ep)[0].n_members. For models with an infinite number of outcomes, n_members is not defined, but n_outcomes is defined and refers to “enough outcomes” (at the user’s discretion) to make estimates of quantities bayes_risk.

Simulation

Both models and simulators allow for simulated data to be drawn from the model distribution using the simulate_experiment() method. This method takes a matrix of model parameters and a vector of experiment parameter records or scalars (depending on the model or simulator), then returns an array of sample data, one sample for each combination of model and experiment parameters.

>>> modelparams = np.linspace(0, 1, 100)
>>> expparams = np.arange(1, 10) * np.pi / 2
>>> D = m.simulate_experiment(modelparams, expparams, repeat=3)
>>> print(isinstance(D, np.ndarray))
True
>>> D.shape == (3, 100, 9)
True

If exactly one datum is requested, simulate_experiment() will return a scalar:

>>> print(m.simulate_experiment(np.array([0.5]), np.array([3.5 * np.pi]), repeat=1).shape)
()

Note that in NumPy, a shape tuple of length zero indicates a scalar value, as such an array has no indices.

Note

For models with fancy outcome datatypes, it is demanded that the outcome data types, [d.dtype for d in m.domain(expparams)], be identical for every experimental parameter expparams being simulated. This can be checked with are_expparam_dtypes_consistent.

Likelihooods

The core functionality of Model, however, is the likelihood() method. This takes vectors of outcomes, model parameters and experiment parameters, then returns for each combination of the three the corresponding probability \(\Pr(d | \vec{x}; \vec{e})\).

>>> modelparams = np.linspace(0, 1, 100)
>>> expparams = np.arange(1, 10) * np.pi / 2
>>> outcomes = np.array([0], dtype=int)
>>> L = m.likelihood(outcomes, modelparams, expparams)

The return value of likelihood() is a three-index array of probabilities whose shape is given by the lengths of outcomes, modelparams and expparams. In particular, likelihood() returns a rank-three tensor \(L_{ijk} := \Pr(d_i | \vec{x}_j; \vec{e}_k)\).

>>> print(isinstance(L, np.ndarray))
True
>>> L.shape == (1, 100, 9)
True

Implementing Custom Simulators and Models

In order to implement a custom simulator or model, one must specify metadata describing the number of outcomes, model parameters, experimental parameters, etc. in addition to implementing the simulation and/or likelihood methods.

Here, we demonstrate how to do so by walking through a simple subclass of FiniteOutcomeModel. For more detail, please see the API Reference.

Suppose we wish to implement the likelihood function

\[\Pr(0 | \omega_1, \omega_2; t_1, t_2) = \cos^2(\omega_1 t_1 / 2) \cos^2(\omega_2 t_2 / 2),\]

as may arise in looking, for instance, at an experiment inspired by 2D NMR. This model has two model parameters, \(\omega_1\) and \(\omega_2\), and so we start by creating a new class and declaring the number of model parameters as a property:

from qinfer import FiniteOutcomeModel
import numpy as np

class MultiCosModel(FiniteOutcomeModel):
    
    @property
    def n_modelparams(self):
        return 2

Next, we proceed to add a property and method indicating that this model always admits two outcomes, irrespective of what measurement is performed. This will also automatically define the domain method.

    @property
    def is_n_outcomes_constant(self):
        return True
    def n_outcomes(self, expparams):
        return 2

We indicate the valid range for model parameters by returning an array of dtype bool for each of an input matrix of model parameters, specifying whether each model vector is valid or not (this is important in resampling, for instance, to make sure particles don’t move to bad locations). Typically, this will look like some typical bounds checking, combined using logical_and and all. Here, we follow that model by insisting that all elements of each model parameter vector must be at least 0, and must not exceed 1.

    def are_models_valid(self, modelparams):
        return np.all(np.logical_and(modelparams > 0, modelparams <= 1), axis=1)

Next, we specify what a measurement looks like by defining expparams_dtype. In this case, we want one field that is an array of two float elements:

    @property
    def expparams_dtype(self):
        return [('ts', 'float', 2)]

Finally, we write the likelihood itself. Since this is a two-outcome model, we can calculate the rank-two tensor \(p_{jk} = \Pr(0 | \vec{x}_j; \vec{e}_k)\) and let pr0_to_likelihood_array() add an index over outcomes for us so \(L_{0jk}=p_{jk}\) and \(L_{1jk}=1-p_{jk}\). To compute \(p_{jk}\) efficiently, it is helpful to do a bit of index gymnastics using NumPy’s powerful broadcasting rules. In this example, we set up the calculation to produce terms of the form \(\cos^2(x_{j,l} e_{k,l} / 2)\) for \(l \in \{0, 1\}\) indicating whether we’re referring to \(\omega_1\) or \(\omega_2\), respectively. Multiplying along this axis then gives us the product of the two cosine functions, and in a way that very nicely generalizes to likelihood functions of the form

\[\Pr(0 | \omega_1, \omega_2; t_1, t_2) = \prod_l \cos^2(\omega_l t_l / 2).\]

Running through the index gymnastics, we can implement the likelihood function as:

    def likelihood(self, outcomes, modelparams, expparams):
        # We first call the superclass method, which basically
        # just makes sure that call count diagnostics are properly
        # logged.
        super(MultiCosModel, self).likelihood(outcomes, modelparams, expparams)
        
        # Next, since we have a two-outcome model, everything is defined by
        # Pr(0 | modelparams; expparams), so we find the probability of 0
        # for each model and each experiment.
        #
        # We do so by taking a product along the modelparam index (len 2,
        # indicating omega_1 or omega_2), then squaring the result.
        pr0 = np.prod(
            np.cos(
                # shape (n_models, 1, 2)
                modelparams[:, np.newaxis, :] *
                # shape (n_experiments, 2)
                expparams['ts']
            ), # <- broadcasts to shape (n_models, n_experiments, 2).
            axis=2 # <- product over the final index (len 2)
        ) ** 2 # square each element
        
        # Now we use pr0_to_likelihood_array to turn this two index array
        # above into the form expected by SMCUpdater and other consumers
        # of likelihood().
        return FiniteOutcomeModel.pr0_to_likelihood_array(outcomes, pr0)

Our new custom model is now ready to use! To simulate data from this model, we set up modelparams and expparams as before, taking care to conform to the expparams_dtype of our model:

>>> mcm = MultiCosModel()
>>> modelparams = np.dstack(np.mgrid[0:1:100j,0:1:100j]).reshape(-1, 2)
>>> expparams = np.empty((81,), dtype=mcm.expparams_dtype)
>>> expparams['ts'] = np.dstack(np.mgrid[1:10,1:10] * np.pi / 2).reshape(-1, 2)
>>> D = mcm.simulate_experiment(modelparams, expparams, repeat=2)
>>> print(isinstance(D, np.ndarray))
True
>>> D.shape == (2, 10000, 81)
True

Finally, we mention a useful tool for doing a set of tests on the custom model, which make sure its pieces are working as expected. These tests look at things like data types and index dimensions of various functions. They also plug the outputs of some methods into the inputs of other methods, and so forth. Although they can’t check the statistical soundness of your model, if they all pass, you can be pretty confident you won’t run into weird indexing bugs in the future.

We just need to pass test_model() an instance of the custom model, a prior that samples valid model parameters, and an array of valid expparams.

>>> from qinfer.tests import test_model
>>> from qinfer import UniformDistribution
>>> prior = UniformDistribution([[0,1],[0,1]])
>>> test_model(mcm, prior, expparams)
.......
----------------------------------------------------------------------
Ran 7 tests in 0.013s

OK

Note

Creating expparams as an empty array and filling it by field name is a straightforward way to make sure it matches expparams_dtype, but it comes with the risk of forgetting to initialize a field, so take care when using this method.

Adding Functionality to Models with Other Models

QInfer also provides model classes which add functionality or otherwise modify other models. For instance, the BinomialModel class accepts instances of two-outcome models and then represents the likelihood for many repeated measurements of that model. This is especially useful in cases where experimental concerns make switching experiments costly, such that repeated measurements make sense.

To use BinomialModel, simply provide an instance of another model class:

>>> from qinfer import SimplePrecessionModel
>>> from qinfer import BinomialModel
>>> bin_model = BinomialModel(SimplePrecessionModel())

Experiments for BinomialModel have an additional field from the underlying models, called n_meas. If the original model used scalar experiment parameters (e.g.: expparams_dtype is float), then the original scalar will be referred to by a field x.

>>> eps = np.array([(12.1, 10)], dtype=bin_model.expparams_dtype)
>>> print(eps['x'], eps['n_meas'])
[ 12.1] [10]

Another model which decorates other models in this way is PoisonedModel, which is discussed in more detail in Performance and Robustness Testing. Roughly, this model causes the likeihood functions calculated by its underlying model to be subject to random noise, so that the robustness of an inference algorithm against such noise can be tested.