The Monte-Carlo engine

library(ambre)
set.seed(2024)

What separates ambre from a spreadsheet is that almost nothing in it is a single number. Concentrations, exposure volumes and log-reductions are all probability distributions, and the package draws from them thousands of times to report a range of risk rather than one fragile point estimate. This vignette opens up that Monte-Carlo core: the two loops it runs, the distributions it understands, and how to reproduce a result. It is the machinery behind the boxplots in vignette("a-get-started", package = "ambre") and the spread you interpret in vignette("d-interpreting-risk", package = "ambre").

Two dimensions: repeats and events

Every simulated quantity is drawn on a grid of two indices, both controlled by config_ambre$exposure:

config_ambre$exposure[, c("name", "type", "value", "min", "max")]
#> # A tibble: 3 × 5
#>   name                 type     value   min   max
#>   <chr>                <chr>    <dbl> <dbl> <dbl>
#> 1 number_of_repeatings value     1000  NA      NA
#> 2 number_of_exposures  value      365  NA      NA
#> 3 volume_perEvent      triangle    NA   0.5     3

So a single draw of, say, exposure volume is really a repeatID x eventID table – 1000 x 365 values – and the whole pipeline carries that shape through to the final DALYs.

The distribution catalog

The third row above, volume_perEvent, is a triangle with min 0.5 and max 3. That type / value / min / max / mode / mean / sd / meanlog / sdlog block is the shared vocabulary of the whole database: the same nine columns describe an inflow concentration, a treatment log-reduction and a cost. create_random_distribution() reads that vocabulary. The supported type values are:

type drawn with parameters used
value (constant) value
uniform runif min, max
log10_uniform 10^runif min, max (as log10)
norm rnorm mean, sd
log10_norm 10^rnorm mean, sd (as log10)
lognorm rlnorm meanlog, sdlog
triangle EnvStats::rtri min, max, mode

Call it directly to see what a draw looks like. It returns an events table (one row per repeatID x eventID) and the paras actually used:

draw <- create_random_distribution(
  type = "uniform",
  number_of_repeatings = 3,
  number_of_events = 5,
  min = 1, max = 10,
  debug = FALSE
)
head(draw$events)
#>   repeatID eventID   values
#> 1        1       1 8.532483
#> 2        1       2 3.887808
#> 3        1       3 7.123270
#> 4        1       4 7.283558
#> 5        1       5 5.113083
#> 6        2       1 7.312783
draw$paras
#>      type repeatings events min max
#> 1 uniform          3      5   1  10

In practice you rarely call it by hand. The pipeline calls generate_random_values(), a thin wrapper that reads one row of a config table and fills in sensible defaults before delegating. Feed it the real volume_perEvent row and you get the exact draw the engine would use for exposure volume:

vol_row <- config_ambre$exposure[config_ambre$exposure$name == "volume_perEvent", ]
gv <- generate_random_values(
  vol_row,
  number_of_repeatings = 2,
  number_of_events = 4,
  debug = FALSE
)
summary(gv$events$volume_perEvent)
#> Length  Class   Mode 
#>      0   NULL   NULL

The distribution schema is described further in vignette("h-config-ambre", package = "ambre").

From min/max to a spread

For the normal and log-normal families you often supply only a min and a max, not a mean and standard deviation. ambre derives the missing parameters so that a given fraction of the mass falls inside [min, max], using get_percentile() (a wrapper around the normal quantile) together with default_min() / default_max(). For example, the z-score that places 90% of a normal distribution within the min-max window is:

ambre:::get_percentile(0.9)
#> [1] 1.644854

You do not normally call these yourself, but knowing they exist explains why a norm row with only min and max still produces a sensible bell curve.

A special case worth knowing

generate_random_values() has one silent behaviour. A triangle whose min equals its max cannot be a triangle, so it is quietly converted to a (degenerate) uniform – effectively a constant:

degenerate <- data.frame(
  type = "triangle", value = NA,
  min = 2, max = 2, mode = NA,
  mean = NA, sd = NA, meanlog = NA, sdlog = NA
)
gv2 <- generate_random_values(degenerate, number_of_repeatings = 1,
                              number_of_events = 3, debug = FALSE)
unique(gv2$events$values)
#> [1] 2

You will see the message “Distribution set from ‘triangle’ to ‘uniform’ because ‘min’ equals ‘max’” in a real run when this happens – it is expected, not an error.

Why barrier credits add up

Log-reductions combine by addition in log space. A train that removes 2 log at the plant and 3 log in the field removes 5 log overall, and 5 log means the surviving concentration is 10^-5 of the inflow. Internally ambre builds a wide repeatID x eventID table of each barrier’s draw and sums the columns, so the additivity holds run by run, preserving the correlation structure across the Monte-Carlo sample. This is why you can compare a treatment strategy and a barrier strategy on the same footing (see vignette("b-initial-vs-new-scenario", package = "ambre")).

Reproducibility

Because every draw is random, set a seed before a run you want to reproduce. Same seed, same numbers:

set.seed(1)
a <- create_random_distribution(type = "uniform", number_of_repeatings = 1,
                                number_of_events = 4, min = 0, max = 1, debug = FALSE)
set.seed(1)
b <- create_random_distribution(type = "uniform", number_of_repeatings = 1,
                                number_of_events = 4, min = 0, max = 1, debug = FALSE)
identical(a$events$values, b$events$values)
#> [1] TRUE

To read the uncertainty a run expresses, summarise a quantity across repeatID with quantiles – exactly what the final DALY boxplots do, and what vignette("d-interpreting-risk", package = "ambre") walks through.