QInfer - Library for Statistical Inference in Quantum Information¶

Related Examples¶

• simple_precession_example
• randomized_benchmarking

Introduction¶

Probability distributions such as prior distributions over model parameters are reprented in QInfer by objects of type Distribution that are responsible for producing samples according to those distributions. This is especially useful, for instance, when drawing initial particles for use with an SMCUpdater.

The approach to representing distributions taken by QInfer is somewhat different to that taken by, for example, scipy.stats, in that a QInfer Distribution is a class that produces samples according to that distribution. This means that QInfer Distribution objects provide much less information than do those represented by objects in scipy.stats, but that they are much easier to write and combine.

Sampling Pre-made Distributions¶

QInfer comes along with several distributions, listed in Specific Distributions. Each of these is a subclass of Distribution, and hence has a method sample() that produces an array of samples.

>>> from qinfer import NormalDistribution
>>> dist = NormalDistribution(0, 1)
>>> samples = dist.sample(n=5)
>>> samples.shape == (5, 1)
True

Combining Distributions¶

Distribution objects can be combined using other distribution objects. For instance, if \(a \sim \mathcal{N}(0, 1)\) and \(b \sim \text{Uni}(0, 1)\), then the product distribution on \((a,b)\) can be produced by using ProductDistribution:

>>> from qinfer import UniformDistribution, ProductDistribution
>>> a = NormalDistribution(0, 1)
>>> b = UniformDistribution([0, 1])
>>> ab = ProductDistribution(a, b)
>>> samples = ab.sample(n=5)
>>> samples.shape == (5, 2)
True

Making Custom Distributions¶

To make a custom distribution, one need only implement sample() and set the property n_rvs to indicate how many random variables the new distribution class represents.

For example, to implement a distribution over \(x\) and \(y\) such that \(\sqrt{x^2 + y^2} \sim \mathcal{N}(1, 0.1)\) and such that the angle between \(x\) and \(y\) is drawn from \(\text{Uni}(0, 2 \pi)\):

from qinfer import Distribution

class RingDistribution(Distribution):
    @property
    def n_rvs(self):
        return 2

    def sample(self, n=1):
        r = np.random.randn(n, 1) * 0.1 + 1
        th = np.random.random((n, 1)) * 2 * np.pi

        x = r * np.cos(th)
        y = r * np.sin(th)

        return np.concatenate([x, y], axis=1)

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¶

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

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

>>> 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

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

>>> 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)

>>> 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

>>> 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],)]

>>> 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

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

>>> 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)

>>> 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

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

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

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.

Introduction¶

Arguably the core of QInfer, the qinfer.smc module implements the sequential Monte Carlo algorithm in a flexible and robust manner. At its most basic, using QInfer’s SMC implementation consists of specifying a model, a prior, and a number of SMC particles to use.

The main component of QInfer’s SMC support is the SMCUpdater class, which performs Bayesian updates on a given prior in response to new data. In doing so, SMCUpdater will also ensure that the posterior particles are properly resampled. For more details on the SMC algorithm as implemented by QInfer, please see [GFWC12].

Using SMCUpdater¶

>>> from qinfer import SMCUpdater, UniformDistribution, SimplePrecessionModel
>>> model = SimplePrecessionModel()
>>> prior = UniformDistribution([0, 1])
>>> updater = SMCUpdater(model, 1000, prior)

>>> true_model = prior.sample()
>>> experiment = np.array([12.1], dtype=model.expparams_dtype)
>>> outcome = model.simulate_experiment(true_model, experiment)
>>> updater.update(outcome, experiment)

>>> posterior_samples = updater.sample(n=100)
>>> posterior_samples.shape == (100, 1)
True

>>> est = updater.est_mean()
>>> print(est)
[ 0.53147953]

Plotting Posterior Distributions¶

The SMCUpdater also provides tools for producing plots to describe the updated posterior. For instance, the plot_posterior_marginal() method uses kernel density estimation to plot the marginal over all but a single parameter over the posterior.

prior = UniformDistribution([0, 1])
model = SimplePrecessionModel()
updater = SMCUpdater(model, 2000, prior)

# Plot according to the initial prior.
updater.plot_posterior_marginal()

# Simulate 50 different measurements and use
# them to update.
true = prior.sample()
heuristic = ExpSparseHeuristic(updater)

for idx_exp in range(25):
    expparams = heuristic()
    datum = model.simulate_experiment(true, expparams)
    updater.update(datum, expparams)

# Plot the posterior.
updater.plot_posterior_marginal()

# Add a legend and show the final plot.
plt.legend(['Prior', 'Posterior'])
plt.show()

(Source code, svg, pdf, hires.png, png)

For multi-parameter models, the plot_covariance() method plots the covariance matrix for the current posterior as a Hinton diagram. That is, positive elements are shown as white squares, while negative elements are shown as black squares. The relative sizes of each square indicate the magnitude, making it easy to quickly identify correlations that impact estimator performance. In the example below, we use the Simple Estimation Functions to quickly analyze Randomized Benchmarking data and show the resulting correlation between the \(p\), \(A\) and \(B\) parameters. For more detail, please see the randomized benchmarking example.

p = 0.995
A = 0.5
B = 0.5

ms = np.linspace(1, 800, 201).astype(int)
signal = A * p ** ms + B

n_shots = 25
counts = np.random.binomial(p=signal, n=n_shots)

data = np.column_stack([counts, ms, n_shots * np.ones_like(counts)])
mean, cov, extra = simple_est_rb(data, return_all=True, n_particles=12000, p_min=0.8)
extra['updater'].plot_covariance()

plt.show()

(Source code, svg, pdf, hires.png, png)

Advanced Usage¶

>>> from qinfer import LiuWestResampler
>>> updater = SMCUpdater(model, 1000, prior, resampler=LiuWestResampler(0.9))

Using Heuristics in Updater Loops¶

During an experiment, the current posterior distribution represented by an SMCUpdater instance can be used to adaptively make decisions about which measurements should be performed. For example, utility functions such as the information gain and negative variance can be used to choose optimal measurements.

On the other hand, this optimization is generally quite computationally expensive, such that many less optimal measurements could have been performed before an optimization step completes. Philosophically, there are no “bad” experiments, in that even suboptimal measurements can still carry useful information about the parameters of interest.

