Correlated sample of spectral accelerations

Draws NS joint realizations of the spectral accelerations at the oscillator periods Tn through a Gaussian copula: the dependence comes from the correlation matrix of the log residuals and each marginal from the spectrum, given either as tabulated quantiles or as an analytic lognormal mixture. The rows that share realization form one paired realization of every ordinate.

Usage

sampleSpectra(spectrum, Tn, correlation, NS, zTrunc)

Arguments

spectrum

A quantile spectrum (Tn, p, Sa in g, with at least three quantile rows at every requested period, each with a probability strictly between 0 and 1 and a finite positive Sa, and optionally a "mean" row with a finite positive Sa), or an analytic mixture of \(\ln S_a\) (Tn, k, muLn, sigmaLn, w, with rows at Tn == 0 and weights summing to one at every period).

Tn

Requested oscillator periods in seconds, finite and non-negative, in any order; the result lists them in this order.

correlation

The symmetric correlation matrix of the log residuals at Tn, with unit diagonal, for example from correlateResiduals(). Correlations given only between a reference period Tn[1] and the others, a vector rho, form the matrix C <- diag(length(Tn)); C[1, -1] <- C[-1, 1] <- rho.

NS

Number of realizations, one positive integer.

zTrunc

Symmetric truncation of the copula marginals in standard deviations, one finite positive number for an analytic mixture, or NA for untruncated sampling.

Value

A data.table of paired realizations with realization (realization, 1 to NS), Tn (s) and Sa (g), one row per realization and requested period, ordered by the position of the period in Tn and then by realization. The rows that share realization form one joint draw.

Mathematics

Standard normal vectors \(Z_s \sim N(0, C)\), \(s = 1, \dots, N_S\), are drawn with mvtnorm::rmvnorm (eigen factorization), after adding \(|\lambda_{min}| + 10^{-8}\) to the diagonal of \(C\) when its smallest eigenvalue is not positive, and mapped to uniforms \(U = \Phi(Z)\). With a truncation \(z_t\), the uniforms are remapped to \(\Phi(-z_t) + U\,(\Phi(z_t) - \Phi(-z_t))\), which keeps their ranks and restricts them to \([\Phi(-z_t), \Phi(z_t)]\). Each column is then inverted through the quantile function \(Q\) of \(\ln S_a\) of its marginal and rescaled to its mean, $$S_{a,s} = e^{Q(U_s)}\, \frac{\bar S_a}{\frac{1}{N_S} \sum_{r=1}^{N_S} e^{Q(U_r)}};$$ this finite-sample recentering multiplies every draw of one period by one positive constant: it keeps their ranks, but no fixed bound on \(S_a\).

For a quantile spectrum, \(Q\) is buildQuantileSpline() of the rows at that exact period and \(\bar S_a\) is its "mean" row; a period without a mean row takes the sample mean, which leaves the draws at \(e^{Q(U_s)}\). A quantile spectrum is sampled untruncated, so \(U\) spans (0, 1) and reaches the end pieces that the spline extrapolates beyond its first and last tabulated probabilities. For a mixture \(F(x) = \sum_k w_k \Phi((x - \mu_k)/\sigma_k)\) of \(\ln S_a\), \(Q\) is the bisection inverse of the untruncated \(F\), so with a truncation the draws of a tabulated period lie, before the recentering, between the mixture quantiles of probability \(\Phi(-z_t)\) and \(\Phi(z_t)\). The mean is \(\sum_k w_k e^{\mu_k + \sigma_k^2/2}\), or with truncation that of the components truncated individually at \(\pm z_t\), $$\sum_k w_k e^{\mu_k + \sigma_k^2/2} \frac{\Phi(z_t - \sigma_k) - \Phi(-z_t - \sigma_k)}{\Phi(z_t) - \Phi(-z_t)},$$ while the draws invert the untruncated mixture over the restricted probabilities. A mixture period that is not tabulated interpolates \(Q\) linearly and the mean geometrically in the coordinate \(\ln(T + 10^{-8})\) of the tabulated periods, the requested period entering as \(\ln T\) (\(T > 0\)) or \(-\infty\) (\(T = 0\)), between the two neighbouring tabulated periods, and takes the end period outside them.

Validation, including the quantile and mean rows of every requested period, happens before any draw, so an invalid request leaves the random stream untouched; a valid one consumes NS multivariate normal draws of the caller's stream.

Examples

Spectrum <- data.frame(Tn = rep(c(0, 1), each = 4),
                       p = rep(c("0.16", "0.50", "0.84", "mean"), 2),
                       Sa = c(0.30, 0.48, 0.72, 0.50, 0.10, 0.19, 0.31, 0.20))
C <- diag(2)
C[1, 2] <- C[2, 1] <- correlateResiduals(Tn = 1, T0 = 0.01)
set.seed(1)
Sample <- sampleSpectra(Spectrum, Tn = c(0, 1), correlation = C, NS = 20, zTrunc = NA)
Sample[, .(mean = mean(Sa)), by = Tn]
#>       Tn  mean
#>    <num> <num>
#> 1:     0   0.5
#> 2:     1   0.2