QInfer: Statistical inference software for quantum applications

Characterizing quantum systems through experimental data is critical to applications as diverse as metrology and quantum computing. Analyzing this experimental data in a robust and reproducible manner is made challenging, however, by the lack of readily-available software for performing principled statistical analysis. We improve the robustness and reproducibility of characterization by introducing an open-source library, QInfer, to address this need. Our library makes it easy to analyze data from tomography, randomized benchmarking, and Hamiltonian learning experiments either in post-processing, or online as data is acquired. QInfer also provides functionality for predicting the performance of proposed experimental protocols from simulated runs. By delivering easy-to-use characterization tools based on principled statistical analysis, QInfer helps address many outstanding challenges facing quantum technology.


Introduction
Statistical modeling and parameter estimation play a critical role in many quantum applications. In quantum information in particular, the pursuit of largescale quantum information processing devices has motivated a range of different characterization protocols, and in turn, new statistical models. For example, quantum state and process tomography are widely used to characterize quantum systems, and are in essence matrix-valued parameter estimation problems [1,2]. Similarly, randomized benchmarking is now a mainstay in assessing quantum devices, motivating the use of rigorous statistical analysis [3] and algorithms [4]. Quantum metrology, meanwhile, is intimately concerned with what parameters can be extracted from measurements of a physical system, immediately necessitating a statistical view [5,6].
The prevalence of statistical modeling in quantum applications should not be surprising: quantum mechanics is an inherently statistical theory, thus inference is an integral part of both experimental and theoretical practice. In the former, experimentalists need to model their systems and infer the value of parameters for the purpose of improving control as well as validating performance. In the latter, numerical experiments utilizing simulated data are now commonplace in theoretical studies, such that the same inference problems are encountered usually as a necessary step to answer questions about optimal data processing protocols or experiment design. In both cases, we lack tools to rapidly prototype and access inference strategies; QInfer addresses this need by providing a modular interface to a Monte Carlo algorithm for performing statistical inference.
Critically, in doing so, QInfer also supports and enables open and reproducible research practices. Parallel to the challenges faced in many other disciplines [7], physics research cannot long survive its own current practices. Open access, open source, and open data provide an indispensable means for research to be reproducible, ensuring that research work is useful to the communities invested in that research [8]. In the particular context of quantum information research, open methods are especially critical given the impact of statistical errors that can undermine the claims of published research [9,10]. Ensuring the reproducibility of research is critical for evaluating the extent to which statistical and methodological errors undermine the credibility of published research [11].
QInfer also constitutes an important step towards a more general framework for quantum verification and validation (QCVV). As quantum information processor prototypes become more complex, the challenge of ensuring that noise processes affecting these devices conform to some agreed-upon standard becomes more difficult. This challenge can be managed, at least in principle, by developing confidence in the truth of certain simplifying assumptions and approximations. The value of randomized benchmarking, for example, depends strongly upon the extent to which noise is approximately Pauli [12]. QInfer provides a valuable framework for the design of automated and efficient noise assessment methods that will enable the compar-ison of actual device performance to the specifications demanded by theory.
To the end of enabling reproducible and accessible research, and hence providing a reliable process for interpreting advances in quantum information processing, we base QInfer using openly-available tools such as the Python programming language, the IPython interpreter, and Jupyter [13,14]. Jupyter in particular has already proven to be an invaluable tool for reproducible research, in that it provides a powerful framework for describing and explaining research software [15]. We provide our library under an open-source license along with examples [16] of how to use QInfer to support reproducible research practices. In this way, our library builds on and supports recent efforts to develop reproducible methods for physics research [17].
QInfer is a mature open-source software library written in the Python programming language which has now been extensively tested in a wide range of inferential problems by various research groups. Recognizing its maturity through its continuing development, we now formally release version 1.0. This maturity has given its developers the opportunity to step back and focus on the accessibility of QInfer such that other researchers can benefit from its utility. This short paper is the culmination of that effort. A full Users' Guide is available in the ancillary files.
We proceed as following. In Section 2, we give a brief introduction to Bayesian inference and particle filtering, the numerical algorithm we use to implement Bayesian updates. In Section 3, we describe applications of QInfer to common tasks in quantum information processing. Next, we describe in Section 4 additional features of QInfer before concluding in Section 5.

