Coastal observation design with fixed scenario banks¶
This example compares 27 hypothetical coastal observation designs. It shows how to keep null and event simulations separate, share random-number banks across designs, summarize replicated detections with Wilson intervals, apply a false-alarm rule, and extract a three-objective Pareto front.
Teaching inputs only
Every signal, measurement-error scale, threshold, cost, and decision rule in this example is invented. The values are not field measurements, instrument specifications, facility quotes, partner commitments, protocol choices, or recommendations for a real monitoring program.
The complete script is
examples/coastal_observation_design.py.
The fast check runs the complete workflow with 60 Monte Carlo trials per design
and scenario on one CPU process:
Omit --check to use the example's prespecified 400 trials per design and
scenario. The --n-jobs option changes CPU parallelism for that full run;
values other than one require the parallel extra, while --check always uses
one process.
Decision question and quantities¶
The teaching question is whether additional sampling events, technical replicates, or a different measurement mode improve event detection and intensity estimation enough to justify their hypothetical cost. The latent event intensity and measurement errors use arbitrary intensity units. Cost is reported in invented cost units, without a currency conversion.
Each design has three factors:
| Factor | Levels | Meaning |
|---|---|---|
| Measurement mode | rapid sensor, reference assay, paired cross-platform | Three invented error and cost settings |
| Sampling events | 4, 7, 13 | Fixed subsets of 13 times in one normalized monitoring window |
| Technical replicates | 1, 2, 3 | Nested measurements averaged at each sampling time |
The full factorial grid therefore contains \(3 \times 3 \times 3 = 27\) designs. The synthetic scorer reports a binary detection indicator, squared intensity error, and deterministic cost. Event intensity root mean squared error (RMSE) is the square root of mean squared error across event-bank trials. For each trial, estimated intensity is the nonnegative maximum of the sampled values, and detection is one when that estimate is at least 0.80. The mean detection indicator estimates detection probability in the event bank and false-alarm probability in the null bank. Event-bank RMSE compares estimated intensity with the latent event amplitude.
def coastal_design_grid() -> list[dict[str, Any]]:
"""Return the 27-design full factorial teaching grid."""
factors = [
Factor(
"measurement_mode",
FactorType.CATEGORICAL,
levels=[
"rapid_sensor",
"reference_assay",
"paired_cross_platform",
],
),
Factor("sampling_events", FactorType.DISCRETE, levels=[4, 7, 13]),
Factor("technical_replicates", FactorType.DISCRETE, levels=[1, 2, 3]),
]
return build_grid(factors, method="full")
def coastal_observables() -> list[Observable]:
"""Return the raw outcomes and their optimization directions."""
return [
Observable("detected", Direction.MAXIMIZE),
Observable("squared_intensity_error", Direction.MINIMIZE),
Observable("cost_units", Direction.MINIMIZE),
]
Prespecify the null and event banks¶
The null bank fixes the latent event amplitude at zero. The event bank draws a
positive Gaussian-shaped event with a random amplitude, center, and width. A
bank contains n_reps independent Monte Carlo trials indexed from zero through
n_reps - 1; the bank names and root seed are fixed before evaluating designs.
@dataclass(frozen=True)
class ScenarioBank:
"""Declare one prespecified synthetic scenario bank."""
name: str
event_present: bool
NULL_BANK = ScenarioBank("null", event_present=False)
EVENT_BANK = ScenarioBank("event", event_present=True)
SCENARIO_BANKS = (NULL_BANK, EVENT_BANK)
Within a bank, the same trial index selects the same latent event, calibration draw, and standardized measurement-error array for every design. This is a common-random-number construction: designs are compared on matched synthetic conditions, while their measurement mode, sampling schedule, and technical replication transform or subset the shared draws. Common random numbers can reduce Monte Carlo noise in design contrasts, but they do not add independent trials or validate the assumed data-generating process.
class SyntheticCoastalWorld:
"""Generate shared latent events and nested hypothetical measurement errors.
Within a scenario bank, a replicate identifier selects the same latent event,
calibration draw, and standardized error array for every design. Designs
transform or subset that common random-number bank according to their
declared settings.
"""
def __init__(self, scenario: ScenarioBank) -> None:
"""Store the fixed scenario declaration."""
self.scenario = scenario
def _truth(self, rep: int) -> SyntheticTruth:
"""Return the shared latent truth for one bank replicate."""
truth_rng = np.random.default_rng(
stable_seed(self.scenario.name, rep, "latent-event")
)
if self.scenario.event_present:
amplitude = max(0.20, float(truth_rng.normal(1.20, 0.20)))
center = float(truth_rng.uniform(0.35, 0.65))
width = float(truth_rng.uniform(0.07, 0.13))
else:
amplitude, center, width = 0.0, 0.50, 0.10
return SyntheticTruth(amplitude, center, width)
def _measurement_bank(self, rep: int) -> tuple[float, NDArray[np.float64]]:
"""Return the shared calibration draw and standardized error bank."""
common_rng = np.random.default_rng(
stable_seed(self.scenario.name, rep, "measurement-bank")
)
return float(common_rng.normal()), common_rng.normal(
size=(3, MASTER_TIMES.size)
)
def generate(
self, config: dict[str, Any], *, rep: int = 0
) -> tuple[SyntheticTruth, dict[str, NDArray[np.float64]]]:
"""Return latent truth and observations for one synthetic replicate.
Returns:
Shared latent truth and the design-specific observed values.
Raises:
ValueError: If the replicate identifier is negative.
"""
if rep < 0:
message = "rep must be non-negative"
raise ValueError(message)
truth = self._truth(rep)
latent = truth.amplitude * np.exp(
-0.5 * ((MASTER_TIMES - truth.center) / truth.width) ** 2
)
calibration_draw, standardized_errors = self._measurement_bank(rep)
mode = str(config["measurement_mode"])
parameters = MODE_PARAMETERS[mode]
technical_replicates = int(config["technical_replicates"])
sampling_events = int(config["sampling_events"])
indices = SAMPLING_INDICES[sampling_events]
measurement_error = standardized_errors[:technical_replicates].mean(axis=0)
observed = (
latent
+ parameters["bias_scale"] * calibration_draw
+ parameters["noise_scale"] * measurement_error
)
return truth, {"values": observed[indices]}
class SyntheticCoastalScorer:
"""Convert hypothetical observations into detection, error, and cost."""
def score(
self,
truth: SyntheticTruth,
observations: dict[str, NDArray[np.float64]],
config: dict[str, Any],
) -> dict[str, float]:
"""Score one replicate in declared intensity and invented cost units.
Returns:
Detection indicator, squared intensity error, and invented cost.
"""
estimated_intensity = max(0.0, float(np.max(observations["values"])))
mode = str(config["measurement_mode"])
parameters = MODE_PARAMETERS[mode]
cost = parameters["fixed_cost"] + parameters["sample_cost"] * int(
config["sampling_events"]
) * int(config["technical_replicates"])
return {
"detected": float(estimated_intensity >= DETECTION_THRESHOLD),
"squared_intensity_error": (estimated_intensity - truth.amplitude) ** 2,
"cost_units": float(cost),
}
Aggregate detections and intensity errors¶
For each bank, run_grid(..., n_reps=n_reps) retains every design and trial
identity. ResultsTable.aggregate_replicates then computes the mean of each
outcome. Passing proportion_observables=["detected"] also records the number
of detections, number of trials, observed proportion, and a two-sided marginal
95% Wilson interval without continuity correction.
The Wilson intervals describe fixed-sample Monte Carlo uncertainty under the example's independent-trial construction. They are marginal intervals: the example does not adjust them for examining 27 designs, using the null intervals for selection, or reporting the selected Pareto front. See Fixed-sample proportions for the API contract and precision-planning boundaries.
Apply the false-alarm rule before Pareto extraction¶
The example calls a design feasible when the upper endpoint of its null-bank 95% Wilson interval is no greater than the hypothetical false-alarm limit of 0.10. Equality passes. This conservative teaching rule uses the entire fixed null bank and is not a regulatory, ecological, or operational standard.
Among feasible evaluated designs, extract_front uses three objectives:
| Objective | Direction | Scenario bank |
|---|---|---|
| Detection probability | Maximize | Event |
| Intensity RMSE | Minimize | Event |
| Invented cost units | Minimize | Deterministic, identical in both banks |
def run_coastal_study(
*, n_reps: int = DEFAULT_REPLICATES, n_jobs: int = 1
) -> list[DesignResult]:
"""Evaluate fixed null/event banks and extract a feasible Pareto front.
The false-alarm gate requires the upper bound of a marginal 95% Wilson
interval to be no greater than 0.10. The three Pareto objectives among
designs that pass that rule are event-detection probability (maximize),
event-intensity root mean squared error (minimize), and hypothetical cost
units (minimize).
Args:
n_reps: Fixed number of Monte Carlo trials per design and scenario.
n_jobs: Number of CPU workers passed to ``run_grid``.
Returns:
All evaluated designs, sorted with feasible Pareto members first.
Raises:
ValueError: If fewer than two trials are requested or ``n_jobs`` is zero.
"""
if n_reps < 2:
message = "n_reps must be at least 2"
raise ValueError(message)
if n_jobs == 0:
message = "n_jobs must not be zero"
raise ValueError(message)
bank_summaries = {
scenario: _summarize_bank(scenario, n_reps=n_reps, n_jobs=n_jobs)
for scenario in SCENARIO_BANKS
}
event = bank_summaries[EVENT_BANK]
null = bank_summaries[NULL_BANK]
results: list[DesignResult] = []
for config in coastal_design_grid():
key = _design_key(config)
event_summary, null_summary = event[key], null[key]
results.append(
DesignResult(
measurement_mode=key[0],
sampling_events=key[1],
technical_replicates=key[2],
detection_probability=event_summary.probability,
detection_low=event_summary.interval_low,
detection_high=event_summary.interval_high,
false_alarm_probability=null_summary.probability,
false_alarm_low=null_summary.interval_low,
false_alarm_high=null_summary.interval_high,
event_detections=event_summary.successes,
false_alarms=null_summary.successes,
event_trials=event_summary.trials,
null_trials=null_summary.trials,
intensity_rmse=event_summary.intensity_rmse,
cost_units=event_summary.cost_units,
interval_confidence=event_summary.interval_confidence,
interval_method=event_summary.interval_method,
interval_scope=event_summary.interval_scope,
root_seed=ROOT_SEED,
passes_false_alarm_rule=(
null_summary.interval_high <= FALSE_ALARM_LIMIT
),
)
)
feasible_positions = [
index for index, result in enumerate(results) if result.passes_false_alarm_rule
]
if feasible_positions:
objectives = np.array(
[
[
results[index].detection_probability,
results[index].intensity_rmse,
results[index].cost_units,
]
for index in feasible_positions
],
dtype=float,
)
local_front = extract_front(
objectives,
[Direction.MAXIMIZE, Direction.MINIMIZE, Direction.MINIMIZE],
)
for local_index in local_front:
position = feasible_positions[int(local_index)]
results[position] = replace(
results[position], pareto_among_evaluated_feasible=True
)
return sorted(
results,
key=lambda result: (
not result.pareto_among_evaluated_feasible,
result.cost_units,
-result.detection_probability,
result.measurement_mode,
result.sampling_events,
result.technical_replicates,
),
)
The returned flag pareto_among_evaluated_feasible has a deliberately narrow
meaning. A flagged design is nondominated only among the feasible members of
this 27-design grid under the finite synthetic banks. The front is not a global
optimum and does not identify a preferred real-world design.
Limitations¶
- The event shape, Gaussian errors, threshold, and mode parameters are hypothetical; real observations may have temporal, spatial, calibration, censoring, failure, and distributional structure absent here.
- The peak-based intensity estimator can be positively biased because it takes a maximum over noisy observations. More sampling events can change both this bias and the detection probability.
- The Wilson intervals are approximate, marginal, and fixed-sample. They do not supply simultaneous coverage, a post-selection guarantee, or permission to stop simulation when a desired result appears.
- The same finite banks support feasibility screening and Pareto extraction. An independent validation bank would be needed to assess a selected design's operating characteristics without this reuse.
- Common random numbers support matched comparisons and can improve their precision when the coupling induces favorable covariance. They do not increase the number of independent event or null trials.
- Invented cost units omit procurement, staffing, logistics, maintenance, failures, data management, permitting, and other program-specific costs.