Impact probability computation of Near-Earth Objects using Monte Carlo Line Sampling and Subset Simulation

Camilla Colombo; Matteo Losacco; Matteo Romano; Pierluigi Di Lizia

REVIEW 1 major objections 5 minor 2 references

Impact probability computation of Near-Earth Objects using Monte Carlo Line Sampling and Subset Simulation

T0 review · 1 major / 5 minor · reviewed 2026-08-14 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗ Review JSON

Pith's one-line read Line sampling and subset simulation compute rare near-Earth-object impact probabilities at the same accuracy as standard Monte Carlo while using one to three orders of magnitude fewer orbital propagations.

desk verdict Line sampling is a genuine new application for NEA impact probabilities, but the reported LS error bars leave out the dominant uncertainty (the reference direction), so the 'same accuracy' claim is not yet supported as written. read the letter →

arxiv 1908.03063 v2 pith:FE2W3I34 submitted 2019-08-08 astro-ph.EP

Matteo Romano , Matteo Losacco , Camilla Colombo , Pierluigi Di Lizia This is my paper · ORCID

classification astro-ph.EP

keywords near-EarthasteroidsimpactprobabilitylinesamplingsubsetsimulationMonteCarlomethodsrareeventestimationuncertaintypropagation

verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

Estimating the probability that an asteroid will hit Earth normally means propagating tens of thousands to millions of sampled orbits through a full Solar System model. This paper adapts two rare-event Monte Carlo schemes, line sampling and subset simulation, to the impact-monitoring problem. Line sampling probes the uncertainty region with parallel lines and integrates the probability along each line, while subset simulation writes the small impact probability as a product of larger conditional probabilities. On three real near-Earth objects the two methods reproduce the standard Monte Carlo estimates while using one to three orders of magnitude fewer propagations when the impact probability is rare, although they do not beat standard Monte Carlo when the probability is high. The paper argues that these methods are a practical, cheaper alternative whenever the usual line-of-variations approach is unreliable or too complex.

What carries the argument

The load-bearing identities are the line-sampling integral of Eq. (14), where each sampled line gives a conditional impact probability $\Phi(c_k^2)-\Phi(c_k^1)$ and the total is their average, and the subset-simulation product rule of Eq. (17), where a small probability is obtained as $P(I_1)\prod_i P(I_{i+1}\mid I_i)$. These are carried by a probability-preserving mapping into independent standard normal coordinates, a Markov-chain Monte Carlo step that finds the impact region and the reference direction, and the fixed conditional probability $p_0$ with its nested sample levels. The machinery converts a rare-event counting problem into a set of cheap one-dimensional or conditional estimates.

What would settle it

Run line sampling with a single reference direction on a synthetic uncertainty set whose impact region is a thin arc or two disconnected lobes, and compare its probability estimate to a brute-force Monte Carlo run with ten million samples; a deviation larger than the Monte Carlo standard deviation would show that the at-most-two-intersections-per-line assumption has been violated.

Watch

Extended reading notes

Core claim

The paper establishes that impact probability can be computed by concentrating the sampling effort where it matters instead of spreading samples uniformly. For line sampling, the orbital uncertainty is mapped to a standard normal space, a reference direction is chosen from a Markov-chain sample lying inside the impact region, and each line parallel to that direction contributes a conditional probability $\Phi(c_k^2) - \Phi(c_k^1)$, so the total estimate is the average of these one-dimensional integrals. For subset simulation, the impact event $I$ is reached through nested intermediate events, and $P(I_n) = P(I_1)\prod_i P(I_{i+1}\mid I_i)$ is approximated as $p_0^{n-1} N_I / N$ with a fixed conditional probability $p_0 = 0.2$. On asteroids 2010 RF12, 2017 RH16, and 99942 Apophis the two methods match the standard Monte Carlo probability estimate while reducing the required number of propagations by one to three orders of magnitude for the rarer events, with line sampling giving particularly low variance at very low probabilities.

Load-bearing premise

The line-sampling estimate assumes each sampling line crosses the impact region at most twice, so the impact region must behave like a flat or slightly curved surface inside the selected time window; fragmented or strongly curved impact regions, or impact regions lying partly outside the sampled uncertainty ellipsoid, would bias the result.

Editorial extensions

If this is right

