Simulating uncertainty

When point estimates of basic-event probabilities are uncertain, you can simulate rates or probabilities, evaluate the top event for each draw with quantify(..., prob = TRUE), and summarize the resulting distribution. This follows the IT security workflow in the tidyfault paper.

Model and bundled data

Use the IT security tree as in quantify()-ing binary outcomes.
it_security_probs holds point estimates at \(t = 10{,}000\) hours.

library(tidyfault)
library(dplyr)
library(tidyr)

data("it_security_nodes")
data("it_security_edges")
data("it_security_probs")

f <- formulate(equate(curate(it_security_nodes, it_security_edges)))

Simulate uncertainty in failure rates

Outline:

Convert each point probability at time \(t\) to an exponential rate \(\lambda = -\log(1-p)/t\).
Perturb \(\lambda\) with a normal draw using an illustrative coefficient of variation.
Map back to \(p(t)\) with pexp().
Call quantify() once on the full matrix of simulated probabilities (vectorized over rows).

set.seed(7)
t_horizon <- 10000
cv <- 0.2
n_draws <- 2500

batch_probs <- it_security_probs |>
  pivot_longer(everything(), names_to = "event", values_to = "prob") |>
  mutate(
    lambda = -log(1 - prob) / t_horizon,
    lambda_se = lambda * cv
  ) |>
  group_by(event) |>
  reframe(
    prob = pexp(
      q = t_horizon,
      rate = rnorm(n = n_draws, mean = lambda, sd = lambda_se)
    ),
    scenario = seq_len(n_draws)
  ) |>
  ungroup() |>
  pivot_wider(id_cols = scenario, names_from = event, values_from = prob)

event_cols <- setdiff(names(batch_probs), "scenario")
sim_probs <- batch_probs
sim_probs$outcome <- quantify(
  f,
  newdata = sim_probs[, event_cols, drop = FALSE],
  prob = TRUE
)

knitr::kable(head(sim_probs, 3), digits = 4)

scenario	DA	EP	IM	LR	MN	PC	PM	PO	VS	WB	outcome
1	0.3038	0.1884	0.0732	0.1237	0.1680	0.1923	0.1579	0.1612	0.0951	0.1665	0.0144
2	0.1722	0.1382	0.0733	0.1105	0.2018	0.2375	0.1557	0.1502	0.1042	0.1418	0.0140
3	0.1926	0.1363	0.0522	0.0849	0.2360	0.1464	0.1621	0.1850	0.1196	0.1662	0.0134

Summarize the top-event distribution

sim_probs |>
  summarize(
    estimate = mean(outcome),
    lower = quantile(outcome, 0.025),
    upper = quantile(outcome, 0.975)
  ) |>
  knitr::kable(digits = 4)

estimate	lower	upper
0.0157	0.0101	0.0223

Interpretation

The interval reflects propagation of rate uncertainty through the Boolean structure, not sampling error from field data.
Tighten or widen the simulation by changing cv, n_draws, or the distributional assumption on \(\lambda\).

For deterministic probability scenarios, see quantify()-ing probabilities.