Inference and Particle Filtering
QInfer is primarily intended to serve as a toolset for implementing Bayesian approaches to statistical inference. In this section, we provide a brief review of the Bayesian formalism for statistical inference. This section is not intended to be comprehensive; our aim is rather to establish the language needed to describe the QInfer codebase.
In the Bayesian paradigm, statistical inference is the process of evaluating data obtained by sampling an unknown member of a family of related probability distributions, then using these samples to assign a relative plausibility to each distribution. Colloquially, we think of this family of distributions as a model parameterized by a vector x of model parameters. We then express the probability that a dataset D was obtained from the model parameters x as Pr(D|x) and read it as "the probability of D given that the model specified by x is the correct model." The function Pr(·|x) is called the likelihood function, and computing it is equivalent to simulating an experiment 1 . For example, the Born rule is a likelihood function, in that it maps a known or hypothetical quantum density matrix x ≡ ρ to a distribution over measurement outcomes of a measurement The problem of estimating model parameters is as follows. Suppose an agent is provided with a dataset D and is tasked with judging the probability that the model specified by a given vector x is in fact the correct one. According to Bayes' rule, where Pr(x) is a probability distribution called the prior distribution and Pr(x|D) is called the posterior distribution. If the agent is provided with a prior distribution, then they can estimate parameters using Bayes' rule. Note that Pr(D) can be computed through marginalization, which is to say that the value can in principle be calculated via the equation For the inference algorithm used by QInfer, Pr(D) is an easily computed normalization constant and there is no need to compute a possibly complicated integral. Importantly, we will demand that the agent's data processing approach works in an iterative manner. Consider the example in which the data D is in fact a set D = {d 1 , . . . , d N } of individual observations. In most if not all classical applications, each individual datum is distributed independently of the rest of the data set, conditioned on the true state. Formally, we write that for all j and k such that This may not hold in quantum models where measurement back-action can alter the state. In such cases, we can simply redefine what the parameters x label, such that this independence property can be taken as a convention, instead of as an assumption. Then, we have that In other words, the agent can process the data sequentially where the prior for each successive datum is the posterior from the last. This Bayes update can be solved analytically in some important special cases, such as frequency estimation [19,20], but is more generally intractable. Thus, to develop a robust and generically useful framework for parameter estimation, the agent relies on numerical algorithms. In particular, QInfer is largely concerned with the particle filtering algorithm [21], also known as the sequential Monte Carlo (SMC) algorithm. In the context of quantum information, SMC was first proposed for learning from continuous measurement records [22], and has since been used to learn from state tomography [23], Hamiltonian learning [24], and randomized benchmarking [4], as well as other applications.
The aim of particle filtering is to replace a continuous probability distribution Pr(x) with a discrete approximation where w = (w k ) is a vector of probabilities. The entry w k is called the weight of the particle, labeled k, and x k is the location of particle k. Of course, the particle filter k w k δ(x − x k ) does not directly approximate Pr(x) as a distribution; the particle filter, if considered as a distribution, is supported only a discrete set of points. Instead, the particle filter is used to approximate expectation values: if f is a function whose domain is the set of model vectors x, we want the particle filter to satisfy The posterior distribution can also be approximated using a particle filter. In fact, a posterior particle filter can be computed directly from a particle filter for the prior distribution as follows. Let {(w k , x k )} be the set of weights and locations for a particle filter for some prior distribution Pr(x). We then compute a particle filter {(w k , x k )} for the posterior distribution by setting x k = x k and where D is the data set used in the Bayesian update. In practice, updating the weights in this fashion causes the particle filter to become unstable as data is collected; by default, QInfer will periodically apply the Liu-West algorithm to restore stability [25]. See Appendix B for details.
At any point during the processing of data, the expectation of any function with respect to the posterior is approximated as In particular, the expected error in x is given by the posterior covariance, Cov(x|D) . This can be used, for instance, to adaptively choose experiments which minimize the posterior variance [24]. This approach has been used to exponentially improve the number of samples required in frequency estimation problems [19,20], and in phase estimation [26,27]. Alternatively, other cost functions can be considered, such as the information gain [23,28]. QInfer allows for quickly computing either the expected posterior variance or the information gain for proposed experiments, making it straightforward to develop adaptive experiment design protocols. The functionality exposed by QInfer follows a simple object model, in which the experiment is described in terms of a model, and background information is described in terms of a prior distribution. Each of these classes is abstract, meaning that they define what behavior a QInfer user must specify in order to fully specify an inference procedure. For convenience, QInfer provides several pre-defined implementations of each, as we will see in the following examples. Concrete implementations of a model and a prior distribution are then used with SMC to update the prior based on data. In summary, the iterative approach described above is formalized in terms of the following Python object model: A complete description of the QInfer object model can be found in the Users' Guide.
Notably qinfer . SMCUpdater relies only on the behavior specified by each of the abstract classes in this object model. Thus, it is straightforward for the user to specify their own prior and likelihood function by either implementing these classes (as in the example of Appendix A), or by using one of the many concrete implementations provided with QInfer.
The concrete implementations provided with QInfer are useful in a range of common applications, as described in the next Section. We will demonstrate how these classes are used in practice with examples drawn from quantum information applications. We will also consider the qinfer . Heuristic class, which is useful in contexts such as online adaptive experiments and simulated experiments.