For impact probabilities around $10^{-3}$ or lower, impact monitoring can use one to three orders of magnitude fewer orbital propagations than standard Monte Carlo while keeping the same accuracy, as shown for 2017 RH16 and Apophis.
Line sampling becomes increasingly attractive as the impact probability decreases, offering much smaller estimate variance and higher figures of merit than both standard Monte Carlo and subset simulation in the very-rare-event regime.
For high impact probabilities, such as the $6.5\times10^{-2}$ case of 2010 RF12, neither advanced method beats standard Monte Carlo, so standard Monte Carlo remains the appropriate tool in that regime.
The reliability of both methods depends on parameter choices: the line-sampling reference direction and the accuracy of the impact-boundary root finding, and the subset-simulation conditional probability, with $p_0$ between about 0.1 and 0.4 recommended.
The methods are positioned as complements to the standard line-of-variations monitoring, intended for short-arc or strongly nonlinear cases where the line of variations is unreliable.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

Because the efficiency gain grows as the target probability shrinks, the same machinery should transfer to planetary-protection and space-debris re-entry problems, where the contamination or impact probabilities of interest are typically far below $10^{-4}$.
The line-sampling assumption of at most two crossings per line means strongly curved or fragmented impact regions will bias the estimate; an adaptive extension that re-orients the reference direction for each connected impact region, which the paper lists as future work, would remove the main limitation.
The Rosenblatt-style mapping used in line sampling is not tied to Gaussian uncertainty: any continuous input distribution can be transformed to standard normal coordinates, so the method could be tested on non-Gaussian orbit-determination posteriors.
A natural testable extension is to run both methods on synthetic impact regions with known probabilities, using the deviation from the known value to calibrate how many Newton iterations and how many lines are needed before the bias becomes negligible.

Share X Bluesky LinkedIn Reddit HN

Signed reviews

No signed human review yet.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Desk Editor's Note

read the letter

Matt,

This paper transfers two reliability methods, line sampling and subset simulation, to NEA impact probability estimation, with three test cases spanning a wide probability range. The LS application is new, and the comparison against SS is useful. They are honest about the high-probability case (2010 RF12) where neither method beats MC, and the sensitivity analysis is a real plus.

The numerical work is solid in its core: the propagator setup is standard, the test cases are real asteroids with realistic covariances, and the trend you'd expect—bigger gains for rarer events—is supported by the three cases. The SS implementation follows Au and Beck faithfully, and the choice of p0=0.2 and N=1000 is reasonable and justified.

The main problem is the LS variance estimate. Equation (16) only captures the spread across sampling lines for a fixed reference direction α, but α itself is estimated from a short MCMC chain. The paper's own sensitivity analysis shows for Apophis that a modest perturbation of α changes the probability estimate from 3.12e-5 to 2.73e-5, a shift of 3.9e-6, while the reported σ is 2.64e-7—about fifteen times smaller. So the error bars in Tables 5 and 6 for LS are not a fair basis for the 'same accuracy' claim, and the FoM values that depend on those error bars are correspondingly optimistic. This is not a formula nitpick; it directly affects the central conclusion. The authors acknowledge α is critical but never propagate its uncertainty into the reported variance.

A second, smaller issue: the configurations labeled σ_MC and N_MC_P are not generated by a reproducible procedure. For Apophis, they say 216 lines achieve MC-level accuracy, but they don't describe how that number was found, leaving the door open to optimistic selection. That is addressable with a simple adaptive rule or by reporting multiple seeds.

The two-intersection assumption is a genuine limitation, but the paper flags it and defers multi-region cases to future work. For the single-event scenarios presented, it is probably acceptable, though untested.

None of this makes the paper a waste. The LS transfer is a genuine contribution, and the sensitivity analysis is useful for practitioners. But the headline claim that LS grants 'the same accuracy' with orders-of-magnitude fewer propagations is not yet established until α uncertainty is folded into the error budget. SS is on firmer ground because its variance estimate comes from the Bayesian post-processor and the sensitivity results there are more stable.

I'd send this to peer review with a request for major revision focused on the LS variance treatment. The paper deserves a serious referee; the authors need to either account for α uncertainty in σ or soften the claim. If they can fix that, it becomes a solid reference for the community.

Best, [Your name]

Referee Report

1 major / 5 minor

