#  Space Oddity? A Statistical Formulation of Conjunction Assessment

S. Elkantassi 111Ph.D candidate, Institute of Mathematics, Ecole Polytechnique
Fédérale de Lausanne (EPFL), Station 8, 1015 Lausanne, Switzerland and A. C.
Davison222Professor, Institute of Mathematics, Ecole Polytechnique Fédérale de
Lausanne (EPFL), Station 8, 1015 Lausanne, Switzerland Institute of
Mathematics, Ecole Polytechnique Fédérale de Lausanne (EPFL), Lausanne, 1015,
Switzerland

###### Abstract

Satellite conjunctions involving ‘near misses’ of space objects are becoming
increasingly likely. One approach to risk analysis for them involves the
computation of the collision probability, but this has been regarded as having
some counter-intuitive properties  and its interpretation has been debated.
This paper formulates an approach to satellite conjunction based on a simple
statistical model and discusses inference on the miss distance between the two
objects, both when the relative velocity can be taken as known and when its
uncertainty must be taken into account. It is pointed out that the usual
collision probability estimate can be badly biased, but that highly accurate
inference on the miss distance is possible. The ideas are illustrated with
case studies and Monte Carlo results that show its excellent performance.

## Nomenclature

(Nomenclature entries should have the units identified)

$f(y;\vartheta)$ = |  probability density function for data $y$, depends on unknown parameter $\vartheta$
---|---
$f(\vartheta)$ = |  prior density for $\vartheta$ in Bayesian formulation of satellite conjunction problem
$y$ = | generic value of data vector
$y^{\mathrm{o}}$ = | observed value of data vector
$\Omega$ = |  dispersion matrix of $y$, considered to be known
$\vartheta$ = | column vector of parameters of statistical model
$\psi$ = | scalar interest parameter for statistical model, usually the unknown miss distance (m or km)
$\psi_{\rm min}$ = | combined hard-body radius HBR(m or km)
$\psi_{0}$ = |  safety threshold (m or km), with $\psi_{0}>\psi_{\rm min}$
$\lambda$ = | column vector of nuisance parameters for statistical model
${\rm Pr}(\cdot)$ = |  probability of event ‘$\cdot$’
$\operatorname*{E}(\cdot)$ = |  expectation with respect to distribution ${\rm Pr}(\cdot)$
${\rm Pr}_{0}(\cdot)$ = |  probability of event ‘$\cdot$’ computed under null hypothesis $\psi=\psi_{0}$
$p_{c}$ = | collision probability
$\widehat{p}_{c}$ = | estimated collision probability
$r(\psi)$ = | likelihood root, an approximate pivot
$q(\psi)$ = | correction to likelihood root
$r^{*}(\psi)$ = |  modified likelihood root $r(\psi)+\log\\{q(\psi)/r(\psi)\\}/r(\psi)$, an approximate pivot
$r^{*}_{B}(\psi)$ = |  Bayesian version of $r^{*}(\psi)$
$w(\psi)$ = | Wald statistic, an approximate pivot
$I_{1-2\alpha}$ = |  Two-sided $(1-2\alpha)$ confidence interval for unknown parameter $\psi$
$L_{\alpha},U_{\alpha}$ = |  lower and upper limits to an equi-tailed $(1-2\alpha)$ confidence interval
$\Phi(\cdot)$ = |  cumulative probability function for standard normal, ${{\mathcal{N}}}(0,1)$, distribution
$z_{p}$ = |  $p$ quantile of the standard normal distribution, satisfying $\Phi(z_{p})=p$ for $0<p<1$
$\mu$ = | true current position vector for secondary space object relative to first (m or km)
$\nu$ = | true current velocity vector for secondary space object relative to first (m/s or km/s)
$\|\cdot\|$ = |  Euclidean norm of the vector ‘$\cdot$’
$\beta$ = |  angle between $\mu$ and $\nu$ (radians)
$\theta_{1},\theta_{2},\phi_{1},\phi_{2}$ = | angles of spherical polar coordinates (radians)
$\xi$ = |  two-dimensional projection of $\mu$ in the direction of $\nu$ onto the encounter plane (m or km)
$x$ = |  two-dimensional projection of $y$ in the direction of $\nu$ onto the encounter plane (m or km)
$D$ = |  ${\rm diag}(d_{1}^{2},d_{2}^{2})$, variance matrix of $x$ (m2 or km2)
$I(\cdot)$ = |  indicator function for event ‘$\cdot$’
$L(\vartheta)$ = |  likelihood function, i.e., $f(y;\vartheta)$ regarded as a function of unknown parameter $\vartheta$
$\ell(\vartheta)$ = |  log-likelihood function, i.e., $\log f(y;\vartheta)$ regarded as a function of unknown parameter $\vartheta$
$\jmath(\vartheta)$ = |  observed information matrix, i.e., $-\partial^{2}\ell(\vartheta)/\partial\vartheta\partial\vartheta^{\mathrm{\scriptscriptstyle T}}$
$\jmath_{\lambda\lambda}(\vartheta)$ = |  $(\lambda,\lambda)$ block of observed information matrix, i.e., $-\partial^{2}\ell(\vartheta)/\partial\lambda\partial\lambda^{\mathrm{\scriptscriptstyle T}}$
$\jmath_{p}(\psi)$ = |  profile information, $|\jmath(\widehat{\vartheta}_{\psi})|/|\jmath_{\lambda\lambda}(\widehat{\vartheta}_{\psi})|$
$\varphi(\vartheta)$ = | constructed parameter used for higher-order approximation
Sub- and super-scripts
${\mathrm{o}}$ = |  indicates quantity evaluated for observed data $y=y^{\mathrm{o}}$
$\widehat{\vartheta}$ = |  indicates maximum likelihood estimator of parameter $\vartheta$
$\widehat{\lambda}_{\psi}$ = |  indicates maximum likelihood estimator of parameter $\lambda$ when parameter $\psi$ is fixed
$\widehat{\vartheta}_{\psi}$ = |  indicates maximum likelihood estimator of parameter $\vartheta$ when parameter $\psi$ is fixed, i.e., $\widehat{\vartheta}_{\psi}=(\psi,\widehat{\lambda}_{\psi})$
$x^{{\mathrm{\scriptscriptstyle T}}}$ = |  transpose of $x$

##  1 Introduction

The expansion of the aerospace industry and the increasing number of space
objects, especially in Low Earth Orbit where most spacecraft operate, means
that risk assessment and collision avoidance manoeuvres are vital to ensure
their safety. A great deal of effort has been put into conjunction assessment
for orbiting objects, generally by attempting to estimate the probability of
collision between them, as in [1, Chapter 5] or [2, Section 11.7]. This is
calculated at the time of the closest approach using the estimated position
and velocity vectors for the two objects and the associated error covariances
and, in a short-term conjunction, can be expressed as an integral of a two-
dimensional Gaussian probability density function over the collision cross-
sectional area. Although unavailable in explicit form, it can readily be
evaluated numerically [3, 4, 5, 6, 7, 8]. However both more precise and less
precise measurement reduce the collision probability, a ‘dilution’ property
that has been seen as paradoxical, and its interpretation is not seen as
clear-cut [9]. Another criterion for risk assessment is the closest approach,
or miss, distance, as a miss distance that is likely to be lower than a
specified safety threshold indicates a situation that requires action.  This
distance is estimated repeatedly over the approach to a conjunction and the
collision probability is updated accordingly. There are many algorithms to
compute the miss distance, often as the root of a polynomial equation [10, 11,
12] that depends on the relative path of the objects, which may be highly
nonlinear when they are far apart.

In this work, we formulate a statistical model for conjunction assessment that
resolves the apparent difficulties with the collision probability and suggest
that the miss distance is a more appropriate focus of interest. We discuss
inference for this distance based on standard likelihood theory [13, chapter
9], and also describe an improved theory that is both highly accurate and
should give results similar to a Bayesian formulation [14]. Our approach is
based on significance functions [15] and provides both point and interval
estimates for the miss distance, with the intervals containing the true miss
distance with a specified probability under the model. We can also test
whether the true distance is higher than the safety threshold, in order to
inform decisions about avoidance manoeuvres [16].

The paper is organized as follows. In Section 2 we formulate the conjunction
assessment problem in statistical terms and discuss the relationship between
the conjunction probability and miss distance.  We point out that the
conjunction probability tends to be biased downwards and elucidate the so-
called dilution paradox of the former. In Section 3, we introduce significance
functions and discuss how they provide calibrated frequentist inference,
present elements of modern likelihood inference and link them to the Bayesian
approach. In Section 4 we apply these ideas to the satellite conjunction
problem and in Section 5 study their properties in Monte Carlo studies and
case studies. Section 6 comments on possible extensions and limitations of the
approach.

##  2 Statistical modeling of conjunction assessment

###  2.1 Problem formulation

Modelling and analysis of the relative motion of two space objects has been
successfully applied to many space missions [17]. Algebraic models for
relative errors were proposed by Hill [18] and Clohessy and Wiltshire [19] and
developed in Chen et al. [1]. Many authors [3, 4, 20, 5, 7] have considered
the two objects as ellipsoids and attempted to estimate the probability that
they will collide. As mentioned above, this probability has been considered by
some to have paradoxical properties [9], but these evaporate when the problem
is formulated using a statistical model, as we shall see. Recall that a
parametric statistical model treats the available data $y$ as the realisation
of a random variable whose probability density function $f(y;\vartheta)$ is
determined by unknown parameters $\vartheta$, and the  usual goal is to use
$y$ for inference about the value of a scalar parameter
$\psi=\psi(\vartheta)$. Below we take $\psi$ to be the miss distance and
suppose that a collision occurs if the two objects pass within a minimum
distance of each other, namely the combined hard-body radius, $\psi_{\rm
min}$.  We assume that the conjunction is sufficiently close that the relative
motion can be taken to be linear; Hall [21] summarises the literature on
probability calculation for nonlinear and repeated conjunctions.

Suppose initially that the current positions $\mu_{s1}$ and $\mu_{s2}$ and
velocities $\nu_{s1}$ and $\nu_{s2}$ of the two objects are known; all these
are $3\times 1$ vectors. Define $\mu=\mu_{s2}-\mu_{s1}$ and
$\nu=\nu_{s2}-\nu_{s1}$, so the second object is considered relative to an
origin at the first. In this frame of reference and under linear relative
motion, the second object traverses the line $\mu+t\nu$, where
$t\in\mathbb{R}$. Its distance from the origin,
$(\mu+t\nu)^{\mathrm{\scriptscriptstyle T}}(\mu+t\nu)$, is minimised by
choosing $t=-\nu^{\mathrm{\scriptscriptstyle
T}}\mu/\nu^{\mathrm{\scriptscriptstyle T}}\nu$, at which point the minimum
squared distance is $\psi^{2}=\mu^{\mathrm{\scriptscriptstyle
T}}\mu-\left(\mu^{\mathrm{\scriptscriptstyle
T}}\nu\right)^{2}/\nu^{\mathrm{\scriptscriptstyle T}}\nu$. In terms of
spherical polar coordinates we have

$\displaystyle\mu$ $\displaystyle=$
$\displaystyle\|\mu\|\left(\sin\theta_{1}\cos\phi_{1},\sin\theta_{1}\sin\phi_{1},\cos\theta_{1}\right)^{\mathrm{\scriptscriptstyle
T}},$ (1) $\displaystyle\nu$ $\displaystyle=$
$\displaystyle\|\nu\|\left(\sin\theta_{2}\cos\phi_{2},\sin\theta_{2}\sin\phi_{2},\cos\theta_{2}\right)^{\mathrm{\scriptscriptstyle
T}},$ (2)

where $\|\cdot\|$ denotes the Euclidean norm and
$0\leq\theta_{1},\theta_{2}\leq\pi$ and $0\leq\phi_{1},\phi_{2}<2\pi$ are the
polar and azimuthal angles for $\mu,\nu$. The minimum distance between the two
objects, the miss distance, is

$\psi=\|\mu\|(1-\cos^{2}\beta)^{1/2}=\|\mu\||\sin\beta|,$ (3)

where $\beta$, the angle between the location and velocity vectors $\mu$ and
$\nu$, satisfies

$\cos\beta=\sin\theta_{1}\cos\phi_{1}\sin\theta_{2}\cos\phi_{2}+\sin\theta_{1}\sin\phi_{1}\sin\theta_{2}\sin\phi_{2}+\cos\theta_{1}\cos\theta_{2}.$
(4)

When $\beta=0$ we distinguish two cases: if $\mu^{\mathrm{\scriptscriptstyle
T}}\nu<0$ the second object will pass through the origin, leading to a
collision, whereas if $\mu^{\mathrm{\scriptscriptstyle T}}\nu>0$ the second
object is heading away from the origin, so its current position is the closest
it will come to the first object. More generally, $\psi<\|\mu\|$ only if
$\cos\beta<0$, i.e., $\pi/2<\beta<3\pi/2$. For $\mu$ and $\nu$ to be collinear
but pointing in opposite directions we need $\phi_{2}=\pi+\phi_{1}$ and
$\theta_{2}=\pi-\theta_{1}$, and then
$\cos\beta=\cos(\theta_{1}+\theta_{2})=-1$, so $\beta=\pi$ and hence $\psi=0$,
as expected. To lighten the notation below we write
$\lambda=(\theta_{1},\phi_{1},\|\nu\|,\theta_{2},\phi_{2})$, and
$\vartheta=(\psi,\lambda)$.

In the above deterministic setting the collision probability $p\equiv
p(\vartheta)$ takes two possible values,

$p(\vartheta)=\begin{cases}0,&\psi>\psi_{\rm min},\\\
1,&0\leq\psi\leq\psi_{\rm min},\end{cases}$ (5)

and a decision-maker can make an ideal decision. In reality, of course, both
$\mu$ and $\nu$ are observed with error, and we follow the literature and
assume that the available observations on the positions and velocities of the
two objects have a multivariate normal distribution with known covariance
matrix. If so, then the vector $y$ containing the observed position and
velocity of the second object relative to the first has a six-dimensional
normal distribution, and we suppose that this has mean vector
$\eta(\vartheta)=\left\\{\mu(\psi,\lambda)^{\mathrm{\scriptscriptstyle
T}},\nu(\lambda)^{\mathrm{\scriptscriptstyle
T}}\right\\}^{{\mathrm{\scriptscriptstyle T}}}$ given by equations (1)–(4) and
known $6\times 6$ covariance matrix $\Omega^{-1}$, whose inverse $\Omega$ is
known as the dispersion matrix.  In reparametrizing the six-dimensional model
to $(\psi,\lambda),$ it is important that the new parameters are variation
independent. If the allowable values of $\psi$ were to depend on $\lambda$,
this would introduce irregularities into the statistical model and the usual
asymptotic theory would not apply. Such problems arise when the miss distance
$\psi$ is expressed in Cartesian coordinates.

###  2.2  Two-dimensional case

A simplified version of the problem treats the relative velocity vector $\nu$
as known. In this case the last three components of $y$ and of
$\eta(\vartheta)$ corresponding to the relative velocity can be dropped, and
only the $3\times 3$ corner of $\Omega^{-1}$ pertaining to the relative
position need be retained. The conjunction can then be visualised in the
encounter plane, also called the conjunction plane, which is normal to $\nu$
and passes through the origin. The observed position $y$ of the second object,
its true position $\mu$, and the density of $y$, can be projected orthogonally
in the direction of $\nu$ into this plane, leading to respective projected
observed and true positions $x$ and $\xi$ and bivariate normal density
$f(x;\xi)$ for $x$ with unknown mean vector $\xi$ and known diagonal
covariance matrix $D={\rm diag}(d_{1}^{2},d_{2}^{2})$; see Appendix A. As
shown in the left-hand panel of Figure 1, $\xi$ and $x$ are two-dimensional
vectors: the true point at which the second object will pass through the
encounter plane, $\xi$, is at a distance $\psi=\|\xi\|$ from the origin, and
the collision probability is zero unless $\psi\leq\psi_{\rm min}$. The right-
hand panel of Figure 1 shows how the collision probability is computed: the
density is assumed to be centered at $x$ and the estimator $\widehat{p}_{c}$
is the integral of this density over the disk of radius $\psi_{\rm min}$
around the origin. Thus we can write $\widehat{p}_{c}=p_{c}(x)$, where

