# Detection Rate of $<$50-meter Interstellar Objects with LSST

Carson Ezell111These authors contributed equally to this work.
<EMAIL_ADDRESS>Abraham Loeb22footnotemark: 2
<EMAIL_ADDRESS>Department of Astronomy, Harvard University, 60 Garden
St, Cambridge, Massachusetts 02138, United States

###### Abstract

The previous decade saw the discovery of the first four known interstellar
objects due to advances in astronomical viewing equipment. Future sky surveys
with greater sensitivity will allow for more frequent detections of such
objects, including increasingly small objects. We consider the capabilities of
the Legacy Survey of Space and Time (LSST) of the Vera C. Rubin Observatory to
detect interstellar objects of small sizes during its period of operation over
the next decade. We use LSST’s detection capabilities and simulate populations
of interstellar objects in the range of 1-50m in diameter to calculate the
expected number of small interstellar objects that will be detected. We use
previous detections of interstellar objects to calibrate our object density
estimates. We also consider the impact of the population’s albedo on detection
rates by considering populations with two separate albedo distributions: a
constant albedo of 0.06 and an albedo distribution that resembles near earth
asteroids. We find that the number of detections increases with the diameter
over the range of diameters we consider. We estimate a detection rate of up to
a small ISO every two years of LSST’s operation with an increase by a factor
of ten for future surveys that extend a magnitude deeper.

###### keywords:

interstellar objects , astronomical instrumentation

abstract

##  1 Introduction

The NASA Authorization Act of 2005 included a mandate—known as the George E.
Brown Jr. mandate—to “detect, track, catalog, and characterize the physical
characteristics of near-Earth objects equal or greater than 140 meters in
diameter in order to assess the threat of such near-Earth objects to the
Earth” [1]. Fulfilling the mandate has served as a guiding principle for
astronomical instruments tracking objects in the solar system since it was
created. As a result, less attention has historically been paid to objects
smaller than $140$ m in diameter as well as rare interstellar objects (ISOs)
when designing equipment, data processing pipelines, or conducting data
analysis. While surveys of the solar system would have catalogued any such
objects that they detected, most previous sky surveys would have been unable
to detect small objects, which we define as objects in the $\leq$50 m diameter
range, because of the small region around the Earth in which they would have
been observable.

Over the past decade, the main contributors to cataloging the population of
near-Earth objects (NEO) and potentially hazardous asteroids (PHAs) have been
the Catalina Sky Survey [2] and Pan-STARRS project [3]. The Legacy Survey of
Space and Time (LSST) on the Vera C. Rubin Observatory in Chile will be more
sensitive than previous telescopes, resulting in significant scientific
implications. Importantly for NASA, LSST will increase the completeness of the
catalog of NEOs. Numerous studies have estimated that LSST will identify
between 62 percent and 75 percent of this PHA population using its baseline
design [4, 5, 6, 7], and the completeness of the PHA catalog would then reach
about 86 percent [7].

The increased sensitivity of LSST also suggests that it will improve our
ability to detect objects of smaller sizes and at greater distances, including
small ISOs. In this paper, we place constraints on LSST’s detection
capabilities for small ISOs and calculate the expected number of such objects
that we expect it to detect over its decade of operation.

The Vera C. Rubin Telescope on which the LSST will be conducted has a diameter
of 8.4 m (with a 6.67 m diameter mirror), providing a 9.6 deg2 field of view
[4] with a 3.2 billion pixel camera that surveys the entire southern sky every
4 days. Taking an inventory of the solar system is one of the four main
science themes of the LSST design [8]. Previous studies have estimated the
rate of detection of interstellar objects by LSST, but these studies have
focused on ‘Oumuamua-like objects [9, 10] or asteroids and comets with a
diameter greater than 1 km [11]. Cook et al. [12] considered a variety of
objects, including interstellar asteroids as small in diameter as $20$ m.
However, these estimates were made prior to the detection of several
interstellar objects in as small as the $\sim$ 1 m range, which have led to
updated estimates of ISO density [13]. Previous studies also do not
characterize the detection capability of astronomical surveys for objects
$\leq$10 m in diameter, which may be observable under particular conditions.
We do not consider objects under $5$ m in diameter in our analysis because
they must pass very close to Earth to be observed, leading to significant
trailing losses that render them undetectable [12].

There would be important scientific benefits to detecting small ISOs,
especially given the scarcity of previous ISO observations altogether. Over
the past decade, the first four interstellar objects (ISOs) have been
discovered. The first interstellar object reported was ‘Oumuamua, which had a
diameter of about $200$m and was detected by the PAN-STAARS survey [14]. The
first interstellar meteor discovered (IM1) had a diameter of about $0.45$m
[15], and the second interstellar meteor (IM2) displayed similar physical
characteristics and was about twice as large as IM1 [13]. However, both
interstellar meteors were detected upon entering the atmosphere of the Earth
and burning up, rather than being observed by sky survey equipment. Finally,
the interstellar comet 2I/Borisov was detected in 2019 and was imaged by the
William Herschel Telescope and Gemini North Telescope [16]. Further detections
would enable an improved characterization of the ISO population, including the
distribution of their sizes, albedos, and densities. Such insights could allow
us to place further constraints on our explanations for the origin of ISO
populations. The observed interstellar meteors have also demonstrated
interesting properties of scientific interest, including unusual material
strength [13]. Furthermore, small ISOs could be evidence of extraterrestrial
equipment, which may have gone undetected by previous sky surveys if the
objects did not remain near the earth for extended periods of time and had
relatively low albedos.

##  2 Methodology

Our calculation proceeds as follows. First, we determine the maximum distance
at which an ISO is detectable by LSST as a function of its diameter, albedo,
and phase angle. Then, we can determine the size of the observable region for
objects with a given diameter and albedo. Next, we calculate the distribution
of object chord lengths through the observable region and the distribution of
ISO velocities relative to the Earth to estimate the proportion of objects
crossing through the observable region that will be detected by LSST. Then, we
estimate the crossing rate of ISOs through the observable region based on
their inferred density to determine the rate at which we expect to detect such
objects. Finally, we show how the estimated detection rate of objects with
LSST varies by object diameter and albedo, where the rate of detection for
smaller objects is lower than larger objects despite the higher number density
of smaller objects. We also find that the detection of small ISOs is
sensitively dependent on our assumed albedo distribution, and the number of
detections drops significantly if most ISOs have an albedo near 0.06.

The distance out to which an object is detectable from Earth depends upon its
brightness, which we can calculate from its diameter, albedo, and phase angle.
The absolute magnitude of objects can be determined as,

$H=-5\log_{10}(\frac{D\sqrt{p_{v}}}{1329}),$ (1)

where $D$ is the diameter (km) and $p_{v}$ is the albedo [17].

The albedo distribution of ISOs is largely unknown because we lack a
significant quantity of previous observations. We consider the case of the ISO
population having a constant albedo of 0.06. We also consider the case of the
albedo distribution of ISOs being similar to near earth asteroids. To create
the latter distribution, for each object in our sample, we draw an albedo from
the bimodal Rayleigh distribution described in Wright et al. [18], where,

$p(p_{v})=f_{D}\frac{p_{v}e^{-p_{v}^{2}/2d^{2}}}{d^{2}}+(1-f_{D})\frac{p_{v}e^{-p_{v}^{2}/2b^{2}}}{b^{2}}.$
(2)

We use the same parameters which were found by Wright et al. [18] where
$f_{D}=0.253$, $d=0.03$, and $b=0.168$. The distribution reflects that the
Wide-field Infrared Survey Explorer (WISE) found two populations of asteroids
classified by albedo where $d$ is the peak albedo of the darker population,
$b$ is the peak albedo of the brighter population, and $f_{D}$ is the
proportion of asteroids in the darker population [18]. We assume that small
interstellar objects are not comets because comets of small sizes entering the
solar system would evaporate before they could make a close approach to Earth
[19]. Thus, cometary tails would be more easily detectable at large distances
and are ignored here.

We can calculate the absolute magnitude of ISOs if we know their albedo and
diameter. To determine the apparent magnitudes of interstellar meteors, we use
a similar process to Seligman and Laughlin [9] and Hoover et al. [10]. The
apparent magnitude can be calculated from the absolute magnitude as,

$m=H+2.5\log_{10}(\frac{d_{BS}^{2}d_{BO}^{2}}{p(\Theta)d_{OS}^{4}}),$ (3)

where $d_{BS}$ is the distance from the body to the Sun, $d_{BO}$ is the
distance from the body to the observer—in this case Earth—and $d_{OS}$ is the
distance from the observer to the Sun. The phase integral, $p(\Theta)$, used
in equation (3) is calculated as,

$p(\Theta)=\frac{2}{3}((1-\frac{\Theta}{\pi})\cos(\Theta)+\frac{1}{\pi}\sin\Theta),$
(4)

where the phase angle, $\Theta$, can be calculated from the previously defined
distances where,

$\cos\Theta=\frac{d_{BO}^{2}+d_{BS}^{2}-d_{OS}^{2}}{2d_{BO}d_{BS}}.$ (5)

For a given phase angle, $\Theta$, we can then calculate the maximum distance
at which an object of a given absolute magnitude would be observable. We can
then approximate the total cross-sectional area of the observable region with
the ecliptic as,

$A=\frac{\pi}{N}\sum\limits_{n}(d_{BO}^{n})^{2},$ (6)

where $N$ is the number of samples, $n\in[0,\pi]$ are the phase angles and are
uniformly spread across the entire interval, and $d_{BO}^{n}$ is the maximum
distance at which an object is observable at phase angle $n$.

Figure 1 shows the dependence of the maximum distance from the Earth at which
an object is detectable on the phase angle for objects of various diameters,
given an albedo of 0.06. Figure 2 shows the cross-section of the detectable
region of objects with the ecliptic for objects with various diameters and an
albedo of 0.06.