Applications in Quantum Information
In this Section, we describe various possible applications of QInfer to existing experimental protocols. In doing so, we highlight both functionality built-in to QInfer and how this functionality can be readily extended with custom models and distributions. We begin with the problems of phase and frequency learning, then describe the use of QInfer for state and process tomography, and conclude with applications to randomized benchmarking.

Phase and Frequency Learning
One of the primary applications for particle filtering is for learning the Hamiltonian H under which a quantum system evolves [24]. For instance, consider the single-qubit Hamiltonian H = ωσ z /2 for an unknown parameter ω. An experiment on this qubit may then consist of preparing a state |+ = (|0 + |1 )/ √ 2, evolving for a time t and then measuring in the σ x basis. This model commonly arises from Ramsey interferometry, and gives a likelihood function Note that this is also the same model for Rabi interferometry as well, with the interpretation of H as drive term rather than the internal Hamiltonian for a system. Similarly, this model forms the basis of Bayesian and maximum likelihood approaches to phase estimation.
In any case, QInfer implements (9) as the SimplePrecessionModel class, making it easy to quickly perform Bayesian inference for Ramsey or Rabi estimation problems. We demonstrate this in Listing 1, using ExpSparseHeuristic to select the kth measurement time t k = (9/8) k , as suggested by analytic arguments [19].
Listing 1: Frequency estimation example using SimplePrecessionModel. updater . update ( datum , experiment ) >>> p r i n t ( updater . est mean ( ) ) More complicated models for learning Hamiltonians with particle filtering have also been considered [29][30][31][32]; these can be readily implemented in QInfer as custom models by deriving from the Model class, as described in Appendix A.