$p_{c}(\xi)=\int_{\\{x^{\prime}:\|x^{\prime}\|\leq\psi_{\rm
min}\\}}f(x^{\prime};\xi)\,\mathrm{d}{x^{\prime}}.$ (6)

One might regard $p_{c}(x)$ as an estimate of $p_{c}(\xi)$, since the unknown
$\xi$ is replaced with the known $x$ in computing the integral. Seen through
the lens of the statistical model, however, it is not obvious why $p_{c}(\xi)$
is of interest, as it is the probability that the noisy observation $x$ will
appear to pass inside the hard-body radius, rather than the probability (6)
that the satellite itself will do so. Moreover the symmetry of $f(x;\xi)$
around $\xi$ shown in the left-hand panel of Figure 1 implies that there is
more probability outside the circle of radius $\psi$ than inside, so $x$ will
tend to be further away from the origin than $\xi$. Another way to see this is
to note that the Euclidean norm of $x$, $\|x\|$, satisfies

$\operatorname*{E}(\|x\|^{2})=\operatorname*{E}(x_{1}^{2}+x_{2}^{2})=\xi_{1}^{2}+\xi_{2}^{2}+d_{1}^{2}+d_{2}^{2}=\psi^{2}\left\\{1+(d_{1}^{2}+d_{2}^{2})/\psi^{2}\right\\},$

and hence
$\operatorname*{E}\\{\|x\|\\}\approx\psi\left\\{1+(d_{1}^{2}+d_{2}^{2})/\psi^{2}\right\\}^{1/2}$
for large $\psi$ , i.e., the mean length of $x$ exceeds $\psi$ by an amount
that depends on $\xi$ and $D$. Hence $\widehat{p}_{c}=p_{c}(x)$ will tend to
be smaller than $p_{c}(\xi)$. Figure 2 shows the effect of this in a test case
described more fully in Section 5.4. Each panel takes
$\xi=(11.84,-1.36)^{\mathrm{\scriptscriptstyle T}}$m, so $\psi=11.92$m, and
compares the values of $p_{c}(\xi)$ with boxplots showing the values of
$p_{c}(x)$ for $20,000$ values of $x$ generated as bivariate normal,
${{\mathcal{N}}}_{2}\\{\xi,c^{2}{\rm diag}(d^{2}_{1},d^{2}_{2})\\}$, with
different values of $c^{2}$. The probability of collision is computed using a
contour integral transformation [7], which is then approximated numerically
using the trapezoidal rule. The left-hand panel for $c^{2}=10^{-2},$ shows
that although $p_{c}(\xi)\approx 10^{-3}$, around 25% of the values of
$p_{c}(x)$ are less than $10^{-4}$. Overall the values of $p_{c}(x)$ can vary
over many orders of magnitude and can be much lower than $p_{c}(\xi)$. Figure
3, which shows the values of $\xi$ and $10^{3}$ simulated values of $x$ for
three values of $c^{2}=10^{-2},1,2$, illustrates how $p_{c}(x)$ varies greatly
when $\psi_{\rm min}=5$m and when $\psi_{\rm min}=20$m. Clearly the bias is
also present when the velocity is not treated as known.

Figure 1: Statistical formulation of satellite conjunction in the encounter
plane. The primary object is at the origin $O$ and the solid circle around it
has radius $\psi_{\rm min}$, the combined hard-body radius. The true position
at which the second object will cross the encounter plane, $\xi$, can be
expressed in terms of the angular coordinates $\psi$ and $\lambda$. The noisy
observation of the second object will cross the encounter plane at $x$. Left
panel: the ellipses indicate the true density of $x$, with mean at $\xi$.
Right panel: the ellipses indicate the assumed density of $x$ when computing
the conjunction probability estimate $\widehat{p}_{c}$.

Figure 2: Behaviour of $\widehat{p}_{c}=p_{c}(x)$ for hard-body radius
$\psi_{\rm min}=5$m (left) and 20m (right) in Case Study C, with $\psi\approx
11.92$m, as a function of the variance of
$x\sim{{\mathcal{N}}}_{2}\\{\xi,c^{2}{\rm diag}(d^{2}_{1},d^{2}_{2})\\}$. For
each $c^{2}$ in $0.005,0.01,\ldots,10$ we computed $p_{c}(\xi)$ (blue
segments) and $p_{c}(x)$ (boxplots) for 20,000 simulated values of $x$. The
average values of $p_{c}(x)$ are shown by the red segments. In the left-hand
panel $\psi>\psi_{\min}$, so $p_{c}(\xi)\to 0$ as $c\to 0$, whereas in the
right-hand panel $\psi<\psi_{\min}$, so $p_{c}(\xi)\to 1$ as $c\to 0$.

Figure 3: Scatter plots showing 1000 observed locations in the conjunction
plane sampled as
$x\sim\mathcal{N}_{2}\left\\{\xi,c^{2}\text{diag}(d_{1}^{2},d_{2}^{2})\right\\}$
with $c^{2}=10^{-2},1$, and $2$ (top to bottom). The true position of the
second object, $\xi$, is shown by the red point and the observed positions,
$x$, are shown by blue points whose sizes and shades depend on the value of
$p_{c}(x)$. The shaded grey disks show the hard-body radius, $\psi_{\rm
min}=5$ m on the left and $\psi_{\rm min}=20$ m on the right.

Despite the comments above, one reason to use $\widehat{p}_{c}$ is that it is
a Bayesian estimator. If a prior density $f(\xi)$ for $\xi$ is placed on the
encounter plane, then the posterior probability that $\xi$ lies within the
hard-body radius is given in terms of the posterior density of $\xi$, i.e.,

$f(\xi\mid x)={f(x;\xi)f(\xi)\over\int f(x;\xi)f(\xi)\,\mathrm{d}{\xi}},$

by

${\rm Pr}(\psi\leq\psi_{\rm min}\mid x)={\rm Pr}(\|\xi\|\leq\psi_{\rm min}\mid
x)=\int_{\\{\xi:\|\xi\|\leq\psi_{\rm min}\\}}f(\xi\mid x)\,\mathrm{d}{\xi},$
(7)

and if $f(\xi)$ is constant and $f(x;\xi)$ is bivariate normal, then
$f(\xi\mid x)=f(x;\xi)$ and (7) equals $\widehat{p}_{c}=p_{c}(x)$. This
explains why $\widehat{p}_{c}$ is a plausible estimator of (5), but does not
alter its downward bias. Moreover a constant prior for $\xi$ is improper, with
the undesirable property that the ratio of the probability inside any disk
around the origin is zero relative to the probability outside that disk, thus
expressing a prior belief that $\xi$ is infinitely far from the origin, i.e.,
the second object will traverse the encounter plane infinitely far from the
first. This gives an alternative explanation of the behaviour seen in Figure
2.

###  2.3 Conjunction probability or miss distance?

Two further points are immediately clear from the above very simple
statistical model for satellite conjunction.

First, the conventional target of inference, the true collision probability
(5),  depends on the unknown relative position and velocity of the two
objects, but not on the covariance matrix $\Omega^{-1}$ or the data $y$. The
collision probability is generally estimated by $\widehat{p}_{c}$, which
depends on both $\Omega^{-1}$ and $y$. The ‘paradoxical’ behaviour whereby
$\widehat{p}_{c}$ is very tiny when the data have a very large variance and
then increases when that variance decreases is the behaviour of an estimator,
not of a parameter of the model. The estimator $\widehat{p}_{c}$ depends on
$\Omega^{-1}$, and as the variance decreases we expect that either
$\widehat{p}_{c}\to 1$, if a collision will occur, or $\widehat{p}_{c}\to 0$,
otherwise; in both cases $\widehat{p}_{c}\to p(\vartheta)$, as we should
expect when the data become noiseless. Thus probability dilution is the
natural behaviour of an estimator in response to changes in the variability of
the underlying data, not a probability paradox.

Second, the fact that $p(\vartheta)$ takes just two values whereas the miss
distance $\psi$ takes values in a continuum suggests that $\psi$ is a better
overall target of inference. For example, the maximum likelihood estimator of
$p(\vartheta)$ is $I(\widehat{\psi}\leq\psi_{\rm min})$, where $I(\cdot)$ and
$\widehat{\psi}$ denote the indicator function and the maximum likelihood
estimator of $\psi$, and clearly $\widehat{\psi}$ is more informative.
Moreover in a Bayesian framework the unknown parameters are regarded as random
variables and the posterior probability of collision given the data is
obtained by integrating over the posterior density of $\psi$ given $y$; c.f.
(7) when $\nu$ is known. In this setting also the collision probability is
based on the available knowledge about the miss distance $\psi$, which is the
more fundamental quantity. In Section 3.4 we discuss Bayesian inference in
more detail, and in Section 4 we show how inference on $\psi$ provides a
significance probability with an interpretation akin to that of
$\widehat{p}_{c}$.

Motivated by these considerations, we turn to inference on $\psi$ based on the
observed value of $y$. If the data suggest that $\psi$ is lower than a safety
threshold $\psi_{0}$, possibly with $\psi_{0}>\psi_{\rm min}$ for a safety
margin, then action to avert a collision should be considered. Our goal below
is therefore inference on the unknown miss distance $\psi$, allowing for the
fact that $\lambda$ is also unknown and must be estimated; in the six-
dimensional case we can write the relative distance vector using (1) but with
$\|\mu\|$ replaced by $\psi/|\sin\beta|$, for $\beta\neq\pi$.  In statistical
terms the scalar $\psi$ is the primary object of inference, the so-called
interest parameter, whereas the $5\times 1$ vector $\lambda$ of nuisance
parameters, while essential for realistic modelling, is of only secondary
concern. In the special case with known velocity vector $\nu$ and considering
only the encounter plane, the miss distance $\psi$ remains the parameter of
interest and its interpretation is unchanged, but the nuisance parameter
$\lambda$ is scalar. This simplifies the problem but the statistical issue
remains the same.

We return to this model in Section 4 after outlining elements of the theory of
inference.

##  3 Inference

###  3.1 Calibration,  decisions and significance functions

Statistical inference involves statements about the properties of a
probability distribution that is assumed to have given rise to observed data
$y^{\mathrm{o}}$.  In the simplest setting the distribution depends only on a
scalar parameter $\psi$, which is the focus of interest, and the likelihood
function $L(\psi)=f(y^{\mathrm{o}};\psi)$ is used to compare the plausibility
of different values of $\psi$ as explanations for $y^{\mathrm{o}}$. The best-
fitting model is provided by the maximum likelihood estimate based on $y^{o}$,
$\widehat{\psi}^{\mathrm{o}}$, and the relative likelihood function
$L(\psi)/L(\widehat{\psi}^{o})$, which has maximum value 1, allows values of
$\psi$ to be compared. A ‘pure likelihood’ approach [e.g., 22, Chapter 3]
treats any $\psi$ for which $L(\psi)\geq u\,L(\widehat{\psi}^{o})$ as
plausible, but with $u$ chosen essentially arbitrarily. In practice further
information is typically  required in order to choose $c$.

Bayesian inference treats $\psi$ as a random variable and  is based on a
density $f(\psi)$ that weights values of $\psi$ according to their
plausibility prior to seeing the data. This is updated in light of the
observed data $y^{\mathrm{o}}$ using Bayes’ formula, resulting in the
posterior distribution function

${\rm Pr}(\psi\leq\psi_{0}\mid
y^{\mathrm{o}})={\int_{-\infty}^{\psi_{0}}L(\psi)f(\psi)\mathrm{d}{\psi}\over\int_{-\infty}^{\infty}L(\psi)f(\psi)\mathrm{d}{\psi}}.$
(8)

Clearly this calculation depends on the prior density $f(\psi)$; if this is
badly chosen then (8) may have poor properties when used repeatedly.  The most
obvious choice of prior in the satellite conjunction setting is uniform on the
position of the secondary space object, but this has the undesirable
properties mentioned in Section 2.2.

When the losses due to possible evasive actions can be specified, the data can
be used to choose the action that minimises the expected posterior loss. For
the simplest possible formulation in the collision avoidance context, suppose
that $\psi$ represents the unknown miss distance, that the two actions $a=0$
and $a=1$ correspond to ‘do nothing’ and ‘take evasive action’, and that the
loss $l_{ae}$ when action $a$ is taken and event $e$ occurs is as given in
Table 2; the loss in doing nothing in the case of no collision is zero. If
$\psi\leq\psi_{0}$ results in a collision, then the posterior expected loss
due to taking action $a\in\\{0,1\\}$ is

$l_{a0}\,{\rm Pr}(\psi>\psi_{0}\mid y^{\mathrm{o}})+l_{a1}\,{\rm
Pr}(\psi\leq\psi_{0}\mid y^{\mathrm{o}})=l_{a0}(1-p)+l_{a1}p,$

say, which is minimised by doing nothing if $l_{01}p<l_{10}(1-p)+l_{11}p$.
Equivalently, evasive action should be taken if ${\rm Pr}(\psi\leq\psi_{0}\mid
y^{\mathrm{o}})\geq l_{10}/(l_{10}+l_{01}-l_{11})$. Thus if the losses are
known explicitly, they provide a threshold for action, and the use of a
decision rule such as ‘take evasive action if the posterior probability of
collision exceeds $10^{-4}$’ corresponds to an implicit ratio of losses. This
decision setup could be made more realistic, and in any case would only be
regarded as guidance in a practical setting; our point is that it provides a
rational basis for considering action when (8) exceeds a threshold, and
explicitly links that threshold to potential losses.

Table 2: Basic decision analysis for satellite conjunction, with losses $l_{ae}$ corresponding to action $a$ and event $e$. Action | Event
---|---
| No collision | Collision
Do nothing, $a=0$ | $l_{00}=0$ | $l_{01}$
Evasive action, $a=1$ | $l_{10}$ | $l_{11}$

Bayesian inference has certain advantages, but other approaches to calibrating
the likelihood are often preferred. Repeated sampling inference invokes
hypothetical repetition of the random experiment that is presumed to have led
to the observed data [23, pp. 33–38]. In the simplest case
$\widehat{\psi}^{\mathrm{o}}$ is regarded as a realization of a random
variable $\widehat{\psi}$ that has a normal distribution,
${{\mathcal{N}}}(\psi,\lambda^{2})$, under repeated sampling, with $\lambda$
known. This implies that

${\rm
Pr}(\widehat{\psi}\leq\widehat{\psi}^{o};\psi)=\Phi\left\\{(\widehat{\psi}^{\mathrm{o}}-\psi)/\lambda\right\\},$
(9)

where $\Phi$ denotes the standard normal cumulative distribution function. The
significance function (9), also called the confidence distribution or P-value
function [24, 25, 26], is then used for inference on $\psi$. For example, the
null hypothesis that $\psi=\psi_{0}$ can be tested against the alternative
that $\psi>\psi_{0}$ by computing the significance probability

$p_{\rm obs}={\rm
Pr}\left(\widehat{\psi}\geq\widehat{\psi}^{\mathrm{o}};\psi_{0}\right),$ (10)

small values of which are regarded as evidence against the null hypothesis in
favour of the alternative; here $p_{\rm
obs}=1-\Phi\left\\{(\widehat{\psi}^{\mathrm{o}}-\psi_{0})/\lambda\right\\}$.
Likewise a two-sided $(1-2\alpha)\times 100\%$ confidence interval
${\mathcal{I}}_{1-2\alpha}$ for the value of $\psi$ underlying the data, the
so-called ‘true value’, has as its upper and lower limits $U_{\alpha}$ and
$L_{\alpha}$ the solutions to the equations

${\rm Pr}(\widehat{\psi}\leq\widehat{\psi}^{o};U_{\alpha})=\alpha,\quad{\rm
Pr}(\widehat{\psi}\leq\widehat{\psi}^{o};L_{\alpha})=1-\alpha,$