Figure 1: The maximum distance at which an ISO with an albedo of 0.06 can be
detected as a function of the phase angle between the Sun, object, and Earth
for various ISO diameters. Figure 2: Cross-sectional area of the detection
region for ISOs with an albedo of 0.06 in the ecliptic plane for various
diameters. The plot is centered on Earth where the radial direction (x-axis)
points away from the Sun and the perpendicular direction (y-axis) is
tangential to the Earth’s position vector relative to the Sun.

Next, we consider the velocity distribution of ISOs with respect to Earth. We
assume that ISOs have the same velocity dispersion as local stars around the
LSR of $(\sigma_{U},\sigma_{V},\sigma_{W})=(33\pm 4,38\pm 4,23\pm 2)$ km s-1,
and we generate a population of 100,000 ISOs according to this distribution.
The velocity distribution of the objects does not depend on the object albedo
or diameter.

The velocity distribution of the objects relative to the Earth depends on the
direction of the Earth’s velocity with respect to the Sun. For each ISO, we
draw a corresponding velocity vector for the Earth by selecting a random
position within its orbit. The velocity of the Earth for each draw is given in
ecliptic coordinates as
$(v_{x},v_{y},v_{z})=(v_{\oplus}\cos\theta,v_{\oplus}\sin\theta,0)$ where
$\theta\in[0,2\pi]$. We also assume the magnitude of the Earth’s velocity is
constant at $v_{\oplus}=30$ km s-1.

We draw the velocity vector for the Earth in ecliptic coordinates relative to
the Sun for simplicity, but the velocity vectors for the population of
interstellar objects are given in galactocentric coordinates relative to the
LSR. As a result, we must convert the velocity of the Earth to the galactic
coordinate system relative to the LSR, where we use the procedure outlined in
McCabe [20]. Further details about this procedure are given in Appendix A.

As the object passes through the solar system, its velocity will increase
because of gravitational focusing from the Sun and the Earth. To account for
gravitational focusing from the Sun, we multiply the velocity vectors of the
interstellar objects by a factor of
$\sqrt{v_{\infty}^{2}+v_{esc}^{2}}/v_{\infty}$, where $v_{esc}$ is the escape
velocity from the solar system at 1 au from the Sun. The escape velocity is
calculated as,

$v_{esc}=\sqrt{2GM_{\odot}/d_{\oplus}},$ (7)

where $M_{\oplus}$ is mass of Sun and $d_{\oplus}$ = 1 au is the Earth-Sun
separation. We repeat the same procedure to calculate the effect of
gravitational focusing from the Earth. We find that effects from gravitational
focusing from the Earth are negligible because most objects in the population
are too far from the Earth for its gravitational focusing effect to be
significant, but we nevertheless include this adjustment in our calculation.
If a similar procedure to ours is used in the future to calculate the
detection rates of very small objects near the Earth by more sensitive
telescopes, the effect of gravitational focusing from the Earth may be more
significant. In Figure 3, we show the differential probability density
function (PDF) of the velocity relative to the Earth of our population of
ISOs.

Figure 3: The probability per logarithmic bin in velocity, where $v$ is the
velocity relative to the Earth for the sample population of ISOs. The minimum
velocity relative to the Earth is equivalent to the escape speed from the
Earth. We assume the velocity distribution is ISOs is similar to local stars
around the LSR and is independent of other variables, such as albedo or
diameter.

LSST will survey the sky about 100 times per year, resulting in a complete
survey of the sky every 3.65 days [4]. Objects must be detected by LSST three
times to form a track for the object to be recorded, so objects need to remain
within a detectable region for about 10.95 days. Since small ISOs must pass
close to Earth to be detectable, they may pass through the detectable region
in less than 10.95 days despite being within the lower tail for velocity
relative to the Earth. Nevertheless, a sufficiently high density of small ISOs
might allow for the detection of many such objects that have a relatively low
velocity and long chord length through the detectable region.

The amount of time that an object remains visible depends upon its chord
length through the observable region. However, since the maximum observable
distance depends upon the phase angle, the observable region is not spherical
(as shown in Figure 2). Furthermore, the size of the observable region is
different for each object since it depends on the object’s albedo. We can
generate a distribution of chord lengths from our population by randomly
assigning each object two points on the boundary of the observable region and
calculating the distance between them.

At each possible phase angle between $0^{\circ}<\Theta<180^{\circ}$, the
maximum distance of the observed object is fully constrained for a given
diameter and albedo. As a result, given the phase angle, the boundary of the
observable region is a circle perpendicular to the ecliptic plane and centered
on the vector passing through the center of the Earth and the Sun. The
boundary circle is a cross-section of the boundary of the observable region
formed by constraining the phase angle. For phase angles
$0^{\circ}<\Theta<90^{\circ}$, the object would be in opposition, so the
cross-sectional circle would be on the side of the Earth opposite the Sun. For
phase angles $90^{\circ}<\Theta<180^{\circ}$, the cross-sectional circle would
be between the Earth and the Sun. At each phase angle, the radius of the
cross-sectional circle also varies, where it has a radius of $0$ when
$\Theta=0^{\circ}$ or $\Theta=180^{\circ}$. Using the law of sines, we can
calculate the radius of the cross-sectional circle given the phase angle as,

$r=\frac{d_{bo}d_{bs}\sin(\Theta)}{d_{os}}.$ (8)

When we draw our random sample of boundary points to generate chords through
the region and the distribution of chord lengths, we select each phase angle
with a probability in proportion to the radius of the cross-sectional circle
which demarcates the boundary of the observable region at that phase angle. We
calculate a list of $N$ radii which each corresponds to a separate phase
angle, $n$, which are uniformly distributed over the interval
$0^{\circ}<\Theta<180^{\circ}$. We let $r_{n}$ be the radius of the circle
demarcating the boundary of the observable region at phase angle $n$. Figure 4
shows the relationship between the radius of the observable region for ISOs
and the phase angle. Then, the probability that an interstellar object enters
or leaves the observable region at a phase angle closest to $\Theta_{n}$ for
any $n$ is,

$P(\Theta_{n})=\frac{r_{n}}{\sum\limits_{n}r_{n}}.$ (9)

We draw phase angles from the resulting probability distribution and assign
one to each boundary point. Given the phase angle, we can calculate the value
$d_{bo}$, which is constrained since the object is at the maximum observable
distance for the phase angle. Each boundary point we draw must also be
randomly located on the cross-sectional circle corresponding to its phase
angle. Thus, for each boundary point, we also randomly draw another angle,
$\phi$, from the uniform distribution $\sim U[0,2\pi]$, which corresponds to
the angle of the object within the cross-sectional circle. Then, we then let
the position of the boundary point in ecliptic coordinates be
$(d_{bo}^{2}-r_{n}^{2},r\cos\phi,r\sin\phi)$. Here, the x-direction is the
Sun-Earth vector along the ecliptic, the y-direction is the perpendicular
vector within the ecliptic plane, and the z-direction is perpendicular to the
ecliptic plane.

Figure 4: Radius of the cross-sectional circle which demarcates the boundary
of the region in which ISOs can be detected. We plot the phase angle on the
x-axis, assume an albedo of 0.06, and plot the curve for three separate
diameters. Figure 5: The probability per logarithmic bin in chord length,
where $l$ is the chord length through the observable region, for the
distribution of chord lengths for ISO populations with various diameters. Two
albedo distributions are depicted: the albedo distribution from equation (2)
on the left and a population with a constant albedo of 0.06 on the right.

We choose pairs of two boundary points and determine the chord length through
the observable region between them. We pair each chord with an ISO in a sample
population and determine its distance covered relative to the Earth over a
period of 11 days. If the total distance covered is less than the chord
length, we consider the object to meet the criteria for detectability. Figure
5 shows the differential PDF of chord lengths for ISO populations with various
diameters and albedo distributions. The albedo affects the distance of points
on the boundary region from the Earth, where objects with higher albedos are
visible at greater distances.

Next, we use past detection rates of ISOs to calibrate our estimate for the
number density of ISOs. ISOs on the order of $\leq 1$m have been detected as
the IM1 and IM2 fireballs burning up in the atmosphere of the Earth. If we
assume that we detect all ISOs which enter the atmosphere of the Earth, we can
estimate the passing rate of ISOs of similar diameters to those objects. Then,
we can calculate the number density of such objects as,

$n\simeq\frac{\Gamma}{v_{\infty}\pi
R_{\oplus}^{2}[1+(v_{esc}/v_{\infty})^{2}]},$ (10)

which was demonstrated in Siraj and Loeb [13]. The implied rates of detection
for IM1 and IM2 are $\Gamma\sim 0.1$ yr-1 [15, 13]. The resulting number
density for IM1-like and IM2-like objects is about $n=2\times 10^{6}$ au-3.
The corresponding local mass density for IM1-like objects is 1.2 $M_{\oplus}$
pc-3, which is an order of magnitude lower than the estimated local mass
density for IM2. We use the given value as the estimate for mass density of
IM1-like objects, and we assume a power law size distribution of ISOs with an
equal total mass of objects per object log mass bin. This assumption is an
approximate estimate that holds across a large distribution of ISO masses
[21]. IM1 has a diameter of about 0.45m, which we use to constrain the
relationship [15]. We assume that ISOs have a similar density independent of
size, so we consider the resulting relation $m\sim D^{3}$. Let $\Gamma_{IM1}$
be our empirical detection rate of IM1-like interstellar meteors determined to
which we calibrate our estimate. Then, our rate of detection of interstellar
meteors of a given diameter $D$ with LSST can be calculated as,

$\Gamma_{D}=\Gamma_{IM1}(\frac{D}{45\text{cm}})^{-3}(\frac{A_{D}}{A_{\oplus}})p_{D},$
(11)

