---
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`)

Rock-level spectral ordinates are amplified to the target site Vs30 using
the Seyhan & Stewart (2014) nonlinear site-factor model (ST17).

- **Mode A** (mean-only input): rock Sa is deterministic; dispersion comes
  from ST17 sigma alone.
- **Mode B** (quantile input): Sa is sampled jointly across periods via
  Gaussian copula; dispersion reflects combined hazard and site-factor
  uncertainty.

See `?fitSaF` 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()`.

## 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: mean plus *p* ∈ {0.16, 0.50, 0.84}.
  Extreme-tail reporting (*p* ∈ {0.05, 0.95}) is requested
  explicitly via the `p` argument to `invertDnDraws()`.

---

## 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.
- Seyhan & Stewart (2014). *Earthquake Spectra* 30(3):1241–1256.