In light of this, it can be substantially less expensive to use a heuristic function of prior information to select experiments without explicit optimization. For example, consider the single-parameter inversion precession model with model parameters \(\vec{x} = (\omega)\), experiment parameters \(e = (\omega_-, t)\) and likelihood function

For a given posterior distribution, the particle guess heuristic (PGH) [WGFC13a] then chooses \(\omega_-\) and \(t\) for the next experiment by first sampling two particles \(\omega_-\) and \(\omega_-'\) from the posterior. The PGH then assigns the time \(t = 1 / |\omega_- \omega_-'|\).

QInfer implements heuristics as subtypes of Heuristic, each of which take an updater and can be called to produce experimental parameters. For example, the PGH is implemented by the class PGH, and can be used in an updater loop to adaptively select experiments.

model = SimpleInversionModel()
prior = UniformDistribution([0, 1])
updater = SMCUpdater(model, 1000, prior)
heuristic = PGH(updater, inv_field='w_', t_field='t')

true_omega = prior.sample()

ts = []
est_omegas = []

for idx_exp in range(100):
    experiment = heuristic()
    datum = model.simulate_experiment(true_omega, experiment)
    updater.update(datum, experiment)

    ts.append(experiment['t'])
    est_omegas.append(updater.est_mean())

ax = plt.subplot(2, 1, 1)
plt.semilogy((est_omegas - true_omega) ** 2)
plt.ylabel('Squared Error')

plt.subplot(2, 1, 2, sharex=ax)
plt.semilogy(ts)
plt.xlabel('# of Measurements')
plt.ylabel('$t$')

plt.show()

(Source code, svg, pdf, hires.png, png)

Changing Heuristic Parameters¶

Essentially, heuristics in QInfer are functions that take SMCUpdater instances and return functions that then yield experiments. This design allows for specializing heuristics by providing other arguments along with the updater. For instance, the ExpSparseHeuristic class implements exponentially-sparse sampling \(t_k = ab^k\) for SimplePrecessionModel. Both \(a\) and \(b\) are parameters of the heuristic, named scale and base, respectively. Thus, it is easy to override the defaults to obtain different heuristics.

model = SimplePrecessionModel()
prior = UniformDistribution([0, 1])
updater = SMCUpdater(model, 1000, prior)
heuristic = ExpSparseHeuristic(updater, scale=0.5)

true_omega = prior.sample()
est_omegas = []

for idx_exp in range(100):
    experiment = heuristic()
    datum = model.simulate_experiment(true_omega, experiment)
    updater.update(datum, experiment)

    est_omegas.append(updater.est_mean())

plt.semilogy((est_omegas - true_omega) ** 2)
plt.xlabel('# of Measurements')
plt.ylabel('Squared Error')

plt.show()

(Source code, svg, pdf, hires.png, png)

In overriding the default parameters of heuristics, the functools.partial() function provided with the Python standard library is especially useful, as it allows for easily making new heuristics with different default parameter values:

>>> from qinfer import ExpSparseHeuristic
>>> from functools import partial
>>> rescaled_heuristic_class = partial(ExpSparseHeuristic, scale=0.01)

Later, once we have an updater, we can then make new instances of our rescaled heuristic.

>>> heuristic = rescaled_heuristic_class(updater) 

This technique is especially useful in Performance Testing, as it makes it easy to modify existing heuristics by changing default parameters through partial application.

Introduction¶

Randomized benchmarking allows for extracting information about the fidelity of a quantum operation by exploiting twirling errors over an approximate implementation of the Clifford group [KL+08]. This provides the advantage that the fidelity can be learned without simulating the dynamics of a quantum system. Instead, benchmarking admits an analytic form for the survival probability for an arbitrary input state in terms of the strength \(p\) of an equivalent depolarizing channel.

QInfer supports randomized benchmarking by implementing this survival probability as a likelihood function. This allows for randomized benchmarking to be used together with Sequential Monte Carlo, such that prior information can be incorporated and robustness to finite sampling can be obtained [GFC14].

Regardless of the order or interleaving mode, each model instance for randomized benchmarking yields 0/1 data, with 1 indicating a survival (measuring the same state after applying a gate sequence as was initially prepared). To use these models with data batched over many sequences, model instances can be augmented by BinomialModel.

Zeroth-Order Model¶

The RandomizedBenchmarkingModel class implements randomized benchmarking as a QInfer model, both in interleaved and non-interleaved modes. For the non-interleaved mode, there are three model parameters, \(\vec{x} = (p, A_0, B_0)\), given by [MGE12] as

where \(E_\psi\) is the measurement, \(\rho_\psi\) is the state preparation, \(\Lambda\) is the average map over timesteps and gateset elements, \(d\) is the dimension of the system, and where \(F_\ave\) is the average gate fidelity, taken over the gateset. The functions p and F convert back and forth between depolarizing parameters and fidelities.

An experiment parameter vector for this model is simply a specification of \(m\), the length of the Clifford sequence used for that datum. Since RandomizedBenchmarkingModel represents 0/1 data, it is common to wrap this model in a BinomialModel:

>>> from qinfer import BinomialModel
>>> from qinfer import RandomizedBenchmarkingModel
>>> model = BinomialModel(RandomizedBenchmarkingModel(order=0, interleaved=False))
>>> expparams = np.array([
...    (100, 1000) # 1000 shots of sequences with length 100.
... ], dtype=model.expparams_dtype)

>>> from qinfer.rb import RandomizedBenchmarkingModel
>>> model = RandomizedBenchmarkingModel(interleaved=True)

Introduction¶

Tomography is the most common quantum statistical problem being the subject of both theoretical and practical studies. The qinfer.tomography module has rich support for many of the common models of tomography including standard distributions and heuristics, and also provides convenient plotting tools.

If the quantum state of a system is \(\rho\) and a measurement is obtained, then the probability of obtaining the outcome associated with effect \(E\) is \(\operatorname{Pr}(E|\rho) = \operatorname{Tr}(\rho E)\), the Born rule. The tomography problem is the inverse problem and is succinctly stated as follows. Given an unknown state \(\rho\), and a set of count statistics \(\{n_k\}\) from measurements, the corresponding operators of which are \(\{E_k\}\), determine \(\rho\).

In the context of Bayesian statistics, priors are distributions of quantum states and models define the Hilbert space dimension and the Born rule as a likelihood function (and perhaps additional complications associated with drift parameters and measurement errors). The tomography module implements these details and has built-in support for common specifications.