and this yields
${\mathcal{I}}_{1-2\alpha}=(L_{\alpha},U_{\alpha})=(\widehat{\psi}^{\mathrm{o}}-\lambda
z_{1-\alpha},\widehat{\psi}^{\mathrm{o}}-\lambda z_{\alpha})$, where $z_{p}$,
the $p$ quantile of the standard normal distribution, satisfies
$\Phi(z_{p})=p$ for $0<p<1$. The limits of ${\mathcal{I}}_{1-2\alpha}$
simplify to the familiar $\widehat{\psi}^{\mathrm{o}}\pm\lambda z_{1-\alpha}$
on recalling that $z_{1-\alpha}=-z_{\alpha}>0$. In this ideal case the
inferences are perfectly calibrated: under repeated sampling with
$\psi=\psi_{0}$ the significance probability $p_{\rm obs}$ has an exact
uniform distribution and ${\mathcal{I}}_{1-2\alpha}$ contains $\psi_{0}$ with
probability exactly $1-2\alpha$, for any $\alpha\in(0,0.5)$. When
$\psi=\psi_{0}$, therefore, there is a probability $p_{\rm obs}$ that a
decision to reject this hypothesis in favor of the alternative based on a
significance probability $p_{\rm obs}$ will be incorrect.

Significance functions are decreasing in the parameter and have the properties
of survivor functions, so the confidence density based on differentiating (9),

$-{\partial{\rm
Pr}(\widehat{\psi}\leq\widehat{\psi}^{o};\psi)\over\partial\psi},$ (11)

has the formal properties of a probability density function for $\psi$. The
confidence density and its integral, the confidence distribution, can be
regarded as frequentist summaries of the information about $\psi$ based on the
observed data; see Schweder and Hjort [26]. Unlike the posterior density
obtained from (8), no prior information is involved, and despite the use of
the word ‘density’ for (11), $\psi$ is regarded as an unknown constant.

###  3.2 Approximate inference

In Section 3.1, $\widehat{\psi}$ was assumed to have an exact normal
distribution, but the argument there often applies approximately when some
measure of precision, typically the sample size $n$, becomes large. Under mild
regularity conditions the maximum likelihood estimator $\widehat{\psi}$ has an
approximate normal distribution centered at the true parameter $\psi$ with
variance $1/\jmath(\widehat{\psi})$, where $\jmath(\psi)=-\partial^{2}\log
L(\psi)/\partial\psi^{2}$ is the observed information. This allows calibration
of $\psi$ using the Wald statistic

$w(\psi)=\jmath(\widehat{\psi})^{1/2}(\widehat{\psi}-\psi),$

corresponding to a quadratic approximation to the log likelihood function, or
using the likelihood root

$r(\psi)={\rm sign}(\widehat{\psi}-\psi)\\{2\log L(\widehat{\psi})-2\log
L(\psi)\\}^{1/2}.$

Both $w(\psi)$ and $r(\psi)$ are approximate pivots: they are functions of the
data and parameter and have approximate standard normal distributions  under
repeated sampling if $\psi$ equals its true value. The approximations
introduce so-called first-order error, of size $O(n^{-1/2})$, in (9). For more
details, see, for example, Chapter 9 of Cox and Hinkley [13].

To illustrate these approximations, consider a single observation from the
Rayleigh distribution, which arises as the length $y$ of a bivariate normal
vector whose components are independent $N(0,\psi^{2})$ variables. The
probability density function for $y$ is

$f(y;\psi)={y\over\psi^{2}}e^{-y^{2}/(2\psi^{2})},\quad y>0,\quad\psi>0,$

and one can readily check that $\widehat{\psi}=y/\surd{2}$,
$\jmath(\widehat{\psi})=4/\widehat{\psi}^{2}$ and

$w(\psi)=2(1-\psi/\widehat{\psi}),\quad r(\psi)={\rm
sign}(\widehat{\psi}-\psi)[2\\{2\log(\psi/\widehat{\psi})+(\widehat{\psi}/\psi)^{2}-1\\}]^{1/2}.$

The left-hand panel of Fig 4 shows the significance functions
$\Phi\left\\{w^{\mathrm{o}}(\psi)\right\\}$ and
$\Phi\left\\{r^{\mathrm{o}}(\psi)\right\\}$ when $y^{\mathrm{o}}=\sqrt{2}$, so
$\widehat{\psi}$ has observed value $\widehat{\psi}^{\mathrm{o}}=1$. The
functions are decreasing in $\psi$ and the maximum likelihood estimate and the
limits of the 90% confidence interval are the values of $\psi$ for which the
functions equal $0.5$ and $0.05,0.95$, respectively. Evidence against the null
hypothesis $\psi=\psi_{0}=0.3$ in favor of the alternative $\psi>\psi_{0}$ is
found from the P-value (10) given by the intersections of the significance
functions with the vertical line $\psi=0.3$; for the Wald statistic we obtain
$1-\Phi\left\\{w^{\mathrm{o}}(\psi_{0})\right\\}=0.0808$ and for the
likelihood root we obtain
$1-\Phi\left\\{r^{\mathrm{o}}(\psi_{0})\right\\}=4.33\times 10^{-5}$.  Thus
these pivots give quite different evidence about $\psi$. The right-hand panel
of Fig 4, which shows how pivots themselves depend on $\psi$, makes it easier
to read off confidence intervals; using $r^{\mathrm{o}}(\psi)$ to obtain one-
sided intervals $(L_{\alpha},\infty)$, for example, we see that $\psi_{0}$
lies below the intervals for $\alpha=10^{-2}$, $10^{-3}$ and $10^{-4}$, and
just inside that for $\alpha=10^{-5}$. A nested set of confidence intervals
such as this is more informative than a single significance probability based
on $r^{\mathrm{o}}(\psi_{0})$.

Figure 4: Left: Significance functions based on likelihood root
$r^{\mathrm{o}}(\psi)$ (solid black),Wald statistic $w(\psi)$ (dotted blue),
exact (wide cyan) and modified likelihood root $r^{*{\mathrm{o}}}(\psi)$ (red
dashes). The dashed vertical line shows $\psi_{0}=0.3$, and the horizontal
lines correspond to probabilities $10^{-6},10^{-5}$, $10^{-4}$, $10^{-3}$,
$10^{-2}$, 0.025, 0.05, 0.5, 0.95, 0.975, 0.999, 0.9999, 0.99999 and 0.999999.
Right: functions from (a) transformed to the standard normal scale.

In the Rayleigh example $\widehat{\psi}/\psi$ has a known distribution, so the
left-hand side of (9) provides an exact significance function. The sample size
here is $n=1$, so it is no surprise that the exact function differs greatly
from the large-sample approximations based on $w^{\mathrm{o}}(\psi)$ or
$r^{\mathrm{o}}(\psi)$; the exact significance probability for testing
$\psi_{0}=3$ is $1.49\times 10^{-5}$. Remarkably, however, an approximation
that treats the modified likelihood root

$r^{*}(\psi)=r(\psi)+{1\over r(\psi)}\log\left\\{{q(\psi)\over
r(\psi)}\right\\}$ (12)

as standard normal is essentially exact; its significance probability is
$1.55\times 10^{-5}$, giving a relative error of $3.9\%$. The approximation
based on $r^{*}(\psi)$ is said to be third-order accurate, i.e., its error is
$O(n^{-3/2})$, and moreover this error is relative, so high accuracy can be
expected even when the true probabilities are small. Here $n=1$, so the
approximation yields quasi-perfect calibration of different values of $\psi$.
In this particular example $q(\psi)=1-\widehat{\psi}^{2}/\psi^{2}$, but in
general it is uglier; see Appendix B.

###  3.3 Several parameters

Models with a single scalar parameter are unusual: other parameters are
usually present and must be estimated. Despite this, near-exact inferences
remain possible in many cases, as we now outline.

Suppose now that the density function $f(y;\vartheta)$ depends on an unknown
$d$-dimensional parameter $\vartheta$, which comprises a scalar interest
parameter $\psi$ and a nuisance parameter $\lambda$; $\psi$ and $\lambda$ are
supposed to be variation independent. As $\lambda$ is unknown, it must be
replaced by an estimate, and we see from (9) that this would introduce errors.

In this situation the maximum likelihood estimate
$\widehat{\vartheta}^{\mathrm{o}}=(\widehat{\psi}^{\mathrm{o}},\widehat{\lambda}^{\mathrm{o}})$
maximises the log-likelihood $\ell(\vartheta)=\log
f(y^{\mathrm{o}};\vartheta)$ with respect to $\vartheta$, and the partial
maximum likelihood estimate
$\widehat{\vartheta}^{\mathrm{o}}_{\psi}=(\psi,\widehat{\lambda}_{\psi}^{\mathrm{o}})$
maximizes $\ell(\vartheta)$ with respect to $\lambda$ for fixed $\psi$. The
large-sample properties of the maximum likelihood estimator
$\widehat{\vartheta}$ under repeated sampling are well-established [13,
Chapter 9]: as the sample size $n\to\infty$ and under mild regularity
conditions, $\widehat{\vartheta}$ has an approximate $d$-dimensional normal
distribution with mean the true parameter $\vartheta$ and covariance matrix
$\jmath(\widehat{\vartheta})^{-1}$, where
$\jmath(\vartheta)=-\partial^{2}\ell(\vartheta)/\partial\vartheta\partial\vartheta^{\mathrm{\scriptscriptstyle
T}}$ is the $d\times d$ observed information matrix and
$\vartheta^{\mathrm{\scriptscriptstyle T}}$ is the transpose of the $d\times
1$ vector $\vartheta$. Under these conditions the error committed by replacing
parameters in (9) by their estimates is $O(n^{-1/2})$, giving first-order
approximations, and the same error is committed by treating the likelihood
root

$r(\psi)={\rm
sign}(\widehat{\psi}-\psi)\left[2\left\\{\ell(\widehat{\vartheta})-\ell(\widehat{\vartheta}_{\psi})\right\\}\right]^{1/2}$
(13)

or the Wald statistic

$w(\psi)=\jmath_{\mathrm{p}}(\widehat{\psi})^{1/2}(\widehat{\psi}-\psi),$ (14)

as standard normal; here

$\jmath_{\mathrm{p}}({\psi})=\dfrac{\left|\jmath\left(\widehat{\vartheta}_{\psi}\right)\right|}{\left|\jmath_{\lambda\lambda}\left(\widehat{\vartheta}_{\psi}\right)\right|},$

$|\cdot|$ indicates the determinant, and
$\displaystyle\jmath_{\lambda\lambda}(\vartheta)$ is the
$\displaystyle(\lambda,\lambda)$ corner of $\jmath(\vartheta)$. The
improvement to third-order accuracy seen in the scalar case is again given by
(12), now using (13) and

$q(\psi)=\dfrac{\left|\varphi(\widehat{\vartheta})-\varphi(\widehat{\vartheta_{\psi}})~{}~{}\varphi_{\lambda}(\widehat{\vartheta_{\psi}})\right|}{\left|\varphi_{\vartheta}(\widehat{\vartheta})\right|}\dfrac{\left|\jmath(\widehat{\vartheta})\right|^{1/2}}{\left|\jmath_{\lambda\lambda}(\widehat{\vartheta}_{\psi})\right|^{1/2}}.$
(15)

Here $\varphi(\vartheta)$ is a $d\times 1$ constructed parameter that can be
viewed as a directional derivative of the log-likelihood, and
$\varphi_{\vartheta}(\vartheta)=\partial\varphi(\vartheta)/\partial\vartheta^{\mathrm{\scriptscriptstyle
T}}$ and
$\varphi_{\lambda}(\vartheta)=\partial\varphi(\vartheta)/\partial\lambda^{\mathrm{\scriptscriptstyle
T}}$ are respectively $d\times d$ and $d\times(d-1)$ matrices; see Appendix B.

In applications, expressions (13)–(15) are evaluated at the observed data
$y^{\mathrm{o}}$ and the corresponding estimates
$\widehat{\vartheta}^{\mathrm{o}}$ and $\widehat{\vartheta}^{o}_{\psi}$,
leading to the approximate significance function

${\rm
Pr}(\widehat{\psi}\leq\widehat{\psi}^{o};\psi)=\Phi\left\\{r^{*{\mathrm{o}}}(\psi)\right\\}$
(16)

which can be used in the same way as (9). The relative error committed in
using this expression, $O(n^{-3/2})$, is appreciably lower than that using a
normal approximation to the likelihood root (13).

###  3.4 Bayesian approximation

In the setting with several parameters, expression (8) generalizes to the
marginal posterior distribution of the scalar interest parameter $\psi$, i.e.,

${\rm Pr}(\psi\leq\psi_{0}\mid y^{\mathrm{o}})={\int_{-\infty}^{\psi_{0}}\int
L(\psi,\lambda)f(\psi,\lambda)\,\mathrm{d}{\lambda}\,\mathrm{d}{\psi}\over\int_{-\infty}^{\infty}\int
L(\psi,\lambda)f(\psi,\lambda)\,\mathrm{d}{\lambda}\,\mathrm{d}{\psi}},$ (17)

where $f(\psi,\lambda)$ is the prior density. On approximating both integrals
using Laplace’s method, it turns out that [14, Section 8.7]

${\rm Pr}(\psi\leq\psi_{0}\mid
y^{\mathrm{o}})=\Phi\left\\{r_{B}^{*{\mathrm{o}}}(\psi_{0})\right\\}\left\\{1+O(n^{-1})\right\\},$
(18)

where $r_{B}^{*{\mathrm{o}}}(\psi)$ is defined by equation (12) but with
$q^{\mathrm{o}}(\psi)$ replaced by

$q_{B}(\psi)=\left.{\mathrm{d}{\ell(\widehat{\vartheta}_{\psi})}\over\mathrm{d}{\psi}}\right.\times{|\jmath_{\lambda\lambda}(\widehat{\vartheta}_{\psi})|^{1/2}\over|\jmath(\widehat{\vartheta})|^{1/2}}\times{f(\widehat{\vartheta})\over
f(\widehat{\vartheta}_{\psi})}$ (19)

evaluated at the observed data $y=y^{\mathrm{o}}$ and the corresponding
estimates $\widehat{\vartheta}=\widehat{\vartheta}^{\mathrm{o}}$ and
$\widehat{\vartheta}_{\psi}=\widehat{\vartheta}^{o}_{\psi_{0}}$. Thus there is
a close parallel between frequentist and the Bayesian approximations and, with
an appropriate choice of the prior, the probabilities that the corresponding
confidence intervals contain $\psi_{0}$ differ by only $O(n^{-1})$. The
Jeffreys prior gives inferences invariant to 1–1 transformation of $\vartheta$
and  thus is often regarded as a natural choice, though it is criticised by
Fraser et al. [27]. In the conjunction problem this prior has the undesirable
property mentioned at the end of Section 2.2 of attributing zero probability
to any sphere around the origin, relative to outside that region.

In a similar context Davison and Sartori [28] compare the performances of
$r_{B}^{*{\mathrm{o}}}(\psi)$ and $r^{*{\mathrm{o}}}(\psi)$ and find that the
former performs rather worse. The Bayesian confidence intervals are slightly
shorter and tend to contain the true parameter less often,  and we would
expect this in the present setting also, owing to the downward bias of
$\widehat{p}_{c}$. Thus we do not implement Bayesian approximations; rather we
find it reassuring that Bayesian and frequentist inferences can be approached
in the same way. If reliable prior information was available on each likely
conjunction, then a Bayesian approach would be indicated.

Further details of the Bayesian approximation can be found in Appendix C.

##  4 Inference for conjunction assessment

###  4.1 Likelihood

In this section we apply the standard theory described in Section 3.3 to the
statistical model described in Section 2. After dropping irrelevant additive
constants, the log-likelihood function is

$\ell(\vartheta)=-{\textstyle{1\over
2}}\\{y-\eta(\vartheta)\\}^{{\mathrm{\scriptscriptstyle
T}}}\Omega\\{y-\eta(\vartheta)\\},$ (20)

