Quantiles and mean of a spectral mixture

Reduces a spectrum represented as a finite lognormal mixture at every oscillator period, such as the rock spectra of the hazard products or the result of integrateSiteMixture(), to the spectral accelerations at the requested probabilities and to the mean, period by period.

Usage

summarizeMixture(mixture, levels)

Arguments

mixture

A mixture as integrateSiteMixture() takes and returns it, validated by the same rules: one row per (Tn, k) with columns Tn (s), k (component label), muLn, sigmaLn (finite, sigmaLn non-negative) and w (positive weights summing to one at every period), including the PGA rows (Tn == 0).

levels

Probabilities in \((0, 1)\).

Value

A data.table ordered by Tn with Tn, level (the probability as text, then "mean") and Sa (g), one row per period and level, unrounded.

Mathematics

At each period the weights are normalized, \(\tilde w_k = w_k / \sum w\), and the mixture distribution function is \(F(x) = \sum_k \tilde w_k \Phi((x - \mu_k)/\sigma_k)\) for \(x = \ln S_a\), a component with \(\sigma_k = 0\) contributing the step \(1\{x \ge \mu_k\}\). The quantile at probability \(p\) is \(\exp(x_p)\) with \(x_p = \inf\{x : F(x) \ge p\}\), the smallest value whose distribution function reaches \(p\) (generalized inverse). When every dispersion is positive, \(F\) is continuous and strictly increasing and \(x_p\) is the root of \(F(x) = p\), found by stats::uniroot() to \(10^{-12}\) on \([\min(\mu - 12 s), \max(\mu + 12 s)]\), \(s = \max(\sigma, 1)\); otherwise \(F\) jumps at the atoms, \(F(x) = p\) may have no root, and \(x_p\) is the upper end of 100 bisections of that interval, each keeping \(F \ge p\) at the upper end, so that a probability inside a jump returns the atom. The mean is \(\sum_k \tilde w_k e^{\mu_k + \sigma_k^2/2}\). These are marginal statistics of the ordinate at each period; they do not describe one spectrum across periods.

Examples

Mixture <- data.table::CJ(Tn = c(0, 0.2, 1), k = c("a", "b"))
Mixture[, `:=`(muLn = log(0.4) + 0.1 * (k == "b") - 0.5 * Tn, sigmaLn = 0.5, w = 0.5)]
#> Key: <Tn, k>
#>       Tn      k       muLn sigmaLn     w
#>    <num> <char>      <num>   <num> <num>
#> 1:   0.0      a -0.9162907     0.5   0.5
#> 2:   0.0      b -0.8162907     0.5   0.5
#> 3:   0.2      a -1.0162907     0.5   0.5
#> 4:   0.2      b -0.9162907     0.5   0.5
#> 5:   1.0      a -1.4162907     0.5   0.5
#> 6:   1.0      b -1.3162907     0.5   0.5
summarizeMixture(Mixture, levels = c(0.16, 0.50, 0.84))
#> Key: <Tn>
#>        Tn  level        Sa
#>     <num> <char>     <num>
#>  1:   0.0   0.16 0.2551234
#>  2:   0.0    0.5 0.4205084
#>  3:   0.0   0.84 0.6931052
#>  4:   0.0   mean 0.4770942
#>  5:   0.2   0.16 0.2308452
#>  6:   0.2    0.5 0.3804918
#>  7:   0.2   0.84 0.6271475
#>  8:   0.2   mean 0.4316927
#>  9:   1.0   0.16 0.1547402
#> 10:   1.0    0.5 0.2550513
#> 11:   1.0   0.84 0.4203895
#> 12:   1.0   mean 0.2893723