---
title: "Simulating uncertainty"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Simulating uncertainty}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(comment = "#>", message = FALSE, warning = FALSE)
```

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)`](https://tidyfault.netlify.app/reference/quantify.html), 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](https://tidyfault.netlify.app/articles/quantify-binary).
- [`it_security_probs`](https://tidyfault.netlify.app/reference/it_security_probs.html) holds point estimates at $t = 10{,}000$ hours.

```{r}
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:

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

```{r}
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)
```

---

## Summarize the top-event distribution

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

---

## 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](https://tidyfault.netlify.app/articles/quantify-probabilities).