where $\vartheta=(\psi,\lambda)$ contains the miss distance $\psi$ for which
inference is required and the five-dimensional nuisance parameter vector
$\lambda=(\theta_{1},\phi_{1},\|\nu\|,\theta_{2},\phi_{2})$. The maximum
likelihood estimator $\widehat{\vartheta}$ based on the observed relative
distance and velocity contained in the $6\times 1$ vector $y$ satisfies
$\eta(\widehat{\vartheta})=y$ and the observed information matrix is

$\jmath(\widehat{\vartheta})=\left.{\partial\eta^{\mathrm{\scriptscriptstyle
T}}(\vartheta)\over\partial\vartheta}\Omega{\partial\eta(\vartheta)\over\partial\vartheta^{\mathrm{\scriptscriptstyle
T}}}\right|_{\vartheta=\widehat{\vartheta}},$ (21)

where ${\partial\eta^{\mathrm{\scriptscriptstyle
T}}(\vartheta)/\partial\vartheta}$ is a $6\times 6$ matrix. Hence
$\ell(\widehat{\vartheta})=0$, $\widehat{\vartheta}$ is simply a
transformation of $y$, and
$\widehat{\vartheta}_{\psi}=(\psi,\widehat{\lambda}_{\psi})$ minimises the
weighted sum of squares in (20) for fixed $\psi$. These quantities allow
inference on $\psi$ based on the likelihood root (13) and the Wald statistic
(14), while the more accurate modified likelihood root (12) also requires
(15). The normal model is a curved exponential family [29, Section 5.2] in
terms of $\vartheta$, so we can take $\varphi(\vartheta)=\eta(\vartheta)$, the
computation of $r^{*}(\psi)$ only involves $\eta(\vartheta)$ and its
derivatives, and expression (15) simplifies to

$q(\psi)=|y-\eta(\widehat{\vartheta}_{\psi})~{}~{}\eta_{\lambda}(\widehat{\vartheta}_{\psi})|\times|\Omega|^{1/2}\times|\jmath_{\lambda\lambda}(\widehat{\vartheta}_{\psi})|^{-1/2}.$
(22)

Appendices B , C and D contain the detailed calculations for this, for the
Bayesian approximation with the Jeffreys prior
$|\eta_{\vartheta}(\vartheta)|$, and implementation details.

Here there are six parameters and a single six-dimensional observation $y$, so
the sample size is $n=1$, and it appears that we cannot expect large-sample
approximations to apply. However the covariance matrix for an average of $n$
independent observations would be $(n\Omega)^{-1}$, so a large sample size $n$
is equivalent to a small variance for the observations or equivalently large
$\Omega$, which is the correct gauge of accuracy.

As we saw in Section 2.2, the model simplifies greatly when the velocity
vector $\nu$ is known. In this case $y$ and $\eta(\vartheta)$ are replaced by
the projections of the observed position vector and the true position of the
second object into the encounter plane, i.e., $x=(x_{1},x_{2})$ and
$\xi=(\xi_{1},\xi_{2})=(\psi\cos\lambda,\psi\sin\lambda)$, and as the variance
matrix of $x$ can be diagonalised, the log likelihood becomes (see Figure 1)

$\ell(\psi,\lambda)=-{1\over 2}\left\\{{(x_{1}-\psi\cos\lambda)^{2}\over
d_{1}^{2}}+{(x_{2}-\psi\sin\lambda)^{2}\over
d_{2}^{2}}\right\\},\quad\psi>0,0\leq\lambda<2\pi,$ (23)

The approach advocated by Carpenter [30] is related to the discussion above,
but the confidence statements therein are based on the marginal quantiles of
the miss distance distribution rather than on likelihood theory. That approach
does not allow for the uncertainty about the nuisance parameter and appears to
be equivalent to basing inference on the Wald statistic, which can perform
very poorly in nonlinear settings.

###  4.2  Evidence and decisions

In Section 3.1 we argued that significance functions such as (10) and (16)
allow inference on the true miss distance $\psi$, either by constructing
confidence intervals or as an assessment of the evidence that $\psi$ equals
some particular value $\psi_{0}$. In the present context $\psi_{0}$ might be a
safety threshold, and then one approach to inference on $\psi$ is to test the
null hypothesis $H_{0}:\psi=\psi_{0}$ against the alternative hypothesis
$H_{+}:\psi>\psi_{0}$, with evasive action to be considered if (16) exceeds
some threshold ${\varepsilon}$, i.e., $H_{0}$ cannot be rejected at level
${\varepsilon}$. In practice ${\varepsilon}$ is often taken to be $10^{-4}$.
If the significance probability is correctly calibrated and the true miss
distance is $\psi_{0}$, then the false positive probability, that of
considering action unnecessarily, would be $1-{\varepsilon}$, whatever
${\varepsilon}$ is chosen. The choice of $H_{+}$ as alternative hypothesis
ensures that $p_{\rm obs}$ is small when the estimated conjunction probability
$\widehat{p}_{c}$ and the Bayesian posterior probability ${\rm
Pr}(\psi\leq\psi_{0}\mid y)$ would also be small, despite their different
interpretations and properties. Hejduk et al. [31] argue that the null
hypothesis $\psi\leq\psi_{0}$ is unnatural, since it implies that the ‘null’
situation is to anticipate a collision, but it appears more important to us
that small values of $p_{\rm obs}$ correspond to small conjunction
probabilities, which is ensured by the above setup. Our approach is supported
by regarding hypothesis testing as attempting ‘proof by stochastic
contradiction’: the null hypothesis represents an assumption that is regarded
as absurd (disproved) when the corresponding significance level is
sufficiently small. If so, it makes sense to take $\psi=\psi_{0}$ as the null
hypothesis, as we hope that this will be contradicted by the data and no
evasive action need be considered.

The interpretation of the threshold ${\varepsilon}$ in terms of an elementary
decision analysis was described in Section 3.1, and since the purpose of
conjunction analysis is to assist decision-making, consideration of losses
seems an appropriate basis for choosing ${\varepsilon}$. Although our earlier
discussion suggested considering evasive action when the posterior probability
${\rm Pr}(\psi\leq\psi_{0}\mid y)$ exceeds ${\varepsilon}$, it seems better to
replace ${\rm Pr}(\psi\leq\psi_{0}\mid y)$ by a significance probability
$p_{\rm obs}$ computed as $\Phi\left\\{r^{{\mathrm{o}}}(\psi_{0})\right\\}$ or
$\Phi\left\\{r^{*{\mathrm{o}}}(\psi_{0})\right\\}$, which are approximately
uniformly distributed under $H_{0}$.

The abuse of hypothesis tests has been much discussed [e.g., 32], and it is
often suggested that they be systematically replaced by confidence intervals.
Plotting $\Phi\left\\{r^{*{\mathrm{o}}}(\psi)\right\\}$ as a function of
$\psi$, as in Figure 4, allows the construction of confidence intervals for
$\psi$, the form of which expresses the uncertainty and suggests what power is
available: a narrow interval corresponds to more precise estimation and hence
higher power for rejecting hypotheses such as $H_{0}$ above. Unlike a two-
sided $(1-2\alpha)$ confidence interval $[L_{\alpha},U_{\alpha}]$, one-sided
intervals such as $[L_{\alpha},+\infty)$ or $[0,U_{\alpha})$ give no
information about the accuracy of the estimate, so a two-sided interval
provides a better basis for risk assessment. A one-sided confidence interval
$(L_{\varepsilon},\infty)$ that does not contain $\psi_{0}$ leads to the same
decision as observing $p_{\rm obs}<{\varepsilon}$, and it seems partly a
matter of taste which is preferred: computing $p_{\rm obs}$ alone is quicker
but is less informative than a plot of
$\Phi\left\\{r^{*{\mathrm{o}}}(\psi)\right\\}$. In navigating the literature
it is enlightening to recognize that hypothesis testing has a variety of
distinct uses [33], one of which is to flag situations that merit more
detailed scrutiny. In conjunction analysis one might therefore plot
significance functions only when $p_{\rm obs}>{\varepsilon}$, thereby focusing
on those conjunctions requiring careful consideration.

##  5 Numerical results

###  5.1 General

Below we investigate the accuracy of the normal approximations to the Wald
statistic, the likelihood root and the modified likelihood root in two
scenarios. We do so in terms of one-sided error rates for confidence intervals
$(L_{\alpha},U_{\alpha})$ for the true miss distance $\psi_{0}$ and use ${\rm
Pr}_{0}$ to indicate probability computed when $\psi=\psi_{0}$. An ideal two-
sided equi-tailed confidence interval with coverage probability $1-2\alpha$
should satisfy ${\rm Pr}_{0}(L_{\alpha}\leq\psi_{0}\leq U_{\alpha})=1-2\alpha$
and have one-sided left-tail and right-tail error rates ${\rm
Pr}_{0}(\psi_{0}<L_{\alpha})$ and ${\rm Pr}_{0}(U_{\alpha}<\psi_{0})$ both
equal to $\alpha$, for any $\alpha\in(0,0.5)$. Departures from this will
indicate deficiencies of the confidence intervals and the corresponding tests,
whereas close agreement will indicate that the inferences are well-calibrated.
As the pivots are decreasing in $\psi$ and should ideally have standard normal
distributions, inaccurate left-tail error rates correspond to the pivot
departing from normality in its upper tail. In the present setting, accurate
left-tail error rates are most important, since they correspond to well-
calibrated significance probabilities and confidence intervals of form
$(L_{\alpha},\infty)$.

The form of the covariance matrix used in practice depends on the type of the
conjunction [1, Chapter 5]. In a short-term conjunction, uncertainty on the
velocity is negligible compared to uncertainty on the position. In a long-term
conjunction, the motion is nonlinear and the computations are more involved.
In both cases, the quality of risk assessment depends heavily on the
covariance matrix. Below we suppose that the error covariance matrix for the
relative distance and velocity of the second satellite relative to the first
is given by

$\Omega^{-1}=\left[\begin{array}[]{cc}{P}_{1}&{P}_{12}\\\
{P}_{12}&{P}_{2}\end{array}\right],$

where $P_{1}$, $P_{2}$, and $P_{12}$ are the position, the velocity and the
cross-correlation covariance matrices,  of units $\mathrm{m}^{2}$,
$\mathrm{m}^{2}\mathrm{~{}s}^{-2}$ and $\mathrm{m}^{2}\mathrm{~{}s}^{-1}$
respectively. The six eigenvalues of $\Omega^{-1}$ are difficult to interpret
physically, and can vary greatly.

In our first two case studies, we assume that $P_{12}=0_{3\times 3}$, and
choose $P_{1}=\tau\sigma^{2}I_{3}$ and $P_{2}=\sigma^{2}I_{3}$. This choice
implies that the standard deviation of position errors along each axis
direction is $\sqrt{\tau}\sigma$ (km) and the standard deviation of velocity
errors is $\sigma$ (km/s). Uncertainty on the position is typically larger
than that on the velocity, and then $\tau>1$.

###  5.2 Case study A

The relative quantities and spherical coordinates of two satellites are as
given in Table 3. The relative distance is around 102 km and the relative
speed around 11.7 km/s, the value of ${\sigma}^{2}$ varies from $10^{-3}$ km2
to $2$ km2, and that of $\tau$ varies from 1 to 3.

Table 3: Conjunction geometry of case studies: A: a simulated example; B, the U.S. and Russian satellite collision event; and C an event with high probability of conjunction. Variable | Case study
---|---
| A | B | C
Miss distance (m) | $10^{3}\times 35.267$ | 698.011 | 11.917
$\Delta X$ (m) | $-10^{3}\times 100$ |  $-$258.909 |  $-$7.678
$\Delta Y$ (m) | $-10^{3}\times 20$ |  $-$635.813 |  $-$9.152
$\Delta Z$ (m) | 0 | 126.229 | 0.564
$\Delta V_{x}$ (km/s) | 10 | 10.580 | 9.926
$\Delta V_{y}$ (km/s) | 6 |  $-$3.733 |  $-$9.653
$\Delta V_{z}$ (km/s) | 1 | 3.126 |  $-$4.110
$\theta_{1}$ | 1.570 | 1.389 | 1.618
$\theta_{2}$ | 1.485 | 1.299 | 1.860
$\phi_{1}$ |  $-$2.944 | 1.957 |  $-$2.269
$\phi_{2}$ |  $-$2.944 |  $-$0.339 |  $-$ 0.772

Table 4 shows the error rates for the Wald statistic, the likelihood root $r$
and the modified likelihood root $r^{*}$ based on $10^{4}$ datasets simulated
for various combinations of values of $\sigma$ and $\tau$. For very small
$\sigma^{2}$ all three sets of error rates are close to the nominal values,
but problems with the Wald statistic and to a lesser extent the likelihood
root start to appear when $\sigma^{2}\geq 10^{-1}$, with the left-tail error
systematically too high and the right-tail error systematically too low. The
modified root behaves much better overall, though its right-tail error also
rises as $\sigma^{2}$ increases. As mentioned in Section 5.1, tests of $H_{+}$
require accuracy in the left tail, so right-tail error is less important.

These remarks are confirmed by the Gaussian QQ-plots of simulated values of
the three quantities in Figure 5. If the distribution is exactly Gaussian,
then the confidence intervals are exactly calibrated, so a departure from the
line of unit slope through the origin implies a lack of calibration. For small
$\sigma^{2}$, all three statistics have standard normal distributions and give
comparable results, but for larger $\sigma^{2}$, the Wald statistic and the
likelihood root are shifted to the right and right-skewed, more strikingly for
larger values of $\tau$.  This reflects the upward bias of the estimated
distance, discussed in Section 2.2, and explains the asymmetric error rates in
Table 4, with lower probabilities for the right than for the left. The
asymmetry increases with larger uncertainties on the relative distance and
velocity and with smaller nominal error rates.

The third column of Figure 5 shows that the modified likelihood root $r^{*}$
corrects the departure from normality in the upper tail even for
$\sigma^{2}=5$, and its error rates are closer to the nominal rates in all
cases considered. For $\sigma^{2}>2$ and for $1\%$ nominal levels, the Wald
statistic and the likelihood root show extreme overcoverage on the right and
undercoverage on the left; although the modified likelihood root provides a
considerable improvement, its right-tail error is somewhat smaller than the
nominal value.

Table 4:  Left and right error rates ($\%$) for two-sided nominal $10\%,5\%,$
and $1\%$ confidence intervals for the true parameter $\psi_{0}$ of Case Study
A, estimated from $10^{4}$ Monte Carlo samples. The standard errors (SE)
appear in the last line.