The tomography module was developed concurrently with new results on Bayesian priors in [GCC16]. Please see the paper for more detailed discussions of these more advanced topics.

Bases¶

Bases define the map from abstract quantum theory to concrete representations. Once a basis is chosen, the tomography problem becomes a special case of a generic parameter estimation problem.

The tomography module is used in the same way as other but requires the specification of a basis via the TomographyBasis class. QInfer comes with the following Built-in bases: the Gell-Mann basis, the Pauli basis and function to combine bases with the tensor product. Thus, the first step in using the tomography module is to define a basis. For example, here we define the 1-qubit Pauli basis:

>>> from qinfer.tomography import pauli_basis
>>> basis = pauli_basis(1)
>>> print(basis) 
<TomographyBasis pauli_basis dims=[2] at ...>

User defined bases must be orthogonal Hermitian operators and have a 0’th component of \(I/\sqrt{d}\), where \(d\) is the dimension of the quantum system and \(I\) is the identity operator. This implies the remaining operators are traceless.

Built-in Distributions¶

QInfer comes with several built-in distributions listed in Specific Distributions. Each of these is a subclass of DensityOperatorDistribution. Distributions of quantum channels can also be subclassed as such with appeal to the Choi-Jamiolkowski isomorphism.

For density matrices, the GinibreDistribution defines a prior over mixed quantum which allows for support also for states of fixed rank [OSZ10]. For example, we can draw a sample from the this prior as follows:

>>> from qinfer.tomography import GinibreDistribution
>>> prior = GinibreDistribution(basis)
>>> print(prior.sample()) 
[[ 0.70710678 -0.17816233  0.45195168 -0.08341437]]

Recall this is this representation of of a qubit in the Pauli basis defined above. Quantum states are in general high dimensional objects which makes visualizing distributions of them challenging. The only 2-dimensional example is that of a rebit, which is usually defined as a qubit in the Pauli (or Bloch) representation with one of the Pauli expectations constrained to zero (usually \(\operatorname{Tr}(\rho \sigma_y)=0\)).

Here we create a distribution of rebits accord to the Ginibre ensemble and use plot_rebit_prior() to depict this distribution through (by default) 2000 random samples. While discussing models below, we will see how to depict the particles of an SMCUpdater directly.

basis = tomography.bases.pauli_basis(1)
prior = tomography.distributions.GinibreDistribution(basis)
tomography.plot_rebit_prior(prior, rebit_axes=[1, 3])
plt.show()

(Source code, svg, pdf, hires.png, png)

Using TomographyModel¶

The core of the tomography module is the TomographyModel. The key assumption in the current version is that of two-outcome measurements. This has the convenience of allowing experiments to be specified by a single vectorized positive operator:

>>> from qinfer.tomography import TomographyModel
>>> model = TomographyModel(basis)
>>> print(model.expparams_dtype) 
[('meas', <type 'float'>, 4)]

Suppose we measure \(\sigma_z\) on a random state. The measurement effects are \(\frac12 (I\pm \sigma_z)\). Since they sum to identity, we need only specify one of them. We can use TomographyModel to calculate the Born rule probability of obtaining one of these outcomes as follows:

>>> expparams = np.zeros((1,), dtype=model.expparams_dtype)
>>> expparams['meas'][0, :] = basis.state_to_modelparams(np.sqrt(2)*(basis[0]+basis[1])/2)
>>> print(model.likelihood(0, prior.sample(), expparams)) 
[[[ 0.62219803]]]

Built-in Heuristics¶

In addition to analyzing given data sets, the tomography module is well suited for testing measurement strategies against standard heuristics. These built-in heuristics are listed at Abstract Heuristics. For qubits, the most commonly used heuristic is the random sampling of Pauli basis measurements, which is implemented by RandomPauliHeuristic.

>>> from qinfer.tomography import RandomPauliHeuristic
>>> from qinfer import SMCUpdater
>>> updater = SMCUpdater(model, 100, prior)
>>> heuristic = RandomPauliHeuristic(updater)
>>> print(model.simulate_experiment(prior.sample(), heuristic())) 
0

Time-Dependent Parameters¶

In addition to learning static parameters, QInfer can be used to learn the values of parameters that change stochastically as a function of time. In this case, the model parameter vector \(\vec{x}\) is interpreted as a time-dependent vector \(\vec{x}(t)\) representing the state of an underlying process. The resulting statistical problem is often referred to as state-space estimation. By using an appropriate resampling algorithm, such as the Liu-West algorithm [LW01], state-space and static parameter estimation can be combined such that a subset of the components of \(\vec{x}\) are allowed to vary with time.

QInfer represents state-space filtering by the use of the Simulatable.update_timestep() method, which samples how a model parameter vector is updated as a function of time. In particular, this method is used by SMCUpdater to draw samples from the distribution \(f\)

