Returns (XS, YS, ZS), with ZS of shape (Nx,Ny) and XS and YS of shape (Nx+1,Ny+1) giving the height and box corners, respectively, of the histogram estimate of the PDF from which pts are drawn. The bin widths and orientations are chosen optimally for the convergence of the histogram to the true PDF under a squared-error metric; the automatic binning is tuned to work particularly well for multivariate Gaussians.
Note: the first index of ZS varies with the X coordinate, while the second varies with the y coordinate. This is consistent with pp.pcolor(), but inconsistent with pp.imshow() and pp.contour().
| Parameters: |
|
|---|
Returns an estimate of the area within the given credible levels for the posterior represented by pts.
The algorithm uses a two-step process (see this document) so that the resulting credible areas will be unbiased.
| Parameters: |
|
|---|---|
| Returns: | If truth is None, areas, an array of the same shape as levels giving the credible areas; if truth is not None then (areas, searched_area, p_value). |
Returns a linearly interpolated quantile value.
| Parameters: |
|
|---|
Load data from the file, using header for column names.
| Parameters: |
|
|---|
Plots the 1D normalized empirical CDF for the given points. Additional arguments are passed to matplotlib’s plot.
Produces a chain plot of the mean values of each variable at each step. The chain should have the shape (nwalkers, nsteps, nvars), and the resulting grid of plots will be as close to square as possible, showing the mean walker position at each step.
| Parameters: |
|
|---|
Plots a single-figure representation of the chain evolution of the given chain. The figure shows the evolution of the mean of each coordinate of the ensemble, normalised to zero-mean, unit-standard deviation.
| Parameters: |
|
|---|
Plot probability interval contours estimated from a histogram PDF of the given points. The number of bins in each dimension is chosen optimally for minimizing the squared error with a Gaussian PDF. Additional arguments passed to pp.contour().
| Parameters: |
|
|---|
Plots the given probability interval contours, using a greedy selection algorithm. Additional arguments passed to pp.contour().
The algorithm uses a two-step process (see this document) so that the resulting credible areas will be unbiased.
| Parameters: |
|
|---|
Plots a histogram estimate of the posterior from which pts are drawn. Extra arguments are passed to pp.hist().
| Parameters: |
|
|---|
Plots a 2D histogram density of the given points.
| Parameters: |
|
|---|
Plot probability intervals corresponding to levels in 1D. Additional args are passed to pp.axvline(). The chosen levels are symmetric, in that they have equal probability mass outside the interval on each side.
| Parameters: |
|
|---|
Plots the a KDE estimate of the posterior from which pts are drawn. Extra arguments are passed to pp.plot().
| Parameters: |
|
|---|---|
| Returns: | (xs, ys), the coordinates of the plotted line. |
Plot a 2D KDE estimated posterior.
| Parameters: |
|
|---|
Like plot_emcee_chains() but for kombine.
Like plot_emcee_chains_one_fig() but for kombine.
Module containing functions to compute the autocorrelation function and estimate the associated autocorrelation length of series.
The estimate of the autocorrelation function is based on the method described at http://www.math.nyu.edu/faculty/goodman/software/acor/ and implemented in the associated acor C++ code, though the implementation details differ.
Returns the autocorrelation function of the given series. The function is normalised so that it is 1 at zero lag.
If series is an N-dimensional array, the ACF will be computed along axis and the result will have the same shape as series.
Returns an estimate of the autocorrelation length of the given series:
The estimate is the smallest \(L\) such that
In words: the ACL is estimated over a window that is at least \(M\) ACLs long, with the constraint that \(ML < N/2\).
Defined in this way, the ACL gives the reduction factor between the number of samples and the “effective” number of samples. In particular, the variance of the estimated mean of the series is given by
Returns nan if there is no such estimate possible (because the series is too short to fit \(2M\) ACLs).
For an N-dimensional array, returns an array of ACLs of the same shape as series, but with the dimension along axis removed.
Returns an array giving the ACL for each parameter in the given emcee chain.
| Parameters: |
|
|---|
Returns the Gelman-Rubin R convergence statistic applied to individual walkers’ trajectories in each parameter.
Returns an array of shape (Ntemp, Nparams) giving the estimated autocorrelation lengths for each parameter across each temperature of the parallel-tempered set of chains. If a particular ACL cannot be estimated, that element of the array will be None. See emcee_chain_autocorrelation_lengths() for a description of the optional arguments.
Returns a thinned, burned-in version of the emcee chain.
| Parameters: |
|
|---|
Returns a thinned, burned in version of the emcee parallel-tempered chains in ptchain, or None if it is not possible to estimate an ACL for some component of the chain.
Just like emcee_chain_autocorrelation_lengths() but for kombine.
Plots a grid of the autocorrelation function (post burnin) for each of the parameters in the given chain.
Returns an estimate of the WAIC from an emcee sampler’s lnlike (should be of shape (nwalkers, nsteps)). The WAIC is defined by
See Gelman, Hwang, and Vehtari (2013) for a motivation for this quantity in terms of an unbiased estimate of the expected log pointwise predictive density.
Bases: object
Runner object for an emcee sampler.
Return the estimate of the current chain’s autocorrelation lengths, using plotutils.autocorr.emcee_chain_autocorrelation_lengths().
Returns a chain with the first 1/6 of the samples removed. This chain corresponds to the samples that are then thinned to produce the thin_chain property. There is, of course, no guarantee that this chain is actually burned in.
Resets the stored sampler and internal state. If no pts argument is given, will use the last position of the sampler as the new starting position; otherwise, the new starting position will be given by pts.
Increase the thinning parameter by a factor of two, modifying the stored chain and lnprob states accordingly.
Run the associated sampler to produce nthinsteps worth of stored ensembles (i.e. the sampler will be run for nthinsteps*self.thin total steps).
Run the sampler, thinning as necessary, until neff effective ensembles are obtained. When savedir is not None, the sampler will be periodically saved to the given directory.
Save the state of the runner stored chain, lnprob, and thin parameter in the given directory. Three files will be created (approximately atomically, so the operation is nearly safe from interruption):
- chain.npy.bz2
- lnprob.npy.bz2
- thin.txt
Storing the current chain, lnprob and the thin parameter.
In addition, the file runner.pkl.bz2 will be created storing a pickled version of the runner object.
Bases: plotutils.runner.EnsembleSamplerRunner
Runner class for the PTSampler; basically behaves as EnsembleSamplerRunner, but tailored to the different chain shape and behaviour of emcee.PTSampler.
Return the estimate of the current chain’s autocorrelation lengths, using plotutils.autocorr.emcee_ptchain_autocorrelation_lengths().
Returns a chain with the first 1/6 of the samples removed. This chain corresponds to the samples that are then thinned to produce the thin_chain property. There is, of course, no guarantee that this chain is actually burned in.
Increase the thinning parameter by a factor of two, modifying the stored chain and lnprob states accordingly.
Run the sampler, thinning as necessary, until neff effective ensembles are obtained. When savedir is not None, the sampler will be periodically saved to the given directory.
Save the state of the runner stored chain, lnprob, lnlike, betas, and thin parameter in the given directory. Three files will be created (approximately atomically, so the operation is nearly safe from interruption):
- chain.npy.bz2
- lnprob.npy.bz2
- lnlike.npy.bz2
- beta_history.npy.bz2
- nswap.txt
- nswap_accepted.txt
- nprop.txt
- nprop_accepted.txt
- betas.txt
- time.txt
- thin.txt
Storing the current chain, lnprob, lnlike, betas (inverse temperatures), and the thin parameter.
In addition, the file runner.pkl.bz2 will be created storing a pickled version of the runner object.
Useful parameterizations for various commonly-found elements of models.
Returns the log of the Jacobian factor
for the bounded parameters p.
Returns the parameters associated with the values x that are bounded between low and high.
Returns the values, each bounded between low and high (one of these can be None) associated with the parameters p.
The parameterisation is
if both lower and upper limits are given, and
or
if only one limit is given.
| Parameters: |
|
|---|
Returns the log of the determinant of the Jacobian of the transformation that produces the covariance matrix from parameters x:
Returns a covariance matrix from the parameters x, which should be of shape ((N+1)*N/2,) for an (N,N) covariance matrix.
The parametersation is taken from the Stan sampler, and is as follows. The parameters x are the lower-triangluar elements of a matrix, \(y\). Let
Then the covariance matrix is
With this parameterization, there are no constraints on the components of x.
Given a covariance matrix, m, returns the parameters associated with it through the cov_matrix() function.
Returns the log of the Jacobian factor
for the parameters p.
Returns the parameters associated with the values x which should be increasing.
Returns the values for parameters p which are constrained to be always increasing.
The parameterisation is
Note that \(-\infty < p < \infty\).
A parameterisation of the roots of a real, ‘stable’ polynomial.
A stable polynomial has roots with a negative real part; it is the characteristic polynomial for a linear ODE with decaying solutions. The parameterisation provides a mapping from \(\mathbb{R}^n\) to the roots that is one-to-one; there are no root-permutation degeneracies. The log-Jacobian function produces a flat distribution on the real and imaginary (if any) parts of the roots.
| Parameters: |
|
|---|
The parameterisation uses a ‘bounded logit’ transformation to map ranges of reals to \(\pm \infty\):
The mapping of the roots proceeds as follows. First, discard the roots with strictly negative imaginary parts (these are the conjugates of corresponding roots with strictly positive imaginary parts, so we lose no information). Then, sort the remaining roots in order of decreasing imaginary part; if there are any strictly real roots, with imaginary part zero, sort these in decreasing order, from least negative to most negative.
Let \(a = -\left( r_\mathrm{max} - r_\mathrm{min} \right)\) and \(b = r_\mathrm{max}\). Then, proceeding by pairs of roots:
If the imaginary part of root i is greater than zero, then
\[\begin{split}p_i = \mathrm{logit}\left(\mathrm{im}\left(r_i\right); a, b\right) \\ p_{i+1} = \mathrm{logit}\left(\mathrm{re}\left(r_i\right); -r_\mathrm{max}, -r_\mathrm{min} \right) \end{split}\]Then set \(b = \mathrm{im} r_i\), and proceed to the next pair.
If the imaginary part of root i is zero, then
\[\begin{split}p_i = \mathrm{logit}\left( r_i + r_\mathrm{min}; a, b \right) \\ p_{i+1} = \mathrm{logit}\left( r_{i+1} + r_\mathrm{min}; a, r_i + r_\mathrm{min} \right)\end{split}\]then set \(b = r_{i+1} + r_\mathrm{min}\).
If the number of roots is odd, then the final root must be real, and is the smallest of all the real roots:
\[p_{N-1} = \mathrm{logit}\left( r_{N-1} + r_\mathrm{min}; a, \mathrm{min}(b, 0) \right).\]
Intuitively, you can imagine constructing the parameterisation using a line that sweeps down the imaginary axis, starting from \(i r_\mathrm{max}\); as it hits each complex root, it records the logit of that root’s imaginary and real parts (which must lie within certain bounds), and then sets the maximum bound for the imaginary part of the next root. Once all the complex roots have been parameterised, a line begins at \(-r_\mathrm{min}\) on the real axis and sweeps left; as it hits each real root, it records the logit of that root between the current bounds, and resets the bound on the maximum value of the next root to the current real root.
The reverse transformation begins by unpacking the first parameter value using the inverse logit transform; if this value is positive, then the root is complex, and the next parameter correponds to the real part. If the value is negative, then this a real root, and the following roots are also real (and negative). Bounds on the subsequent root values are set accordingly in either case for the next inverse logit transformation.
The parameterisation maps the allowed, sorted root space onto the entire \(\mathbb{R}^N\) real space in a one-to-one way.
Given N-1 parameters, p, returns N positive values that sum to one. The transformation comes from the Stan manual. Imagine a stick that begins with unit length; the parameters are logit-transformed fractions of the amount of the stick remaining that is broken off in each of N-1 steps to produce the N lengths.
Represents a one-dimensional Gaussian kernel density estimator conducted in log-space, according to the formula
The KDE occurs on the left-hand-side, but the reported distribution is \(dN/dx(x)\).
A contour calculator. Takes as input a set of points drawn from some distribution, and a set of credible levels. Then Contour.in_contours() computes for any set of points which ones like within the given credible contours and which are outside.
Returns the lower bound on the posterior estimate for inclusion in each of the contours.
Returns arrays indicating whether the given points are within each of the contours of the object.
| Parameters: | pts – A (Npts, Ndim) array of points to test against each contour. |
|---|---|
| Returns: | A sequence of the same shape as Contours.contours containing boolean arrays indicating which points lie within the contour. |