Left tail $(\%)$ Right tail $(\%)$ Uncertainty Statistic 5 2.5 0.5 5 2.5 0.5
$(\sigma^{2},\tau)=(10^{-3},1)$ $w$ $5.22$ $2.58$ $0.54$ $5.03$ $2.49$ $0.49$
$r$ $5.20$ $2.56$ $0.53$ $5.04$ $2.51$ $0.50$ $r^{*}$ $5.15$ $2.54$ $0.53$
$5.13$ $2.54$ $0.51$ $(\sigma^{2},\tau)=(10^{-3},2)$ $w$ $5.43$ $2.61$ $0.43$
$4.96$ $2.61$ $0.59$ $r$ $5.43$ $2.61$ $0.43$ $4.98$ $2.61$ $0.61$ $r^{*}$
$5.40$ $2.58$ $0.43$ $5.02$ $2.62$ $0.62$ $(\sigma^{2},\tau)=(10^{-3},3)$ $w$
$4.93$ $2.48$ $0.53$ $4.83$ $2.46$ $0.54$ $r$ $4.91$ $2.47$ $0.52$ $4.88$
$2.46$ $0.55$ $r^{*}$ $4.85$ $2.40$ $0.50$ $4.90$ $2.49$ $0.55$
$(\sigma^{2},\tau)=(10^{-1},1)$ $w$ 5.45 2.88 0.64 4.15 2.09 0.37 $r$ 5.26
2.78 0.59 4.27 2.15 0.37 $r^{*}$ 4.89 2.53 0.54 4.65 2.37 0.42
$(\sigma^{2},\tau)=(10^{-1},2)$ $w$ $5.62$ $2.98$ $0.67$ $4.60$ $2.28$ $0.50$
$r$ $5.48$ $2.75$ $0.61$ $4.70$ $2.32$ $0.54$ $r^{*}$ $5.06$ $2.48$ $0.54$
$5.06$ $2.56$ $0.64$ $(\sigma^{2},\tau)=(10^{-1},3)$ $w$ $5.66$ $2.86$ $0.65$
$4.49$ $2.26$ $0.49$ $r$ $5.54$ $2.70$ $0.63$ $4.58$ $2.35$ $0.50$ $r^{*}$
$5.06$ $2.50$ $0.59$ $4.92$ $2.57$ $0.53$ $(\sigma^{2},\tau)=(1,1)$ $w$ $7.21$
$4.10$ $1.04$ $3.18$ $1.37$ $0.16$ $r$ $6.45$ $3.36$ $0.68$ $3.35$ $1.53$
$0.21$ $r^{*}$ $5.39$ $2.64$ $0.55$ $4.56$ $2.37$ $0.31$
$(\sigma^{2},\tau)=(1,2)$ $w$ $6.58$ $3.72$ $0.92$ $3.64$ $1.50$ $0.24$ $r$
$5.92$ $3.16$ $0.57$ $3.80$ $1.62$ $0.29$ $r^{*}$ $4.81$ $2.40$ $0.41$ $5.32$
$2.59$ $0.48$ $(\sigma^{2},\tau)=(1,3)$ $w$ $6.40$ $3.50$ $0.96$ $3.73$ $1.64$
$0.25$ $r$ $5.74$ $3.03$ $0.64$ $3.94$ $1.76$ $0.25$ $r^{*}$ $4.74$ $2.49$
$0.54$ $5.19$ $2.72$ $0.46$ $(\sigma^{2},\tau)=(2,1)$ $w$ $7.78$ $4.57$ $1.41$
$2.12$ $0.69$ $0.03$ $r$ $6.52$ $3.56$ $0.78$ $2.27$ $0.75$ $0.03$ $r^{*}$
$5.19$ $2.64$ $0.58$ $4.46$ $1.85$ $0.10$ $(\sigma^{2},\tau)=(2,2)$ $w$ $7.94$
$4.59$ $1.33$ $2.30$ $0.71$ $0.05$ $r$ $6.77$ $3.49$ $0.71$ $2.43$ $0.88$
$0.05$ $r^{*}$ $5.18$ $2.60$ $0.54$ $4.62$ $2.09$ $0.18$
$(\sigma^{2},\tau)=(2,3)$ $w$ $7.67$ $4.07$ $1.15$ $2.19$ $0.73$ $0$ $r$
$6.49$ $3.24$ $0.65$ $2.39$ $0.86$ $0$ $r^{*}$ $4.71$ $2.34$ $0.49$ $4.48$
$1.98$ $0.11$ SE 0.22 0.16 0.07 0.22 0.16 0.07

Figure 5: Normal QQ-plots of $w(\psi^{o})$ (left), $r(\psi^{o})$ (middle) and
$r^{*}(\psi^{o})$ (right) based on $10^{4}$ Monte Carlo sample quantiles, with
$(\sigma^{2},\tau)$ equal to $(10^{-3},1),(1,1),(5,2),(5,5)$ (top to bottom).

###  5.3 Case study B

Our second example is the 10 February 2009 collision of the U.S. operational
communications satellite Iridium 33 and the decommissioned Russian
communications satellite Cosmos $2251$, whose relative configuration is given
in Table 3. Figure 6 shows the significance functions for this conjunction.
With a safety threshold of 20m, corresponding to $\psi_{0}=0.02$ (the dashed
vertical line in both panels), the significance functions all suggest that
evasive action is essential except that the probability of collision is
$1.14687\times 10^{-5}$ smaller than typically-employed remediation threshold
of $10^{-4}$, meaning that a conjunction remediation action is not needed. The
Wald statistic and the likelihood root are indistinguishable, but for the
modified likelihood root is shifted slightly to the left, reducing the upward
bias and slightly weakening any evidence that the true miss distance is
greater than $\psi_{0}$. The significance probabilities for testing
$\psi=\psi_{0}=20$ m against $H_{+}:\psi>\psi_{0}$ are $1.2\times 10^{-3}$ for
for the Wald statistic and the likelihood root, and $7.2\times 10^{-3}$ for
the modified likelihood root, so the same conclusions would be drawn for both
approaches if the threshold is ${\varepsilon}=10^{-4}$, but clearly there is
potential for different conclusions if we consider different uncertainties,
and other values of ${\varepsilon}$ and $\psi_{0}$. The posterior probability
${\rm Pr}(\psi\leq\psi_{0}\mid y^{\mathrm{o}})$ would be similar to the
significance probability for the modified likelihood root.

Figure 6: Summaries for US and Russian satellite collision event. Left:
significance functions based on likelihood root $r^{\mathrm{o}}(\psi)$ (solid
black),Wald statistic $w(\psi)$ (dotted blue), and modified likelihood root
$r^{*{\mathrm{o}}}(\psi)$ (red dashes). Right: the likelihood root
$r^{\mathrm{o}}(\psi)$ (solid black),Wald statistic $w(\psi)$ (dotted blue),
and modified likelihood root $r^{*{\mathrm{o}}}(\psi)$ themselves.

Table 6 gives left and right error rates for different values of $\sigma^{2}$
and $\tau$ for this conjunction. Its first row corresponds to the case where
the standard deviation of the position error along each axis is $10\mathrm{m}$
and the standard deviation of the velocity errors is
$10\mathrm{m}/\mathrm{s}$, giving $(\sigma,\tau)=(10^{-2},1)$. In lower rows
we first increase uncertainty on the position while keeping uncertainty on the
velocity fixed by increasing $\tau$, and then increase both velocity and
position errors by increasing $\sigma^{2}$. For these simulations the true
miss distance and relative velocity are taken to have $\psi=698$ m and
$\|\nu\|=11.648\times 10^{3}$ m/s, so the errors we consider vary from $10$ to
$200$ for both velocity (m/s) and distance (m).

The error rates for the Wald statistic and the likelihood root are almost
identical, implying that the corresponding pivots are indistinguishable. For
small velocity errors, with $\sigma^{2}<10^{-3}$, there is no significant
difference in the error rates for the three statistics. For larger velocity
variance, although the overall error rates found by summing the left and the
right error rates equal the nominal values, left error always dominates the
sum. In these cases the modified likelihood root is more symmetric and shows
fewer extreme values, especially for large $\tau$, and thus interval estimates
based on $r^{*}$ are more reliable.

Table 5:  Left and right error rates ($\%$) for two-sided nominal $10\%,5\%,$
and $1\%$ confidence intervals for the parameter $\psi$ of the U.S. and
Russian satellite collision event, Case Study B, based on $10^{4}$ Monte Carlo
replications. The standard errors (SE) appear in the last line.

Left tail $(\%)$ Right tail $(\%)$ Uncertainty Statistic 5 2.5 0.5 5 2.5 0.5
$(\sigma^{2},\tau)=(10^{-4},1)$ $w$ $5.06$ $2.54$ $0.42$ $5.00$ $2.53$ $0.64$
$r$ $5.06$ $2.54$ $0.42$ $5.00$ $2.53$ $0.64$ $r^{*}$ $4.98$ $2.46$ $0.41$
$5.05$ $2.62$ $0.65$ $(\sigma^{2},\tau)=(10^{-4},2)$ $w$ $5.27$ $2.73$ $0.60$
$4.66$ $2.29$ $0.52$ $r$ $5.27$ $2.73$ $0.60$ $4.66$ $2.290$ $0.52$ $r^{*}$
$5.17$ $2.67$ $0.56$ $4.78$ $2.33$ $0.53$ $(\sigma^{2},\tau)=(10^{-4},4)$ $w$
$4.67$ $2.29$ $0.49$ $4.53$ $2.34$ $0.43$ $r$ $4.67$ $2.29$ $0.49$ $4.53$
$2.34$ $0.43$ $r^{*}$ $4.58$ $2.21$ $0.47$ $4.67$ $2.40$ $0.47$
$(\sigma^{2},\tau)=(10^{-4},10^{2})$ $w$ $5.89$ $2.98$ $0.74$ $4.33$ $2.05$
$0.47$ $r$ $5.89$ $2.98$ $0.74$ $4.33$ $2.05$ $0.47$ $r^{*}$ $5.18$ $2.73$
$0.65$ $5.14$ $2.59$ $0.61$ $(\sigma^{2},\tau)=(10^{-3},1)$ $w$ $5.17$ $2.66$
$0.50$ $4.60$ $2.16$ $0.40$ $r$ $5.17$ $2.66$ $0.50$ $4.60$ $2.16$ $0.40$
$r^{*}$ $4.96$ $2.55$ $0.45$ $4.88$ $2.37$ $0.44$
$(\sigma^{2},\tau)=(10^{-3},2)$ $w$ $5.00$ $2.50$ $0.48$ $4.93$ $2.35$ $0.37$
$r$ $5.00$ $2.50$ $0.48$ $4.93$ $2.35$ $0.37$ $r^{*}$ $4.76$ $2.34$ $0.45$
$5.31$ $2.54$ $0.42$ $(\sigma^{2},\tau)=(10^{-3},4)$ $w$ $5.25$ $2.73$ $0.61$
$4.54$ $2.23$ $0.38$ $r$ $5.25$ $2.73$ $0.61$ $4.54$ $2.23$ $0.38$ $r^{*}$
$4.74$ $2.52$ $0.52$ $5.00$ $2.45$ $0.43$ $(\sigma^{2},\tau)=(10^{-2},1)$ $w$
$5.75$ $2.96$ $0.54$ $4.53$ $2.06$ $0.31$ $r$ $5.75$ $2.96$ $0.54$ $4.54$
$2.07$ $0.32$ $r^{*}$ $5.09$ $2.56$ $0.44$ $5.33$ $2.70$ $0.43$
$(\sigma^{2},\tau)=(10^{-2},2)$ $w$ $5.96$ $2.94$ $0.51$ $3.69$ $1.92$ $0.30$
$r$ $5.96$ $2.94$ $0.51$ $3.69$ $1.92$ $0.30$ $r*$ $4.89$ $2.39$ $0.38$ $4.87$
$2.53$ $0.49$ $(\sigma^{2},\tau)=(10^{-2},4)$ $w$ $6.20$ $3.26$ $0.57$ $2.73$
$1.14$ $0$ $r$ $6.20$ $3.26$ $0.57$ $2.73$ $1.14$ $0$ $r*$ $4.78$ $2.45$
$0.41$ $4.96$ $2.30$ $0.31$ SE 0.22 0.16 0.07 0.22 0.16 0.07

###  5.4 Case study C

Our third case study is based on a test case posted by NASA333
https://github.com/nasa/CARA_Analysis_Tools/tree/master/two-
dimension_Pc/UnitTest/InputFiles as one of the limited datasets for their
publicly-released Conjunction Assessment Risk Analysis (CARA) tool. Each test
case contains data that can be readily converted to the Earth-centered
inertial (ECI) reference frame for both the primary and secondary, $\eta_{s1}$
and $\eta_{s2}$, respectively. The test cases also contain state covariance
data, originally expressed in the Cartesian radial, in-track and cross-track
(RIC) reference frame (also called the UVW frame). The RIC covariance matrices
can be transformed to the ECI state covariances for the primary and secondary,
respectively denoted by $\Omega_{s1}^{-1}$ and $\Omega_{s1}^{-1}$, which are
$6\times 6$ matrices. These transformed quantities can then be used to define
the relative ECI state at close approach, $\eta=\eta_{s1}-\eta_{s2}$, and the
associated combined covariance,
$\Omega^{-1}=\Omega_{s1}^{-1}+\Omega_{s2}^{-1}$.

Column C of Table 3 shows the conjunction elements for this example. Since the
motion is linear and uncertainty on the velocity can be ignored, we use the
reduced two-dimensional model representing conjunction in the encounter plane.
The projected quantities in this plane, as descried in Appendix C, are given
by

$A=(CV,v/\|v\|)=\begin{bmatrix}0.72&0.11&-0.69\\\ 0.68&-0.30&0.67\\\
0.13&0.95&0.28\\\ \end{bmatrix},\quad
x^{\mathrm{o}}=A^{\mathrm{T}}y^{\mathrm{o}}=\begin{bmatrix}11.84\\\ -1.36\\\
\end{bmatrix},\quad
D=\operatorname{diag}\left(d_{1}^{2},d_{2}^{2}\right)=\begin{bmatrix}25.1^{2}&0.00\\\
0.00&11.61^{2}\\\ \end{bmatrix},$

so $\widehat{\psi}=\|x\|=11.92$ m. Not surprisingly, the significance
functions plotted in Figure 7 show that with $\psi_{0}=10$ m the significance
probabilities $p_{\rm obs}$ are all of order 0.5, though that for the modified
likelihood root is closer to 0.6; clearly evasive action would be essential in
this case.

Table 6 shows empirical error rates for the three pivots based on $10^{6}$
values of $x$ simulated from the bivariate normal distribution with
$\xi=x^{\mathrm{o}}$ and covariance matrix $c^{2}D$, with $c^{2}$ varying from
$0.005$ to $4$, as in Figure 2. In order to study the error rates for very
rare events, we took values of $\alpha$ ranging from 2.5% down to 0.005%, the
latter corresponding to two-sided 99.99% confidence intervals; the one-sided
error rates for $\alpha=0.05\%$ and $\alpha=0.005\%$ span the level
${\varepsilon}=0.01\%$ above which evasive action might be considered. Under
the conditions of the CARA event, $c^{2}=1$, the left-tail error is
systematically high and the right-tail error is exactly zero. This is
unsurprising, because in this setting the pivots can be close to zero, and
with increasing uncertainty the upper confidence limit is almost invariably
larger than $\psi_{0}$. As $c^{2}$ decreases, the uncertainty becomes
unrealistically small and the properties of all three pivots improve. When
$c^{2}$ increases, the error rates for $r$ and $w$ remain poor, and we see
that $r^{*}$ has better error rate properties, particularly in the left tail,
which is of most interest. Although the errors for $r^{*}$ are closer to the
nominal level, it should be jointly interpreted with the other pivots,
especially when the uncertainties are large relative to the miss distance. In
a second experiment we fixed $D$ to the real uncertainty matrix but increased
$\xi$ to $c^{\prime}\xi$, for $c^{\prime}>1$. A situation in which we expect
$\operatorname*{E}(\|x\|)\approx\psi_{0}$ for large $c^{\prime}$, hence
smaller bias in the probability of collision. Table 7 shows the resulting
error rates, again based on $10^{6}$ simulated values of $x$. All the error
rates approach their nominal values as $c^{\prime}$ increases, but $r^{*}$
again performs best overall, particularly in the left tail.

Figure 7: Summaries of the CARA event. Left: significance functions based on
likelihood root $r^{\mathrm{o}}(\psi)$ (solid black),Wald statistic $w(\psi)$
(dotted blue), and modified likelihood root $r^{*{\mathrm{o}}}(\psi)$ (red
dashes), Right: the same quantities transformed to the standard normal scale.
Table 6: Empirical left- and right-tail error rates ($\%$) at nominal levels
$\alpha=2.5\%,0.5\%,0.05\%$ and $0.005\%$ for the CARA event, Case Study C,
with variance matrix $c^{2}D$, based on $10^{6}$ Monte Carlo samples. The
standard errors (SE) appear in the last line.

