Fit a period-wise lognormal approximation to one spectrum

Uses OpenQuake's reported mean as the authoritative spectral ordinate and the seven reported quantiles only to calibrate the natural-log standard deviation at each oscillator period. The quantiles are not returned and do not remain part of the fitted representation.

Usage

fitSaLognormal(uhs)

Arguments

uhs

A data.table representing one spectrum. It must contain columns Tn, p, and Sa. Every period must contain exactly one "mean" row and exactly the seven probability levels 0.05, 0.10, 0.16, 0.50, 0.84, 0.90, and 0.95 (equivalent text forms such as "0.5" and "0.50" are accepted). A Tn == 0 PGA period must be present. Sa values must be finite and strictly positive, and the quantiles must be non-decreasing with probability.

Value

A data.table, ordered by Tn, with exactly the columns Tn, Sa, muLn, and sigmaLn. Sa is copied from the input mean; sigmaLn is the period-specific fitted value. No input quantile is returned. The function is deterministic and does not modify uhs.

Details

For each period, sigmaLn is obtained by a mean-anchored Huber fit. The natural-log location is then

$$\mu_{\ln}(T_n)=\ln\{\operatorname{mean} Sa(T_n)\} -\sigma_{\ln}(T_n)^2/2,$$

which preserves every supplied mean exactly at the population level.

Examples

p <- c(0.05, 0.10, 0.16, 0.50, 0.84, 0.90, 0.95)
sigmaLn <- 0.4
UHS <- data.table::CJ(Tn = c(0, 0.2, 1), p = sprintf("%.2f", p))
UHS[, Sa := c(0.5, 0.8, 0.2)[match(Tn, c(0, 0.2, 1))] *
      exp(sigmaLn * stats::qnorm(as.numeric(p)) - sigmaLn^2 / 2)]
#> Key: <Tn, p>
#>        Tn      p         Sa
#>     <num> <char>      <num>
#>  1:   0.0   0.05 0.23904854
#>  2:   0.0   0.10 0.27643825
#>  3:   0.0   0.16 0.31007833
#>  4:   0.0   0.50 0.46155817
#>  5:   0.0   0.84 0.68703913
#>  6:   0.0   0.90 0.77064569
#>  7:   0.0   0.95 0.89118278
#>  8:   0.2   0.05 0.38247767
#>  9:   0.2   0.10 0.44230120
#> 10:   0.2   0.16 0.49612533
#> 11:   0.2   0.50 0.73849308
#> 12:   0.2   0.84 1.09926261
#> 13:   0.2   0.90 1.23303311
#> 14:   0.2   0.95 1.42589245
#> 15:   1.0   0.05 0.09561942
#> 16:   1.0   0.10 0.11057530
#> 17:   1.0   0.16 0.12403133
#> 18:   1.0   0.50 0.18462327
#> 19:   1.0   0.84 0.27481565
#> 20:   1.0   0.90 0.30825828
#> 21:   1.0   0.95 0.35647311
#>        Tn      p         Sa
#>     <num> <char>      <num>
UHS <- data.table::rbindlist(list(
  UHS,
  data.table::data.table(
    Tn = c(0, 0.2, 1), p = "mean", Sa = c(0.5, 0.8, 0.2)
  )
))
fitSaLognormal(uhs = UHS)
#>       Tn    Sa       muLn sigmaLn
#>    <num> <num>      <num>   <num>
#> 1:   0.0   0.5 -0.7731472     0.4
#> 2:   0.2   0.8 -0.3031436     0.4
#> 3:   1.0   0.2 -1.6894379     0.4