Summary. This paper adapts two variance-reduction techniques from structural reliability, line sampling (LS) and subset simulation (SS), to the computation of near-Earth object impact probabilities. The LS method maps the initial-state uncertainty to a standard normal space, determines a reference direction from an MCMC sample of the impact region, integrates the Gaussian density along sampling lines between the two intersections of each line with the impact region, and averages the resulting conditional probabilities. The SS method estimates the rare impact probability as a product of larger conditional probabilities over nested intermediate event regions defined by decreasing thresholds on the minimum planetocentric distance. Both methods are tested on three NEAs with decreasing impact probability: 2010 RF12 (~6.5e-2), 2017 RH16 (~1.4e-3), and 99942 Apophis (~3e-5). The results are compared against standard Monte Carlo in terms of number of propagations, estimated probability, standard deviation, coefficient of variation, and figure of merit. A sensitivity analysis examines the effect of the LS reference direction, the LS root-finding accuracy, and the SS conditional probability. The authors conclude that for rare and very rare events both methods reduce the number of required propagations while granting the same accuracy as standard Monte Carlo.

Significance. If the quantitative claims are supported, this is a useful demonstration that rare-event simulation techniques, well established in reliability engineering, can be brought to bear on asteroid impact risk assessment and can reduce the computational cost of impact probability estimation by orders of magnitude for low-probability scenarios. The three test cases with decreasing probabilities support the expected trend of growing gains as the event becomes rarer, and the sensitivity analysis provides practical parameter guidance (e.g., p0 between 0.1 and 0.4 for SS, and a warning about the LS reference direction and Newton iterations). The paper is clearly written, the algorithms are described in enough detail to be reproduced, and the initial conditions and covariance matrices are tabulated. The qualitative conclusion—that LS and SS become increasingly attractive as the impact probability decreases—is credible and supported by the three demonstrations.

major comments (1)

[§2.4 and §5.1.1, Eqs. (15)–(16), Tables 6 and 8] The variance estimate in Eq. (16) treats the reference direction alpha as fixed, but alpha is itself an estimated quantity obtained from a finite MCMC chain (Eq. (7)). The sensitivity analysis in Tables 7 and 8 shows that perturbing alpha for the Apophis case changes the LS probability estimate from 3.12e-5 to 2.73e-5, a shift of about 3.9e-6, while Table 6 reports sigma = 2.64e-7 for the nominal configuration. The omitted alpha-uncertainty is therefore about an order of magnitude larger than the reported standard deviation. Because the FoM values and the 'same level of accuracy' comparison in Tables 5 and 6 rely directly on this sigma, the quantitative savings claim for LS is not yet supported. I recommend that the authors either (i) quantify and propagate the uncertainty in alpha (e.g., by repeating the MCMC estimation several times or by bootstrapping the MCMC output), (ii) report the sigma in Tables 5 and 6 explicitly as conditional on alpha and provide a separate total-variance estimate, or (iii) perform repeated independent LS runs that include the alpha-estimation step, so that the run-to-run variability is measured empirically. Without one of these, the reader cannot assess whether the LS accuracy is as high as claimed.

minor comments (5)

[Table 8] In Table 8, the LS_nom (sigma_MC) row reports sigma = 5.48e-7 and delta = 1.72e-2, whereas the same configuration in Table 6 gives sigma = 5.48e-6 and delta = 1.72e-1; the FoM values also differ (1.32e9 versus 1.33e7). Please correct this typo and ensure that the same configuration is reported consistently across tables.
[§5.2, discussion of Fig. 9] The sentence 'This trend confirms the results reported in Sect. 4 for the same test case, where MC outperformed SS with a value of conditional probability equal to 0.1' appears to refer to a p0 value of 0.2, which is the value used in Sect. 4; please clarify or correct the stated value.
[§6 (Conclusions)] There is a typographical error in the Conclusions: 'Y et, some improvements are still possible' should read 'Yet, some improvements are still possible.'
[§2.3] The two-intersection assumption is stated and later acknowledged as a limitation in the Conclusions and in Fig. 10. The paper would benefit from explicitly stating in Sect. 4 that all three test cases satisfy this assumption (a single connected impact region within the selected time window), so that the reader can judge the scope of the demonstrated efficiency.
[§2.3, Eq. (13)] The finite-difference step Delta_c used for the numerical derivative is introduced as 'arbitrarily small' but no value or selection criterion is given in the paper; please report the value used in the simulations or cite a reference for its choice.

Circularity Check

0 steps flagged · score 1.0 of 10

No circular derivation: LS and SS probabilities are independently sampled and benchmarked against standard MC; minor self-citations are not load-bearing.

full rationale