$\begin{array}[]{ll cccc c
cccc}\hline\cr\hline\cr&&\lx@intercol\hfil\text{Left tail
}(\%)\hfil\lx@intercol&&\lx@intercol\hfil\text{Right tail
}(\%)\hfil\lx@intercol\\\ \cline{3-6}\cr\cline{8-11}\cr
c^{2}&\text{Statistic}&2.5&0.5&0.05&0.005&&2.5&0.5&0.05&0.005\\\
\hline\cr&w&2.5986&0.5229&0.0504&0.0062&&2.4616&0.4985&0.0542&0.0062\\\
0.005&r&2.6115&0.5275&0.0507&0.0062&&2.4274&0.4817&0.0494&0.0047\\\
&r^{*}&2.5206&0.5059&0.0490&0.0060&&2.5256&0.5031&0.0524&0.0048\\\
\hline\cr&w&2.6242&0.5278&0.0531&0.0054&&2.4530&0.5194&0.0665&0.0117\\\
0.01&r&2.6411&0.5347&0.0548&0.0057&&2.3701&0.4644&0.0433&0.0039\\\
&r^{*}&2.5064&0.5044&0.0500&0.0053&&2.5087&0.5008&0.0476&0.0042\\\
\hline\cr&w&2.7484&0.5611&0.0538&0.0051&&3.1974&1.3743&0.4281&0.0616\\\
0.05&r&2.8128&0.5804&0.0570&0.0055&&0.8226&0.0000&0.0000&0.0000\\\
&r^{*}&2.5127&0.5070&0.0474&0.0047&&1.3042&0.0000&0.0000&0.0000\\\
\hline\cr&w&2.7974&0.5681&0.0556&0.0057&&2.0243&0.3167&0.0000&0.0000\\\
0.1&r&2.9171&0.6048&0.0623&0.0062&&0.0000&0.0000&0.0000&0.0000\\\
&r^{*}&2.4808&0.5020&0.0481&0.0048&&0.0041&0.0004&0.0001&0.0000\\\
\hline\cr&w&2.8968&0.5855&0.0602&0.0061&&0.1192&0.0000&0.0000&0.0000\\\
0.2&r&3.1142&0.6499&0.0680&0.0071&&0.0000&0.0000&0.0000&0.0000\\\
&r^{*}&2.4780&0.4964&0.0497&0.0050&&0.0414&0.0079&0.0013&0.0004\\\
\hline\cr&w&3.1192&0.6139&0.0650&0.0082&&0.0000&0.0000&0.0000&0.0000\\\
0.5&r&3.6289&0.7466&0.0807&0.0104&&0.0000&0.0000&0.0000&0.0000\\\
&r^{*}&2.4578&0.4864&0.0510&0.0070&&0.3016&0.0862&0.0212&0.0080\\\
\hline\cr&w&3.3729&0.6352&0.0666&0.0088&&0.0000&0.0000&0.0000&0.0000\\\
0.8&r&4.1558&0.8421&0.0891&0.0110&&0.0000&0.0000&0.0000&0.0000\\\
&r^{*}&2.4974&0.4877&0.0508&0.0066&&0.6050&0.2070&0.0618&0.0264\\\
\hline\cr&w&3.5192&0.6687&0.0619&0.0050&&0.0000&0.0000&0.0000&0.0000\\\
1&r&4.4715&0.9130&0.0920&0.0075&&0.0000&0.0000&0.0000&0.0000\\\
&r^{*}&2.5252&0.5021&0.0475&0.0036&&0.7748&0.2829&0.0899&0.0380\\\
\hline\cr&w&4.1574&0.7724&0.0737&0.0081&&0.0000&0.0000&0.0000&0.0000\\\
2&r&5.7200&1.1922&0.1211&0.0129&&0.0000&0.0000&0.0000&0.0000\\\
&r^{*}&2.6649&0.5265&0.0512&0.0052&&1.5817&0.6505&0.2498&0.1222\\\
\hline\cr&w&4.6041&0.8318&0.0761&0.0070&&0.0000&0.0000&0.0000&0.0000\\\
3&r&6.5660&1.3730&0.1428&0.0131&&0.0000&0.0000&0.0000&0.0000\\\
&r^{*}&2.7466&0.5273&0.0516&0.0052&&2.2104&0.9820&0.4040&0.2065\\\
\hline\cr&w&4.9401&0.8807&0.0772&0.0070&&0.0000&0.0000&0.0000&0.0000\\\
4&r&7.2084&1.5009&0.1495&0.0147&&0.0000&0.0000&0.0000&0.0000\\\
&r^{*}&2.8042&0.5438&0.0520&0.0057&&2.8894&1.3420&0.5913&0.3171\\\
\hline\cr\hline\cr\mbox{SE }(\times
10^{-3})&&15.6&7.05&2.24&0.71&&15.6&7.05&2.24&0.71\\\
\hline\cr\hline\cr\end{array}$

Table 7: Empirical left- and right-tail error rates ($\%$) at nominal levels
$\alpha=2.5\%,0.5\%,0.05\%$ and $0.005\%$ for the CARA event, Case Study C,
with true position vector $c^{\prime}\xi$ in the encounter plane, based on
$10^{6}$ Monte Carlo samples. The standard errors (SE) appear in the last
line.

$\begin{array}[]{lc cccc c
cccc}\hline\cr\hline\cr&&\lx@intercol\hfil\text{Left tail
}(\%)\hfil\lx@intercol&&\lx@intercol\hfil\text{Right tail
}(\%)\hfil\lx@intercol\\\ \cline{3-6}\cr\cline{8-11}\cr
c^{\prime}&\text{Statistic}&2.5&0.5&0.05&0.005&&2.5&0.5&0.05&0.005\\\
\hline\cr&\text{Wald}&2.9354&0.5837&0.0570&0.0066&&0.0045&0.0000&0.0000&0.0000\\\
2&r&3.2085&0.6631&0.0677&0.0079&&0.0000&0.0000&0.0000&0.0000\\\
&r^{*}&2.4625&0.4821&0.0478&0.0055&&0.0749&0.0157&0.0027&0.0007\\\
\hline\cr&\text{Wald}&2.8360&0.5736&0.0620&0.0063&&1.6491&0.1615&0.0000&0.0000\\\
3&r&2.9632&0.6124&0.0675&0.0070&&0.0000&0.0000&0.0000&0.0000\\\
&r^{*}&2.5088&0.5057&0.0530&0.0048&&0.0052&0.0005&0.0000&0.0000\\\
\hline\cr&\text{Wald}&2.7575&0.5586&0.0580&0.0053&&3.0591&1.2037&0.1975&0.0077\\\
4&r&2.8315&0.5846&0.0615&0.0058&&0.0000&0.0000&0.0000&0.0000\\\
&r^{*}&2.4990&0.4999&0.0506&0.0046&&0.0052&0.0000&0.0000&0.0000\\\
\hline\cr&\text{Wald}&2.7127&0.5559&0.0523&0.0048&&3.1485&1.3268&0.5772&0.2385\\\
5&r&2.7627&0.5703&0.0554&0.0052&&2.1960&0.4215&0.0380&0.0027\\\
&r^{*}&2.5079&0.5045&0.0475&0.0040&&2.5059&0.4962&0.0480&0.0037\\\
\hline\cr&\text{Wald}&2.6819&0.5358&0.0523&0.0047&&2.7548&0.9488&0.4112&0.2252\\\
6&r&2.7249&0.5494&0.0537&0.0049&&2.2379&0.4364&0.0428&0.0051\\\
&r^{*}&2.4966&0.4969&0.0483&0.0044&&2.4794&0.4978&0.0508&0.0060\\\
\hline\cr\hline\cr\text{SE }(\times
10^{-3})&&15.6&7.05&2.24&0.71&&15.6&7.05&2.24&0.71\\\
\hline\cr\hline\cr\end{array}$

##  6 Further remarks

In some cases successive six-dimensional observation vectors
$y_{1},\ldots,y_{n}$ and corresponding $6\times 6$ variance matrices
$\Omega^{-1}_{1},\ldots,\Omega^{-1}_{n}$ are available, with the variance
matrices increasingly concentrated as information accrues on a conjunction. If
the observations can be regarded as independent, then the corresponding log
likelihood is

$\ell(\psi,\lambda_{1},\ldots,\lambda_{n})=-{\textstyle{1\over
2}}\sum_{j=1}^{n}\\{y_{j}-\eta(\vartheta_{j})\\}^{\mathrm{\scriptscriptstyle
T}}\Omega_{j}\\{y_{j}-\eta(\vartheta_{j})\\},$

where $\vartheta_{j}=(\psi,\lambda_{j})$, with $\psi$ representing the miss
distance common to all the observations and $\lambda_{1},\ldots,\lambda_{n}$
representing the $n$ $5\times 1$ vectors of nuisance parameters corresponding
to $y_{1},\ldots,y_{n}$. The more precise $y_{j}$ are automatically given
higher weight, since the corresponding dispersion matrices $\Omega_{j}$ are
larger. In this case the overall parameter vector is
$\vartheta=(\psi,\lambda_{1},\ldots,\lambda_{n})$, and the approach of Section
4 can be applied with minor changes.  A similar but more complex
generalisation should be feasible when the observations are dependent due to
batch updates.

More complicated geometry leading to a different form for
$\eta(\vartheta_{j})$ would be needed if the relative motions could not be
considered to be rectilinear.

Likelihood methods rest on distributional assumptions, and one might query
whether it is appropriate to assume normality of the observation vector.
Provided suitable data are available, standard methods could be used to check
model adequacy and used to elaborate the model if this was found to be
necessary. For example, the multivariate Student $t$ or Laplace distributions
might be used, though the numerical details would be more complex.

The incorporation of reliable prior information would be valuable, and might
for example be based on the output of a filtering approach to tracking. One
might then treat conjunction analysis as a prediction problem rather than an
estimation problem, and then the Bayesian formalism would be attractive.

##  7 Conclusion

In this paper we formulated conjunction assessment in statistical terms,
discussed approximate likelihood inference on the miss distance, and studied
the repeated-sampling properties of likelihood statistics with reference to
conjunction assessment.  Viewed in our framework the usual conjunction
probability estimator has a downward bias that seems not to have been noticed
previously, and the so-called probability dilution paradox vanishes, since it
refers to the properties of an estimator of the conjunction probability,
rather than to the probability itself. Examples illustrate inference on the
miss distance, suggest that standard likelihood confidence intervals may need
improvement when uncertainty on the relative distance and velocity is large,
and numerical work show how an improved approximation gives appreciably better
inferences. In our setup the estimated conjunction probability is replaced by
a significance probability for testing that the true miss distance is larger
than a specified safety threshold. If the model is correctly specified then
this significance probability is calibrated in a repeated-sampling sense and
thus provides a statistically well-founded basis for avoidance decisions.

##  Appendix A: Known velocity

The statistical model simplifies when the velocity $\nu$ of the second object
relative to the first can be assumed known. The random data then consist of a
$3\times 1$ vector $y$ containing the estimated position of the second object
relative to the first, and it is assumed that
$y\sim{{\mathcal{N}}}_{3}(\mu,\Omega^{-1})$, where $\mu$ and $\Omega$ are of
respective dimensions $3\times 1$ and $3\times 3$. Now let $A$ be a $3\times
3$ orthogonal matrix whose final column is $\nu/\|\nu\|$ and whose other
columns are chosen so that the upper left $2\times 2$ corner of
$A^{\mathrm{\scriptscriptstyle T}}\Omega^{-1}A$ is diagonal, say, $D={\rm
diag}(d_{1}^{2},d_{2}^{2})$. If so, then $x=A^{\mathrm{\scriptscriptstyle
T}}y\sim{{\mathcal{N}}}_{3}(A^{\mathrm{\scriptscriptstyle
T}}\mu,A^{\mathrm{\scriptscriptstyle T}}\Omega A)$, and the first two
components of $x$ are independent, with distribution
${{\mathcal{N}}}_{2}(\xi,D)$, say. Hence $\xi$ is the projection of the true
position of the second object along its velocity vector onto the encounter
plane, and thus $\psi=\|\xi\|$ is the miss distance. In terms of the
coordinates in the encounter plane defined by $A$ we can write
$\xi=(\psi\cos\lambda,\psi\sin\lambda)$ and $x=(x_{1},x_{2})$, where $\psi\geq
0$ is the parameter of interest and $\lambda\in[0,2\pi)$ is the nuisance
parameter. We describe how to obtain $A$ at the end of this appendix.

The transformation from $y$ to $x$ is invertible and depends only on the known
quantities $\nu$ and $\Omega$, so no statistical information is lost in using
$x$ instead of $y$. In satellite conjunction analysis it is customary to
ignore the coordinate $x_{3}$ of $y$ in the direction orthogonal to the
encounter plane, which in statistical terms amounts to basing inference on the
joint density of $x_{1}$ and $x_{2}$, leading to the log likelihood (23).

To obtain a suitable $3\times 3$ projection matrix $A$, note that if we define

$B=(b_{1},b_{2},\nu)=\begin{pmatrix}0&\nu_{2}^{2}+\nu_{3}^{2}&\nu_{1}\\\
\nu_{3}&-\nu_{1}\nu_{2}&\nu_{2}\\\
-\nu_{2}&-\nu_{1}\nu_{3}&\nu_{3}\end{pmatrix},\quad
N=\begin{pmatrix}\|b_{1}\|^{-1}&0&0\\\ 0&\|b_{2}\|^{-1}&0\\\
0&0&\|\nu\|^{-1}\end{pmatrix},$

then the first two columns of the orthonormal matrix $BN$ span the encounter
plane. If $C$ denotes the $3\times 2$ matrix containing these columns, then
$C^{\mathrm{\scriptscriptstyle T}}y$ is the orthogonal projection of $y$ onto
the encounter plane. If $VDV^{\mathrm{\scriptscriptstyle T}}$ is the spectral
decomposition of $C^{\mathrm{\scriptscriptstyle T}}\Omega^{-1}C$, then
$(CV)^{\mathrm{\scriptscriptstyle
T}}\Omega^{-1}(CV)=V^{\mathrm{\scriptscriptstyle
T}}C^{\mathrm{\scriptscriptstyle T}}\Omega^{-1}CV=D={\rm
diag}(d^{2}_{1},d^{2}_{2})$. Thus $(CV)^{\mathrm{\scriptscriptstyle T}}y$ is
bivariate normal with mean $\xi=(CV)^{\mathrm{\scriptscriptstyle T}}\mu$ and
covariance matrix $D$. Hence if we let $A=(CV,\nu/\|\nu\|)$, then the first
two elements of $x=A^{\mathrm{\scriptscriptstyle T}}y$ are independent, and
satisfy $x_{1}\sim{{\mathcal{N}}}(\xi_{1},d_{1}^{2})$ and
$x_{2}\sim{{\mathcal{N}}}(\xi_{2},d_{2}^{2})$.

## Appendix B: Theoretical background

The key to improving first-order approximations is to account for the presence
of nuisance parameters. One important improvement, the modified likelihood
root [34, 35], is given by equation (12). If the response distribution is
continuous and has true parameter $\psi_{0}$, then $r^{*}(\psi_{0})$ is
asymptotically standard normal to order $\mathrm{O}\left(n^{-3/2}\right)$.

The tangent exponential model proposed by Fraser and co-authors [e.g., 35] may
be used to construct (12). This is based on a local exponential family
approximation to the log-likelihood function for $n$ independent observations
$y_{1},\ldots,y_{n}$ that has $d\times 1$ canonical parameter

$\varphi^{\mathrm{\scriptscriptstyle
T}}(\vartheta)=\ell_{;V}(\vartheta;y^{\mathrm{o}})=V^{\mathrm{\scriptscriptstyle
T}}{\partial\ell(\vartheta)\over\partial
y}=\sum_{i=1}^{n}\dfrac{\partial\ell(\vartheta;y)}{\partial
y_{i}}\Bigr{\rvert}_{y_{i}=y_{i}^{\mathrm{o}}}V_{i},$ (24)

and involves a directional derivative of the log-likelihood function evaluated
at the observed data point
$y^{\mathrm{o}}=(y_{1}^{\mathrm{o}},\ldots,y_{n}^{\mathrm{o}})^{\mathrm{\scriptscriptstyle
T}}$, as $\ell_{;V}(\vartheta;y^{\mathrm{o}})$ is the derivative of
$\ell(\vartheta)$ in the directions given by the $n$ rows $V_{1},\ldots,V_{n}$
of the $n\times d$ matrix $V$. This matrix can be constructed using a vector
of pivotal quantities $z\equiv
z(y;\vartheta)=\\{z_{1}(y_{1},\vartheta),\ldots,z_{n}(y_{n},\vartheta)\\}^{\mathrm{\scriptscriptstyle
T}}$ through