where $A_{D}$ is the cross-sectional area of the observable region of
interstellar objects of diameter $D$ with LSST, $A_{\oplus}$ is the cross-
sectional area of the Earth, and $p_{D}$ is the detection rate of interstellar
objects with diameter $D$. We can use this relationship to estimate the LSST
detection rate of interstellar objects with various diameters. However, the
estimate becomes less accurate for larger objects, such as those which are
‘Oumuamua-sized ($200$m), because we are extrapolating farther from our
observed number density of smaller objects. For these objects, the detection
rate of an object larger than IM1, such as ‘Oumuamua, can be used instead of
$\Gamma_{IM1}$. However, we do not concern ourselves with ’Oumuamua-sized
objects for the purpose of this analysis of small ISOs.

There are also further limitations of the viewing capabilities and survey
strategies of LSST that will lead to the detection of less interstellar
objects. For example, LSST will only detect interstellar objects within its
field-of-view (FoV), which is determined by its survey strategy. The Wide Fast
Deep survey (WFD), which will take up about 90 percent of the observing time
for LSST, will include declinations from $-60^{\circ}$ to $5^{\circ}$,
including all right ascension values [7]. Hence, LSST will primarily detect
objects that fall within this range, which is about $18,000\text{ deg}^{2}$
[22]. We update our estimated number of detections by the proportion of the
sky within the FoV of the WFD survey.

Above, we outlined a procedure to calculate the estimated number of ISO
detections with LSST given an object diameter, velocity distribution, and
albedo distribution. We repeat this calculation using various diameters and
our two albedo distributions for ISOs. We consider objects between 1m and 50m
in diameter, and we find that the estimated number of objects detected grows
exponentially with diameter in this range. This is consistent with our
expectation that smaller objects will not be detectable by LSST because they
will not remain within the observable region for a sufficient amount of time.
For an albedo distribution based on equation (2), we estimate the detection of
about $5-6$ ISOs with a diameter $\leq 50$m in the decade of LSST’s operation.
Given a constant albedo of 0.06 throughout the population, we estimate the
detection of $1-2$ ISOs with a diameter $\leq 50$m in the same period of
operation.

Our estimates for the number of ISOs detected by LSST of various diameters,
given an albedo distribution, are shown in Table 1. Our estimates are
relatively aligned with those of Hoover et al. [10], which estimated about 15
’Oumuamua-sized objects ($\sim 200$-m diameter) in LSST’s decade of operation.
Our estimate is about twice the estimate of Cook et al. [12] for a constant
albedo population, which may reflect increased ISO density estimates following
the detection of the first four ISOs.

Diameter (m) Detection Rate (Wright et al. [18] albedo distribution) Detection
Rate (albedo = 0.06) $1$ $0$ $0$ $5$ $0.1$ $0$ $10$ $0.42$ $0$ $20$ $3.40$
$0.57$ $30$ $7.63$ $2.03$ $40$ $12.45$ $3.94$ $50$ $15.53$ $6.17$

Table 1: The estimated detection rate (objects per decade) for small
interstellar objects by LSST based on our simulated results.

Some of our estimates for smaller ISOs might be overestimates because of the
trailing loss—the objects spend most of their time in the observable region in
very close proximity to the Earth. This effect is less substantial for larger
objects because their angular velocity through the sky is small in most of the
observable region, which is sufficiently far from the Earth. However, such
objects may become undetectable because of trailing loss when very close to
the Earth—if one of the three observations occurs at such a time, a track of
the object may not be formed.

Figure 6: Probability density functions of albedo distributions of entire
sample interstellar object population (left) and detected interstellar object
populations for various diameters (right). The populations are constructed
from albedo distributions which resemble that of near-earth asteroids.

In Figure 6, we show the differential probability for the albedos of detected
ISOs within the sample population that follows the albedo distribution from
equation (2), and we also show the differential probability of the original
sample population. We find that the high albedo population allows for
significantly more detections, especially for smaller ISOs. If the true albedo
distribution of ISOs were similar to that of near earth asteroids as suggested
by this sample population, we would likely have seen more ISO detections with
previous telescopes, suggesting this is unlikely to be the albedo distribution
of interstellar objects. However, interstellar objects are likely to originate
from different populations, and we cannot rule out that some populations may
contain objects with higher albedos like the near earth asteroids population.
Further detections of ISOs can allow us to produce better estimates of their
albedo distribution, including distinguishing between populations.

We also consider the extent to which the detection rate is impacted by changes
to the sensitivity of the telescope. In Figure 7, we show the impact of
improvements in telescope sensitivity on the detection rate for objects of
various diameters. A telescope that can detect all objects where the limiting
magnitude $m\leq 25$ can detect approximately an order of magnitude more
objects in the $50$ m-diameter range than LSST, which can only detect objects
where $m\leq 24$.

Figure 7: Detection rates of ISOs of various diameters given the limiting
magnitude (higher numbers are fainter) at which an object can be detected.
LSST has a minimum magnitude of $m=24$ (shown in green). Our plot shows
estimates for the sample population where all objects have an albedo of 0.06.

##  3 Conclusion

Our calculations forecast that LSST will be able to detect a small ISO in the
diameter range of 1-50m every 1-2 years. Extending the sensitivity of the
survey by one magnitude will increase this number by a factor of ten.

While solar system surveys have historically been of scientific and political
interest because of the desire to catalog the NEO and PHA
populations—including for planetary protection—the emerging ability to detect
ISOs and smaller objects will lead to new scientific research opportunities.
Furthermore, we expect future telescopes to continue the trend of enabling the
detection of even smaller objects.

In addition to ISOs, very sensitive telescopes may be able to detect other
objects of economic or commercial interest, including tracking probes in orbit
around celestial bodies and in transit through the solar system. Telescopes
may also be more equipped to detect rovers and probes on commercial or
scientific missions throughout the solar system, small asteroids of interest
for a future asteroid mining industry, and weather events or future human
space settlement activities on celestial bodies. In the very long-term, space-
based astronomical instrumentation may transform geopolitics and scientific
research through the ability to make observations from extremely small probes
to detailed accounts of rare events [23]. When assessing the ability of our
instruments to detect small objects besides ISOs, the same assumptions cannot
be made about the albedos and position and velocity distributions of such
objects, but a similar methodology to our own can be used to determine whether
the viewing capabilities of an astronomical survey are sufficient for
detection.

##  Appendix A Velocity Conversion

To compare the velocity vectors of the ISOs to the velocity vector of the
Earth, we need to convert the velocity vector of the Earth into galactic
coordinates. In addition, the velocity vector of the Earth is given relative
to the Sun, and the velocity vector of the ISOs are given relative to the LSR.
Thus, we also subtract the peculiar velocity of the Sun from the velocity
vector of the ISOs so that we can compare the object’s velocity to that of
Earth. We use the procedure outlined in McCabe [20].

First, we convert the velocity vector from ecliptic to equatorial coordinates
using the rotation matrix,

$R=\begin{bmatrix}1&0&0\\\ 0&\cos\epsilon&-\sin\epsilon\\\
0&\sin\epsilon&\cos\epsilon\end{bmatrix}$

where $\vec{v}_{\oplus,eq}=R\vec{v}_{\oplus,ec}$ and is the obliquity of the
ecliptic plane. Then, we convert the velocity vector from equatorial to
galactic coordinates through the relation
$\vec{v}_{\oplus,gal}=M\vec{v}_{\oplus,eq}$ where,

$M_{11}=-\sin l_{CP}\sin\alpha_{GP}-\cos l_{CP}\cos\alpha_{GP}\sin\delta_{GP}$
$M_{12}=\sin l_{CP}\cos\alpha_{GP}-\cos l_{CP}\sin\alpha_{GP}\sin\delta_{GP}$
$M_{13}=\cos l_{CP}\cos\delta_{GP}$ $M_{21}=\cos l_{CP}\sin\alpha_{GP}-\sin
l_{CP}\cos\alpha_{GP}\sin\delta_{GP}$ $M_{22}=-\cos l_{CP}\cos\alpha_{GP}-\sin
l_{CP}\sin\alpha_{GP}\sin\delta_{GP}$ $M_{23}=\sin l_{CP}\cos\delta_{GP}$
$M_{31}=\cos\alpha_{GP}\cos\delta_{GP}$
$M_{32}=\sin\alpha_{GP}\cos\delta_{GP}$ $M_{33}=\sin\delta_{GP}$

Since we have the velocity of the Earth relative to the Sun, we then subtract
the peculiar vector of the Sun relative to the LSR
$\vec{v}_{pec}=(v_{U},v_{V},v_{W})=(10\pm 1,11\pm 2,7\pm 0.5)$ km s-1 from the
velocity vectors of the interstellar meteors to find the objects’ velocity
vectors relative to the Sun in the galactic coordinate system,

$\vec{v}_{o,Sun}=\vec{v}_{o,LSR}-\vec{v}_{pec},$ (12)

where $v_{o,Sun}$ is the object’s velocity relative to the Sun, $v_{o,LSR}$ is
the object’s velocity relative to the LSR, and $v_{pec}$ is the velocity of
the Sun relative to the LSR from the velocity vectors of the ISO population
drawn from the LSR. Then, we can calculate the velocities of the ISOs relative
to the Earth as,

$v_{e}=||\vec{v}_{o,Sun}-\vec{v}_{\oplus,gal}||.$ (13)