State and Process Tomography
Though originally conceived of as a algebraic inverse problem, quantum tomography is also a problem of parameter estimation. Many have also considered the problem in a Bayesian framework [33,34] and the sequential Monte Carlo algorithm has been used in both theoretical and experimental studies [23,28,[35][36][37].
To define the model, we start with a basis for traceless Hermitian operators {B j } d 2 −1 j=1 . In the case of a qubit, this could be the basis of Pauli matrices, for ex-ample. Then, any state ρ can be written for some vector of parameters θ. These parameters must be constrained such that ρ ≥ 0.
In the simplest case, we can consider two-outcome measurements represented by the pair {E, 1 − E}. The Born rule defines the likelihood function For multiple measurements, we simply iterate. For many trials of the same measurement, we can use a derived model as discussed below. QInfer's TomographyModel abstracts many of the implementation details of this problem, exposing tomographic models and estimates in terms of QuTiP's Qobj class [38]. This allows for readily integrating QInfer functionality with that of QuTiP, such as fidelity metrics, diamond norm calculation, and other such manipulations.
Tomography support in QInfer requires one of the bases mentioned above in order to parameterize the state. Many common choices of basis are included as TomographyBasis objects, such as the Pauli or Gell-Mann bases. Many of the most commonly used priors are already implemented as a QInfer Distribution.
Listing 2: Rebit state tomography example using TomographyModel. For simulations, common randomized measurement choices are already implemented. For example, RandomPauliHeuristic chooses random Pauli measurements for qubit tomography.
In Listing 2, we demonstrate QInfer's tomography support for a rebit. By analogy to the Bloch sphere, a rebit may be represented by a point in the unit disk, making rebit tomography useful for plotting examples. More generally, with different choices of basis, QInfer can be used for qubits or higher-dimensional states. For example, recent work has demonstrated the use of QInfer for tomography procedures on sevendimensional states [39]. Critically, QInfer provides a region estimate for this example, describing a region that has a 95% probability of containing the true state. We will explore region estimation further in Section 4.1.
Finally, we note that process tomography is a special case of state tomography [36], such that the same functionality described above can also be used to analyze process tomography experiments. In particular, the qinfer . ProcessTomographyHeuristic class represents the experiment design constraints imposed by process tomography, while BCSZChoiDistribution uses the distribution over completely positive trace-preserving maps proposed by Bruzda et al. [40] to represent a prior distribution over the Choi states of random channels.

Randomized Benchmarking
In recent years, randomized benchmarking (RB) has reached a critical role in evaluating candidate quantum information processing systems. By using random sequences of gates drawn from the Clifford group, RB provides a likelihood function that depends on the fidelity with which each Clifford group element is implemented, allowing for estimates of that fidelity to be drawn from experimental data [4].
In particular, suppose that each gate is implemented with fidelity F , and consider a fixed initial state and measurement. Then, the survival probability over sequences of length m is given by [41] where p : an estimate of the fidelity of interest F . The likelihood function for randomized benchmarking is extremely simple, and requires only scalar arithmetic to compute, making it especially useful for avoiding the computational overhead typically required to characterize large quantum systems with classical resources. Multiple generalizations of RB have been recently developed which extend these benefits to estimating crosstalk [42], coherence [43], and to estimating fidelities of non-Clifford gates [44,45]. RB has also been extended to provide tomographic information as well [46]. The estimates provided by randomized benchmarking have also been applied to design improved control sequences [47,48].
QInfer supports RB experiments through the qinfer . RandomizedBenchmarkingModel class. For common priors, QInfer also provides a simplified interface, qinfer . simple est rb, that reports the mean and covariance over an RB model given experimental data. We provide an example in Listing 3.

Additional Functionality
Having introduced common applications for QInfer, in this Section we describe additional functionality which can be used with each of these applications, or with custom models.

Region Estimation and Error Bars
As an alternative to specifying the entire posterior distribution approximated by qinfer . SMCUpdater, we provide methods for reporting credible regions over the posterior, based on covariance ellipsoids, convex hulls, and minimum volume enclosing ellipsoids [35]. These region estimators provide a rigorous way of summarizing one's uncertainty following an experiment (colloquially referred to as "error bars"), and owing to the Bayesian approach, do so in a manner consistent with experimental experience.
Posterior credible regions can be found by using the SMCUdater.est credible region method. This method returns a set of particles such that the sum of their weights corresponding weights is at least a specified ratio of the total weight. For example, a 95% credible regions is represented as a collection of particles whose weight sums to at least 0.95. This does not necessarily admit a very compact description since many of the particles would be interior to the regions. In such cases, it is useful to find region estimators containing all of the particles describing a credible region. The SMCUpdater.region est hull method does this by finding a convex hull of the credible particles. Such a hull is depicted in Figure 2.
The convex hull of an otherwise random set of points is also not necessarily easy to describe or intuit. In such cases, SMCUdpater.region est ellipsoid finds the minimumvolume enclosing ellipse (MVEE) of the convex hull region estimator. As the name suggests, this is the smallest ellipsoid containing the credible particles. It is strictly larger than the hull and thus maintains credibility. Ellipsoids are specified by their center and covariance matrix. Visualizing the covariance matrix can also usually provide important diagnostic information, as in Figure 3. In that example, we can quickly see that the p and A parameters estimated from a randomized benchmarking experiment are anti- correlated, such that we can explain more preparation and measurement errors by assuming better gates and vice versa.

Derived Models
QInfer allows for the notion of a model chain, where the likelihood of a given model in the chain is a function of the likelihoods of models below it, and possibly new model or experiment parameters. This abstraction is useful for a couple of reasons. It encourages more robust programs, since models in the chain will often be debugged independently. It also often makes writing new models easier since part of the chain may be included by default in the QInfer library, or may overlap with other similar models the user is implementing. Finally, in quantum systems, it is common to have a likelihood function which is most naturally expressed as a hierarchical probability distribution, with base models describing quantum physics, and overlying models describing measurement processes. Model chains are typically implemented through the use of the abstract class DerivedModel. Since this class itself inherits from the Model class, subclass instances must provide standard model properties and methods such as likelihood, n outcomes, and modelparam names. Additionally, DerivedModel accepts an argument model, referring to the underlying model directly below it in the model chain. Class properties exist for referencing models at arbitrary depths in the chain, all the way down to the base model. As an example, consider a base model which is the precession model discussed in Section 3.1. This is a two-outcome model whose outcomes correspond to measuring the state |+ or the orthogonal state |− , which can be viewed as flipping a biased coin. Perhaps an actual experiment of this system consists of flipping the coin N times with identical settings, where the individual results are not recorded, only the total number n + of |+ results. In this case, we can concatenate this base model with the built-in DerivedModel called BinomialModel. This model adds an additional experiment parameter n meas specifying how many times the underlying model's coin is flipped in a single experiment. Note that parallelization, discussed in Section 4.5, is implemented as a DerivedModel whose likelihood batches the underlying model's likelihood function across processors.

Time-Dependent Models
So far, we have only considered time-independent (parameter estimation) models, but particle filtering is useful more generally for estimating time-dependent (state-space) models. Following the work of Isard and Blake [49], when performing a Bayes update, we may also incorporate state-space dynamics by adding a time-step update. For example, to follow a Wiener process, we move each particle x i (t k ) at time t k to its new position with η ∼ N(0, Σ) for a covariance matrix Σ. Importantly, we need not assume that timedependence in x follows specifically a Wiener process.
For instance, one may consider timestep increments describing stochastic evolution of a system undergoing weak measurement [22], such as an atomic ensemble undergoing probing by an optical interferometer [50]. In each case, QInfer uses the timestep increment implemented by the Model.update timestep method, which specifies the time step that SMCUpdater should perform after each datum. This design allows for the specification of more complicated time step updates than the representative example of (13). For instance, the qinfer . RandomWalkModel class adds diffusive steps to existing models and can be used to quickly learn time-dependent properties, such as shown in Listing 5. Moreover, QInfer provides the DiffusiveTomographyModel for including time-dependence in tomography by truncating time step updates to lie within the space of valid states [36]. A video example of time-dependent tomography can be found on YouTube [51].
In this way, by following the Isard and Blake [49] algorithm, we obtain a very general solution for timedependence. Importantly, other approaches exist that may be better suited for individual problems, including modifying resampling procedures to incorporate additional noise [52,53], or adding hyperparameters to describe deterministic time-dependence [36].

Performance and Robustness Testing
One important application of QInfer is predicting how well a particular parameter estimation experiment will work in practice. This can be formalized by considering the risk R( In both cases, QInfer automates the process of performing many independent estimation trials through the perf test multiple function. This function will run an updater loop for a given model, prior, and experiment design heuristic, returning the errors incurred after each measurement in each trial. Taking an expectation value with numpy.mean returns the risk or Bayes risk, depending if the true model keyword argument is set. For example, Listing 6 finds the Bayes risk for a frequency estimation experiment (Section 3.1) as a function of the number of measurements performed.
Performance evaluation can also easily be paral- lelized over trials, as discussed in Section 4.5, allowing for efficient use of computational resources. This is especially important when comparing performance for a range of different parameters. For instance, one might want to consider how the risk and Bayes risk of an estimation procedure scale with errors in a faulty simulator; QInfer supports this usecase with the qinfer . PoisonedModel derived model, which adds errors to an underlying "valid" model. In this way, QInfer enables quickly reasoning about how much approximation error can be tolerated by an estimation procedure.

Parallelization
At each step of the SMC algorithm, the likelihood Pr(d n |x) of an experimental datum d n is computed for every particle x k in the distribution. Typically, the total running time of the algorithm is overwhelmingly spent calculating these likelihoods. However, individual likelihood computations are independent of each other and therefore may be performed in parallel. On a single computational node with multiple cores, limited parallelization is performed automatically by relying on NumPy's vectorization primitives [54].
More generally, however, if the running time of Pr(d n |x) is largely independent of x, we may divide our particles into L disjoint groups, and send each group along with d n to a separate processor to be computed in parallel. In QInfer, this is handled by the derived model (Sec-tion 4.2) qinfer . DirectViewParallelizedModel which uses the Python library ipyparallel [55]. This library supports everything from simple parallelization over the cores of a single processor, to make-shift clusters set up over SSH, to professional clusters using standard job schedulers. Passing the model of interest as well as an ipyparallel . DirectView of the processing engines is all that is necessary to parallelize a model. In Figure 7, a roughly 12× speed-up is demonstrated by parallelizing a model over the 12 cores of a single computer. This model was contrived to demonstrate the parallelization potential of a generic Hamiltonian learning problem which uses dense operators and states. A single likelihood call generates a random 16 × 16 anti-hermitian matrix (representing the generator of a four qubit system), exponentiates it, and returns overlap with the |0000 state. Implementation details can be found in the QInfer examples repository [16], or in the ancillary files.
So far, we have discussed parallelization from the perspective of traditional processors (CPUs), which typically have a small number of processing cores on each chip. By contrast, moderately-priced desktop graphical processing units (GPUs) will often contain thousands of cores, while GPU hosts tailored for scientific use can have tens of thousands. This massive parallelization makes GPUs attractive for particle filtering [56]. Using libraries such as PyOpenCL and Py-CUDA [57] or Numba [58], custom models can be written which take advantage of GPUs within QInfer [52]. For example, qinfer . AcceleratedPrecessionModel offloads its computation of cos 2 to GPUs using PyOpenCL.

Other Features
In addition to the functionality described above, QInfer has a wide range of other features that we describe more briefly here. A complete description can be found in the provided Users' Guide (see ancillary files or docs.qinfer.org).

Plotting and Visualization
QInfer provides plotting and visualization support based on matplotlib [59] and mpltools [60]. In particular, qinfer . SMCUpdater provides methods for plotting posterior distributions and covariance matrices. These methods make it straightforward to visually diagnose the operation of and results obtained from particle filtering.
Similarly, the qinfer . tomography module provides several functions for producing plots of states and distributions over rebits (qubits restricted to real numbers). Rebit visualization is in particular useful for demonstrating the conceptual operation of particle filterbased tomography in a clear and attractive manner.

Fisher Information Calculation
In evaluating estimation protocols, it is important to establish a baseline of how accurately one can estimate a model even in principle. Similarly, such a baseline can be used to compare between protocols by informing as to how much information can be extracted from a proposed experiment. The Cramér-Rao bound and its Bayesian analog, the van Trees inequality (a.k.a. the Bayesian Cramér-Rao bound), formalize this notion in terms of the Fisher information matrix [61,62]. For any model which specifies its derivative in terms of a score, QInfer will calculate each of these bounds, providing useful information about proposed experimental and estimation protocols. The qinfer . ScoreMixin class builds on this by calculating the score of an arbitrary model using numerical differentiation.

Model Selection and Averaging
Statistical inference does not require asserting a priori the correctness of a particular model (that is, likelihood function), but allows a model to be taken as a hypothesis and compared to other models. This is made formal by model selection. From a Bayesian perspective, the ratio of the posterior normalizations for two different models gives a natural and principled model selection criterion, known as the Bayes factor [63]. The Bayes factor provides a model selection rule that is significantly more robust to outlying data than conventional hypothesis testing approaches [64]. For example, in quantum applications, the Bayes factor is particularly useful in tomography, and can be used to decide the rank or dimension of a state [35]. QInfer implements this criterion as the SMCUpdater.normalization record property, allowing for model selection and averaging to be performed in a straightforward manner.

Approximate Maximum-Likelihood Estimation
As opposed to the Bayesian approach, one may also consider maximum likelihood estimation (MLE), in which a model is estimated asx MLE := arg max x Pr(D|x). MLE can be approximated as the mean of an artificially tight posterior distribution obtained by performing Bayesian inference with a likelihood function Pr (D|x) related to the true likelihood by Pr (D|x) = (Pr(D|x)) γ (15) for a quality parameter γ > 1 [65]. Similarly, taking γ < 1 with appropriate resampling parameters allows the user to anneal updates [66], avoiding the dangers posed by strongly multimodal likelihood functions. In this way, taking γ < 1 is roughly analogous to the use of "reset rule" techniques employed in other filtering algorithms [53]. In QInfer, both cases are implemented by the class qinfer . MLEModel, which decorates another model in the manner of Section 4.2.

Likelihood-Free Estimation
For some models, explicitly calculating the likelihood function Pr(D|x) is intractable, but good approaches may exist for drawing new data sets consistent with a hypothesis. This is the case, for instance, if a quantum simulator is used in place of a classical algorithm, as recently proposed for learning in large quantum systems [29]. In the absence of an explicit likelihood function, Bayesian inference must be implemented in a likelihood-free manner, using hypothetical data sets consistent to form a approximate likelihood instead [67]. This introduces an estimation error which can be modeled in QInfer by using the qinfer . PoisonedModel class discussed in Section 4.4.

Simplified Estimation
For the frequency estimation and randomized benchmarking examples described in Section 3, QInfer provides functions to perform estimation using a "standard" updater loop, making it easy to load data from NumPy-, MATLABor CSV-formatted files.
Jupyter Integration Several QInfer classes, including qinfer . Model and qinfer . SMCUpdater, integrate with Jupyter Notebook to provide additional information formatted using HTML. Moreover, the qinfer . IPythonProgressBar class provides a progress bar as a Jupyter Notebook widget with a QuTiP-compatible interface, making it easy to report on performance testing progress.
MATLAB/Julia Interoperability Finally, QInfer functionality is also compatible with MATLAB 2016a and later, and with Julia (using the PyCall. jl package [68]), enabling integration both with legacy code and with new developments in scientific computing.

Conclusions
In this work, we have presented QInfer, our opensource library for statistical inference in quantum information processing. QInfer is useful for a range of different applications, and can be readily used for custom problems due to its modular and extensible design, addressing a pressing need in both quantum information theory and in experimental practice. Importantly, our library is also accessible, in part due to the extensive documentation that we provide (see ancillary files or docs.qinfer.org). In this way, QInfer supports the goal of reproducible research by providing open-source tools for data analysis in a clear and understandable manner.
Listing 8: Example of a custom FiniteOutcomeModel subclass implementing the multi-cos likelihood (16).

B Resampling
The purpose of this appendix is to offer a brief discussion of resampling for particle filters. In QInfer, the standard resampler is the one proposed by Liu and West [25]. We begin by motivating the development of resamplers by explaining the problem of impoverishment in a particle filter. We then describe resamplers by explaining that their effect is to produce a particle filter that approximates a smoothed version of the underlying probability distribution. Finally, we explain the Liu-West resampling algorithm.
Particle filters are intended to allow for the approximation of the expectation values of functions, but admit an ambiguity between using the locations and the weights to do so. Assuming that the particle locations are primarily responsible for representing the particle filtering approximation, and assuming that we want to approximate the expectation value of a function that is not pathological in some way, the number of particles then serves as a reasonable proxy for the quality of the particle filter. This follows from exactly the same argument as for Monte Carlo integration, as in this case, the particle locations can be seen to directly correspond to samples of an integrand. On the other hand, if either of these assumptions are violated, one cannot trust the numerical answers obtained using the particle filter method without an additional argument. Since the weights will in general become less even as Bayesian inference proceeds, we will rely on resampling to provide us with precisely such an argument.
More precisely, the purpose of resampling is to mitigate against the loss in numerical stability caused by having a large number of low-weight particles. If a particle has a small weight, we could neglect it from the computation of an expectation value without introducing a large difference in the result, such that it no longer contributes to the approximation quality of (5). That is, the particle's effectiveness at contributing to the stability of the algorithm decreases as its weight decreases. This observation then motivates using the effective sample size as a criterion to ensure the numerical stability of particle filtering [69].
In the case of roughly equal weights, the effective sample size n ess is roughly equal to the actual number of particles n. If the weights are distributed rather unevenly, however, n ess n. In the latter case, we say that the particle filter is impoverished. Notably, if a particle filter becomes extremely impoverished, it may not be feasible to effectively recover numerical stability, such that the minimum observed value of n ess serves as a diagnostic criteria for when the numerical approximations used by particle filtering have failed [39]. For this reason, QInfer will by default call the resampler when n ess falls below half of its initial value, and will warn if n ess ≤ 10 is ever observed.
Impoverishment is the result of choosing particle locations according to prior information alone; with near certainty, the posterior distribution will be tightly centered away from all initial particle locations such that some re-discretization will be needed to represent the final posterior. The goal of resampling is therefore to modify the choice of particle locations not based on the interpretation of new data, but rather to ensure that the particles are concentrated so as to accurately represent the probability distribution of interest. A resampling algorithm is then any algorithm designed to replace an input particle filter that may be impoverished with a new particle filter that is not impoverished but approximates the same probability distribution. Since such a procedure cannot exist in full generality, each resampling procedure works by enforcing that a particular set of invariants remains true before and after the resampling step.
Returning to the particular case of the Liu-West resampler [25] employed by default in QInfer, we note that the Liu-West algorithm is particularly simple to understand from the perspective of kernel density estimation [70]. We will not provide a complete algorithm here; such explicit algorithms can be found many places in the literature, including Algorithm 2 of Granade [52] and Algorithm 2.5 of Sanders [71]. Rather, we explain that the Liu-West algorithm acts by preserving the first two moments (the expectation and the variance) of its input distribution.
The Liu-West algorithm starts by using that, as in the celebrated method of kernel density estimation, each sample from a continuous distribution can be used to infer properties about the neighborhood around that sample, given moderate assumptions on the smoothness of a distribution. Thus, the particle filter might be used to define a direct approximation to the true probability distribution by first defining some function K, called a kernel, such that the estimated distribution is a good approximation of the true probability distribution. This then leaves open the question of what kernel function K should be chosen. Owing to its generality, choosing K to be a normal distribution works well under only mild assumptions, leaving us to choose the variance of the kernel. Specializing to the single-parameter case for simplicity of notation, we denote Gaussian kernels as K(x; σ 2 ), where σ 2 is the variance to be chosen. The multi-parameter case follows by replacing this variance by a covariance matrix.
The key insight at the core of the Liu-West algorithm is the observation that the posterior distribution will narrow over time, so that we can choose the variance of the kernel to narrow over time in proportion to the variance of the distribution being approximated. In particular, let a ∈ [0, 1] be a parameter and let h such that a 2 + h 2 = 1. Then choose the kernel to be K(x; h 2 V[x]). Since this alone would increase the variance of the distribution under approximation to 1 + h 2 , the Liu-West resampler instead draws new particle locations x from the distribution