$\displaystyle V=\left.\dfrac{\partial
y}{\partial\vartheta^{{\mathrm{\scriptscriptstyle
T}}}}\right|_{y=y^{\mathrm{o}},\vartheta=\widehat{\vartheta}^{\mathrm{o}}}=-\left.\left(\dfrac{\partial
z}{\partial y^{{\mathrm{\scriptscriptstyle
T}}}}\right)^{-1}\right|_{y=y^{\mathrm{o}},\vartheta=\widehat{\vartheta}^{\mathrm{o}}}\left.\times\dfrac{\partial
z}{\partial\vartheta^{{\mathrm{\scriptscriptstyle
T}}}}\right|_{y=y^{\mathrm{o}},\vartheta=\widehat{\vartheta}^{\mathrm{o}}},$
(25)

where $\widehat{\vartheta}^{\mathrm{o}}$ is the maximum likelihood estimate at
$y^{\mathrm{o}}$. The resulting correction term can be written in the form
(15). The numerator of the first term of $q$ is the determinant of a $d\times
d$ matrix whose first column is
$\varphi(\widehat{\vartheta})-\varphi(\widehat{\vartheta_{\psi}})$ and whose
remaining columns are $\varphi_{\lambda}(\widehat{\vartheta_{\psi}})$. Further
details can be found in Brazzale et al. [14, Section 8.5] and a guide to the
literature may be found in Davison and Reid [36].

In order to obtain the matrix $V$ for use when the velocity vector $\nu$ is
unknown, we define the vector of pivots as
$z(y,\vartheta)=\Omega^{1/2}\\{y-\eta(\vartheta)\\}$; these are independent
and standard normal under the model. Partial differentiation yields

$V={\partial\eta(\vartheta)\over\partial\vartheta^{\mathrm{\scriptscriptstyle
T}}}=\eta_{\vartheta}(\vartheta),$

evaluated at the maximum likelihood estimate
$\widehat{\vartheta}^{\mathrm{o}}$ corresponding to $y^{\mathrm{o}}$. The log-
likelihood, $-{\textstyle{1\over 2}}(y-\eta)^{\mathrm{\scriptscriptstyle
T}}\Omega(y-\eta)$, has derivative $\Omega(\eta-y)$ with respect to $y$, so
from (24) the canonical parameter of the tangent exponential model may be
written in the form $\varphi(\vartheta)=G\eta(\vartheta)+a$, where
$G=-\eta^{\mathrm{\scriptscriptstyle
T}}_{\vartheta}(\widehat{\vartheta}^{\mathrm{o}})\Omega$ and
$a=\eta^{\mathrm{\scriptscriptstyle
T}}_{\vartheta}(\widehat{\vartheta}^{\mathrm{o}})\Omega y^{\mathrm{o}}$ are
both constant with respect to $\vartheta$ and $G$ is full-rank. Any canonical
parameter that is an affine transformation of $\eta$ gives the same expression
for $q(\psi)$, because

$\displaystyle\dfrac{\left|\varphi(\widehat{\vartheta})-\varphi\left(\widehat{\vartheta}_{\psi}\right)\quad\varphi_{\lambda}\left(\widehat{\vartheta}_{\psi}\right)\right|}{\left|\varphi_{\vartheta}(\widehat{\vartheta})\right|}=$
$\displaystyle\dfrac{\left|G\left\\{\eta(\widehat{\vartheta})-\eta(\widehat{\vartheta}_{\psi})\right\\}\quad
G\eta_{\lambda}(\widehat{\vartheta}_{\psi})\right|}{\left|G\eta_{\vartheta}(\widehat{\vartheta})\right|}=\dfrac{\left|\eta(\widehat{\vartheta})-\eta(\widehat{\vartheta}_{\psi})\quad\eta_{\lambda}(\widehat{\vartheta}_{\psi})\right|}{\left|\eta_{\vartheta}(\widehat{\vartheta})\right|},$
(26)

where
$\eta_{\lambda}(\vartheta)=\partial{\eta}(\vartheta)/\partial\lambda^{{\mathrm{\scriptscriptstyle
T}}}$. Hence we can take the constructed parameter (24) to be
$\varphi(\vartheta)=\eta(\vartheta)$. To compute (26), we need the $6\times 6$
Jacobian $\eta_{\vartheta}$ and the second derivatives
$\eta_{\vartheta\vartheta}$, a $6\times 6\times 6$ tensor containing the
second derivatives of $\eta$, which is needed to compute
$\jmath_{\lambda\lambda}(\widehat{\vartheta}_{\psi})$. The score equation and
the observed Fisher information can respectively be given as

${\partial\eta^{\mathrm{\scriptscriptstyle
T}}(\vartheta)\over\partial\vartheta}\Omega\\{y-\eta(\vartheta)\\}=0,\quad{\partial\eta^{\mathrm{\scriptscriptstyle
T}}(\vartheta)\over\partial\vartheta}\Omega{\partial\eta(\vartheta)\over\partial\vartheta_{r}}-{\partial^{2}\eta^{\mathrm{\scriptscriptstyle
T}}(\vartheta)\over\partial\vartheta\partial\vartheta_{r}}\Omega\\{y-\eta(\vartheta)\\},\quad
r=1,\ldots,6.$ (27)

The observed information matrix evaluated at the maximum likelihood estimate
for the full model equals
$\jmath(\widehat{\vartheta})=\eta^{\mathrm{\scriptscriptstyle
T}}_{\vartheta}(\widehat{\vartheta})\Omega\eta_{\vartheta}(\widehat{\vartheta})$,
because the score equation for the full model implies that
$\eta(\widehat{\vartheta})=y$.

A fuller version of a similar computation is in Fraser et al. [37].

In the simplified problem with known $\nu$ there are just two parameters
$\psi\geq 0$ and $\lambda\in[0,2\pi)$, the log likelihood reduces to (23) and
$\varphi(\psi,\lambda)=\xi=(\psi\cos\lambda,\psi\sin\lambda)^{\mathrm{\scriptscriptstyle
T}}$. The overall maximum likelihood estimates are

$\widehat{\psi}=(x_{1}^{2}+x_{2}^{2})^{1/2},\quad\widehat{\lambda}=\arctan(x_{2}/x_{1}).$
(28)

If $\psi$ is fixed, then $\widehat{\lambda}_{\psi}$ is readily found as the
unique minimum of the sum of squares in (23); a good starting value should be
$\widehat{\lambda}$. Then (13, 15) reduce to

$\displaystyle r(\psi)$
$\displaystyle=\text{sign}(\widehat{\psi}-\psi)\dfrac{1}{d_{1}d_{2}}\left\\{d_{2}^{2}\left(x_{1}-\psi\cos\lambda\right)^{2}+d_{1}^{2}\left(x_{2}-\psi\sin\lambda\right)^{2}\right\\}^{1/2},$
$\displaystyle q(\psi)$
$\displaystyle=\psi^{1/2}{x_{1}\cos\widehat{\lambda}_{\psi}+x_{2}\sin\widehat{\lambda}_{\psi}-\psi\over\left\\{d_{2}^{2}(x_{1}\cos\widehat{\lambda}_{\psi}-\psi\cos
2\widehat{\lambda}_{\psi})+d_{1}^{2}(x_{2}\sin\widehat{\lambda}_{\psi}+\psi\cos
2\widehat{\lambda}_{\psi})\right\\}^{1/2}}.$

When $d_{1}=d_{2}=d$, say, then
$\widehat{\lambda}_{\psi}\equiv\widehat{\lambda}$ does not depend on $\psi$
and after simplification we obtain

$r(\psi)=w(\psi)=(\widehat{\psi}-\psi)/d,\quad
q(\psi)=r(\psi)(\psi/\widehat{\psi})^{1/2},\quad
r^{*}(\psi)={\widehat{\psi}-\psi\over d}+{d\over
2(\widehat{\psi}-\psi)}\log(\psi/\widehat{\psi}),\quad\psi>0;$

note that $r^{*}(\psi)\to-d/(2\widehat{\psi})$ when $\psi\to\widehat{\psi}$.
The fact that $r^{*}(\psi)<r(\psi)$ for all $\psi$ implies that confidence
intervals for $\psi$ based on $r^{*}(\psi)$ will be closer to the origin than
those based on $r(\psi)$, and significance levels for fixed $\psi$ will be
higher, leading to more conservative inferences. We should consider evasive
action based on $r(\psi)$ when
$1-\Phi\\{r^{\mathrm{o}}(\psi_{0})\\}<{\varepsilon}$, i.e., when
$\widehat{\psi}^{\mathrm{o}}-dz_{\epsilon}$ is smaller than the hard-body
radius $\psi_{0}$; if $\epsilon=10^{-4}$ for example, then
$z_{\epsilon}=-3.72$. Notice that this rule relates the observed distance of
the second object from the origin, $\widehat{\psi}^{\mathrm{o}}$, to the
measurement uncertainty, $d$. Use of $r^{*{\mathrm{o}}}(\psi)$ will lead to
very similar conclusions in most cases.

## Appendix C: Bayesian approximation

The derivation of equation (19) is outlined in Section 8.7 of Brazzale et al.
[14]. As $\widehat{\vartheta}^{\mathrm{\scriptscriptstyle
T}}=(\psi,\widehat{\lambda}_{\psi}^{\mathrm{\scriptscriptstyle T}})$,

${\mathrm{d}{\ell(\widehat{\vartheta}_{\psi})}\over\mathrm{d}{\psi}}=\ell_{\psi}(\widehat{\vartheta}_{\psi})+{\partial\widehat{\lambda}^{\mathrm{\scriptscriptstyle
T}}_{\psi}\over\partial\psi}\ell_{\lambda}(\widehat{\vartheta}_{\psi})=\ell_{\psi}(\widehat{\vartheta}_{\psi})-\ell_{\psi\lambda}(\widehat{\vartheta}_{\psi})\ell^{-1}_{\lambda\lambda}(\widehat{\vartheta}_{\psi})\ell_{\lambda}(\widehat{\vartheta}_{\psi}),$

where subscripts on $\ell$ denote partial differentiation and the second
equality follows from noting that differentiating the equation
$\ell_{\lambda}(\widehat{\vartheta}_{\psi})=0$ that defines
$\widehat{\lambda}_{\psi}$ yields

$\ell_{\psi\lambda}(\widehat{\vartheta}_{\psi})+{\partial\widehat{\lambda}^{\mathrm{\scriptscriptstyle
T}}_{\psi}\over\partial\psi}\,\ell_{\lambda\lambda}(\widehat{\vartheta}_{\psi})=0,$

with the matrix $\ell_{\lambda\lambda}(\widehat{\vartheta}_{\psi})$ invertible
because it is the Hessian corresponding to the maximum of $\ell$ in the
$\lambda$ direction for fixed $\psi$. A standard identity for the determinant
of a partitioned matrix gives

$\left|\begin{matrix}\ell_{\vartheta}(\widehat{\vartheta}_{\psi})&\ell_{\vartheta\lambda}(\widehat{\vartheta}_{\psi})\end{matrix}\right|=\left|\begin{matrix}\ell_{\psi}(\widehat{\vartheta}_{\psi})&\ell_{\psi\lambda}(\widehat{\vartheta}_{\psi})\\\
\ell_{\lambda}(\widehat{\vartheta}_{\psi})&\ell_{\lambda\lambda}(\widehat{\vartheta}_{\psi})\end{matrix}\right|=\left\\{\ell_{\psi}(\widehat{\vartheta}_{\psi})-\ell_{\psi\lambda}(\widehat{\vartheta}_{\psi})\ell^{-1}_{\lambda\lambda}(\widehat{\vartheta}_{\psi})\ell_{\lambda}(\widehat{\vartheta}_{\psi})\right\\}\times\left|\ell_{\lambda\lambda}(\widehat{\vartheta}_{\psi})\right|,$
(29)

and the expression for the observed information in (27) implies that we can
write

$-\ell_{\vartheta\lambda_{r}}(\vartheta)=\eta^{\mathrm{\scriptscriptstyle
T}}_{\vartheta}(\vartheta)\Omega\eta_{\lambda_{r}}(\vartheta)+A_{r}(\vartheta)\\{y-\eta(\vartheta)\\},\quad
r=1,\ldots,d-1,$

where $A(\vartheta)$ involves second derivatives of $\eta$. Now
$y=\eta(\widehat{\vartheta})$, so if $\widehat{\psi}-\psi$ is of order
$n^{-1/2}$, then
$y-\eta(\widehat{\vartheta}_{\psi})=\eta(\widehat{\vartheta})-\eta(\widehat{\vartheta}_{\psi})$
is also $O(n^{-1/2})$, and hence so too is the second term of
$-\ell_{\vartheta\lambda}(\vartheta)$. Thus (29) equals

$\displaystyle\left|\begin{matrix}\ell_{\vartheta}(\widehat{\vartheta}_{\psi})&\ell_{\vartheta\lambda}(\widehat{\vartheta}_{\psi})\end{matrix}\right|$
$\displaystyle=$
$\displaystyle(-1)^{d-1}\left|\begin{matrix}\ell_{\vartheta}(\widehat{\vartheta}_{\psi})&-\ell_{\vartheta\lambda}(\widehat{\vartheta}_{\psi})\end{matrix}\right|,$
$\displaystyle=$
$\displaystyle(-1)^{d-1}\left|\begin{matrix}\eta_{\vartheta}(\widehat{\vartheta}_{\psi})\Omega\\{y-\eta(\widehat{\vartheta}_{\psi})\\}&\eta^{\mathrm{\scriptscriptstyle
T}}_{\vartheta}(\widehat{\vartheta}_{\psi})\Omega\eta_{\lambda}(\widehat{\vartheta}_{\psi})\end{matrix}\right|+O(n^{-1/2}),$
$\displaystyle=$
$\displaystyle(-1)^{d-1}\left|\eta_{\vartheta}(\widehat{\vartheta}_{\psi})\right|\left|\Omega\right|\left|\begin{matrix}\eta(\widehat{\vartheta})-\eta(\widehat{\vartheta}_{\psi})&\eta_{\lambda}(\widehat{\vartheta}_{\psi})\end{matrix}\right|+O(n^{-1/2}).$

As $\vartheta$ is of dimension $d$ and
$|\jmath(\widehat{\vartheta})|=\left|\eta^{\mathrm{\scriptscriptstyle
T}}_{\vartheta}(\widehat{\vartheta})\right|^{2}|\Omega|$, equation (19) equals

$\displaystyle q_{B}(\psi)$ $\displaystyle=$
$\displaystyle{\left|\begin{matrix}\ell_{\vartheta}(\widehat{\vartheta}_{\psi})&\ell_{\vartheta\lambda}(\widehat{\vartheta}_{\psi})\end{matrix}\right|\over{|\ell_{\lambda\lambda}(\widehat{\vartheta}_{\psi})|}}\times{|\jmath_{\lambda\lambda}(\widehat{\vartheta}_{\psi})|^{1/2}\over|\jmath(\widehat{\vartheta})|^{1/2}}{f(\widehat{\vartheta})\over
f(\widehat{\vartheta}_{\psi})},$ (30) $\displaystyle=$
$\displaystyle{\left|\begin{matrix}\ell_{\vartheta}(\widehat{\vartheta}_{\psi})&-\ell_{\vartheta\lambda}(\widehat{\vartheta}_{\psi})\end{matrix}\right|\over{|\jmath_{\lambda\lambda}(\widehat{\vartheta}_{\psi})|^{1/2}|\jmath(\widehat{\vartheta})|^{1/2}}}\times{f(\widehat{\vartheta})\over
f(\widehat{\vartheta}_{\psi})},$ $\displaystyle=$
$\displaystyle{\left|\eta_{\vartheta}(\widehat{\vartheta}_{\psi})\right|\left|\Omega\right|\left|\begin{matrix}\eta(\widehat{\vartheta})-\eta(\widehat{\vartheta}_{\psi})&\eta_{\lambda}(\widehat{\vartheta}_{\psi})\end{matrix}\right|\over|\jmath_{\lambda\lambda}(\widehat{\vartheta}_{\psi})|^{1/2}\left|\eta^{\mathrm{\scriptscriptstyle
T}}_{\vartheta}(\widehat{\vartheta})\right||\Omega|^{1/2}}\times{f(\widehat{\vartheta})\over
f(\widehat{\vartheta}_{\psi})}+O(n^{-1/2}),$ $\displaystyle=$ $\displaystyle
q(\psi)\times{\left|\eta^{\mathrm{\scriptscriptstyle
T}}_{\vartheta}(\widehat{\vartheta}_{\psi})\right|\over{\left|\eta^{\mathrm{\scriptscriptstyle
T}}_{\vartheta}(\widehat{\vartheta})\right|}}\times{f(\widehat{\vartheta})\over
f(\widehat{\vartheta}_{\psi})}+O(n^{-1/2}).$