The derivation chain is self-contained. The line sampling estimate uses Eq. (15) to average Gaussian line integrals from Eq. (14), with the reference direction alpha from Eq. (7) being an input estimated by MCMC, not a parameter fitted to the Monte Carlo probability. The subset simulation estimate uses Eq. (20) as the product of conditional probabilities p0^(n-1) * N_I / N, with p0 = 0.2 and N = 1000 fixed a priori rather than tuned to reproduce the standard MC outputs. The comparisons in Tables 4-6 and the FoM values are computed from these independent estimators against standard MC; no constant is fitted to the target probabilities, and the sensitivity analysis in Sect. 5 explicitly quantifies the effects of perturbing alpha, reducing Newton iterations, and changing p0 instead of hiding them. The self-citations to Losacco et al. (2018), Morselli et al. (2015), and Romano et al. (2020) provide background, methodology provenance, or future-work context, and none is load-bearing for the probability estimates. The skeptic's concern that the variance in Eq. (16) conditions on a fixed alpha and omits reference-direction uncertainty is a legitimate accuracy and robustness issue, and it is visible in Tables 7-8, where a perturbed alpha for Apophis changes the estimate from 3.12e-5 to 2.73e-5; however, that is not circularity under the enumerated patterns, because the LS estimate is not defined in terms of the standard MC probability and no fitted parameter is renamed as a prediction. Similarly, the two-intersections assumption in Sect. 2.3 is a stated modeling limitation, acknowledged by the paper and deferred to future work for disconnected regions, not a circular reduction. No circular step can be quoted because the core probability estimates are independently benchmarked against standard Monte Carlo rather than derived from it.

Assumptions & free parameters 7 free parameters · 7 assumptions · 0 invented entities

The methods rest on standard probability transformations (Rosenblatt, Gaussian CDFs), a standard dynamical model, and a set of user-chosen hyperparameters (p0, N, Delta_t_e, proposal scales, iteration counts). No new physical entities, forces, or conserved quantities are introduced; the paper is purely computational.

free parameters (7)

SS conditional probability p0 = 0.2
Chosen a priori for all test cases (Sect. 3 and 5.2). The sensitivity analysis in Fig. 9 suggests 0.1 to 0.4 is reasonable; it is not fitted to the target impact probabilities, but it controls the hierarchy of conditional levels and the total sample count via Eqs. (19) and (20).
SS samples per conditional level N = 1000
Fixed to 1000 for all test cases (Sect. 3). It directly enters Eq. (19) for total sample count and governs the accuracy of each conditional level, as discussed in Sect. 3.
LS MCMC proposal distribution scale = 1/10 of the initial uncertainty distribution
User-defined proposal scaling in Sect. 2.2. It controls the Markov chain used to estimate the reference direction alpha, which Sect. 5.1.1 identifies as the main accuracy driver of the line sampling method.
LS event time window Delta_t_e = 200 days centered on the event epoch
Arbitrary window chosen around each close approach (Sect. 2). It defines the performance function Y(c) in Eq. (10) and therefore the shape of the impact region under analysis.
LS finite-difference step Delta_c = Arbitrarily small increment, value not reported
Used in Eq. (13) for numerical differentiation in the Newton iterations that locate the impact region boundaries along each sampling line.
LS Newton iteration count and root-finding tolerance = 4 to 5 iterations in nominal runs; tolerance unspecified
Tables 9 and 10 show that reducing the iteration count changes the probability estimate significantly (2 iterations for 2010 RF12 gives 1.35e-1 instead of 6.57e-2), so this is a load-bearing numerical control.
LS MCMC chain length NS = Not reported
The number of Markov chain samples used to compute the reference direction alpha in Eq. (7) is not specified in the numerical sections, making exact reproduction impossible.

assumptions (7)