## References

  * United States Code [2022] United States Code, 51 USC Ch. 711: Near-Earth Objects, 2022. URL: https://uscode.house.gov/view.xhtml?path=/prelim@title51/subtitle7/chapter711&edition=prelim.
  * Larson [2007] S. Larson,  Current NEO surveys 236 (2007) 323–328. URL: https://ui.adsabs.harvard.edu/abs/2007IAUS..236..323L. doi:10.1017/S1743921307003390, conference Name: Near Earth Objects, our Celestial Neighbors: Opportunity and Risk ADS Bibcode: 2007IAUS..236..323L.
  * Denneau et al. [2013] L. Denneau, R. Jedicke, T. Grav, M. Granvik, et al.,  The Pan-STARRS Moving Object Processing System,  Publications of the Astronomical Society of the Pacific 125 (2013) 357. URL: https://dx.doi.org/10.1086/670337. doi:10.1086/670337, publisher: University of Chicago Press.
  * Ivezic et al. [2006] Z. Ivezic, J. A. Tyson, M. Juric, J. Kubica, A. Connolly, F. Pierfederici, A. W. Harris, E. Bowell, t. L. Collaboration,  LSST: Comprehensive NEO Detection, Characterization, and Orbits,  Proceedings of the International Astronomical Union 2 (2006) 353–362. URL: http://arxiv.org/abs/astro-ph/0701506. doi:10.1017/S1743921307003420, arXiv:astro-ph/0701506.
  * Grav et al. [2016] T. Grav, A. K. Mainzer, T. Spahr,  Modeling the Performance of the LSST in Surveying the Near-Earth Object Population,  The Astronomical Journal 151 (2016) 172. URL: https://dx.doi.org/10.3847/0004-6256/151/6/172. doi:10.3847/0004-6256/151/6/172, publisher: The American Astronomical Society.
  * Vereš and Chesley [2017] P. Vereš, S. R. Chesley,  High-fidelity Simulations of the Near-Earth Object Search Performance of the Large Synoptic Survey Telescope,  The Astronomical Journal 154 (2017) 12. URL: https://dx.doi.org/10.3847/1538-3881/aa73d1. doi:10.3847/1538-3881/aa73d1, publisher: The American Astronomical Society.
  * Jones et al. [2018] R. L. Jones, C. T. Slater, J. Moeyens, L. Allen, T. Axelrod, K. Cook, Z. Ivezic, M. Juric, J. Myers, C. E. Petry,  The Large Synoptic Survey Telescope as a Near-Earth Object discovery machine,  Icarus 303 (2018) 181–202. URL: https://www.sciencedirect.com/science/article/pii/S0019103517301434. doi:10.1016/j.icarus.2017.11.033.
  * Ivezic et al. [2019] Z. Ivezic, S. M. Kahn, J. A. Tyson, B. Abel, et al.,  LSST: From Science Drivers to Reference Design and Anticipated Data Products,  The Astrophysical Journal 873 (2019) 111. URL: https://ui.adsabs.harvard.edu/abs/2019ApJ...873..111I. doi:10.3847/1538-4357/ab042c, aDS Bibcode: 2019ApJ…873..111I.
  * Seligman and Laughlin [2018] D. Seligman, G. Laughlin,  The Feasibility and Benefits of In Situ Exploration of ‘Oumuamua-like Objects,  The Astronomical Journal 155 (2018) 217. URL: https://dx.doi.org/10.3847/1538-3881/aabd37. doi:10.3847/1538-3881/aabd37, publisher: The American Astronomical Society.
  * Hoover et al. [2022] D. J. Hoover, D. Z. Seligman, M. J. Payne,  The Population of Interstellar Objects Detectable with the LSST and Accessible for In Situ Rendezvous with Various Mission Designs,  The Planetary Science Journal 3 (2022) 71. URL: https://iopscience.iop.org/article/10.3847/PSJ/ac58fe/meta. doi:10.3847/PSJ/ac58fe, publisher: IOP Publishing.
  * Engelhardt et al. [2017] T. Engelhardt, R. Jedicke, P. Vereš, A. Fitzsimmons, L. Denneau, E. Beshore, B. Meinke,  An Observational Upper Limit on the Interstellar Number Density of Asteroids and Comets,  The Astronomical Journal 153 (2017) 133. URL: https://dx.doi.org/10.3847/1538-3881/aa5c8a. doi:10.3847/1538-3881/aa5c8a, publisher: The American Astronomical Society.
  * Cook et al. [2016] N. V. Cook, D. Ragozzine, M. Granvik, D. C. Stephens,  Realistic Detectability of Close Interstellar Comets,  The Astrophysical Journal 825 (2016) 51. URL: https://dx.doi.org/10.3847/0004-637X/825/1/51. doi:10.3847/0004-637X/825/1/51, publisher: The American Astronomical Society.
  * Siraj and Loeb [2022] A. Siraj, A. Loeb, Interstellar Meteors are Outliers in Material Strength, 2022. URL: http://arxiv.org/abs/2209.09905, arXiv:2209.09905 [astro-ph].
  * Meech et al. [2017] K. J. Meech, R. Weryk, M. Micheli, J. T. Kleyna, et al.,  A brief visit from a red and extremely elongated interstellar asteroid,  Nature 552 (2017) 378–381. URL: https://www.nature.com/articles/nature25020. doi:10.1038/nature25020, number: 7685 Publisher: Nature Publishing Group.
  * Siraj and Loeb [2022] A. Siraj, A. Loeb, The 2019 Discovery of a Meteor of Interstellar Origin, 2022. URL: http://arxiv.org/abs/1904.07224, arXiv:1904.07224 [astro-ph].
  * Guzik et al. [2020] P. Guzik, M. Drahus, K. Rusek, W. Waniak, G. Cannizzaro, I. Pastor-Marazuela,  Initial characterization of interstellar comet 2I/Borisov,  Nature Astronomy 4 (2020) 53–57. URL: https://www.nature.com/articles/s41550-019-0931-8. doi:10.1038/s41550-019-0931-8, number: 1 Publisher: Nature Publishing Group.
  * Bowell et al. [1989] E. Bowell, B. Hapke, D. Domingue, K. Lumme, J. Peltoniemi, A. Harris,  Application of photometric models to asteroids (1989).
  * Wright et al. [2016] E. L. Wright, A. Mainzer, J. Masiero, T. Grav, J. Bauer,  The Albedo Distribution of Near Earth Asteroids,  The Astronomical Journal 152 (2016) 79. URL: https://dx.doi.org/10.3847/0004-6256/152/4/79. doi:10.3847/0004-6256/152/4/79, publisher: The American Astronomical Society.
  * Levison et al. [2002] H. F. Levison, A. Morbidelli, L. Dones, R. Jedicke, P. A. Wiegert, W. F. Bottke,  The Mass Disruption of Oort Cloud Comets,  Science 296 (2002) 2212–2215. URL: https://ui.adsabs.harvard.edu/abs/2002Sci...296.2212L. doi:10.1126/science.1070226, aDS Bibcode: 2002Sci…296.2212L.
  * McCabe [2014] C. McCabe,  The Earth’s velocity for direct detection experiments,  Journal of Cosmology and Astroparticle Physics 2014 (2014) 027. URL: https://ui.adsabs.harvard.edu/abs/2014JCAP...02..027M. doi:10.1088/1475-7516/2014/02/027, aDS Bibcode: 2014JCAP…02..027M.
  * Siraj and Loeb [2021] A. Siraj, A. Loeb,  Interstellar objects outnumber Solar system objects in the Oort cloud,  Monthly Notices of the Royal Astronomical Society: Letters 507 (2021) L16–L18. URL: https://doi.org/10.1093/mnrasl/slab084. doi:10.1093/mnrasl/slab084.
  * Bianco et al. [2021] F. B. Bianco, Z. Ivezic, R. L. Jones, M. L. Graham, et al.,  Optimization of the Observing Cadence for the Rubin Observatory Legacy Survey of Space and Time: A Pioneering Process of Community-focused Experimental Design,  The Astrophysical Journal Supplement Series 258 (2021) 1. URL: https://dx.doi.org/10.3847/1538-4365/ac3e72. doi:10.3847/1538-4365/ac3e72, publisher: The American Astronomical Society.
  * Scherrer and Loeb [2023] R. J. Scherrer, A. Loeb,  Ultra long-term cosmology and astrophysics,  New Astronomy 99 (2023) 101940. URL: https://ui.adsabs.harvard.edu/abs/2023NewA...9901940S. doi:10.1016/j.newast.2022.101940, aDS Bibcode: 2023NewA…9901940S.


Non-Bayesian particle filters

Alexandre J. Chorin and Xuemin Tu

Department of Mathematics, University of California at Berkeley

and

Lawrence Berkeley National Laboratory

Berkeley, CA, 94720

###### Abstract

Particle filters for data assimilation in nonlinear problems use “particles”
(replicas of the underlying system) to generate a sequence of probability
density functions (pdfs) through a Bayesian process. This can be expensive
because a significant number of particles has to be used to maintain accuracy.
We offer here an alternative, in which the relevant pdfs are sampled directly
by an iteration. An example is discussed in detail.

Keywords particle filter, chainless sampling, normalization factor, iteration,
non-Bayesian

##  1 Introduction.

There are many problems in science in which the state of a system must be
identified from an uncertain equation supplemented by a stream of noisy data
(see e.g. [1]). A natural model of this situation consists of a stochastic
differential equation (SDE):

$d{\bf x}={\bf f}({\bf x},t)\,dt+g({\bf x},t)\,d{\bf w},$ (1)

where ${\bf x}=(x_{1},x_{2},\dots,x_{m})$ is an $m$-dimensional vector, $d{\bf
w}$ is $m$-dimensional Brownian motion, ${\bf f}$ is an $m$-dimensional vector
function, and $g$ is a scalar (i.e., an $m$ by $m$ diagonal matrix of the form
$gI$, where $g$ is a scalar and $I$ is the identity matrix). The Brownian
motion encapsulates all the uncertainty in this equation. The initial state
${\bf x}(0)$ is assumed given and may be random as well.

