`MCMC`#

class MCMC(filter_obj: FilterOD | SolutionOD, verbose: bool = True, *, walkers_per_dimension: int = 4, min_walkers: int = 32, walker_spread: float = 1.0, max_initial_redraws: int = 100, stretch_scale: float = 2.0, de_weights: tuple = (0.7, 0.3), burn_in_factor: float = 2.0, thin_factor: float = 0.5, tau_multiple: float = 50, min_acceptance: float = 0.2)#

Bases: FilterOD

Markov chain Monte Carlo sampling of a converged orbit determination.

A filter (LSB, SRIFB, LKF, SRIF, …) is first fitted to convergence. Its last pass fixes the reference trajectory, the prefit residuals \(y_k\), the measurement partials \(H_k\) and the STMs \(\Phi(t_k, t_0)\). About that reference the measurements are linear in the deviation \(x\) at the reference epoch,

\[y_k = A_k x + \epsilon_k, \qquad A_k = H_k \Phi(t_k, t_0), \qquad \epsilon_k \sim \mathcal{N}(0, R_k),\]

and the posterior of \(x\) is sampled with the affine-invariant ensemble sampler of emcee:

\[\log p(x \mid y) = -\frac{1}{2} \sum_k (y_k - A_k x)^T R_k^{-1} (y_k - A_k x) + \log p(x) + \mathrm{const}.\]

No trajectory is propagated and no measurement is recomputed while sampling: every step costs two small matrix products. The data enter through the normal matrix and vector

\[D = \sum_k A_k^T R_k^{-1} A_k, \qquad d = \sum_k A_k^T R_k^{-1} y_k.\]

Whitening

The walkers move in \(\theta\), with \(x = \hat{x} + L\theta\), where \(\hat{x} = M^{-1} b\) and \(L L^T = M^{-1}\) are the batch least-squares solution and covariance (\(M = \bar{P}_0^{-1} + D\), \(b = \bar{P}_0^{-1} \bar{x}_0 + d\)). With the filter’s Gaussian a-priori the log-posterior is then \(-\tfrac{1}{2}\theta^T\theta\), whatever the conditioning of \(M\). A user prior \(\log p\) replaces the Gaussian a-priori, and the log-posterior becomes

\[-\frac{1}{2} \theta^T L^T D L\, \theta + \theta^T L^T (d - D \hat{x}) + \log p(\hat{x} + L\theta - \bar{x}_0),\]

up to a constant. The prior is evaluated on the deviation from the a-priori state, so a Gaussian \(\log p(\delta) = -\tfrac{1}{2}\delta^T \bar{P}_0^{-1}\delta\) reproduces the default.

Sequences

For a MissionSequence every leg is mapped to one global vector \(z\): the first leg’s state at its epoch, followed by each block a later leg adds at its node, with its a-priori sigma. Leg \(j\) is \(x_j = T_j z\), with \(T_j = Q_j \Phi_{j-1}(t_{node}) T_{j-1}\) plus an identity on the blocks the node adds, so a measurement of leg \(j\) has \(A_k = H_k \Phi_j(t_k) T_j\). This is the batch of LSB, with every leg mapped back through the nodes, as long as no considered state is dropped at a node: a dropped one stays considered in \(z\), whereas the filter folds its spread into the next leg’s a priori and lets that leg’s data reduce it (ConsiderDroppedWarning).

Consider parameters

Considered states are not sampled. The walkers sample the estimated states with the considered ones at their reference, and the consider spread is added linearly as in the Schmidt batch: each particle draws \(c \sim \mathcal{N}(0, P_{cc})\) and is shifted by \(S_{xc} c\), with \(S_{xc} = -M_{xx}^{-1} M_{xc}\).

Solution

sample() returns a SolutionOD whose deviation and covariance histories are the chain mean and covariance (plus the consider part) mapped to every epoch with the STMs, and whose filter carries particles (the posterior samples mapped the same way), so SolutionOD.uq_metrics is filled. The samples at the reference epoch are also stored as samples, ready for plot_corner().

Parameters:

filter_objFilterOD or SolutionOD: A filter fitted to convergence with FilterOD.fit(), or the solution it returned. It is copied, not modified.
verbosebool, optional: Print the initialization banner. Defaults to True.
walkers_per_dimensionint, optional: Walkers per estimated state when sample() is not given n_walkers. The stretch move needs more than twice the dimension; four times keeps the ensemble from collapsing onto a subspace. Defaults to 4.
min_walkersint, optional: Fewest walkers when sample() is not given n_walkers. Defaults to 32.
walker_spreadfloat, optional: Spread of the initial walkers, in units of the filter’s 1-sigma. The walkers are drawn from the converged filter’s Gaussian scaled by this factor, so with 1 the chain starts close to stationarity on the linearized problem and the burn-in is short. Defaults to 1.0.
max_initial_redrawsint, optional: Times an initial walker is redrawn when the log-posterior is not finite at its position (a log_prior that rules out part of the Gaussian) before sample() gives up. Defaults to 100.
stretch_scalefloat, optional: Scale parameter \(a\) of the affine-invariant stretch move of moves="stretch", the default of Goodman & Weare (2010) and of emcee. Defaults to 2.0.
de_weightstuple of float, optional: Relative weights of the differential-evolution and snooker moves of moves="de"; emcee normalizes them. Defaults to (0.7, 0.3).
burn_in_factorfloat, optional: Burn-in, in multiples of the largest integrated autocorrelation time of the chain, and at most half the chain. Defaults to 2.0.
thin_factorfloat, optional: Thinning, in multiples of the smallest integrated autocorrelation time of the chain, and at least one step. Defaults to 0.5.
tau_multiplefloat, optional: Chain length, in largest integrated autocorrelation times, below which the autocorrelation estimate (and with it burn-in and thinning) is not trusted, and ShortChainWarning is issued. The tol default of emcee. Defaults to 50.
min_acceptancefloat, optional: Mean acceptance fraction below which the ensemble mixes poorly and AcceptanceFractionWarning is issued. There is no upper bound: the stretch move scales its steps with the spread of the ensemble, so a high acceptance does not mean short steps. On the whitened, near-Gaussian posterior of a few states it is about 0.5 or more. Defaults to 0.2.

Attributes:

labels: Name of each entry of the global state samples holds.

Methods

`fit`(args, *kwargs)	Not supported: the source filter is fitted, then `sample()` runs.
`sample`([n_steps, n_walkers, moves, ...])	Sample the linearized posterior and pack it as a solution.

Raises:

FilterTypeError: If filter_obj is not a filter (or the solution of one), or is itself an MCMC.
CovarianceAnalysisError: If the filter ran a covariance analysis: there are no residuals to sample from.
UnfittedFilterError: If the filter was never fitted.
FilterNotConvergedError: If the last fit did not converge: the linearization would not be about the solution.
ProcessNoiseError: If the filter has process noise, which the batch posterior of the reference-epoch state cannot represent.

MCMC#

`MCMC`#