domain assumption Initial state uncertainty is Gaussian in the chosen equinoctial parameters.
Stated in Sect. 2.1: 'the distribution of the initial conditions (position and velocity) is assumed to be Gaussian.' The Rosenblatt transformation and the standard-normal mapping in Eqs. (3)-(6) rely on this.
domain assumption Within the selected time interval Delta_t_e, the performance function Y(c) has at most two zeros and is smooth near the impact region.
Assumed in Sect. 2.3, justified by approximating the impact region as flat or slightly curved. Required for the line integral formula Eq. (14). The paper notes that larger intervals or disconnected impact regions break this assumption.
domain assumption The dynamical model (Sun, major planets, Moon, relativistic corrections) and JPL Horizons ephemerides are accurate enough that propagated close-approach distances are unbiased.
Sect. 4 states that all propagations use this model with an RK78 integrator and tolerances of 1e-12. The comparison to standard MC uses the same model, so model error cancels in the relative comparison, but it is a premise for the absolute probability values.
domain assumption Standard Monte Carlo estimates are converged enough to serve as reference probabilities.
LS and SS results are benchmarked against MC runs of 1e4, 5e4, and 1e6 samples. For Apophis at P = 3e-5 with 1e6 samples, the MC coefficient of variation is 0.18, so the benchmark itself has substantial uncertainty that is not propagated into the comparison.
standard math The Rosenblatt transformation and unit-Gaussian CDF computations are correct under the Gaussian assumption.
Standard probability results invoked without proof in Eqs. (3)-(6) and Eq. (14). Under the Gaussian assumption the transformation is linear, as cited from Zio and Pedroni (2009) and Rosenblatt (1952).
standard math The Bayesian SS+ moment formulas from Zuev et al. (2012) correctly describe the subset simulation estimator.
Eqs. (22)-(24) are borrowed from the cited literature and used to compute variances for the SS estimates in Sect. 4.
domain assumption The MCMC and fmincon procedure in LS finds a reference direction that points toward the true impact region.
The reference direction alpha in Eq. (7) is computed from a Markov chain seeded by an fmincon optimization. If the optimization misses the impact region, alpha is wrong and the LS estimate is biased, as the perturbed-direction results in Sect. 5.1.1 demonstrate.

how reviews work

0 comments

Cite this review

Pith. "Pith review of Impact probability computation of Near-Earth Objects using Monte Carlo Line Sampling and Subset Simulation." pith.science (2026). https://pith.science/paper/FE2W3I34

@misc{pith2026190803063,
  author       = {Pith},
  title        = {Pith review of: Impact probability computation of Near-Earth Objects using Monte Carlo Line Sampling and Subset Simulation},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/FE2W3I34}},
  note         = {Machine review of arXiv:1908.03063}
}

read the original abstract

This work introduces two Monte Carlo (MC)-based sampling methods, known as line sampling and subset simulation, to improve the performance of standard MC analyses in the context of asteroid impact risk assessment. Both techniques sample the initial uncertainty region in different ways, with the result of either providing a more accurate estimate of the impact probability or reducing the number of required samples during the simulation with respect to standard MC techniques. The two methods are first described and then applied to some test cases, providing evidence of the increased accuracy or the reduced computational burden with respect to a standard MC simulation. Finally, a sensitivity analysis is carried out to show how parameter setting affects the accuracy of the results and the numerical efficiency of the two methods.

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

2 extracted references · 2 canonical work pages

[1]

Armellin, R., Di Lizia, P ., Bernelli-Zazzera, F., Berz, M.: Asteroid close encounters characterization using differential algebra: the case of Apophis. Cel. Mech. Dyn. Astron. 107(4), 451–470 (2010) Au, S.-K., Beck, J.L.: Estimation of small failure probabilities in high dimensions by subset simulation. Probab. Eng. Mech. 16, 263–277 (2001) Au, S.-K., Wa...

work page 2010
[241]

In: Pham, H

University of Arizona Press, Tucson, Arizona, USA (1994) Zio, E.: System Reliability and Risk Analysis. In: Pham, H. (ed.) The Monte Carlo Simulation Method for System Reliability and Risk Analysis, pp. 7–17. Springer, London (2013) Zio, E., Pedroni, N.: Subset simulation and line sampling for advanced Monte Carlo reliability analysis. In: Proceedings of ...

work page 1994

Pith tools

Reviewed August 14, 2026 · model on record in the stance chip above.

[1] [1]

Armellin, R., Di Lizia, P ., Bernelli-Zazzera, F., Berz, M.: Asteroid close encounters characterization using differential algebra: the case of Apophis. Cel. Mech. Dyn. Astron. 107(4), 451–470 (2010) Au, S.-K., Beck, J.L.: Estimation of small failure probabilities in high dimensions by subset simulation. Probab. Eng. Mech. 16, 263–277 (2001) Au, S.-K., Wa...

work page 2010

[2] [241]

In: Pham, H

University of Arizona Press, Tucson, Arizona, USA (1994) Zio, E.: System Reliability and Risk Analysis. In: Pham, H. (ed.) The Monte Carlo Simulation Method for System Reliability and Risk Analysis, pp. 7–17. Springer, London (2013) Zio, E., Pedroni, N.: Subset simulation and line sampling for advanced Monte Carlo reliability analysis. In: Proceedings of ...

work page 1994