Paired realizations of a site spectrum from a rock quantile spectrum

Draws, at every oscillator period, NS paired realizations of the rock PGA, the rock ordinate and the site ordinate from a rock spectrum given by its mean and quantiles and from the site factor of evaluateSiteFactor(). The realizations are returned before any summary; reduce them per period with summarizeProfiles() (see the example). integrateSiteSpectrum() evaluates the same model without random numbers.

Usage

simulateSiteSpectrum(spectrum, VS30, Vref, model, NS)

Arguments

spectrum: The rock quantile spectrum, see integrateSiteSpectrum().
VS30: One target velocity in m/s.
Vref: One reference velocity in m/s; 760 or 3000 unless it equals VS30.
model: Site factor model of evaluateSiteFactor(), "ST20" or "ST17".
NS: One positive integer: realizations per period.

Value

A data.table with one row per period and realization: Tn (s), realization (1 to NS), PGA, Sa and SaF (g; rock PGA, rock ordinate and site ordinate of the realization), in the order drawn. Advances the random number stream.

Mathematics

At PGA and at each period the rock quantile function $Q(u) = \ln S_a$ is that of integrateSiteSpectrum(): the Hyman monotone spline through the tabulated knots $(p_k, \ln S_{a,k})$, continued below the first and above the last knot by the adjacent secant slopes (floored at zero) and clamped to $u \in [0, 1]$. It is non-decreasing, so every copula uniform maps to the rock quantile of its rank. At period $T_n$, $(Z_{0j}, Z_{nj})$, $j = 1, \ldots, n_s$, are drawn from the bivariate normal with unit variances and correlation $\rho_{BJ}(T_n, 0)$ of correlateResiduals() (mvtnorm::rmvnorm()), a Gaussian copula between PGA and the ordinate, and each marginal sample is re-centred on its tabulated mean: $$PGA_j = e^{Q_0(\Phi(Z_{0j}))} \frac{\overline{PGA}}{\langle e^{Q_0(\Phi(Z_0))} \rangle}, \quad S_{a,j} = e^{Q_n(\Phi(Z_{nj}))} \frac{\bar S_a}{\langle e^{Q_n(\Phi(Z_n))} \rangle},$$ with $\langle \cdot \rangle$ the sample mean. Then $\ln F_j \sim N(\mu_F(PGA_j), \sigma_F(PGA_j))$ independently and $S_{aF,j} = S_{a,j} F_j$. When VS30 equals Vref, for any reference velocity, $F$ is the unit factor ($\mu_F = 0$, $\sigma_F = 0$) and $S_{aF,j} = S_{a,j}$.

The triple $(PGA_j, S_{a,j}, S_{aF,j})$ is paired within one period only: every period draws its own copula sample, so rows with the same realization at different periods do not form one spectrum. Periods are drawn in increasing order, each consuming from the caller's random stream the NS copula draws and then, where the factor dispersion is positive, the NS factor residuals (stats::rnorm draws nothing where the dispersion is zero, as ST20 gives from 3000 m/s against a 760 m/s reference); when VS30 equals Vref no factor residual is drawn. Sample statistics carry Monte Carlo error of order $n_s^{-1/2}$.

Examples

P <- c("0.05", "0.10", "0.16", "0.50", "0.84", "0.90", "0.95", "mean")
Factor <- c(0.45, 0.60, 0.72, 0.95, 1.30, 1.48, 1.78, 1.00)
Spectrum <- data.table::data.table(Tn = rep(c(0, 0.2, 1), each = 8), p = rep(P, 3),
                                   Sa = rep(c(0.30, 0.45, 0.24), each = 8) * rep(Factor, 3))
set.seed(919)
Draws <- simulateSiteSpectrum(Spectrum, VS30 = 270, Vref = 760, model = "ST20", NS = 300)
Draws[, summarizeProfiles(.SD, levels = c(0.16, 0.5, 0.84, "mean"), properties = "SaF"), by = Tn]
#>        Tn  level     n       SaF
#>     <num> <char> <int>     <num>
#>  1:   0.0   0.16   300 0.1741885
#>  2:   0.0    0.5   300 0.2736075
#>  3:   0.0   0.84   300 0.4255491
#>  4:   0.0   mean   300 0.2971330
#>  5:   0.2   0.16   300 0.3580639
#>  6:   0.2    0.5   300 0.5169956
#>  7:   0.2   0.84   300 0.7403857
#>  8:   0.2   mean   300 0.5536243
#>  9:   1.0   0.16   300 0.3110117
#> 10:   1.0    0.5   300 0.4665861
#> 11:   1.0   0.84   300 0.7124517
#> 12:   1.0   mean   300 0.5091320