As the experiment unfolds, it is observed, and the values ${\bf b}^{n}$ of a
measurement process are recorded at times $t^{n}$; for simplicity assume
$t^{n}=n\delta$, where $\delta$ is a fixed time interval and $n$ is an
integer. The measurements are related to the evolving state ${\bf x}(t)$ by

${\bf b}^{n}={\bf h}({\bf x}^{n})+G{\bf W}^{n},$ (2)

where ${\bf h}$ is a $k$-dimensional, generally nonlinear, vector function
with $k\leq m$, $G$ is a diagonal matrix, ${\bf x}^{n}={\bf x}(n\delta)$, and
${\bf W}^{n}$ is a vector whose components are independent Gaussian variables
of mean 0 and variance 1, independent also of the Brownian motion in equation
(1). The task is to estimate ${\bf x}$ on the basis of equation (1) and the
observations (2).

If the system (1) is linear and the data are Gaussian, the solution can be
found via the Kalman-Bucy filter. In the general case, it is natural to try to
estimate ${\bf x}$ as the mean of its evolving probability density. The
initial state ${\bf x}$ is known and so is its probability density; all one
has to do is evaluate sequentially the density $P_{n+1}$ of ${\bf x}^{n+1}$
given the probability density $P_{n}$ of ${\bf x}^{n}$ and the data ${\bf
b}^{n+1}$. This can be done by following “particles” (replicas of the system)
whose empirical distribution approximates $P_{n}$. In a Bayesian filter (see
e.g [2, 3, 4, 5, 6, 7, 8, 9], one uses the pdf $P_{n}$ and equation (1) to
generate a prior density, and then one uses the new data ${\bf b}^{n+1}$ to
generate a posterior density $P_{n+1}$. In addition, one may have to sample
backward to take into account the information each measurement provides about
the past and avoid having too many identical particles. Evolving particles is
typically expensive, and the backward sampling, usually done by Markov chain
Monte Carlo (MCMC), can be expensive as well, because the number of particles
needed can grow catastrophically (see e.g. [10]).

In this paper we offer an alternative to the standard approach, in which
$P_{n+1}$ is sampled directly without recourse to Bayes’ theorem and backward
sampling, if needed, is done by chainless Monte Carlo [11]. Our direct
sampling is based on a representation of a variable with density $P_{n+1}$ by
a collection of functions of Gaussian variables parametrized by the support of
$P_{n}$, with parameters found by iteration. The construction is related to
chainless sampling as described in [11]. The idea in chainless sampling is to
produce a sample of a large set of variables by sequentially sampling a
growing sequence of nested conditionally independent subsets. As observed in
[12, 13], chainless sampling for a SDE reduces to interpolatory sampling, as
explained below. Our construction will be explained in the following sections
through an example where the position of a ship is deduced from the
measurements of an azimuth, already used as a test bed in [6, 14, 15].

##  2 Sampling by interpolation and iteration.

First we explain how to sample via interpolation and iteration in a simple
example, related to the example and the construction in [12]. Consider the
scalar SDE

$dx=f(x,t)dt+\sqrt{\sigma}\,dw;$ (3)

we want to find sample paths $x=x(t),0\leq t\leq 1$, subject to the conditions
$x(0)=0,x(1)=X$.

Let $N(a,v)$ denote a Gaussian variable with mean $a$ and variance $v$. We
first discretize equation (3) on a regular mesh $t^{0},t^{1},\dots,t^{N}$,
where $t^{n}=n\delta$, $\delta=1/N$, $0\leq n\leq N$, with $x^{n}=x(t^{n})$,
and, following [12], use a balanced implicit discretization [16, 17]:

$x^{n+1}=x^{n}+f(x^{n},t^{n})\delta+(x^{n+1}-x^{n})f^{\prime}(x^{n})\delta+W^{n+1},$

where $f^{\prime}(x^{n},t^{n})=\frac{\partial f}{\partial x^{n}}(x^{n},t^{n})$
and $W^{n+1}$ is $N(0,\sigma/N)$. The joint probability density of the
variables $x^{1},\dots,x^{N-1}$ is $Z^{-1}\exp(-\sum_{0}^{N}V^{i})$, where $Z$
is the normalization constant and

$\displaystyle V_{i}$ $\displaystyle=\frac{\left((1-\delta
f^{\prime})(x^{n+1}-x^{n})-\delta f\right)^{2}}{2\sigma\delta}$
$\displaystyle=\frac{\left(x^{n+1}-x^{n}-\delta f/(1-\delta
f^{\prime})\right)^{2}}{2\sigma_{n}},$

where $f,f^{\prime}$ are functions of the $x^{j}$, and
$\sigma_{n}=\sigma\delta/(1-\delta f^{\prime})^{2}$ (see [18]). One can obtain
sample solutions by sampling this density, e.g. by MCMC, or one can obtain
them by interpolation (chainless sampling), as follows.

Consider first the special case $f(x,t)=f(t)$, so that in particular
$f^{\prime}=0$. Each increment $x^{n+1}-x^{n}$ is now a $N(a_{n},\sigma/N)$
variable, with the $a_{n}=f(t^{n})\delta$ known explicitly. Let $N$ be a power
of $2$. Consider the variable $x^{N/2}$. On one hand,

$x^{N/2}=\sum_{1}^{N/2}(x^{n}-x^{n-1})=N(A_{1},V_{1}),$

where $A_{1}=\sum_{1}^{N/2}a_{n},V_{1}=\sigma/2$. On the other hand,

$X=x^{N/2}+\sum_{N/2+1}^{N}(x^{n}-x^{n-1}),$

so that

$x^{N/2}=N(A_{2},V_{2}),$

with

$A_{2}=X-\sum_{N/2+1}^{N-1}a_{n},\quad V_{2}=V_{1}.$

The pdf of $x^{N/2}$ is the product of the two pdfs; one can check that

$\displaystyle\exp\left(-\frac{(x-A_{1})^{2}}{2V_{1}}\right)\exp\left(-\frac{(x-A_{2})^{2}}{2V_{2}}\right)$
$\displaystyle=$
$\displaystyle\exp\left(-\frac{(x-\bar{a})^{2}}{2\bar{v}}\right)\exp(-\phi),$

where $\bar{v}=\frac{V_{1}V_{2}}{V_{1}+V_{2}}$,
$\bar{a}=\frac{V_{1}A_{1}+V_{2}A_{2}}{V_{1}+V_{2}}$, and
$\phi=\frac{(A_{2}-A_{1})^{2}}{2(V_{1}+V_{2})}$; $e^{-\phi}$ is the
probability of getting from the origin to $X$, up to a normalization constant.

Pick a sample $\xi_{1}$ from the $N(0,1)$ density; one obtains a sample of
$x^{N/2}$ by setting $x^{N/2}=\bar{a}+\sqrt{\bar{v}}\xi_{1}$. Given a sample
of $x^{N/2}$ one can similarly sample $x^{N/4},x^{3N/4}$, then $x^{N/8}$,
$x^{3N/8}$, etc., until all the $x^{j}$ have been sampled. If we define
${\bf\xi}=(\xi_{1},\xi_{2},\dots,\xi_{N-1})$, then for each choice of
${\bf\xi}$ we find a sample $(x^{1},\dots,x^{N-1})$ such that

$\displaystyle\exp\left(-\frac{\xi_{1}^{2}+\cdots+\xi_{N-1}^{2}}{2}\right)\exp\left(-\frac{(X-\sum_{n}a_{n})^{2}}{2\sigma}\right)$
$\displaystyle=$
$\displaystyle\exp\left(-\frac{(x^{1}-x^{0}-a_{0})^{2}}{2\sigma/N}-\frac{(x^{2}-x^{1}-a_{1})^{2}}{2\sigma/N}\right.$
$\displaystyle\left.-\dots-\frac{(x^{N}-x^{N-1}-a_{N-1})^{2}}{2\sigma/N}\right),$
(4)

where the factor $\exp\left(-\frac{(X-\sum_{n}a_{n})^{2}}{2\sigma}\right)$ on
the left is the probability of the fixed end value $X$ up to a normalization
constant. In this linear problem, this factor is the same for all the samples
and therefore harmless. One can repeat this sampling process for multiple
choices of the variables $\xi_{j}$; each sample of the corresponding set of
$x^{n}$ is independent of any previous samples of this set.

Now return to the general case. The functions $f$, $f^{\prime}$ are now
functions of the $x^{j}$. We obtain a sample of the probability density we
want by iteration. First pick $\Xi=(\xi_{1},\xi_{2},\dots,\xi_{N-1})$, where
each $\xi_{j}$ is drawn independently from the $N(0,1)$ density (this vector
remains fixed during the iteration). Make a first guess ${\bf
x}^{0}=(x^{1}_{0},x^{2}_{0},\dots,x^{N-1}_{0})$ (for example, if $X\neq 0$,
pick ${\bf x}=0$). Evaluate the functions $f,f^{\prime}$ at ${\bf x}^{j}$
(note that now $f^{\prime}\neq 0$, and therefore the variances of the various
increments are no longer constants). We are back in previous case, and can
find values of the increments $x^{n+1}_{j+1}-x^{n}_{j+1}$ corresponding to the
values of $f,f^{\prime}$ we have. Repeat the process starting with the new
iterate. If the vectors ${\bf x}^{j}$ converge to a vector ${\bf
x}=(x^{1},\dots,x^{N-1})$, we obtain, in the limit, equation (4), where now on
the right side $\sigma$ depends on $n$ so that $\sigma=\sigma_{n}$, and both
$a_{n},\sigma_{n}$ are functions of the final ${\bf x}$. The left hand side of
(4) becomes:

$\exp\left(-\frac{\xi_{1}^{2}+\cdots+\xi_{N-1}^{2}}{2}\right)\exp\left(-\frac{(X-\sum_{n}a_{n})^{2}}{2\sum_{n}\sigma_{n}}\right).$

Note that now the factor
$\exp\left(-\frac{(X-\sum_{n}a_{n})^{2}}{2\sum_{n}\sigma_{n}}\right)$ is
different from sample to sample, and changes the relative weights of the
different samples. In averaging, one should take this factor as weight, or
resample as described at the end of the following section. In order to obtain
more uniform weights, one also can use the strategies in [11, 12].

One can readily see that the iteration converges if $KTM<1$, where $K$ is the
Lipshitz constant of $f$, $T$ is the length of the interval on which one works
(here $T=1$), and $M$ is the maximum norm of the vectors ${\bf x}^{j+1}-{\bf
x}^{j}$. If this inequality is not satisfied for the iteration above, it can
be re-established by a suitable underrelaxation. One should course choose $N$
large enough so that the results are converged in $N$. We do not provide more
details here because they are extraneous to our purpose, which is to explain
chainless/interpolatory sampling and the use of reference variables in a
simple context.

##  3 The ship azimuth problem.

The problem we focus on is discussed in [6, 14, 15], where it is used to
demonstrate the capabilities of particular Bayesian filters. A ship sets out
from a point $(x_{0},y_{0})$ in the plane and undergoes a random walk,

$\displaystyle x^{n+1}$ $\displaystyle=x^{n}+dx^{n+1},$ $\displaystyle
y^{n+1}$ $\displaystyle=y^{n}+dy^{n+1},$ (5)

for $n\geq 0$, and with $x^{0}=y^{0}$ given, and $dx^{n+1}=N(dx^{n},\sigma)$,
$dy^{n+1}=N(dy^{n},\sigma)$, i.e., each displacement is a sample of a Gaussian
random variable whose variance $\sigma$ does not change from step to step and
whose mean is the value of the previous displacement. An observer makes noisy
measurements of the azimuth $\arctan(y^{n}/x^{n})$, recording

$b^{n}=\arctan\frac{y^{n}}{x^{n}}+N(0,s).$ (6)

where the variance $s$ is also fixed; here the observed quantity $b$ is scalar
and is not be denoted by a boldfaced letter. The problem is to reconstruct the
positions ${\bf x}^{n}=(x^{n},y^{n})$ from equations (5,6). We take the same
parameters as [6]: $x_{0}=0.01,y_{0}=20$, $dx^{1}=0.002$, $dy^{1}=-0.06$,
$\sigma=1\cdot 10^{-6},s=25\cdot 10^{-6}$. We follow numerically $M$
particles, all starting from $X_{i}^{0}=x_{0},Y_{i}^{0}=y_{0}$, as described
in the following sections, and we estimate the ship’s position at time
$n\delta$ as the mean of the locations ${\bf
X}^{n}_{i}=(X^{n}_{i},Y^{n}_{i}),i=1,\dots,M$ of the particles at that time.
The authors of [6] also show numerical results for runs with varying data and
constants; we discuss those refinements in section 6 below.

##  4 Forward step.

Assume we have a collection of $M$ particles ${\bf X}^{n}$ at time
$t^{n}=n\delta$ whose empirical density approximates $P_{n}$; now we find
increments $d{\bf X}^{n+1}$ such that the empirical density of ${\bf
X}^{n+1}={\bf X}^{n}+d{\bf X}^{n+1}$ approximates $P_{n+1}$. $P_{n+1}$ is
known implicitly: it is the product of the density that can be deduced from
the SDE and the one that comes from the observations, with the appropriate
normalization. If the increments were known, their probability $p$ (the
density $P_{n+1}$ evaluated at the resulting positions ${\bf X}^{n+1}$) would
be known, so $p$ is a function of $d{\bf X}^{n+1}$, $p=p(d{\bf X}^{n+1})$. For
each particle $i$, we are going to sample a Gaussian reference density, obtain
a sample of probability $\rho$, then solve (by iteration) the equation

$\rho=p(d{\bf X}^{n+1}_{i})$ (7)

to obtain $d{\bf X}_{i}^{n+1}$.

Define $f(x,y)=\arctan(y/x)$ and $f^{n}=f(X^{n},Y^{n})$. We are working on one
particle at a time, so the index $i$ can be temporarily suppressed. Pick two
independent samples $\xi_{x}$, $\xi_{y}$ from a $N(0,1)$ density (the
reference density in the present calculation), and set
$\rho=\frac{1}{2\pi}\exp\left(-\frac{\xi_{x}^{2}}{2}-\frac{\xi_{y}^{2}}{2}\right)$;
the variables $\xi_{x}$, $\xi_{y}$ remain unchanged until the end of the
iteration. We are looking for displacements $dX^{n+1}$, $dY^{n+1}$, and
parameters $a_{x},a_{y},v_{x},v_{y},\phi$, such that:

$\displaystyle 2\pi\rho=$
$\displaystyle\exp\left(-\frac{(dX^{n+1}-dX^{n})^{2}}{2\sigma}-\frac{(dY^{n+1}-dY^{n})^{2}}{2\sigma}\right.$
$\displaystyle\left.-\frac{(f^{n+1}-b^{n+1})^{2}}{2s}\right)\exp(\phi)$
$\displaystyle=$
$\displaystyle\exp\left(-\frac{(dX^{n+1}-a_{x})^{2}}{2v_{x}}-\frac{(dY^{n+1}-a_{y})^{2}}{2v_{y}}\right)$
(8)

The first equality states what we wish to accomplish: find increments
$dX^{n+1}$, $dY^{n+1}$, functions respectively of $\xi_{x},\xi_{y}$, whose
probability with respect to $P_{n+1}$ is $\rho$. The factor $e^{\phi}$ is
needed to normalize this term ($\phi$ is called below a “phase”). The second
equality says how the goal is reached: we are looking for parameters
$a_{x},a_{y},v_{x},v_{y},$ (all functions of ${\bf X}^{n}$) such that the
increments are samples of Gaussian variables with these parameters, with the
assumed probability. One should remember that in our example the mean of
$dX^{n+1}$ is $dX^{n}$, and similarly for $dY^{n+1}$. We are not representing
$P_{n+1}$ as a function of a single Gaussian- there is a different Gaussian
for every value of ${\bf X}^{n}$.

To satisfy the second equality we set up an iteration for vectors $d{\bf
X}^{n+1,j}(=d{\bf X}^{j}$ for brevity) that converges to $d{\bf X}^{n+1}$.
Start with $d{\bf X}^{0}=0$. We now explain how to compute $d{\bf X}^{j+1}$
given $d{\bf X}^{j}$.

Approximate the observation equation (6) by

$f({\bf
X}^{j})+f_{x}\cdot(dX^{j+1}-dX^{j})+f_{y}\cdot(dY^{j+1}-dY^{j})=b^{n+1}+N(0,s),$
(9)

where the derivatives $f_{x},f_{y}$ are, like $f$, evaluated at ${\bf
X}^{j}={\bf X}^{n}+d{\bf X}^{j}$, i.e., approximate the observation equation
by its Taylor series expansion around the previous iterate. Define a variable
$\eta^{j+1}=(f_{x}\cdot dX^{j+1}+f_{y}\cdot
dY^{j+1})/\sqrt{f_{x}^{2}+f_{y}^{2}}$. The approximate observation equation
says that $\eta^{j+1}$ is a $N(a_{1},v_{1})$ variable, with

$\displaystyle a_{1}$ $\displaystyle=-\frac{f-f_{x}\cdot dX^{j}-f_{y}\cdot
dY^{j}-b^{n+1}}{\sqrt{f_{x}^{2}+f_{y}^{2}}}\,,$ $\displaystyle v_{1}$
$\displaystyle=\frac{s}{f_{x}^{2}+f_{y}^{2}}.$ (10)

On the other hand, from the equations of motion one finds that $\eta^{j+1}$ is
$N(a_{2},v_{2})$, with $a_{2}=(f_{x}\cdot dX^{n}+f_{y}\cdot
dY^{n})/\sqrt{f_{x}^{2}+f_{y}^{2}}$ and $v_{2}=\sigma$. Hence the pdf of
$\eta^{j+1}$ is, up to normalization factors,

$\exp\left(-\frac{(x-a_{1})^{2}}{2v_{1}}-\frac{(x-a_{2})^{2}}{2v_{2}}\right)=\exp\left(-\frac{(x-\bar{a})^{2}}{2\bar{v}}\right)\exp(-\phi),$

where $\bar{v}=\frac{v_{1}v_{2}}{v_{1}+v_{2}}$,
$\bar{a}=\frac{a_{1}v_{1}+a_{2}v_{2}}{v_{1}+v_{2}}$,
$\phi=\frac{(a_{1}-a_{2})^{2}}{2(v_{1}+v_{2})}=\phi^{j+1}$.

We can also define a variable $\eta_{+}^{j+1}$ that is a linear combination of
$dX^{j+1}$, $dY^{j+1}$ and is uncorrelated with $\eta^{j+1}$:

$\eta_{+}^{j+1}=\frac{-f_{y}\cdot dY^{j+1}+f_{x}\cdot
dX^{j+1}}{\sqrt{f_{x}^{2}+f_{y}^{2}}}.$

The observations do not affect $\eta_{+}^{j+1}$, so its mean and variance are
known. Given the means and variances of $\eta^{j+1}$, $\eta^{j+1}_{+}$ one can
easily invert the orthogonal matrix that connects them to $dX^{j+1}$,
$dY^{j+1}$ and find the means and variances $a_{x},v_{x}$ of $dX^{j+1}$ and
$a_{y},v_{y}$ of $dY^{j+1}$ after their modification by the observation (the
subscripts on $a,v$ are labels, not differentiations). Now one can produce
values for $dX^{j+1},dY^{j+1}$:

$dX^{j+1}=a_{x}+\sqrt{v_{x}}\xi_{x},\quad dY^{j+1}=a_{y}+\sqrt{v_{y}}\xi_{y},$

where $\xi_{x}$, $\xi_{y}$ are the samples from $N(0,1)$ chosen at the
beginning of the iteration. This completes the iteration.

This iteration converges to ${\bf X}^{n+1}$ such that $f({\bf
X}^{n+1})=b^{n+1}+N(0,s)$, and the phases $\phi^{j}$ converge to a limit
$\phi=\phi_{i}$, where the particle index $i$ has been restored. The time
interval over which the solution is updated in each step is short, and we do
not expect any problem with convergence, either here or in the next section,
and indeed there is none; in all cases the iteration converges in a small
number of steps. Note that after the iteration the variables
$X^{n+1}_{i},Y^{n+1}_{i}$ are no longer independent- the observation creates a
relation between them.

Do this for all the particles. The particles are now samples of $P_{n+1}$, but
they have been obtained by sampling different densities (remember that the
parameters in the Gaussians in equation (8) vary). One can get rid of this
heterogeneity by viewing the factors $\exp(-\phi)$ as weights and resampling,
i.e., for each of $M$ random numbers $\theta_{k},k=1,\dots,M$ drawn from the
uniform distribution on $[0,1]$, choose a new ${\bf\hat{X}}^{n+1}_{k}={\bf
X}^{n+1}_{i}$ such that $Z^{-1}\sum_{j=1}^{i-1}\exp(-\phi_{j})<\theta_{k}\leq
Z^{-1}\sum_{j=1}^{i}\exp(-\phi_{j})$ (where
$Z=\sum_{j=1}^{M}\exp(-\phi_{j})$), and then suppress the hat. We have traded
the resampling of Bayesian filters for a resampling based on the normalizing
factors of the several Gaussian densities; this is a worthwhile trade because
in a Bayesian filter one gets a set of samples many of which may have low
probability with respect to $P_{n+1}$, and here we have a set of samples each
one of which has high probability with respect to a pdf close to $P_{n+1}$.

Note also that the resampling does not have to be done at every step- for
example, one can add up the phases for a given particle and resample only when
the ratio of the largest cumulative weight $\exp(-\sum\phi_{i})$ to the
smallest such weight exceeds some limit $L$ (the summation is over the weights
accrued to a particular particle $i$ since the last resampling). If one is
worried by too many particles being close to each other (”depletion” in the
Bayesian terminology), one can divide the set of particles into subsets of
small size and resample only inside those subsets, creating a greater
diversity. As will be seen in section 6, none of these strategies will be used
here and we will resample fully at every step.

##  5 Backward sampling.

The algorithm of the previous section is sufficient to create a filter, but
accuracy may require an additional refinement. Every observation provides
information not only about the future but also about the past- it may, for
example, tag as improbable earlier states that had seemed probable before the
observation was made; one may have to go back and correct the past after every
observation (this backward sampling is often misleadingly motivated solely by
the need to create greater diversity among the particles in a Bayesian
filter). As will be seen below, this backward sampling does not provide a
significant boost to accuracy in the present problem, but it is described here
for the sake a completeness.

Given a set of particles at time $(n+1)\delta$, after a forward step and maybe
a subsequent resampling, one can figure out where each particle $i$ was in the
previous two steps, and have a partial history for each particle $i$: ${\bf
X}_{i}^{n-1},{\bf X}_{i}^{n},{\bf X}_{i}^{n+1}$ (if resamples had occurred,
some parts of that history may be shared among several current particles).
Knowing the first and the last member of this sequence, one can interpolate
for the middle term as in section 2, thus projecting information backward.
This requires that one recompute $d{\bf X}^{n}$.

Let $d{\bf X}^{\text{tot}}=d{\bf X}^{n}+d{\bf X}^{n+1}$; in the present
section this quantity is assumed known and remains fixed. In the azimuth
problem discussed here, one has to deal with the slight complication due to
the fact that the mean of each increment is the value of the previous one, so
that two successive increments are related in a slightly more complicated way
than usual. The displacement $dX^{n}$ is a $N(dX^{n-1},\sigma)$ variable, and
$dX^{n+1}$ is a $N(dX^{n},\sigma)$ variable, so that one goes from $X^{n-1}$
to $X^{n+1}$ by sampling first a $(2dX^{n-1},4\sigma)$ variable that takes us
from ${X}^{n-1}$ to an intermediate point $P$, with a correction by the
observation half way up this first leg, and then one samples a
$N(dX^{\text{tot}},\sigma)$ variable to reach $X^{n+1}$, and similarly for
$Y$. Let the variable that connects ${\bf X}^{n-1}$ to $P$ be $d{\bf
X}^{\text{new}}$, so that what replaces $d{\bf X}^{n}$ is $d{\bf
X}^{\text{new}}/2$. Accordingly, we are looking for a new displacement $d{\bf
X}^{\text{new}}=(dX^{\text{new}},dY^{\text{new}})$, and for parameters
$a_{x}^{\text{new}},a_{y}^{\text{new}},v_{x}^{\text{new}},v_{y}^{\text{new}}$
such that

$\displaystyle\exp\left(-\frac{\xi_{x}^{2}+\xi_{y}^{2}}{2}\right)$
$\displaystyle=$
$\displaystyle\exp\left(-\frac{(dX^{\text{new}}-2dX^{n-1})^{2}}{8\sigma}-\frac{(dY^{\text{new}}-2dY^{n-1})^{2}}{8\sigma}\right)$
$\displaystyle\times\exp\left(-\frac{(f^{n}-b^{n})^{2}}{2s}\right)$
$\displaystyle\times\exp\left(-\frac{(dX^{\text{new}}-dX^{\text{tot}})^{2}}{2\sigma}-\frac{(dY^{\text{new}}-dX^{\text{tot}})^{2}}{2\sigma}\right)\exp(\phi)$
$\displaystyle=$
$\displaystyle\exp\left(-\frac{(dX^{\text{new}}-\bar{a}_{x})^{2}}{2v^{\text{new}}_{x}}-\frac{(dY^{\text{new}}-\bar{a}_{y})^{2}}{2v^{\text{new}}_{y}}\right),$

where $f^{n}=f(X^{n-1}+dX^{\text{new}}/2,Y^{n-1}+dY^{\text{new}}/2)$ and
$\xi_{x}$, $\xi_{y}$ are independent $N(0,1)$ Gaussian variables. As in
equation (8), the first equality embodies what we wish to accomplish- find
increments, functions of the reference variables, that sample the new pdf at
time $n\delta$ defined by the forward motion, the constraint imposed by the
observation, and by knowledge of the position at time $(n+1)\delta t$. The
second equality states that this is done by finding particle-dependent
parameters for a Gaussian density.

We again find these parameters as well as the increments by iteration. Much of
the work is separate for the $X$ and $Y$ components of the equations of
motion, so we write some of the equations for the $X$ component only. Again
set up an iteration for variables $dX^{\text{new},j}=dX^{j}$ which converge to
$dX^{\text{new}}$. Start with $dX^{0}=0$. To find $dX^{j+1}$ given $dX^{j}$,
approximate the observation equation (6), as before, by equation (9); define
again variables $\eta^{j+1},\eta^{j+1}_{+}$, one in the direction of the
approximate constraint and one orthogonal to it; in the direction of the
constraint multiply the pdfs as in the previous section; construct new means
$a^{1}_{x},a^{1}_{y}$ and new variances $v^{1}_{x},v^{1}_{y}$ for $dX,dY$ at
time $n$, taking into account the observation at time $n$, again as before.
This also produces a phase $\phi=\phi_{0}$.

Now take into account that the location of the boat at time $n+1$ is known;
this creates a new mean $\bar{a}_{x}$, a new variance $\bar{v}_{x}$, and a new
phase $\phi_{x}$, by $\bar{v}=\frac{v_{1}v_{2}}{v_{1}+v_{2}}$,
$\bar{a}_{x}=\frac{a_{1}v_{1}+a_{2}v_{2}}{v_{1}+v_{2}}$,
$\phi_{x}=\frac{(a_{1}-a_{2})^{2}}{v_{1}+v_{2}}$, where
$a_{1}=2a^{1},v_{1}=4v^{1}_{x},a_{2}=X^{\text{tot}},v_{2}=\sigma$. Finally,
find a new interpolated position
$dX^{j+1}=a_{x}^{\text{new}}/2+\sqrt{v_{x}^{\text{new}}}\xi_{x}$ (the
calculation for $dY^{j+1}$ is similar, with a phase $\phi_{y}$), and we are
done. The total phase for in this iteration is
$\phi=\phi_{0}+\phi_{x}+\phi_{y}$. As the iterates $dX^{j}$ converge to
$dX^{\text{new}}$, the phases converge to a limit $\phi=\phi_{i}$. The
probability of a particle arriving at the given position at time $(n+1)\delta
t$ having been determined in the forward step, there is no need to resample
before comparing samples. Once one has the values of $\bf{X}^{\text{new}}$, a
forward step gives corrected values of ${\bf X}^{n+1}$; one can use this
interpolation process to correct estimates of ${\bf X}^{k}$ by subsequent
observations for $k=n-1,k=n-2,\dots$, as many as are useful.

##  6 Numerical results.

Before presenting examples of numerical results for the azimuth problem, we
discuss the accuracy one can expect. A single set of observations for our
problem relies on 160 samples of a $N(0,\sigma)$ variable. The maximum
likelihood estimate of $\sigma$ given these samples is a random variable with
mean $\sigma$ and standard deviation $.11\sigma$. We estimate the uncertainty
in the position of the boat by picking a set of observations, then making
multiple runs of the boat where the random components of the motion in the
direction of the constraint are frozen while the ones orthogonal to it are
sampled over and over from the suitable Gaussian density, then computing the
distances to the fixed observations, estimating the standard deviation of
these differences, and accepting the trajectory if the estimated standard
deviation is within one standard deviation of the nominal value of $s$. This
process generates a family of boat trajectories compatible with the given
observations. In Table I we display the standard deviations of the differences
between the resulting paths and the original path that produced the
observations after the number of steps indicated there (the means of these
differences are statistically indistinguishable from zero). This Table
provides an estimate of the accuracy we can expect. It is fair to assume that
these standard deviations are underestimates of the uncertainty- a variation
of a single standard deviation in $s$ is a strict constraint, and we allowed
no variation in $\sigma$.

Table I

Intrinsic uncertainty in the azimuth problem

step |  $x$ component |  $y$ component
---|---|---
40 | .0005 | .21
80 | .004 | .58
120 | .010 | .88
160 | .017 | .95

If one wants reliable information about the performance of the filter, it is
not sufficient to run the boat once, record observations, and then use the
filter to reconstruct the boat’s path, because the difference between the true
path and the reconstruction is a random variable which may be accidentally
atypically small or atypically large. We have therefore run a large number of
such reconstructions and computed the means and standard deviations of the
discrepancies between path and reconstruction as a function of the number of
steps and of other parameters. In Tables II and III we display the means and
standard deviations of these discrepancies (not of their mean!) in the the x
and y components of the paths with 2000 runs, at the steps and numbers of
particles indicated, with no backward sampling. (Ref. [6] used 100 particles).
On the average the error is zero, and the error that can be expected in any
one run is of the order of magnitude of the unavoidable error. The standard
deviation of the discrepancy is not significantly smaller with 2 particles
that with 100- the main source of the discrepancy is the uncertainty in the
data. Most of time one single particle (no resampling) is enough; however, a
single particle may temporarily stray into low-probability areas and creates
large arguments and numerical difficulties in the various functions used in
the program. Two particles with resampling keep each other within bounds,
because if one of them strays it gets replaced by a replica of the other. The
various more sophisticated resampling strategies at the end of section 4 make
no discernible difference here, and backward sampling does not help much
either, because they too are unable to remedy the limitation of the data set.

Table IIa

Mean and standard variation of the discrepancy between synthetic data and
their reconstruction, 2000 runs, no back step, 100 particles

n. of steps | x component | y component
---|---|---
| mean | s.d. | mean | s.d.
40 | .0004 | .04 | .0001 | .17
80 | -.001 | .04 | -.01 | .54
120 | -.0008 | .07 | -.03 | 1.02
160 | -.002 | .18 | -.05 | 1.56

Table IIb

Mean and standard variation of the discrepancy between synthetic data and
their reconstruction, 2000 runs, no back step, 2 particles

n. of steps | x component | y component
---|---|---
| mean | s.d. | mean | s.d.
40 | .002 | .17 | -.0004 | .20
80 | .01 | .43 | -.0006 | .58
120 | .01 | .57 | .009 | 1.08
160 | .006 | .54 | .01 | 1.67

In Figure 1 we plot a sample boat path, its reconstruction, and the
reconstructions obtained (i) when the initial data for the reconstruction are
strongly perturbed (here, the initial data for $x,y$ were perturbed initially
by, respectively, $.1$ and $.4$), and (ii) when the value of $\sigma$ assumed
in the reconstruction is random: $\sigma=N(\sigma_{0},\epsilon\sigma_{0})$,
where $\sigma_{0}$ is the constant value used until now and $\epsilon=0.4$ but
the calculation is otherwise identical. This produces variations in $\sigma$
of the order of $40\%$; any larger variance in the perturbations produced
negative value of $\sigma$. The differences between the reconstructions and
the true path remain within the acceptable range of errors. These graphs show
that the filter has little sensitivity to perturbations (we did not calculate
statistics here because the insensitivity holds for each individual run).

Figure 1: Some boat trajectories (explained in the text)

We now estimate the parameter $\sigma$ from data. The filter needs an estimate
of $\sigma$ to function, call this estimate $\sigma_{\text{assumed}}$. If
$\sigma_{\text{assumed}}\neq\sigma$, the other assumptions used to produce the
data set (e.g. independence of the displacements and of the observations) are
also false, and all one has to do is detect the fallacy. We do it by picking a
trajectory of a particle and computing the quantity

$D=\frac{(\sum_{2}^{J}(dX^{j+1}-dX^{j}))^{2}+(\sum_{2}^{J}(dY^{j+1}-dY^{j}))^{2}}{\sum_{2}^{J}(dX^{j+1}-dX^{j})^{2}+\sum_{2}^{J}(dY^{j+1}-dY^{j})^{2}}.$

If the increments are independent then on the average $D=1$; we will try to
find the real $\sigma$ by finding a value of $\sigma_{\text{assumed}}$ for
which this happens. We chose $J=40$ (the early part of a trajectory is less
noisy than the later parts).

As we already know, a single run cannot provide an accurate estimate of
$\sigma$, and accuracy in the reconstruction depends on how many runs are
used. In Table III we display some values of $D$ averaged over 200 and over
5000 runs as a function of the ratio of $\sigma_{\text{assumed}}$ to the value
of $\sigma$ used to generate the data. From the longer computation one can
find the correct value of $\sigma$ with an error of about $3\%$, while with
200 runs the uncertainty is about $10\%$.

Table III

The mean of the discriminant D as a function of
$\sigma_{\text{assumed}}/\sigma$, 30 particles

$\sigma_{\text{assumed}}/\sigma$ | 5000 runs | 200 runs
---|---|---
.5 | 1.14 $\pm$ .01 | 1.21 $\pm$ .08
.6 | 1.08 $\pm$ .01 | 1.14 $\pm$ .07
.7 | 1.05 $\pm$ .01 | 1.10 $\pm$ .07
.8 | 1.04 $\pm$ .01 | 1.14 $\pm$ .07
.9 | 1.00 $\pm$ .01 | 1.01 $\pm$ .07
1.0 | 1.00 $\pm$ .01 | .96 $\pm$ .07
1.1 | .97 $\pm$ .01 | 1.01 $\pm$ .07
1.2 | .94 $\pm$ .01 | .99 $\pm$ .07
1.3 | .93 $\pm$ .01 | 1.02 $\pm$ .07
1.4 | .90 $\pm$ .01 | .85 $\pm$ .06
1.5 | .89 $\pm$ .01 | .93 $\pm$ .07
2.0 | .86 $\pm$ .01 | .78 $\pm$ .05

##  7 Conclusions.

We have exhibited a non-Bayesian filtering method, related to recent work on
chainless sampling, designed to focus particle paths more sharply and thus
require fewer of them, at the cost of an added complexity in the evaluation of
each path. The main features of the algorithm are a representation of a new
pdf by means of a set of functions of Gaussian variables and a resampling
based on normalization factors. The construction was demonstrated on a
standard ill-conditioned test problem. Further applications will be published
elsewhere.

##  8 Acknowledgments.

We would like to thank Prof. R. Kupferman, Prof. R. Miller, and Dr. J. Weare
for asking searching questions and providing good advice. This work was
supported in part by the Director, Office of Science, Computational and
Technology Research, U.S. Department of Energy under Contract No. DE-
AC02-05CH11231, and by the National Science Foundation under grant
DMS-0705910.

## References

  * [1] Doucet, A., de Freitas, N., and Gordon, N. (eds) (2001), Sequential Monte Carlo Methods in Practice, Springer, New York.
  * [2] Maceachern, S., Clyde, M., and Liu, J. (1999), Can. J. Stat. 27:251–267.
  * [3] Liu, J. and Sabati, C. (2000), Biometrika 87:353–369.
  * [4] Doucet, A., Godsill, S., and Andrieu, C. (2000), Stat. Comp. 10:197–208.
  * [5] Arulampalam, M., Maskell, S., Gordon, N., and Clapp, T. (2002), IEEE Trans. Sig. Proc. 50:174–188.
  * [6] Gilks, W. and Berzuini, C. (2001), J. Roy. Statist. Soc. B 63:127–146.
  * [7] Chorin, A.J. and Krause, P. (2004), Proc. Nat. Acad. Sci. USA 101:15013–15017.
  * [8] Dowd, M. (2006), Environmetrics 17:435–455.
  * [9] Doucet, A. and Johansen, A., Particle Filtering and Smoothing: Fifteen years Later, Handbook of Nonlinear Filtering (eds. D. Crisan et B. Rozovsky), Oxford University Press, to appear.
  * [10] Snyder, C., Bengtsson, T., Bickel, P, and Anderson, J. (2008), Mon. Wea. Rev. 136:4629–4640.
  * [11] Chorin, A.J. (2008), Comm. Appl. Math. Comp. Sc. 3:77–93.
  * [12] Weare, J. (2007), Proc. Nat. Acad. Sc. USA 104:12657–12662.
  * [13] Weare, J., Particle filtering with path sampling and an application to a bimodal ocean current model, in press, J. Comput. Phys.
  * [14] Gordon, N., Salmond, D., and Smith A. (1993), IEEE Proceedings-F 140: 107–113.
  * [15] Carpenter, J., Clifford, P., and Fearnhead, P. (1999) IEEE Proceedings-Radar Sonar and Navigation 146:2–7.
  * [16] Milstein, G., Platen, E., and Schurz H. (1998), SIAM J. Num. Anal. 35:1010–1019.
  * [17] Kloeden, P. and Platen, E. (1992), Numerical Solution of Stochastic Differential Equations, Springer, Berlin.
  * [18] Stuart, A., Voss, J., and Wilberg, P. (2004), Comm. Math. Sc.  4:685–697.