where \(t_{\ell}\) is the time at which the experiment \(e_{\ell}\) is measured, and where \(t_{\ell+1}\) is the time step immediately following \(t_{\ell}\). As this distribution is in general dependent on the experiment being performed, update_timestep() is vectorized in a manner similar to likelihood() (see Designing and Using Models for details). That is, given a tensor \(X_{i,j}\) of model parameter vectors and a vector \(e_k\) of experiments, update_timestep() returns a tensor \(X_{i,j,k}'\) of sampled model parameters at the next time step.

Random Walk Models¶

As an example, RandomWalkModel implements update_timestep() by taking as an input a Distribution specifying steps \(\Delta \vec{x} = \vec{x}(t + \delta t) - \vec{x}(t)\). An instance of RandomWalkModel decorates another model in a similar fashion to other derived models. For instance, the following code declares a precession model in which the unknown frequency \(\omega\) changes by a normal random variable with mean 0 and standard deviation 0.005 after each measurement.

>>> from qinfer import (
...     SimplePrecessionModel, RandomWalkModel, NormalDistribution
... )
>>> model = RandomWalkModel(
...     underlying_model=SimplePrecessionModel(),
...     step_distribution=NormalDistribution(0, 0.005 ** 2)
... )

We can then draw simulated trajectories for the true and estimated value of \(\omega\) using a minor modification to the updater loop discussed in Sequential Monte Carlo.

model = RandomWalkModel(
    # Note that we set a minimum frequency of negative
    # infinity to prevent problems if the random walk
    # causes omega to cross zero.
    underlying_model=SimplePrecessionModel(min_freq=-np.inf),
    step_distribution=NormalDistribution(0, 0.005 ** 2)
)
prior = UniformDistribution([0, 1])
updater = SMCUpdater(model, 2000, prior)

expparams = np.empty((1, ), dtype=model.expparams_dtype)

true_trajectory = []
est_trajectory = []

true_params = prior.sample()

for idx_exp in range(400):
    # We don't want to grow the evolution time to be arbitrarily
    # large, so we'll instead choose a random time up to some
    # maximum.
    expparams[0] = np.random.random() * 10 * np.pi
    datum = model.simulate_experiment(true_params, expparams)
    updater.update(datum, expparams)

    # We index by [:, :, 0] to pull off the index corresponding
    # to experiment parameters.
    true_params = model.update_timestep(true_params, expparams)[:, :, 0]

    true_trajectory.append(true_params[0])
    est_trajectory.append(updater.est_mean())

plt.plot(true_trajectory, label='True')
plt.plot(est_trajectory, label='Est.')
plt.legend()
plt.xlabel('# of Measurements')
plt.ylabel(r'$\omega$')

plt.show()

(Source code, svg, pdf, hires.png, png)

Specifying Custom Time-Step Updates¶

The RandomWalkModel example above is somewhat unrealistic, however, in that the step distribution is independent of the evolution time. For a more reasonable noise process, we would expect that \(\mathbb{V}[\omega(t + \Delta t) - \omega(t)] \propto \Delta t\). We can subclass SimplePrecessionModel to add this behavior with a custom update_timestep() implementation.

class DiffusivePrecessionModel(SimplePrecessionModel):
    diffusion_rate = 0.0005 # We'll multiply this by
                             # sqrt(time) below.

    def update_timestep(self, modelparams, expparams):
        step_std_dev = self.diffusion_rate * np.sqrt(expparams)
        steps = step_std_dev * np.random.randn(
            # We want the shape of the steps in omega
            # to match the input model parameter and experiment
            # parameter shapes.
            # The axis of length 1 represents that this model
            # has only one model parameter (omega).
            modelparams.shape[0], 1, expparams.shape[0]
        )
        # Finally, we add a new axis to the input model parameters
        # to match the experiment parameters.
        return modelparams[:, :, np.newaxis] + steps

model = DiffusivePrecessionModel()
prior = UniformDistribution([0, 1])
updater = SMCUpdater(model, 2000, prior)

expparams = np.empty((1, ), dtype=model.expparams_dtype)

true_trajectory = []
est_trajectory = []

true_params = prior.sample()

for idx_exp in range(400):
    expparams[0] = np.random.random() * 10 * np.pi
    datum = model.simulate_experiment(true_params, expparams)
    updater.update(datum, expparams)

    true_params = model.update_timestep(true_params, expparams)[:, :, 0]

    true_trajectory.append(true_params[0])
    est_trajectory.append(updater.est_mean())

plt.plot(true_trajectory, label='True')
plt.plot(est_trajectory, label='Est.')
plt.legend()
plt.xlabel('# of Measurements')
plt.ylabel(r'$\omega$')

plt.show()

(Source code, svg, pdf, hires.png, png)

Introduction¶

In developing statistical inference applications, it is essential to test the robustness of one’s software to errors and noise of various kinds. Thus, QInfer provides tools to do so by repeatedly running estimation tasks to measure performance, and by corrupting likelihood calculations in various realistic ways.

Testing Estimation Performance¶

Given an estimator, a common question of interest is what risk that estimator incurrs; that is, what is the expected error of the estimates reported by the estimator? This is formalized by defining the loss \(L(\hat{\vec{x}}, \vec{x})\) of an estimate \(\hat{\vec{x}}\) given a true value \(\vec{x}\). In QInfer, we adopt the quadratic loss \(L_Q\) as a default, defined as

where \(\matr{Q}\) is a positive-semidefinite matrix that establishes the relative scale of each model parameter.

We can now define the risk for an estimator as

As the risk is defined by an expectation value, we can estimate it by again Monte Carlo sampling over trials. That is, for a given true model \(\vec{x}\), we can draw many different data sets and then find the corresponding estimate for each in order to determine the risk.

Similarly, the Bayes risk \(r(\hat{\vec{x}}(\cdot), \pi)\) is defined by also taking the expectation over a prior distribution \(\pi\),

The Bayes risk can thus be estimated by drawing a new set of true model parameters with each trial. QInfer implements risk and Bayes risk estimation by providing a function perf_test() which simulates a single trial with a given model, an experiment design heuristic, and either a true model parameter vector or a prior distribution. The perf_test_multiple() function then collects the results of perf_test() for many trials, reporting an array of performance metrics that can be used to quickly compute the risk and Bayes risk.

For example, we can use performance testing to evaluate the Bayes risk as a function of the number of particles used in order to determine quality parameters in experimental practice. Consider the simple precession model (SimplePrecessionModel) under an exponentially sparse sampling heuristic (ExpSparseHeuristic). Then, we can test the performance of estimating the precession frequency for several different values of n_particles:

model = SimplePrecessionModel()
prior = UniformDistribution([0, 1])
heuristic_class = ExpSparseHeuristic

for n_particles in [100, 200, 400]:
    perf = perf_test_multiple(
        n_trials=50,
        model=model, n_particles=n_particles, prior=prior,
        n_exp=50, heuristic_class=heuristic_class
    )
    # The array returned by perf_test_multiple has
    # shape (n_trials, n_exp), so to take the mean over
    # trials (Bayes risk), we need to average over the
    # zeroth axis.
    bayes_risk = perf['loss'].mean(axis=0)

    plt.semilogy(bayes_risk, label='{} particles'.format(n_particles))

plt.xlabel('# of Measurements')
plt.ylabel('Bayes Risk')
plt.legend()
plt.show()

(Source code, svg, pdf, hires.png, png)

We can also pass fixed model parameter vectors to evaluate the risk (that is, not the Bayes risk) as a function of the true model:

model = SimplePrecessionModel()
prior = UniformDistribution([0, 1])
heuristic_class = ExpSparseHeuristic
n_exp = 50

omegas = np.linspace(0.1, 0.9, 6)
risks = np.empty_like(omegas)

for idx_omega, omega in enumerate(omegas[:, np.newaxis]):
    perf = perf_test_multiple(
        n_trials=100,
        model=model, n_particles=400, prior=prior,
        n_exp=n_exp, heuristic_class=heuristic_class,
        true_mps=omega
    )
    # We now only take the loss after the last
    # measurement (indexed by -1 along axis 1).
    risks[idx_omega] = perf['loss'][:, -1].mean(axis=0)

plt.semilogy(omegas, risks)

plt.xlabel(r'$\omega$')
plt.ylabel('Risk')
plt.show()

(Source code, svg, pdf, hires.png, png)

Robustness Testing¶

Modeling Faulty or Noisy Simulators¶

In addition, QInfer allows for testing robustness against errors in the model itself by using PoisonedModel. This derived model adds noise to a model’s likelihood() method in such a way as to simulate sampling errors incurred in likelihood-hood free parameter estimation (LFPE) approaches [FG13]. The noise that PoisonedModel adds can be specified as the tolerance of an adaptive likelihood estimation (ALE) step [FG13], or as the number of samples and hedging used for a hedged maximum likelihood estimator of the likelihood [FB12]. In either case, the requested noise is added to the likelihood reported by the underlying model, such that

\[\widehat{\Pr}(d | \vec{x}; \vec{e}) = \Pr(d | \vec{x}; \vec{e}) + \epsilon,\]

where \(\widehat{\Pr}\) is the reported estimate of the true likelihood.

For example, to simulate using adaptive likelihood estimation to reach a threshold tolerance of 0.01:

>>> from qinfer import SimplePrecessionModel, PoisonedModel
>>> model = PoisonedModel(SimplePrecessionModel(), tol=0.01)

We can then use perf_test_multiple() as above to quickly test the effect of noise in the likelihood function on the Bayes risk.

models = [
    SimplePrecessionModel(),
    PoisonedModel(SimplePrecessionModel(), tol=0.25)
]
prior = UniformDistribution([0, 1])
heuristic_class = ExpSparseHeuristic

for model in models:
    perf = perf_test_multiple(
        n_trials=50,
        model=model, n_particles=400, prior=prior,
        true_model=models[0],
        n_exp=50, heuristic_class=heuristic_class
    )
    bayes_risk = perf['loss'].mean(axis=0)

    plt.semilogy(bayes_risk, label=type(model).__name__)

plt.xlabel('# of Measurements')
plt.ylabel('Bayes Risk')
plt.legend()
plt.show()

(Source code, svg, pdf, hires.png, png)

Introduction¶

QInfer provides tools to expedite simulation by distributing computation across multiple nodes using standard parallelization tools.

Distributed Computation with IPython¶

The ipyparallel package (previously IPython.parallel) provides facilities for parallelizing computation across multiple cores and/or nodes. ipyparallel separates computation into a controller that is responsible for one or more engines, and a client that sends commands to these engines via the controller. QInfer can use a client to send likelihood evaluation calls to engines, via the DirectViewParallelizedModel class.

This class takes a DirectView onto one or more engines, typically obtained with an expression similar to client[:], and splits calls to likelihood() across the engines accessible from the DirectView.

>>> from ipyparallel import Client 
>>> from qinfer import SimplePrecessionModel 
>>> from qinfer import DirectViewParallelizedModel 
>>> c = Client() 
>>> serial_model = SimplePrecessionModel() 
>>> parallel_model = DirectViewParallelizedModel(serial_model, c[:]) 

The newly decorated model will now distribute likelihood calls, such that each engine computes the likelihood for an equal number of particles. As a consequence, information shared per-experiment or per-outcome is local to each engine, and is not distributed. Therefore, this approach works best at quickly parallelizing where the per-model cost is significantly larger than the per-experiment or per-outcome cost.

>>> lbview = client.load_balanced_view() 
>>> performance = qi.perf_test_multiple(
...     100, serial_model, 6000, prior, 200, heuristic_class,
...     apply=lbview.apply
... ) 

GPGPU-based Likelihood Computation with PyOpenCL¶

Though QInfer does not yet have built-in support for GPU-based parallelization, PyOpenCL can be used to effectively distribute models as well. Here, the Cartesian product over outcomes, models and experiments matches closely the OpenCL concept of a global ID, as this example demonstrates. Once a kernel is developed in this way, PyOpenCL will allow for it to be used with any available OpenCL-compliant device.

Note that for sufficiently fast models, the overhead of copying data between the CPU and GPU may overwhelm any speed benefits obtained by this parallelization.

Introduction¶

QInfer can be used in conjunction with software written in scientific software platforms other than Python, thanks to the ready availability of interoperability libraries for interacting with Python. In this section, we provide brief examples of using these libraries with QInfer.

MATLAB Interoperability¶

As of version 2016a, MATLAB includes built-in functions for calling Python-language software. In particular, these functions can be used to use QInfer from within MATLAB. For example, the following MATLAB snippet will generate and analyze frequency estimation data using the Simple Estimation functions provided by QInfer.

>> true_omega = 70.3;
>> n_shots = 400;
>>
>> ts = pi * (1:1:100) / (2 * 100);
>>
>> signal = sin(true_omega * ts / 2) .^ 2;
>> counts = binornd(n_shots, signal);
>>
>> setenv MKL_NUM_THREADS 1
>> data = py.numpy.column_stack({counts ts ...
n_shots * ones(1, size(ts, 2))});
>> est = py.qinfer.simple_est_prec(data, ...
pyargs('freq_min', 0, 'freq_max', 100));

Importantly, the setenv command is required to work around a bug internal to the MATLAB interpreter.

Julia Interoperability¶

In Julia, interoperability can be achieved using the PyCall.jl package, which provides macros for making Python modules available as Julia variables. To install PyCall.jl, use Julia’s built-in package installer:

julia> Pkg.add("PyCall")

After installing PyCall.jl, the example above proceeds in a very similar fashion:

julia> true_omega = 70.3
julia> n_shots = 100
julia>
julia> ts = pi * (1:1:100) / (2 * 100)
julia>
julia> signal = sin(true_omega * ts / 2) .^ 2
julia> counts = map(p -> rand(Binomial(n_shots, p)), signal);
julia> @pyimport numpy as np
julia> @pyimport qinfer as qi
julia>
julia> data = [counts'; ts'; n_shots * ones(length(ts))']'
julia> est_mean, est_cov = qi.simple_est_prec(data, freq_min=0, freq_max=100)

Introduction¶

A model in QInfer is a class that describes the probabilities of observing data, given a particular experiment and given a particular set of model parameters. The observation probabilities may be given implicitly or explicitly, in that the class may only allow for sampling observations, rather than finding the a distribution explicitly. In the former case, a model is represented by a subclass of Simulatable, while in the latter, the model is represented by a subclass of Model.

Simulatable - Base Class for Implicit (Simulatable) Models¶

Model - Base Class for Explicit (Likelihood) Models¶

If a model supports explicit calculation of the likelihood function, then this is represented by subclassing from Model.

FiniteOutcomeModel - Base Class for Models with a Finite Number of Outcomes¶

The likelihood function provided by a subclass is used to implement Simulatable.simulate_experiment(), which is possible because the likelihood of all possible outcomes can be computed. This class also concretely implements the domain method by looking at the definition of n_outcomes.

DifferentiableModel - Base Class for Explicit Models with Differentiable Likelihoods¶

Introduction¶

QInfer provides several models which decorate other models, providing additional functionality or changing the behaviors of underlying models.

PoisonedModel - Model corrupted by likelihood errors¶

class qinfer.PoisonedModel(underlying_model, tol=None, n_samples=None, hedge=None)[source]¶

Bases: qinfer.derived_models.DerivedModel

Model that simulates sampling error incurred by the MLE or ALE methods of reconstructing likelihoods from sample data. The true likelihood given by an underlying model is perturbed by a normally distributed random variable \(\epsilon\), and then truncated to the interval \([0, 1]\).

The variance of \(\epsilon\) can be specified either as a constant, to simulate ALE (in which samples are collected until a given threshold is met), or as proportional to the variance of a possibly-hedged binomial estimator, to simulate MLE.

Parameters:	underlying_model (Model) – The “true” model to be poisoned. tol (float) – For ALE, specifies the given error tolerance to simulate. n_samples (int) – For MLE, specifies the number of samples collected. hedge (float) – For MLE, specifies the hedging used in estimating the true likelihood.

likelihood(outcomes, modelparams, expparams)[source]¶

simulate_experiment(modelparams, expparams, repeat=1)[source]¶

Simulates experimental data according to the original (unpoisoned) model. Note that this explicitly causes the simulated data and the likelihood function to disagree. This is, strictly speaking, a violation of the assumptions made about Model subclasses. This violation is by intention, and allows for testing the robustness of inference algorithms against errors in that assumption.

BinomialModel - Model over batches of two-outcome experiments¶

class qinfer.BinomialModel(underlying_model)[source]¶

Bases: qinfer.derived_models.DerivedModel

Model representing finite numbers of iid samples from another model, using the binomial distribution to calculate the new likelihood function.

Parameters:	underlying_model (qinfer.abstract_model.Model) – An instance of a two- outcome model to be decorated by the binomial distribution.

Note that a new experimental parameter field n_meas is added by this model. This parameter field represents how many times a measurement should be made at a given set of experimental parameters. To ensure the correct operation of this model, it is important that the decorated model does not also admit a field with the name n_meas.

decorated_model¶

expparams_dtype¶

is_n_outcomes_constant¶

Returns True if and only if the number of outcomes for each experiment is independent of the experiment being performed.

This property is assumed by inference engines to be constant for the lifetime of a Model instance.

n_outcomes(expparams)[source]¶

Returns an array of dtype uint describing the number of outcomes for each experiment specified by expparams.

Parameters:	expparams (numpy.ndarray) – Array of experimental parameters. This array must be of dtype agreeing with the expparams_dtype property.

domain(expparams)[source]¶

Returns a list of ``Domain``s, one for each input expparam.

Parameters:	expparams (numpy.ndarray) – Array of experimental parameters. This array must be of dtype agreeing with the expparams_dtype property, or, in the case where n_outcomes_constant is True, None should be a valid input.
Return type:	list of Domain

are_expparam_dtypes_consistent(expparams)[source]¶

Returns True iff all of the given expparams correspond to outcome domains with the same dtype. For efficiency, concrete subclasses should override this method if the result is always True.

Parameters:	expparams (np.ndarray) – Array of expparamms of type expparams_dtype
Return type:	bool

likelihood(outcomes, modelparams, expparams)[source]¶

simulate_experiment(modelparams, expparams, repeat=1)[source]¶

update_timestep(modelparams, expparams)[source]¶

MultinomialModel - Model over batches of D-outcome experiments¶

class qinfer.MultinomialModel(underlying_model)[source]¶

Bases: qinfer.derived_models.DerivedModel

Model representing finite numbers of iid samples from another model with a fixed and finite number of outcomes, using the multinomial distribution to calculate the new likelihood function.

Parameters:	underlying_model (qinfer.abstract_model.FiniteOutcomeModel) – An instance of a D-outcome model to be decorated by the multinomial distribution. This underlying model must have is_n_outcomes_constant as True.

Note that a new experimental parameter field n_meas is added by this model. This parameter field represents how many times a measurement should be made at a given set of experimental parameters. To ensure the correct operation of this model, it is important that the decorated model does not also admit a field with the name n_meas.

decorated_model¶

expparams_dtype¶

is_n_outcomes_constant¶

Returns True if and only if the number of outcomes for each experiment is independent of the experiment being performed.

This property is assumed by inference engines to be constant for the lifetime of a Model instance.

n_sides¶

Returns the number of possible outcomes of the underlying model.

underlying_domain¶

Returns the Domain of the underlying model.

n_outcomes(expparams)[source]¶

Returns an array of dtype uint describing the number of outcomes for each experiment specified by expparams.

Parameters:	expparams (numpy.ndarray) – Array of experimental parameters. This array must be of dtype agreeing with the expparams_dtype property.

domain(expparams)[source]¶

Returns a list of Domain objects, one for each input expparam. :param numpy.ndarray expparams: Array of experimental parameters. This

array must be of dtype agreeing with the expparams_dtype property.

Return type:	list of Domain

are_expparam_dtypes_consistent(expparams)[source]¶

Returns True iff all of the given expparams correspond to outcome domains with the same dtype. For efficiency, concrete subclasses should override this method if the result is always True.

Parameters:	expparams (np.ndarray) – Array of expparamms of type expparams_dtype
Return type:	bool

likelihood(outcomes, modelparams, expparams)[source]¶

simulate_experiment(modelparams, expparams, repeat=1)[source]¶

MLEModel - Model for approximating maximum-likelihood estimation¶

class qinfer.MLEModel(underlying_model, likelihood_power)[source]¶

Bases: qinfer.derived_models.DerivedModel

Uses the method of [JDD08] to approximate the maximum likelihood estimator as the mean of a fictional posterior formed by amplifying the Bayes update by a given power \(\gamma\). As \(\gamma \to \infty\), this approximation to the MLE improves, but at the cost of numerical stability.

Parameters:	likelihood_power (float) – Power to which the likelihood calls should be rasied in order to amplify the Bayes update.

simulate_experiment(modelparams, expparams, repeat=1)[source]¶

likelihood(outcomes, modelparams, expparams)[source]¶

Distribution - Abstract Base Class for Probability Distributions¶

class qinfer.Distribution[source]¶

Bases: object

Abstract base class for probability distributions on one or more random variables.

n_rvs¶

The number of random variables that this distribution is over.

Type:	int

sample(n=1)[source]¶

Returns one or more samples from this probability distribution.

Parameters:	n (int) – Number of samples to return.
Return type:	numpy.ndarray
Returns:	An array containing samples from the distribution of shape (n, d), where d is the number of random variables.

Specific Distributions¶

class qinfer.UniformDistribution(ranges=array([[0, 1]]))[source]¶

Bases: qinfer.distributions.Distribution

Uniform distribution on a given rectangular region.

Parameters:	ranges (numpy.ndarray) – Array of shape (n_rvs, 2), where n_rvs is the number of random variables, specifying the upper and lower limits for each variable.

n_rvs¶

sample(n=1)[source]¶

grad_log_pdf(var)[source]¶

class qinfer.DiscreteUniformDistribution(num_bits)[source]¶

Bases: qinfer.distributions.Distribution

Discrete uniform distribution over the integers between 0 and 2**num_bits-1 inclusive.

Parameters:	num_bits (int) – non-negative integer specifying how big to make the interval.

n_rvs¶

sample(n=1)[source]¶

class qinfer.MVUniformDistribution(dim=6)[source]¶

Bases: qinfer.distributions.Distribution

Uniform distribution over the rectangle \([0,1]^{\text{dim}}\) with the restriction that vector must sum to 1. Equivalently, a uniform distribution over the dim-1 simplex whose vertices are the canonical unit vectors of \(\mathbb{R}^\text{dim}\).

Parameters:	dim (int) – Number of dimensions; n_rvs.

n_rvs¶

sample(n=1)[source]¶

class qinfer.NormalDistribution(mean, var, trunc=None)[source]¶

Bases: qinfer.distributions.Distribution

Normal or truncated normal distribution over a single random variable.

Parameters:	mean (float) – Mean of the represented random variable. var (float) – Variance of the represented random variable. trunc (tuple) – Limits at which the PDF of this distribution should be truncated, or None if the distribution is to have infinite support.

n_rvs¶

sample(n=1)[source]¶

grad_log_pdf(x)[source]¶

class qinfer.MultivariateNormalDistribution(mean, cov)[source]¶

Bases: qinfer.distributions.Distribution

Multivariate (vector-valued) normal distribution.

Parameters:	mean (np.ndarray) – Array of shape (n_rvs, ) representing the mean of the distribution. cov (np.ndarray) – Array of shape (n_rvs, n_rvs) representing the covariance matrix of the distribution.

n_rvs¶

sample(n=1)[source]¶

grad_log_pdf(x)[source]¶

class qinfer.SlantedNormalDistribution(ranges=array([[0, 1]]), weight=0.01)[source]¶

Bases: qinfer.distributions.Distribution

Uniform distribution on a given rectangular region with additive noise. Random variates from this distribution follow \(X+Y\) where \(X\) is drawn uniformly with respect to the rectangular region defined by ranges, and \(Y\) is normally distributed about 0 with variance weight**2.

Parameters:	ranges (numpy.ndarray) – Array of shape (n_rvs, 2), where n_rvs is the number of random variables, specifying the upper and lower limits for each variable. weight (float) – Number specifying the inverse variance of the additive noise term.

n_rvs¶

sample(n=1)[source]¶

class qinfer.LogNormalDistribution(mu=0, sigma=1)[source]¶

Bases: qinfer.distributions.Distribution

Log-normal distribution.

Parameters:	mu – Location parameter (numeric), set to 0 by default. sigma – Scale parameter (numeric), set to 1 by default. Must be strictly greater than zero.

n_rvs¶

sample(n=1)[source]¶

class qinfer.ConstantDistribution(values)[source]¶

Bases: qinfer.distributions.Distribution

Represents a determinstic variable; useful for combining with other distributions, marginalizing, etc.

Parameters:	values – Shape (n,) array or list of values \(X_0\) such that \(\Pr(X) = \delta(X - X_0)\).

n_rvs¶

sample(n=1)[source]¶

class qinfer.BetaDistribution(alpha=None, beta=None, mean=None, var=None)[source]¶

Bases: qinfer.distributions.Distribution

The beta distribution, whose pdf at \(x\) is proportional to \(x^{\alpha-1}(1-x)^{\beta-1}\). Note that either alpha and beta, or mean and var, must be specified as inputs; either case uniquely determines the distribution.

Parameters:	alpha (float) – The alpha shape parameter of the beta distribution. beta (float) – The beta shape parameter of the beta distribution. mean (float) – The desired mean value of the beta distribution. var (float) – The desired variance of the beta distribution.

n_rvs¶

sample(n=1)[source]¶

class qinfer.BetaBinomialDistribution(n, alpha=None, beta=None, mean=None, var=None)[source]¶

Bases: qinfer.distributions.Distribution

The beta-binomial distribution, whose pmf at the non-negative integer \(k\) is equal to \(\binom{n}{k}\frac{B(k+\alpha,n-k+\beta)}{B(\alpha,\beta)}\) with \(B(\cdot,\cdot)\) the beta function. This is the compound distribution whose variates are binomial distributed with a bias chosen from a beta distribution. Note that either alpha and beta, or mean and var, must be specified as inputs; either case uniquely determines the distribution.

Parameters:	n (int) – The \(n\) parameter of the beta-binomial distribution. alpha (float) – The alpha shape parameter of the beta-binomial distribution. beta (float) – The beta shape parameter of the beta-binomial distribution. mean (float) – The desired mean value of the beta-binomial distribution. var (float) – The desired variance of the beta-binomial distribution.

n_rvs¶

sample(n=1)[source]¶

class qinfer.GammaDistribution(alpha=None, beta=None, mean=None, var=None)[source]¶

Bases: qinfer.distributions.Distribution

The gamma distribution, whose pdf at \(x\) is proportional to \(x^{-\alpha-1}e^{-x\beta}\). Note that either alpha and beta, or mean and var, must be specified as inputs; either case uniquely determines the distribution.

Parameters:	alpha (float) – The alpha shape parameter of the gamma distribution. beta (float) – The beta shape parameter of the gamma distribution. mean (float) – The desired mean value of the gamma distribution. var (float) – The desired variance of the gamma distribution.

n_rvs¶

sample(n=1)[source]¶

class qinfer.InterpolatedUnivariateDistribution(pdf, compactification_scale=1, n_interp_points=1500)[source]¶

Bases: qinfer.distributions.Distribution

Samples from a single-variable distribution specified by its PDF. The samples are drawn by first drawing uniform samples over the interval [0, 1], and then using an interpolation of the inverse-CDF corresponding to the given PDF to transform these samples into the desired distribution.

Parameters:	pdf (callable) – Vectorized single-argument function that evaluates the PDF of the desired distribution. compactification_scale (float) – Scale of the compactified coordinates used to interpolate the given PDF. n_interp_points (int) – The number of points at which to sample the given PDF.

n_rvs¶

sample(n=1)[source]¶

class qinfer.HilbertSchmidtUniform(dim=2)[source]¶

Bases: qinfer.distributions.SingleSampleMixin, qinfer.distributions.Distribution

Creates a new Hilber-Schmidt uniform prior on state space of dimension dim. See e.g. [Mez06] and [Mis12].

Parameters:	dim (int) – Dimension of the state space.

n_rvs¶

sample()[source]¶

make_Paulis(paulis, d)[source]¶

class qinfer.HaarUniform(dim=2)[source]¶

Bases: qinfer.distributions.SingleSampleMixin, qinfer.distributions.Distribution

Haar uniform distribution of pure states of dimension dim, parameterized as coefficients of the Pauli basis.

Parameters:	dim (int) – Dimension of the state space.

Note

This distribution presently only works for dim==2 and the Pauli basis.

n_rvs¶

class qinfer.GinibreUniform(dim=2, k=2)[source]¶

Bases: qinfer.distributions.SingleSampleMixin, qinfer.distributions.Distribution

Creates a prior on state space of dimension dim according to the Ginibre ensemble with parameter k. See e.g. [Mis12].

Parameters:	dim (int) – Dimension of the state space.

n_rvs¶

Combining Distributions¶

QInfer also offers classes for combining distributions together to produce new ones.

class qinfer.ProductDistribution(*factors)[source]¶

Bases: qinfer.distributions.Distribution

Takes a non-zero number of QInfer distributions \(D_k\) as input and returns their Cartesian product.

In other words, the returned distribution is \(\Pr(D_1, \dots, D_N) = \prod_k \Pr(D_k)\).

Parameters:	factors (Distribution) – Distribution objects representing \(D_k\). Alternatively, one iterable argument can be given, in which case the factors are the values drawn from that iterator.

n_rvs¶

sample(n=1)[source]¶

class qinfer.PostselectedDistribution(distribution, model, maxiters=100)[source]¶

Bases: qinfer.distributions.Distribution

Postselects a distribution based on validity within a given model.

n_rvs¶

sample(n=1)[source]¶

Returns one or more samples from this probability distribution.

Parameters:	n (int) – Number of samples to return.
Return numpy.ndarray:
	An array containing samples from the distribution of shape (n, d), where d is the number of random variables.

grad_log_pdf(x)[source]¶

class qinfer.MixtureDistribution(weights, dist, dist_args=None, dist_kw_args=None, shuffle=True)[source]¶

Bases: qinfer.distributions.Distribution

Samples from a weighted list of distributions.

Parameters:

weights – Length n_dist list or np.ndarray of probabilites summing to 1. dist – Either a length n_dist list of Distribution instances, or a Distribution class, for example, NormalDistribution. It is assumed that a list of Distribution``s all have the same ``n_rvs. dist_args – If dist is a class, an array of shape (n_dist, n_rvs) where dist_args[k,:] defines the arguments of the k’th distribution. Use None if the distribution has no arguments. dist_kw_args – If dist is a class, a dictionary where each key’s value is an array of shape (n_dist, n_rvs) where dist_kw_args[key][k,:] defines the keyword argument corresponding to key of the k’th distribution. Use None if the distribution needs no keyword arguments. shuffle (bool) – Whether or not to shuffle result after sampling. Not shuffling will result in variates being in the same order as the distributions. Default is True.

n_rvs¶

n_dist¶

The number of distributions in the mixture distribution.

sample(n=1)[source]¶

class qinfer.ConstrainedSumDistribution(underlying_distribution, desired_total=1)[source]¶

Bases: qinfer.distributions.Distribution

Samples from an underlying distribution and then enforces that all samples must sum to some given value by normalizing each sample.

Parameters:	underlying_distribution (Distribution) – Underlying probability distribution. desired_total (float) – Desired sum of each sample.

underlying_distribution¶

n_rvs¶

sample(n=1)[source]¶

Mixins for Distribution Development¶

class qinfer.SingleSampleMixin[source]¶

Bases: object

Mixin class that extends a class so as to generate multiple samples correctly, given a method _sample that generates one sample at a time.

sample(n=1)[source]¶

Introduction¶

A Domain represents a collection of objects. They are used by Simulatable (and subclasses like Model and FiniteOutcomeModel) to store relevant information about the possible outcomes of a given experiment. This includes properties like whether or not there are a finite number of possibilities, if so how many, and what their data types are.

Domain - Base Class for Domains¶

All domains should inherit from this base class.

RealDomain - (A subset of) Real Numbers¶

IntegerDomain - (A subset of) Integers¶

This is the default domain for FiniteOutcomeModel.

MultinomialDomain - Tuples of Integers with a Constant Sum¶

This domain is used by MultinomialModel.