---
title: "Newmark ensemble formulation"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Newmark ensemble formulation}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", eval = FALSE)
```

## Overview

`newmark` answers one question: given a site-specific seismic hazard and a
slope with known dynamic properties, what horizontal seismic coefficient
*k*_max is required to limit permanent co-seismic displacement to a
specified target *d*\*?

The performance-based design problem is formulated as: for a given
target tolerable displacement *d*\* > 0, target return period *T*_R,
and target exceedance probability *p* ∈ (0, 1), choose the smallest
yield coefficient *k*_y such that the displacement induced at *T*_R
remains below *d*\* with probability at least 1 − *p*:

$$
k_\text{max}(p \mid d^*, T_R) \;=\; \inf\!\Bigl\{\,k_y \,:\, P\!\bigl[\,D_N(k_y; T_R) > d^*\,\bigr] \,\leq\, p\,\Bigr\}.
$$

The right-hand side has no closed form — neither the joint
distribution of intensity measures nor the convolution with the
empirical regression residual admits one. The framework evaluates
the inverse mapping numerically via Monte Carlo, propagating two
sources of uncertainty:

1. **Hazard uncertainty** — spectral accelerations are drawn from the
   quantile curves of the uniform-hazard spectrum via a Gaussian copula that
   preserves the **full** inter-period correlation matrix (Baker &
   Jayaram 2008); no star/hub simplification.
2. **Model uncertainty** — six empirical Newmark models are combined in a
   logic-tree ensemble with user-assigned epistemic weights *w*_i, ∑*w*_i = 1.

For the full mathematical derivation and calibration discussion, see:

> Verri Kozlowski, A. (2026). *Probabilistic estimation of Newmark
> displacements and seismic coefficients under hazard uncertainty.*
> Working paper.

---

## Three-stage structure

### Stage 1 — Site amplification (`fitSaF`, `fitSaFMixture`)

Rock-level spectral ordinates are amplified to the target site Vs30
using the NGA-East ergodic site-response models: `"ST20"` (default;
Stewart et al. 2020 linear model, Hashash et al. 2020 nonlinear model)
or `"ST17"` (the 2017 PEER-report generation), selected via the
`models` argument.

- **Mode A** (mean-only input): rock Sa is deterministic; dispersion
  comes from the site-model sigma alone.
- **Mode B** (quantile input): Sa is sampled per period from the
  empirical quantile function, with the PGA–Sa(Tn) dependence
  represented by a per-period two-variable Gaussian copula (Baker &
  Jayaram 2008 correlation); dispersion reflects combined hazard and
  site-factor uncertainty. The full inter-period correlation matrix is
  used downstream in `fitDnCurve()`, not in this per-period stage.
- **Mixture input** (`fitSaFMixture`): when the rock hazard is an
  analytic lognormal mixture (DSHA scenario engine), the composition
  with the site-factor lognormal is evaluated exactly by quadrature —
  no sampling, byte-stable output.

See `?fitSaF` and `?fitSaFMixture` for the function-level interface.

### Stage 2 — Displacement curve (`fitDnCurve`)

For each Monte Carlo realisation, one ground-motion scenario is drawn
(a consistent tuple of PGA, Sa(1.3Ts), Sa(1.5Ts)) and passed through all
active displacement models simultaneously using a shared aleatory residual.
The result is a family of coherent displacement curves *D*_N(*k*_y).

### Stage 3 — Coefficient inversion (`invertDnDraws`)

Each curve is projected monotone and inverted in log-log space to find
the yield acceleration *k*_max such that *D*_N(*k*_max) = *d*\*. Mean and
quantiles of *k*_max over all realisations constitute the output.

---

## Displacement models

| ID | Authors | Type | Spectral reference |
|---|---|---|---|
| AM88 | Ambraseys & Menu (1988) | rigid block | PGA |
| JB07 | Jibson (2007) | rigid block | PGA, Arias intensity |
| SR08 | Saygili & Rathje (2008) | rigid block | PGA, Arias intensity |
| BT07 | Bray & Travasarou (2007) | flexible block | Sa(1.5 Ts) |
| BM17 | Bray, Macedo & Travasarou (2018) | flexible block, subduction | Sa(1.5 Ts) |
| BM19 | Bray & Macedo (2019, corr. 2023) | flexible block, crustal | Sa(1.3 Ts) |

Activate/deactivate models via the `weights` argument (0 = inactive).

For BM19, near-fault pulse motions are flagged when PGV > 115 cm s⁻¹
and the equation switches to the Bray-Macedo (2023) D100 / D50 form
controlled by the `NFC` argument (`"D100"` is the maximum-component
default, conservative for slopes within ±45° of fault-normal;
`"D50"` is the median-component case for other orientations). A
sub-regime split at PGV = 150 cm s⁻¹ inside the pulse equation
captures the saturation of seismic displacement at very high PGV.

---

## Key assumptions

- Shared aleatory residual across models per realisation: a single
  *z*^(*n*) ∼ 𝒩(0,1) is drawn and applied to all six models scaled
  by each model's own σ_lnD. This corresponds to ρ = 1
  cross-model residual correlation. Independent residuals per model
  (ρ = 0) underestimates the joint ensemble variability;
  intermediate cross-model correlations are an open refinement
  direction.
- Independent site-factor and Sa uncertainty within each realisation.
- Log-log linear inversion with boundary extrapolation. Per-model
  calibration ranges (`getKyLimits()`) are reported as a diagnostic
  but not enforced as a clamp inside `invertDnDraws()`.
- Zero-displacement screening: the flexible-block papers define a
  discrete probability mass P(D = 0) (negligible displacement) ahead
  of the lognormal branch. The package intentionally omits it and
  returns the unconditional lognormal — a conservative convention
  (small displacements are over- rather than under-predicted).
- Near-fault pulse behaviour (BM19) is a declared attribute of the
  record/scenario (`pulse` argument, default ordinary), not an
  automatic PGV-threshold switch.

## Operational defaults

- Monte Carlo sample size *N*_S: routine practice is
  *N*_S ∈ \[10³, 10⁴\]. Convergence diagnostics (band-width
  stability under *N*_S × 2) accompany each application.
- *k*_y grid points: *N*_k = 30, log-spaced over
  \[0.01, max(PGA, 0.80)\] g (default of `getDnKy()`).
- Reported quantiles: `invertDnDraws()` defaults to the mean plus
  *p* ∈ {0.05, 0.10, 0.16, 0.84, 0.90, 0.95}; the median and any
  other levels are requested via its `p` argument.

---

## References

- Baker & Jayaram (2008). *Earthquake Spectra* 24(1):299–317.
- Bray & Travasarou (2007). *J. Geotech. Geoenviron. Eng.* 133(4):381–392.
- Bray, Macedo & Travasarou (2018). *J. Geotech. Geoenviron. Eng.* 144(3):04017124.
- Bray & Macedo (2019, corr. 2023). *J. Geotech. Geoenviron. Eng.* 145(12); *Soil Dyn. Earthq. Eng.* 168:107835.
- Jibson (2007). *Engineering Geology* 91(2–4):209–218.
- Saygili & Rathje (2008). *J. Geotech. Geoenviron. Eng.* 134(6):790–803.
- Stewart et al. (2020). *Earthquake Spectra* 36(1). doi:10.1177/8755293019878185.
- Hashash et al. (2020). *Earthquake Spectra* 36(1). doi:10.1177/8755293019878193.
- Hashash et al. (2017). PEER Report 2017/05.