where $q(\psi)$ is given in equation (15).

The Jeffreys prior is the root of the determinant of the Fisher information
matrix,

$f(\vartheta)\propto\left|\eta^{\mathrm{\scriptscriptstyle
T}}_{\vartheta}(\vartheta)\Omega\eta_{\vartheta}(\vartheta)\right|^{1/2}\propto\left|\eta_{\vartheta}(\vartheta)\right|,$

as the constant $|\Omega|$ can be ignored. If this prior is used then (30)
simplifies to equation (22), plus a term of $O(n^{-1/2})$. Hence
$r_{B}(\psi)=r^{*}(\psi)+O(n^{-1})$, and inferences from (16) and (18) will be
the same to this order of error.

## Appendix D: Implementation details

The main steps for the numerical implementation of higher-order quantities are
given in table 9. Numerical experiments have been obtained with Intel® Xeon(R)
CPU E5-1650 v3 @3.50GHz × 12 processor with 64GB ram. The R code can be
designed to use parallel loops for speed-up purposes. In terms computational
cost (CPU time), it takes roughly 8 seconds to run $10^{3}$ simulations in the
two-dimensional settings.

Numerical details
---
Input: $6\times 1$ state vector $y$ containing the observed $3\times 1$
relative position and velocity vectors $\widehat{\mu}$ and $\widehat{\nu}$,
$6\times 6$ covariance matrix $\Omega^{-1}$ for $y$, and a nominal confidence
level $\alpha$.
1. Compute the overall maximum likelihood estimates $\widehat{\theta}$: (i) In the six-dimensional setting, $\widehat{\psi}$ and $\widehat{\lambda}=(\widehat{\theta}_{1},\widehat{\phi}_{1},\widehat{\|}\nu\|,\widehat{\theta}_{2},\widehat{\phi}_{2})$ are obtained from $(\widehat{\mu},\widehat{\nu})$ using equations (1)–(4). (ii) In the simplified two-dimensional setting, we first need to project relative quantities into the encounter plane using the matrix $A$ (see Appendix A), then $\widehat{\psi}$ and $\widehat{\lambda}$ are given in equation (28). 2. Compute the observed information matrix $\jmath(\widehat{\vartheta})$ given by (21) (this is a $6\times 6$ or $2\times$ 2 matrix), the estimated variance of each $\widehat{\vartheta}_{r}$ is the $(r,r)$ element of $\jmath(\widehat{\vartheta})^{-1}$. Let $\text{se}(\widehat{\psi})$ denote the square root of the estimated variance for $\widehat{\psi}$. 3. Define a grid ${\mathcal{G}}=\\{\psi^{\prime}_{1},\ldots,\psi^{\prime}_{n_{\psi}}\\}$ of values of $\psi$ that includes the maximum likelihood estimate $\widehat{\psi}$. We use a non-uniform grid in the interval $\left[\max\left(\widehat{\psi}-z_{1-\alpha}\text{se}(\widehat{\psi}),0\right),\widehat{\psi}+z_{\alpha}\text{se}(\widehat{\psi})\right]$, such that the mesh is finer for small $\psi$ and coarser for larger values. 4. For each $\psi\in{\mathcal{G}}$, obtain the constrained estimates $\widehat{\lambda}_{\psi}$ subject to $0\leq\theta_{1},\theta_{2}\leq\pi$, $-\pi\leq\phi_{1},\phi_{2}\leq\pi$ and $0<||\nu||$, by minimising the sum of squares in (20), and store the corresponding values of $\widehat{\vartheta}_{\psi}=(\psi,\widehat{\lambda}_{\psi})$. We have found that • it may help to transform the components of $\vartheta$ to take values in the real line, for example removing the restrictions by maximizing in terms of $\log||\nu||,\log\left\\{\theta_{i}/(\theta_{i}-\pi)\right\\}$, and $\tan(\phi_{i}/2)$, for $i=1,2$, and • constrained optimisation using the ‘Rvmmin’ or ‘nlminb’ solvers in the optimx() function of the R package optimx leads to an algorithm that is overall robust, fast, and generally insensitive to perturbations in initial values. 5. Use partial derivatives of the log-likelihood and the mean vector $\eta$ to evaluate $\ell(\widehat{\vartheta}_{\psi})$, $\jmath_{\lambda\lambda}(\widehat{\vartheta}_{\psi})$, $\eta(\widehat{\vartheta}_{\psi})$ and $\eta_{\lambda}(\widehat{\vartheta})$ and then use expressions (13), (14) and (15) to obtain the values of $r(\psi)$, $w(\psi)$, $q(\psi)$ and $r^{*}(\psi)$ on ${\mathcal{G}}$. 6. Interpolate $r(\psi)$ and $r^{*}(\psi)$ on ${\mathcal{G}}$ by (for example) a cubic smoothing spline in which the values of $\psi$ are treated as functions of those of $r(\psi)$ and $r^{*}(\psi)$. Very large values of $r^{*}$ arising in a few cases due to numerical instabilities when $|r|<0.1$ are excluded. 7. If required, obtain the point estimate $\widehat{\psi}^{*}$ of $\psi$ by using the interpolating spline for $r^{*}(\psi)=0$. This is not needed for the Wald or the likelihood root, as $\widehat{\psi}$ is already known. 8. Obtain a $(1-2\alpha)$ confidence interval $(\psi_{\alpha},\psi_{1-\alpha})$ based on $r(\psi)$ as the solutions of $r(\psi)=\pm z_{\alpha}$. If the equation $r(\psi)=-z_{\alpha}$ cannot be solved, then the lower limit of the confidence interval is $\psi_{\alpha}=0$. Confidence intervals based on $r^{*}(\psi)$ are obtained likewise.
Output: Point estimates $\widehat{\psi}$ and $\widehat{\psi}^{*}$ of the miss
distance and corresponding two-sided $(1-2\alpha)$ confidence intervals.
Table 9: Miss distance interval estimation

## Acknowledgement

We thank Alessandra Brazzale, Nancy Reid and very constructive reviewers for
helpful remarks that led to major improvements in our work. The work was
supported by the Swiss National Science Foundation.

## References

  * Chen et al. [2017] Chen, L., Bai, X. Z., Liang, Y. G., and Li, K. B., _Orbital Data Applications for Space Objects: Conjunction Assessment and Situation Analysis_ , Springer, Singapore, 2017.  10.1007/978-981-10-2963-9.
  * Vallado [2013] Vallado, D. A., _Fundamentals of Astrodynamics and Applications_ , 4th ed., Microcosm Press, Hawthorne, CA, 2013.
  * Foster and Estes [1992] Foster, J. L., and Estes, H. S., _A parametric Analysis of Orbital Debris Collision Probability and Maneuver Rate for Space Vehicles_ , NASA, National Aeronautics and Space Administration, Lyndon B. Johnson Space Center, Houston, 1992.
  * Chan [1997] Chan, K., “Collision Probability Analyses for Earth Orbiting Satellites,” _Pacific Basin Societies, Advances in the Astronautical Sciences, 7th International Space Conference_ , Vol. 96, 1997, pp. 1033–1048.
  * Alfano [2005] Alfano, S., “A Numerical Implementation of Spherical Object Collision Probability,” _The Journal of the Astronautical Sciences_ , Vol. 53, No. 1, 2005, pp. 103–109.  10.1007/BF03546397.
  * Alfano [2006] Alfano, S., “Satellite Collision Probability Enhancements,” _Journal of Guidance Control and Dynamics_ , Vol. 29, 2006, pp. 588–592.  10.2514/1.15523.
  * Patera [2005] Patera, R., “Calculating Collision Probability for Arbitrary Space Vehicle Shapes via Numerical Quadrature,” _Journal of Guidance Control and Dynamics_ , Vol. 28, 2005, pp. 1326–1328.  10.2514/1.14526.
  * Garcia-Pelayo and Hernando-Ayuso [2016] Garcia-Pelayo, R., and Hernando-Ayuso, J., “Series for Collision Probability in Short-Encounter Model,” _Journal of Guidance Control and Dynamics_ , Vol. 39, No. 8, 2016, pp. 1904–1912.  10.2514/1.G001754.
  * Balch et al. [2019] Balch, M. S., Martin, R., and Ferson, S., “Satellite Conjunction Analysis and the False Confidence Theorem,” _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , Vol. 475, No. 2227, 2019.  10.1098/rspa.2018.0565.
  * Alfano and Negron [1993] Alfano, S., and Negron, D. J., “Determining Satellite Close Approaches,” _Journal of the Astronautical Sciences_ , Vol. 41, No. 2, 1993, pp. 217–225.  URL https://ui.adsabs.harvard.edu/abs/1993JAnSc..41..217A.
  * Gronchi [2005] Gronchi, G. F., “An Algebraic Method to Compute the Critical Points of the Distance Function Between Two Keplerian Orbits,” _Celestial Mechanics and Dynamical Astronomy_ , Vol. 93, No. 1, 2005, pp. 295–329.  10.1007/s10569-005-1623-5.
  * Armellin et al. [2010] Armellin, R., Lizia, P. D., Berz, M., and Makino, K., “Computing the Critical Points of the Distance function Between Two Keplerian Orbits Via Rigorous Global Optimization,” _Celestial Mechanics and Dynamical Astronomy_ , Vol. 107, No. 3, 2010, pp. 377–395.  10.1007/s10569-010-9281-7.
  * Cox and Hinkley [1974] Cox, D. R., and Hinkley, D. V., _Theoretical Statistics_ , Chapman & Hall, London, 1974.  10.1201/b14832.
  * Brazzale et al. [2007] Brazzale, A. R., Davison, A. C., and Reid, N., _Applied Asymptotics: Case Studies in Small Sample Statistics_ , Cambridge University Press, Cambridge, 2007\.  10.1017/CBO9780511611131.
  * Cunen et al. [2020] Cunen, C., Hjort, N. L., and Schweder, T., “Confidence in Confidence Distributions!” _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , Vol. 476, No. 2237, 2020, p. 20190781.  10.1098/rspa.2019.0781.
  * Neyman [1937] Neyman, J., “Outline of a Theory of Statistical Estimation Based on the Classical Theory of Probability,” _Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences_ , Vol. 236, No. 767, 1937, pp. 333–380.  URL http://www.jstor.org/stable/91337.
  * Klinkrad [2006] Klinkrad, H., _Space Debris: Models and Risk Analysis_ , Astronautical Engineering, Springer, Berlin, Heidelberg, 2006.  10.1007/3-540-37674-7.
  * Hill [1878] Hill, G. W., “Researches in the Lunar Theory,” _American Journal of Mathematics_ , Vol. 1, No. 1, 1878, pp. 5–26.  10.2307/2369430.
  * Clohessy and Wiltshire [1960] Clohessy, W. H., and Wiltshire, R. S., “Terminal Guidance System for Satellite Rendezvous,” _Journal of the Aerospace Sciences_ , Vol. 27, No. 9, 1960, pp. 653–658.  10.2514/8.8704.
  * Alfriend et al. [1999] Alfriend, K. T., Akella, M. R., Frisbee, J., Foster, J. L., Lee, D. J., and Wilkins, M., “Probability of Collision Error Analysis,” _Space Debris_ , Vol. 1, No. 1, 1999, pp. 21–35.  10.1023/A:1010056509803.
  * Hall [2021] Hall, D. T., “Expected collision rates for tracked satellites,” _Journal of Spacecraft and Rockets_ , Vol. 28, No. 3, 2021, pp. 715–728.  10.2514/1.A34919.
  * Edwards [1972] Edwards, A. W. F., _Likelihood_ , Cambridge University Press, Cambridge, 1972\.
  * Fisher [1973] Fisher, R. A., _Statistical Methods and Scientific Inference_ , Third ed., Collier–Macmillan, New York, 1973.
  * Fraser [2017] Fraser, D. A. S., “$p$-values: The insight to modern statistical inference,” _Annual Review of Statistics and its Application_ , Vol. 4, 2017, pp. 1–14.  10.1146/annurev-statistics-060116-054139.
  * Fraser [2019] Fraser, D. A. S., “The $p$-value function and statistical inference,” _American Statistician_ , Vol. 73, 2019, pp. 135–147.  10.1080/00031305.2018.1556735.
  * Schweder and Hjort [2016] Schweder, T., and Hjort, N. L., _Confidence, Likelihood, Probability: Statistical Inference with Confidence Distributions_ , Cambridge University Press, Cambridge, 2016.  10.1017/CBO9781139046671.
  * Fraser et al. [2016] Fraser, D. A. S., Bédard, M., Wong, A., Lin, W., and Fraser, A. M., “Bayes, reproducibility and the quest for truth,” _Statistical Science_ , Vol. 31, 2016, pp. 578–590.  10.1214/16-STS573.
  * Davison and Sartori [2008] Davison, A. C., and Sartori, N., “The Banff Challenge: Statistical detection of a noisy signal,” _Statistical Science_ , Vol. 23, 2008, pp. 354–364.  10.1214/08-STS260.
  * Davison [2003] Davison, A. C., _Statistical Models_ , Cambridge University Press, Cambridge, 2003.  10.1017/CBO9780511815850.
  * Carpenter [2019] Carpenter, J., “Covariance Realism is not enough,” _AAS Advances in the Astronautical Sciences_ , 2019.
  * Hejduk et al. [2019] Hejduk, M. D., Snow, D. E., and Newman, L. K., “Satellite Conjunction Assessment Risk Analysis for ”Dilution Region" Events: Issues and Operational Approaches,” _Space Traffic Management Conference_ , Vol. 28, 2019\.  URL https://commons.erau.edu/stm/2019/presentations/28.
  * Carpenter et al. [2017] Carpenter, J., Alfano, S., Hall, D. T., Hejduk, M., and Gaebler, J., “Relevance of the American Statistical Society’s Warning on P-Values for Conjunction Assessment,” _AAS/AIAA Astrodynamics Specialist Conference_ , Stevenson WA, 2017.
  * Cox [2020] Cox, D., “Statistical Significance,” _Annual Review of Statistics and Its Application_ , Vol. 7, No. 1, 2020, pp. 1–10.  10.1146/annurev-statistics-031219-041051.
  * Barndorff-Nielsen and Cox [1994] Barndorff-Nielsen, O. E., and Cox, D. R., _Inference and Asymptotics_ , Chapman & Hall, London, 1994.  10.1201/9780203750940.
  * Fraser et al. [1999a] Fraser, D. A. S., Reid, N., and Wu, J., “A simple general formula for tail probabilities for frequentist and Bayesian inference,” _Biometrika_ , Vol. 86, 1999a, pp. 249–264.  10.1093/biomet/86.2.249.
  * Davison and Reid [2022] Davison, A. C., and Reid, N., “The tangent exponential model,” _Handbook of Bayesian, Fiducial and Frequentist Inference_ , edited by J. O. Berger, X.-L. Meng, N. Reid, and M. Xie, Chapman & Hall/CRC, Boca Raton, FL, 2022, p. to appear.
  * Fraser et al. [1999b] Fraser, D. A. S., Wong, A., and Wu, J., “Regression analysis, nonlinear or nonnormal: Simple and accurate calculation of $p$-values from likelihood analysis,” _Journal of the American Statistical Association_ , Vol. 94, 1999b, pp. 1265–1295.  10.2307/2669942.