From a373184e8d1e9af0a0690545e5b98592c5e05c14 Mon Sep 17 00:00:00 2001 From: Ozan Adiguzel Date: Sun, 10 Oct 2021 05:06:09 -0400 Subject: [PATCH 1/3] extend gpdfit to return marginal posterior of k #186 --- DESCRIPTION | 2 +- R/gpdfit.R | 23 +++++++++-------------- R/psis.R | 21 ++++++++++----------- R/psislw.R | 4 ++-- man/gpdfit.Rd | 2 +- tests/testthat/test_gpdfit.R | 6 +++--- 6 files changed, 26 insertions(+), 32 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index d0cd9d35..072c5155 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -51,5 +51,5 @@ Suggests: VignetteBuilder: knitr Encoding: UTF-8 SystemRequirements: pandoc (>= 1.12.3), pandoc-citeproc -RoxygenNote: 7.1.1 +RoxygenNote: 7.1.2 Roxygen: list(markdown = TRUE) diff --git a/R/gpdfit.R b/R/gpdfit.R index 7cd46ea3..6ecc757c 100644 --- a/R/gpdfit.R +++ b/R/gpdfit.R @@ -17,7 +17,7 @@ #' algorithm is to sort the elements of `x`. If `x` is already #' sorted in ascending order then `sort_x` can be set to `FALSE` to #' skip the initial sorting step. -#' @return A named list with components `k` and `sigma`. +#' @return A named list with components `k_hat`, `sigma_hat`, `k` and `w_theta`. #' #' @details Here the parameter \eqn{k} is the negative of \eqn{k} in Zhang & #' Stephens (2009). @@ -39,32 +39,27 @@ gpdfit <- function(x, wip = TRUE, min_grid_pts = 30, sort_x = TRUE) { jj <- seq_len(M) xstar <- x[floor(N / 4 + 0.5)] # first quartile of sample theta <- 1 / x[N] + (1 - sqrt(M / (jj - 0.5))) / prior / xstar - l_theta <- N * lx(theta, x) # profile log-lik + k <- matrixStats::rowMeans2(log1p(-theta %o% x)) + l_theta <- N * (log(-theta / k) - k - 1) # profile log-lik w_theta <- exp(l_theta - matrixStats::logSumExp(l_theta)) # normalize theta_hat <- sum(theta * w_theta) - k <- mean.default(log1p(-theta_hat * x)) - sigma <- -k / theta_hat + k_hat <- mean.default(log1p(-theta_hat * x)) + sigma_hat <- -k_hat / theta_hat if (wip) { - k <- adjust_k_wip(k, n = N) + k_hat <- adjust_k_wip(k_hat, n = N) } - if (is.nan(k)) { - k <- Inf + if (is.nan(k_hat)) { + k_hat <- Inf } - nlist(k, sigma) + nlist(k_hat, sigma_hat, k, w_theta) } # internal ---------------------------------------------------------------- -lx <- function(a,x) { - a <- -a - k <- vapply(a, FUN = function(a_i) mean(log1p(a_i * x)), FUN.VALUE = numeric(1)) - log(a / k) - k - 1 -} - #' Adjust k based on weakly informative prior, Gaussian centered on 0.5. This #' will stabilize estimates for very small Monte Carlo sample sizes and low neff #' cases. diff --git a/R/psis.R b/R/psis.R index e82a1ef9..013ef12e 100644 --- a/R/psis.R +++ b/R/psis.R @@ -97,10 +97,10 @@ psis.array <- function(log_ratios, ..., r_eff = NULL, cores = getOption("mc.cores", 1)) { - importance_sampling.array(log_ratios = log_ratios, ..., - r_eff = r_eff, - cores = cores, - method = "psis") + importance_sampling.array(log_ratios = log_ratios, ..., + r_eff = r_eff, + cores = cores, + method = "psis") } @@ -245,14 +245,14 @@ psis_smooth_tail <- function(x, cutoff) { # save time not sorting since x already sorted fit <- gpdfit(exp(x) - exp_cutoff, sort_x = FALSE) - k <- fit$k - sigma <- fit$sigma + k <- fit$k_hat + sigma <- fit$sigma_hat if (is.finite(k)) { - p <- (seq_len(len) - 0.5) / len - qq <- qgpd(p, k, sigma) + exp_cutoff - tail <- log(qq) + p <- (seq_len(len) - 0.5) / len + qq <- qgpd(p, k, sigma) + exp_cutoff + tail <- log(qq) } else { - tail <- x + tail <- x } list(tail = tail, k = k) } @@ -385,4 +385,3 @@ throw_psis_r_eff_warning <- function() { call. = FALSE ) } - diff --git a/R/psislw.R b/R/psislw.R index 9ed5ab15..b3d7753f 100644 --- a/R/psislw.R +++ b/R/psislw.R @@ -73,8 +73,8 @@ psislw <- function(lw, wcp = 0.2, wtrunc = 3/4, tail_ord <- order(x_tail) exp_cutoff <- exp(cutoff) fit <- gpdfit(exp(x_tail) - exp_cutoff, wip=FALSE, min_grid_pts = 80) - k <- fit$k - sigma <- fit$sigma + k <- fit$k_hat + sigma <- fit$sigma_hat prb <- (seq_len(tail_len) - 0.5) / tail_len qq <- qgpd(prb, k, sigma) + exp_cutoff smoothed_tail <- rep.int(0, tail_len) diff --git a/man/gpdfit.Rd b/man/gpdfit.Rd index 8cb61d12..d2eee365 100644 --- a/man/gpdfit.Rd +++ b/man/gpdfit.Rd @@ -21,7 +21,7 @@ sorted in ascending order then \code{sort_x} can be set to \code{FALSE} to skip the initial sorting step.} } \value{ -A named list with components \code{k} and \code{sigma}. +A named list with components \code{k_hat}, \code{sigma_hat}, \code{k} and \code{w_theta}. } \description{ Given a sample \eqn{x}, Estimate the parameters \eqn{k} and \eqn{\sigma} of diff --git a/tests/testthat/test_gpdfit.R b/tests/testthat/test_gpdfit.R index 97c57555..1397b886 100644 --- a/tests/testthat/test_gpdfit.R +++ b/tests/testthat/test_gpdfit.R @@ -5,13 +5,13 @@ context("generalized pareto") test_that("gpdfit returns correct result", { set.seed(123) x <- rexp(100) - gpdfit_val_old <- unlist(gpdfit(x, wip=FALSE, min_grid_pts = 80)) + gpdfit_val_old <- unlist(gpdfit(x, wip=FALSE, min_grid_pts = 80)[1:2]) expect_equal_to_reference(gpdfit_val_old, "reference-results/gpdfit_old.rds") - gpdfit_val_wip <- unlist(gpdfit(x, wip=TRUE, min_grid_pts = 80)) + gpdfit_val_wip <- unlist(gpdfit(x, wip=TRUE, min_grid_pts = 80)[1:2]) expect_equal_to_reference(gpdfit_val_wip, "reference-results/gpdfit.rds") - gpdfit_val_wip_default_grid <- unlist(gpdfit(x, wip=TRUE)) + gpdfit_val_wip_default_grid <- unlist(gpdfit(x, wip=TRUE)[1:2]) expect_equal_to_reference(gpdfit_val_wip_default_grid, "reference-results/gpdfit_default_grid.rds") }) From cdae01d8c4b04f8c52a20051725eee68ee6c19e7 Mon Sep 17 00:00:00 2001 From: Ozan Adiguzel Date: Wed, 13 Oct 2021 10:18:44 -0400 Subject: [PATCH 2/3] return theta densities --- R/gpdfit.R | 12 +++++++++--- man/gpdfit.Rd | 2 +- tests/testthat/reference-results/gpdfit.rds | Bin 91 -> 2033 bytes .../reference-results/gpdfit_default_grid.rds | Bin 91 -> 1009 bytes tests/testthat/reference-results/gpdfit_old.rds | Bin 90 -> 2032 bytes tests/testthat/test_gpdfit.R | 6 +++--- 6 files changed, 13 insertions(+), 7 deletions(-) diff --git a/R/gpdfit.R b/R/gpdfit.R index 6ecc757c..ace718c1 100644 --- a/R/gpdfit.R +++ b/R/gpdfit.R @@ -17,7 +17,7 @@ #' algorithm is to sort the elements of `x`. If `x` is already #' sorted in ascending order then `sort_x` can be set to `FALSE` to #' skip the initial sorting step. -#' @return A named list with components `k_hat`, `sigma_hat`, `k` and `w_theta`. +#' @return A named list with components `k_hat`, `sigma_hat`, `k`, and `d_theta`. #' #' @details Here the parameter \eqn{k} is the negative of \eqn{k} in Zhang & #' Stephens (2009). @@ -50,11 +50,17 @@ gpdfit <- function(x, wip = TRUE, min_grid_pts = 30, sort_x = TRUE) { k_hat <- adjust_k_wip(k_hat, n = N) } - if (is.nan(k_hat)) { + if (is.na(k_hat)) { k_hat <- Inf + sigma_hat <- NaN } - nlist(k_hat, sigma_hat, k, w_theta) + # estimate the normalization term with the trapezoidal rule + c <- 2 / sum((w_theta[-M] + w_theta[-1]) * (k[-M] - k[-1])) + # convert normalized quadrature weights to normalized densities + d_theta <- c * w_theta + + nlist(k_hat, sigma_hat, k, d_theta) } diff --git a/man/gpdfit.Rd b/man/gpdfit.Rd index d2eee365..952815c4 100644 --- a/man/gpdfit.Rd +++ b/man/gpdfit.Rd @@ -21,7 +21,7 @@ sorted in ascending order then \code{sort_x} can be set to \code{FALSE} to skip the initial sorting step.} } \value{ -A named list with components \code{k_hat}, \code{sigma_hat}, \code{k} and \code{w_theta}. +A named list with components \code{k_hat}, \code{sigma_hat}, \code{k}, \code{d_theta}. } \description{ Given a sample \eqn{x}, Estimate the parameters \eqn{k} and \eqn{\sigma} of diff --git a/tests/testthat/reference-results/gpdfit.rds b/tests/testthat/reference-results/gpdfit.rds index febd170f63c02c28d9f65f60b572df576d9583f6..2a89b57b0b9cbf243df9dc2692c7eddbd118035f 100644 GIT binary patch literal 2033 zcmZ9;c{~$}0|0PR_9WAb^sJ$mo}-O&D@GwlOPHatQRb+*_B_LKmgJM@nb{L7vgVj8 zk=AnzdFEW1(OfA+Q|LuMM`f5D13|iXWn{M6)$p`&kjZ77Ka!`DJwQ(s*fcQcvP1e9$6Tk= zQC$WhnXF*gUIOdA^)lb@lfv^r!#s3J>2po;e9{xcU4nkJ=iMtM7F*MrG6ORXZw@uy zm?gvp)&_)Ul^q}U$<1y#I?=p=@)fdGCq{R{A(RPVig#vJVAqSYXMKwr2&6PGIlmc) zQ|7jQ$?TdFZ`FLi^E#%DPCnZfrUQ3dkFRvHtUBG6ty0vnB0vyt`WeKJKqyYwA60-K zw}4oJeT8SnZyOFBIB#jF31+L^`DURN(Z6JNKNt|^uy{}|dCP7&B!3Bfvx2JgT@`6d zC2_Z^FewkHSf_oBgHqF*Gq`-^Ij00Pu%UPuS5*ZVV zt)>#$rgzpW3bsiTYXS7{X#})}k*JV3xe9=Li`&0g1&4If)k0|YW)GmDE-YFYT(juXR zr*6o{rD$(5ZEiP`7N`x)StWAa3`hddi=)Abi#%~8d;5twA?dOfnk9N$Nw7KjkBe!= zjUeDfqzliWutNoZK8!h*F@Izsm)cm3-`8zQWb%dU85_&(a`oxvE`mO?E1{cvO}8QK z#d|c~CY|3DTJQLdbmY)@d@Vp3vG#zxXaQWvtv+q?RULj6e(POsyjtE_U6_^+N)@RY z7zN~uO3yo@XcGgvzG{f*4PNiWNH86Zac&LZFhkAay@pWL8NM@m)F?ZO@CFy&(DKT(rZ;Hk zdXA4*G==k$P*m{SnHf{!xkW6I|G@OAv)V+P;$P(%SP|2DG}hsgGp103|Fk&rbU6#G z;aRpeUeQ#vHyJa;np6mE$Gn;g-N-s=X^f;O&;88(+JGR0{Q9bi@we13p^Kis)N73}%o-HMR=uaTryh=6AAws_(IX0ftYZrk z;v$Uk^6IcHBf`B7E84}2oxkFB@?9gVedP6#QKJzS!L7^`3Mrzw;<$~ zd5_rwPP7sec0R)O4%eP?S{u-6YkNF=msoM=*M-OPX*R`1_ga97OJgbr-2`holWd|; z(9`!L;9Iz3%@p{fnZ|wPb$3dIBH%57@Ue2%{^#(o#ww_L@1<@na^(TNJw+#W>i)hp z>fsvSYO*VQpWr%~L(}(OS3o*Ilq%R)l1Cv|l>lKflMj*>+ZHyL`ssOPrK4nk9)n|) zkha`jx9MimmuwWZG^9UR?{2#st5o_t2ME4BCEa`$Ynv-4#d!o1jJcb@yU5+MrQw=TU&fD>f6ueAP zj1>JYW6069vhpwB zB=Y*RJ7-1Y+-UhRhEc_(+C+ewRX2jTRNDyk<^vmIJNNjtzD__${IiuiY=ErBwfN5~ zcRB&m8Z|j@5eMPA(Bl8(7PR#Hk(E8EI!*4@prVx2zu|@h1^@s6 literal 91 zcmb2|=3oE==I#ec2?+^F32BKb2}x{*KOI{Bsy*Y|oSRqw#=fFftmx7@pnM`zm$A_X rg94_^YkQ8|xss&d%&BJ1)$-3u!&ik?${bOXfGiTk^SfH5z0GcCX diff --git a/tests/testthat/reference-results/gpdfit_default_grid.rds b/tests/testthat/reference-results/gpdfit_default_grid.rds index d2cf8816396a3e4f238a1413c564d68536fc659a..7b507215a47485857612bd7f6326e1e3c242a9dc 100644 GIT binary patch literal 1009 zcmV&xv;H+ zt|roVdH-F{{pY%VJn#E33?pL10s$lPa3SN75U-PK7)B^@9owKmR2KhvQY*~8`}A1- zJ8|Ij$PBA_HVyNyj_vc6^}!zpRZs0RHn6|!k+iNZ3TDnMl$?9`7ECoB=suxy!ta|q zdvgC=fJvn=Hvi3SFrHX;%J!=t{MHzR@~uC@C@$R{9DE%{dYc`u1c$ppQqywcs|<9TaJCTH)UrRUz2ZG zbiSeBuO~!0o~>)*Zw&s3$r%L z3h?Dx&4QNuxHhG#-=Xa~YxVUSSF@|i3&rutp=`>gj;}rA^=x@`d0J)e9Lom<$PLY9 z?5IO_bK2($2>bR<4K9|0th24gb8;D|i{gx2GY^X&tbfj)DASNxG*b z42~99N~3d5!5Q7_Hq$?`@TEQKf+BJp8oDajBC`Y9LRcLe;|txNHS}+<4}qWjhxQ|WM5hPE&nLM6!`xj->| z=g*p(nv=2L1us2+)}=9fmm$^o;_%*|wI?U_e}gbv=_ceR>?W0)!rc_%CbgS1ZrbpI zMA^#q1eF94K}=9ZFq~inK{Y`QL1IJ1hKLOj8zMGDY>3zpu_0nZ#D<6s6B{NrOl+9g zFtK4`!^DP(4HFwCHr0~NYxgdrA;Xxeq;5prn7UQe9Zua5)UBp&4RtT2c`x;>rw1iH zAbP;`prVIxdWfJ0H9ctPfi6XKDWXdeU5e;ZM3*AE6w#%KE=6=HqDwJdis@2JmtwjU f)1{a$#dImAOEF!F=~CBwGR)$CNmh8AJO=;(l)Uda literal 91 zcmb2|=3oE==I#ec2?+^F32BK*2}x{*KOI{Bsy*XW3|`cCtNz2UEBh`QcQi6;oIh~j t#DPNxIG03z7Ag+zZscLykk+A)Z057g!r0ba@Wc;xhJCqf7sLY11OR%HB`p8| diff --git a/tests/testthat/reference-results/gpdfit_old.rds b/tests/testthat/reference-results/gpdfit_old.rds index 337903700ca4c71125d25040f7cf5bf91a9a5f77..98cd765eae5d130b9c251615d25bcc845bee47c2 100644 GIT binary patch literal 2032 zcmZ9;c{~$}0|4-(A4-~XAJ37(hP)?hIYOjNlnt{ZhDlhC9IZARp?q?rTINp1a)n3k zYkKmqha5vx7|EVoGcxCXpU>~V{l0&FHIt#yGXyGSAG;=C^N_lethb^uJSF#GC@SGt?d(t)l}8Fbxq8Wa@S z3tmXE+G-V&$WuZ~%^329l=NOdv`MgxOcrgR7&vn8eoeKkpGiM2UqeNg1%Qz(YA#%VVnjW9QF=s`XAe z^@3?^g@^2lMAaU5_NWSHkPv>m_(q4b02jmD&)mY=R23}SQC02@=2f-Uvo5GN-yczF7Us}C1`Cb)8s%Z{o31LFD6h?BUcSIUlAFX zWSnCO`7CdPa=g-K++^euz)Vp8F?B`wS0oAmB zh@7zf(fsyJuXynDDIcxfw;`HErSTDqj@S=1X zzug^6uCT5S!2#opM|Z%{{a~n<1ANAddte2Mhv(d6Nb-H$FuI76;Ns5PyQUS=G6)N zQIA2!@%9+xE=gVDtXI2mqxpW~Q-&8w4Mx|rx;C8VvvtJ6?A#$1Wv~Q{S*pFi&-u$}!qRL(O}M_pfv*>O zBuL+i`;xN0Kq(1xfGJe$OOzGLKXsqewOw{kgxZaXt+A_qx%?tR0uVQc5i9-iV{97A zN(V5rz0;%WKX|hF7uaaH7Iu01h)yCd1+an5KLAHR|1bj(jwYAa#WfJxb9w1xK?Pm41dn(R;W-m@HqY-ra zg6+pNdP7MlI={r{wSj6vF(Mel(N@(Ug)s$-$A<#=Srt!Kj?$|s1X?HZ&XC|%=V-tG z%Yd)B;de*BUmj|V`NzAZ9Bn#~PpmJQrw6o}?C(EJtER)aSM+utDG#lMbm!uZw=D0N9zWRYv({{bkV@_gH=nPrM^0DnFyr4g%ya2i zl`b`B4s|^{7ZS-UysqpiEkvWIXwp5kXQ-#Cx>RlIS*i|I&sB@c)`u55q@D&R{a995huVN8oY} z&RXC{Y>@ Date: Thu, 27 Jan 2022 04:32:25 -0500 Subject: [PATCH 3/3] updated gpdfit and created dgpd/pgpd/qgpd/rgpd --- NAMESPACE | 4 + R/gdp.R | 102 ++++++++++++++++++ R/gpdfit.R | 60 +++-------- R/psis.R | 2 +- R/psislw.R | 2 +- man/GPD.Rd | 42 ++++++++ man/gpdfit.Rd | 2 +- tests/testthat/reference-results/gpdfit.rds | Bin 2033 -> 2998 bytes .../reference-results/gpdfit_default_grid.rds | Bin 1009 -> 1419 bytes .../testthat/reference-results/gpdfit_old.rds | Bin 2032 -> 2996 bytes tests/testthat/test_gpdfit.R | 4 +- 11 files changed, 168 insertions(+), 50 deletions(-) create mode 100644 R/gdp.R create mode 100644 man/GPD.Rd diff --git a/NAMESPACE b/NAMESPACE index 180b3e22..ed538194 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -75,6 +75,7 @@ export(.ndraws) export(.thin_draws) export(E_loo) export(compare) +export(dgpd) export(elpd) export(example_loglik_array) export(example_loglik_matrix) @@ -115,12 +116,15 @@ export(pareto_k_ids) export(pareto_k_influence_values) export(pareto_k_table) export(pareto_k_values) +export(pgpd) export(print_dims) export(pseudobma_weights) export(psis) export(psis_n_eff_values) export(psislw) +export(qgpd) export(relative_eff) +export(rgpd) export(sis) export(stacking_weights) export(tis) diff --git a/R/gdp.R b/R/gdp.R new file mode 100644 index 00000000..a89313af --- /dev/null +++ b/R/gdp.R @@ -0,0 +1,102 @@ +#' The Generalized Pareto Distribution +#' +#' Density, distribution function, quantile function and random generation +#' for the generalized Pareto distribution with location equal to \code{mu}, +#' scale equal to \code{sigma}, and shape equal to \code{k}. +#' +#' @param x,q vector of quantiles. +#' @param p vector of probabilities. +#' @param n number of observations. If \code{length(n) > 1}, the length is taken +#' to be the number required. +#' @param mu scalar location parameter +#' @param sigma scalar, positive scale parameter +#' @param k scalar shape parameter +#' @param log,log.p logical; if TRUE, probabilities p are given as log(p). +#' @param lower.tail logical; if TRUE (default), probabilities are \eqn{P[X \le x]} +#' otherwise, \eqn{P[X > x]}. +#' +#' @name GPD + + +#' @rdname GPD +#' @export +dgpd <- function(x, mu = 0, sigma = 1, k = 0, log = FALSE) { + stopifnot(length(mu) == 1 && length(sigma) == 1 && length(k) == 1) + if (is.na(sigma) || sigma <= 0) { + return(rep(NaN, length(x))) + } + d <- (x - mu) / sigma + ind <- (d > 0) & ((1 + k * d) > 0) + ind[is.na(ind)] <- FALSE + if (k == 0) { + d[ind] <- -d[ind] - log(sigma) + } else { + d[ind] <- log1p(k * d[ind]) * -(1 / k + 1) - log(sigma) + } + d[!ind] <- -Inf + if (!log) { + d <- exp(d) + } + d +} + + +#' @rdname GPD +#' @export +pgpd <- function(q, mu = 0, sigma = 1, k = 0, lower.tail = TRUE, log.p = FALSE) { + stopifnot(length(mu) == 1 && length(sigma) == 1 && length(k) == 1) + if (is.na(sigma) || sigma <= 0) { + return(rep(NaN, length(q))) + } + q <- pmax(q - mu, 0) / sigma + if (k == 0) { + p <- 1 - exp(-q) + } else { + p <- -expm1(log(pmax(1 + k * q, 0)) * -(1 / k)) + } + if (!lower.tail) { + p <- 1 - p + } + if (log.p) { + p <- log(p) + } + p +} + +#' @rdname GPD +#' @export +qgpd <- function(p, mu = 0, sigma = 1, k = 0, lower.tail = TRUE, log.p = FALSE) { + stopifnot(length(mu) == 1 && length(sigma) == 1 && length(k) == 1) + if (is.na(sigma) || sigma <= 0) { + return(rep(NaN, length(p))) + } + if (log.p) { + p <- exp(p) + } + if (!lower.tail) { + p <- 1 - p + } + if (k == 0) { + q <- mu - sigma * log1p(-p) + } else { + q <- mu + sigma * expm1(-k * log1p(-p)) / k + } + q +} + +#' @rdname GPD +#' @export +rgpd <- function(n, mu = 0, sigma = 1, k = 0) { + stopifnot( + length(n) == 1 && length(mu) == 1 && length(sigma) == 1 && length(k) == 1 + ) + if (is.na(sigma) || sigma <= 0) { + return(rep(NaN, n)) + } + if (k == 0) { + r <- mu + sigma * rexp(n) + } else { + r <- mu + sigma * expm1(-k * log(runif(n))) / k + } + r +} diff --git a/R/gpdfit.R b/R/gpdfit.R index ace718c1..1a2f7805 100644 --- a/R/gpdfit.R +++ b/R/gpdfit.R @@ -17,7 +17,7 @@ #' algorithm is to sort the elements of `x`. If `x` is already #' sorted in ascending order then `sort_x` can be set to `FALSE` to #' skip the initial sorting step. -#' @return A named list with components `k_hat`, `sigma_hat`, `k`, and `d_theta`. +#' @return A named list with components `k_hat`, `sigma_hat`, `k`, `k_w`, and `k_d` #' #' @details Here the parameter \eqn{k} is the negative of \eqn{k} in Zhang & #' Stephens (2009). @@ -29,7 +29,7 @@ #' for the generalized Pareto distribution. *Technometrics* **51**, 316-325. #' gpdfit <- function(x, wip = TRUE, min_grid_pts = 30, sort_x = TRUE) { - # See section 4 of Zhang and Stephens (2009) + # see section 4 of Zhang and Stephens (2009) if (sort_x) { x <- sort.int(x) } @@ -46,8 +46,19 @@ gpdfit <- function(x, wip = TRUE, min_grid_pts = 30, sort_x = TRUE) { k_hat <- mean.default(log1p(-theta_hat * x)) sigma_hat <- -k_hat / theta_hat + # quadrature weights for k are same as for theta + k_w <- w_theta + # quadrature weights are just the normalized likelihoods + # we get the unnormalized posterior by multiplying these by the prior + k_d <- k_w * dgpd(-theta, mu = -1 / x[N], sigma = 1 / prior / xstar, k = 0.5) + # normalize using the trapezoidal rule + Z <- sum((k_d[-M] + k_d[-1]) * (k[-M] - k[-1])) / 2 + k_d <- k_d / Z + + # adjust k_hat based on weakly informative prior, Gaussian centered on 0.5. + # this stabilizes estimates for very small Monte Carlo sample sizes and low neff if (wip) { - k_hat <- adjust_k_wip(k_hat, n = N) + k_hat <- (k_hat * N + 0.5 * 10) / (N + 10) } if (is.na(k_hat)) { @@ -55,46 +66,5 @@ gpdfit <- function(x, wip = TRUE, min_grid_pts = 30, sort_x = TRUE) { sigma_hat <- NaN } - # estimate the normalization term with the trapezoidal rule - c <- 2 / sum((w_theta[-M] + w_theta[-1]) * (k[-M] - k[-1])) - # convert normalized quadrature weights to normalized densities - d_theta <- c * w_theta - - nlist(k_hat, sigma_hat, k, d_theta) -} - - -# internal ---------------------------------------------------------------- - -#' Adjust k based on weakly informative prior, Gaussian centered on 0.5. This -#' will stabilize estimates for very small Monte Carlo sample sizes and low neff -#' cases. -#' -#' @noRd -#' @param k Scalar khat estimate. -#' @param n Integer number of tail samples used to fit GPD. -#' @return Scalar adjusted khat estimate. -#' -adjust_k_wip <- function(k, n) { - a <- 10 - n_plus_a <- n + a - k * n / n_plus_a + a * 0.5 / n_plus_a -} - - -#' Inverse CDF of generalized pareto distribution -#' (assuming location parameter is 0) -#' -#' @noRd -#' @param p Vector of probabilities. -#' @param k Scalar shape parameter. -#' @param sigma Scalar scale parameter. -#' @return Vector of quantiles. -#' -qgpd <- function(p, k, sigma) { - if (is.nan(sigma) || sigma <= 0) { - return(rep(NaN, length(p))) - } - - sigma * expm1(-k * log1p(-p)) / k + nlist(k_hat, sigma_hat, k, k_w, k_d) } diff --git a/R/psis.R b/R/psis.R index 013ef12e..34c2fa97 100644 --- a/R/psis.R +++ b/R/psis.R @@ -249,7 +249,7 @@ psis_smooth_tail <- function(x, cutoff) { sigma <- fit$sigma_hat if (is.finite(k)) { p <- (seq_len(len) - 0.5) / len - qq <- qgpd(p, k, sigma) + exp_cutoff + qq <- qgpd(p, sigma = sigma, k = k) + exp_cutoff tail <- log(qq) } else { tail <- x diff --git a/R/psislw.R b/R/psislw.R index b3d7753f..636d69ce 100644 --- a/R/psislw.R +++ b/R/psislw.R @@ -76,7 +76,7 @@ psislw <- function(lw, wcp = 0.2, wtrunc = 3/4, k <- fit$k_hat sigma <- fit$sigma_hat prb <- (seq_len(tail_len) - 0.5) / tail_len - qq <- qgpd(prb, k, sigma) + exp_cutoff + qq <- qgpd(prb, sigma = sigma, k = k) + exp_cutoff smoothed_tail <- rep.int(0, tail_len) smoothed_tail[tail_ord] <- log(qq) x_new <- x diff --git a/man/GPD.Rd b/man/GPD.Rd new file mode 100644 index 00000000..23b28a03 --- /dev/null +++ b/man/GPD.Rd @@ -0,0 +1,42 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/gdp.R +\name{GPD} +\alias{GPD} +\alias{dgpd} +\alias{pgpd} +\alias{qgpd} +\alias{rgpd} +\title{The Generalized Pareto Distribution} +\usage{ +dgpd(x, mu = 0, sigma = 1, k = 0, log = FALSE) + +pgpd(q, mu = 0, sigma = 1, k = 0, lower.tail = TRUE, log.p = FALSE) + +qgpd(p, mu = 0, sigma = 1, k = 0, lower.tail = TRUE, log.p = FALSE) + +rgpd(n, mu = 0, sigma = 1, k = 0) +} +\arguments{ +\item{x, q}{vector of quantiles.} + +\item{mu}{scalar location parameter} + +\item{sigma}{scalar, positive scale parameter} + +\item{k}{scalar shape parameter} + +\item{log, log.p}{logical; if TRUE, probabilities p are given as log(p).} + +\item{lower.tail}{logical; if TRUE (default), probabilities are \eqn{P[X \le x]} +otherwise, \eqn{P[X > x]}.} + +\item{p}{vector of probabilities.} + +\item{n}{number of observations. If \code{length(n) > 1}, the length is taken +to be the number required.} +} +\description{ +Density, distribution function, quantile function and random generation +for the generalized Pareto distribution with location equal to \code{mu}, +scale equal to \code{sigma}, and shape equal to \code{k}. +} diff --git a/man/gpdfit.Rd b/man/gpdfit.Rd index 952815c4..9398b44c 100644 --- a/man/gpdfit.Rd +++ b/man/gpdfit.Rd @@ -21,7 +21,7 @@ sorted in ascending order then \code{sort_x} can be set to \code{FALSE} to skip the initial sorting step.} } \value{ -A named list with components \code{k_hat}, \code{sigma_hat}, \code{k}, \code{d_theta}. +A named list with components \code{k_hat}, \code{sigma_hat}, \code{k}, \code{k_w}, and \code{k_d} } \description{ Given a sample \eqn{x}, Estimate the parameters \eqn{k} and \eqn{\sigma} of diff --git a/tests/testthat/reference-results/gpdfit.rds b/tests/testthat/reference-results/gpdfit.rds index 2a89b57b0b9cbf243df9dc2692c7eddbd118035f..443899850db2960db1b7b4c4eab08c30ba99f731 100644 GIT binary patch literal 2998 zcmdVP`#%$kAIEV;iz0Uw%dL}3?spMlE@e1#J>-%crDg7ALul@IYQxMWgoRua+1y4h z!(4OEJ-1A`%*_4!{dRsjk4OK&_w&Q+@qRo%Lh;NDzmCx~4XU2HCe)eyQJ==w_C~q= zNydEE%Iy`NoLS4OgMwm@ZWi0W#@HKT;`=!*Vl$` zdiDivWS|%k0jKcYmH5%t6;C8{uxt)-%}?m59c4GEa$ zX=)ee>UNbwko+)<6XKfzd2sw71r* zx(xyo06Pr}`aS@ux$(CJO7V*hzk{!8?fLHq0&?zxY>^trtwCDN?sqn%>#};|eIXy9 zQ|SUr@H~xq;wRH}t=y6htIZZS5vyYIrQ0Zb06T^0=ClR~vX84O-RWlf6pHZ9xxB)( zcQQt)|3VY{dQ0fq8>fK_ysC+bX)|$<(S_n!ek;a8BBRjAqTYNhNfhq&_9Mgs+9Y3T zaWO93rfiLjzobU)MV|4M zCY_eDG61AI{2;ivbj=1nwGXbF)$QjFPLG-s{Yg z#{QnT?_f<wASS=l;jv-{N~ibBVX`wPd&}&F&~^- z#k1@S=odp}<91?gfgdcyQ!X~m!UOzuZ^9v!-mA;$y0Wr9kx1H0o%GRs-`jh3D(Ph^ z6T`Y(?YsvCxY23_hx*~g?bVHJ2UFkV`8}@{Uvm9%4?oS$r~iSQR5BR*`FRgkud_^k z0f3*}d*Air{65-pVk>IOIiezJKbyXk+c*`y>ygCR_{pA|S8?y9k350)-bfl#lpWjT z2W3A(XEle8o-DOIlghCDSP;H6D#=Y*06H6Pkiz}_J)IBh*71;2PfT;A?DWX&W#A=c zrrQ+qnU#I9>uDL3aNe6q_mPGybE5>M+kbBv*=pu1%z$M0M|JLzblkK74Q^k(utlY8 zO%LBJ1o?G_H-mOF64<%M)0UeUMw6qQ3Ghh3wZ2OCZ^T|yNWAc~{Owril%}w>uh*!? z&&TR>h#7bLH@v`Z^PpRAGb7uw=gU0TK6=ACK4rmiW_na>nI41aJrJeADfravy82A; zUCEweMy>_$w%ka%?ht(;wc@W`-R+LD4*>3{Hj*+ZeX_DA#dpfev9ylYjJINf>WC06 z5cG4bL zYc}C_JI82^mC=o{tNuwIP%zo%-02{8Lv}B|_DIE3wjZ0of1-Vb1W&RTIz|gl#cC)B# zi-d!o_wX^T(A455PqUStKO6b8Ty%YzwO?B(4}+@s4)3C6cPBP08hPT$zDa8J3)n%< z&oYgg<~A%U$ z=TOB|wD7IRm<5GQH*bvsX%cjT!}ZdrKvbf9I0a0>$lNUS-fd~@jds#=NA>x_ziPS2gMf#w!aYY)bevh=99Ut? z0S2v*(rL_+btWQU?ueo@e$6Y^hM2&>KbCe%G(%YX){Pg9 zXs(W!wfgL_ETE@6W#^oeJ8}7WTQhvF*`@cnj2hI@S>+IlJXo*W9X$|GZkH(wI#2S*e`qqxttQ%}_%Xwcvi8%R@ zYmnSdkB_?61U%_q`fNWXh9rKfQ zoW5DE1Z023Z##3!5I$Q0+IiDY3Qcw@U;kQr6!eZB7zld?V&8V+n@!uC;tIC!t-Enj zJMIPE_PEzv20fKMCwI|LLuQp{NHQOlXK5l1#zmM(W~0Rlltl_zBu(!^^^CbKuqi_k z|4^)isXkQHnA;p1Go+6e%U9+s{6(yV3g`|9KM4A-Ay%ZUP#7y|Y65jJ=C;IQhCZWz zQ5Cwwf*%Z(F8tTPMhyLd7Rys+DZD3Pst4sVMnV3eLS@Oqdy=NcPzz&}1r|M2fEN27 znI0w=tuB4{Le_G4a89yLy42X!+;{>WJYg@TviAUm_~;x9pa^z=kjsloddeIqjvXfa z8KC$3v;st7!oNw2Gb2Die*7g7z^dw8c_y6g3hWofR3Wbj{4SN_FWUjqs?Jqs;@NCr zEiR0F0>6k0CQJnUE*3E1cEInVuFn)1lgq)$BO0cs4~b#Ez^N-DDxuGWWVysi6cPRJ lJ33~LgA*tk2GaMAVg7%{t}zZ4rts7Aazb2X?Ys;O{{bQ)y6pe} literal 2033 zcmZ9;c{~$}0|0PR_9WAb^sJ$mo}-O&D@GwlOPHatQRb+*_B_LKmgJM@nb{L7vgVj8 zk=AnzdFEW1(OfA+Q|LuMM`f5D13|iXWn{M6)$p`&kjZ77Ka!`DJwQ(s*fcQcvP1e9$6Tk= zQC$WhnXF*gUIOdA^)lb@lfv^r!#s3J>2po;e9{xcU4nkJ=iMtM7F*MrG6ORXZw@uy zm?gvp)&_)Ul^q}U$<1y#I?=p=@)fdGCq{R{A(RPVig#vJVAqSYXMKwr2&6PGIlmc) zQ|7jQ$?TdFZ`FLi^E#%DPCnZfrUQ3dkFRvHtUBG6ty0vnB0vyt`WeKJKqyYwA60-K zw}4oJeT8SnZyOFBIB#jF31+L^`DURN(Z6JNKNt|^uy{}|dCP7&B!3Bfvx2JgT@`6d zC2_Z^FewkHSf_oBgHqF*Gq`-^Ij00Pu%UPuS5*ZVV zt)>#$rgzpW3bsiTYXS7{X#})}k*JV3xe9=Li`&0g1&4If)k0|YW)GmDE-YFYT(juXR zr*6o{rD$(5ZEiP`7N`x)StWAa3`hddi=)Abi#%~8d;5twA?dOfnk9N$Nw7KjkBe!= zjUeDfqzliWutNoZK8!h*F@Izsm)cm3-`8zQWb%dU85_&(a`oxvE`mO?E1{cvO}8QK z#d|c~CY|3DTJQLdbmY)@d@Vp3vG#zxXaQWvtv+q?RULj6e(POsyjtE_U6_^+N)@RY z7zN~uO3yo@XcGgvzG{f*4PNiWNH86Zac&LZFhkAay@pWL8NM@m)F?ZO@CFy&(DKT(rZ;Hk zdXA4*G==k$P*m{SnHf{!xkW6I|G@OAv)V+P;$P(%SP|2DG}hsgGp103|Fk&rbU6#G z;aRpeUeQ#vHyJa;np6mE$Gn;g-N-s=X^f;O&;88(+JGR0{Q9bi@we13p^Kis)N73}%o-HMR=uaTryh=6AAws_(IX0ftYZrk z;v$Uk^6IcHBf`B7E84}2oxkFB@?9gVedP6#QKJzS!L7^`3Mrzw;<$~ zd5_rwPP7sec0R)O4%eP?S{u-6YkNF=msoM=*M-OPX*R`1_ga97OJgbr-2`holWd|; z(9`!L;9Iz3%@p{fnZ|wPb$3dIBH%57@Ue2%{^#(o#ww_L@1<@na^(TNJw+#W>i)hp z>fsvSYO*VQpWr%~L(}(OS3o*Ilq%R)l1Cv|l>lKflMj*>+ZHyL`ssOPrK4nk9)n|) zkha`jx9MimmuwWZG^9UR?{2#st5o_t2ME4BCEa`$Ynv-4#d!o1jJcb@yU5+MrQw=TU&fD>f6ueAP zj1>JYW6069vhpwB zB=Y*RJ7-1Y+-UhRhEc_(+C+ewRX2jTRNDyk<^vmIJNNjtzD__${IiuiY=ErBwfN5~ zcRB&m8Z|j@5eMPA(Bl8(7PR#Hk(E8EI!*4@prVx2zu|@h1^@s6 diff --git a/tests/testthat/reference-results/gpdfit_default_grid.rds b/tests/testthat/reference-results/gpdfit_default_grid.rds index 7b507215a47485857612bd7f6326e1e3c242a9dc..65497c17a97ec3bb3c8abd48ab67045ebad85282 100644 GIT binary patch literal 1419 zcmbWrdpr{e0Kjoh6vmdv8S|)>ypAh23{ixTFe=O_5wanaN28BPn8s&DLz8)b&RaW$ zIh{w|aVoStn)f`yCYv4G-RJK9?XUa&@%!Vem32()-y~eg4`F!;;Qj!URBzwtuxaIx zFow%+XJ6&W!BcdtsS`cdCA?BddYL!l>_x!1z_FiVuMs?>6?};<bFeh)ua zLAPn^$w*|4Y<3#O^ReNdRPOBU+#%nO9IL&{7C6uL>yJtMR|2&)F^A$7vgM82v6?B_ z3HkR+g?2lB11NO_KO}Vbnedv>`{;`w)y7GPIF>X$3k+Gc8-By>%IrHsgORGvTY5+V zhbfXluKQwq#gC$abJlXiE?nyLlihiL_Tyl)XiZyTwy}FB{c(G4T>F9iO#1YxU0IYu zmEG4tPBqYeGrg#3opR`!Mf<#Hmi`NhOLHFD9;Bn?Dc8Xsr!oWG*}Z^xT zHO)d37t%IIc(Y_8;cZ(1R@;jAv}~!v?Giesc%f{*?gqopUteF6ZafRLM*1NrKJKb- z?rtrY-Xw7*#nWH$Wq*<>yvo@&3zr}6R0jo|$?4zVIr;3FyzSuylj@S%P&EyvIog#* z5?Nn)&6X8*vH`6vbRAb=IR7CMlRr`KeTqL=$DOh?iB7d*b{-uXB&=D`mf=;i&^qmGjT#EPZa_gXk^CCM(b?7Z~WX;=g{X{{amm<4QnpVOtO)ZWH+TlH}xBP<6lx3 zmqM23QZ{#Q|Hg1_(D*@)KXKmxSC!KC0sc-e%;{{1TmS80-Q$;Evkmj-2GYa~Jim6V zBxW>>fmO~|h+B?eC#U1;(*UkhGNnm+OH`S-8A!3@YR8B#{ff9z8Z`Q+;k!;>{El(bnyQ8?WB zLRiw_W~u0;_JalE&0_!|x6L6v zOz>!5l<5=tIV$lVqKhz!M-J>|kAKKe?VZf9&%lDFp|9)?eQHJHaEO*I?vMSOJGY3y zv1Ni^WlKkGxZ7vF=huqwg@m$5c658RBWgT6ci-GuEMQRD_ i@W`r;xlNZqdeki{`Tx*hXfi=fq&0u?TjaE~nAqQcyvQK{ literal 1009 zcmV&xv;H+ zt|roVdH-F{{pY%VJn#E33?pL10s$lPa3SN75U-PK7)B^@9owKmR2KhvQY*~8`}A1- zJ8|Ij$PBA_HVyNyj_vc6^}!zpRZs0RHn6|!k+iNZ3TDnMl$?9`7ECoB=suxy!ta|q zdvgC=fJvn=Hvi3SFrHX;%J!=t{MHzR@~uC@C@$R{9DE%{dYc`u1c$ppQqywcs|<9TaJCTH)UrRUz2ZG zbiSeBuO~!0o~>)*Zw&s3$r%L z3h?Dx&4QNuxHhG#-=Xa~YxVUSSF@|i3&rutp=`>gj;}rA^=x@`d0J)e9Lom<$PLY9 z?5IO_bK2($2>bR<4K9|0th24gb8;D|i{gx2GY^X&tbfj)DASNxG*b z42~99N~3d5!5Q7_Hq$?`@TEQKf+BJp8oDajBC`Y9LRcLe;|txNHS}+<4}qWjhxQ|WM5hPE&nLM6!`xj->| z=g*p(nv=2L1us2+)}=9fmm$^o;_%*|wI?U_e}gbv=_ceR>?W0)!rc_%CbgS1ZrbpI zMA^#q1eF94K}=9ZFq~inK{Y`QL1IJ1hKLOj8zMGDY>3zpu_0nZ#D<6s6B{NrOl+9g zFtK4`!^DP(4HFwCHr0~NYxgdrA;Xxeq;5prn7UQe9Zua5)UBp&4RtT2c`x;>rw1iH zAbP;`prVIxdWfJ0H9ctPfi6XKDWXdeU5e;ZM3*AE6w#%KE=6=HqDwJdis@2JmtwjU f)1{a$#dImAOEF!F=~CBwGR)$CNmh8AJO=;(l)Uda diff --git a/tests/testthat/reference-results/gpdfit_old.rds b/tests/testthat/reference-results/gpdfit_old.rds index 98cd765eae5d130b9c251615d25bcc845bee47c2..50c6b41709d3ef394dbc0088636f5a4aec6ae0db 100644 GIT binary patch literal 2996 zcmdVP_dgYiAII^G$X+EYdqz}_!?CgwM@Av@*wnE_I5{Zm9FFZ|@6y4+F^|2Hx*6dl zD;YVij*MfEV|>52`_uRF=pXoget11TkLOzmb?L(IqcP5es$`P*5sqH*YZ84$fc_^r zi{V8fGU}7&ykmmk9Pp~Vq4AsSG(D7I=0j(18_w9$bRbOPdgEUXpK{Y5m70;}1wJNH zj`!O)6jqK%>Kn_iL)8lfOda&&q%h?0?aP^}Gxsydq4kyOg^8--<=50(Eb1BOMox=C zUI{`{zbIX2-dl8sgmq^-!sZ1Upj5`&Wer5iFxA9T$Yb4k9tr$H43DPUiM>F<;sp&& zDltEO9Ku(jx8T;%7+KSM&gkMiiY#CIxEHD5x>y|tY1(-2777{q`*SlY&WeL6I7dK( z5#qqj&$_Jg&ZyaFEbuF<*2<&PmO(=Z+1O zr+In)2j2HxJe*}Mtd)Vbp#!=$sCEZ0EtXjZ0UnDawzzK}P?KSNzAp&ihP!94_{=zK zJs^G?&sw0tZR6(h^uFr2?mK(0uV0pMth9+(VZF?q-OP|qEbh7fclmIyxqsRCG4;*> ztNgSw<$RRFyj_UhiNMZZT%iee#NmF%*mvGDFItO#SyjM+5aO|q+WiJ;8woPeuI%z< zaaF$*;Af|!;lZavJi3!F1xeLO&7Vtvp|m^;_7Ie8TK%>IT%u|jF^M2%y5E|K5I0arfmLPOuDIJW~ZvL z<;(=prJLFK-nR#~Lw^qnbxarwD54>#d3E~cUVrlWNmu%Se{|Z>Yci`b(R$4!8Mx+9Z+XYX)eV6O~!g};#2WV(kwJZt6h+Xh4^s_XlS>iZ8;IKrd=*t!cwQ}?c_Re=P_O};QEj+>>u~G4Iv&&=dACp7Xum&Cy?;w4W7-#LPCm4ZThp?l zoLRz7h4X@MD~vfwC2l`bY|m@q1RTKPJMxyamCq)`Tk}tyE=$j(u5;5%@?u;!9H!X1 z-*Cqav)VD-Z;q;SHV?D{*3i%JQ2@k_t3y{) zYX5K33yzf4FE&N*YZYBjb*2E^)IRf9Ng*n|b7P37CWmry3x)|FF`KI|Qvaf9<=4m) znO&^<_*v6P1^!_9>yR&p0EyY>*>F{CdCrnT*BuzOOO2z)1H18os-3(Hm zPTe$W;JQ78M@8OFj!vv?{5Cedn3|Sed{Y1OmYer1Uwim8r-2_J`2ZRioe|B-T}Km< z8oaZ+K$5EQe|H;D?%t{N#(P>a!C(I6QuD;V>3}^$0Cr|uf%}n;#of}uteZ+mqW-z> zPx+YCx68KRBa_<4_C|Kk3vVALTl$4Wm@#h{hcylj5_>94Cg7oXln;I zjngyvp42(u?@t!F2TnUFX2_}y+E2HTU)RZ=cuy-{6RM7Vpk@uV!o!20s&=35zH1rN z_YPe2VQXsGSj+n$`pSJ)Zhyx#Fo^Di_T~@WoQgzY)<)edZ)}=cf+bO0UH#@$B<&@~x?W7i^p`zO7{kD>Oo0q19VJ$Mfl2_I5 zB`K*rCtQTD;==jlz|mSto#=!{;Y$voi=~?wU!(!3+1Kw4mq0saxzn##7P+aRc+?D; z7gh?6fur{VRjycM{%HBZPh|eeACq`Hwu4EuU>zK?;CXnO|Hd`JY}}!??Ut>MJu*E?C~j!vwl&Z+-h1u^3)QC zdo_Kvv1xjhVJUw@y-gt{qv*@7=RUu#$*p+5AIh8eOrYrab~`W?rh~H5V+7%y{aH90 z*f1LsuxhTecMlHw2$h?~2M|;CXG0Qm@Gy3aQrJqa?8S#BcKC?)#CRFZHBBe65KlSG zsrX7)MI|5Rdm)a)eoq!t;5KQ!;o1w8wTu&(bsmE9uIKa?*!Dq#Duq6$oPVgm@|cNm|oqLt|z60@hx2J zqqO%+vJ2WZMowE6QH4xI{Ltrz>TGmG%>~4AiA)vaBA#Y3o>K(lP?U*a;$ zeL4C-^2rUl};)F#1a1IF*bnefVm##rV?V%l5YSaWehU^(t}(i zK@Tli3g9S0LFIJF|hII3Dc;@D0k@%9lJz4?fQ8w<9Y26**Td_Q~V{l0&FHIt#yGXyGSAG;=C^N_lethb^uJSF#GC@SGt?d(t)l}8Fbxq8Wa@S z3tmXE+G-V&$WuZ~%^329l=NOdv`MgxOcrgR7&vn8eoeKkpGiM2UqeNg1%Qz(YA#%VVnjW9QF=s`XAe z^@3?^g@^2lMAaU5_NWSHkPv>m_(q4b02jmD&)mY=R23}SQC02@=2f-Uvo5GN-yczF7Us}C1`Cb)8s%Z{o31LFD6h?BUcSIUlAFX zWSnCO`7CdPa=g-K++^euz)Vp8F?B`wS0oAmB zh@7zf(fsyJuXynDDIcxfw;`HErSTDqj@S=1X zzug^6uCT5S!2#opM|Z%{{a~n<1ANAddte2Mhv(d6Nb-H$FuI76;Ns5PyQUS=G6)N zQIA2!@%9+xE=gVDtXI2mqxpW~Q-&8w4Mx|rx;C8VvvtJ6?A#$1Wv~Q{S*pFi&-u$}!qRL(O}M_pfv*>O zBuL+i`;xN0Kq(1xfGJe$OOzGLKXsqewOw{kgxZaXt+A_qx%?tR0uVQc5i9-iV{97A zN(V5rz0;%WKX|hF7uaaH7Iu01h)yCd1+an5KLAHR|1bj(jwYAa#WfJxb9w1xK?Pm41dn(R;W-m@HqY-ra zg6+pNdP7MlI={r{wSj6vF(Mel(N@(Ug)s$-$A<#=Srt!Kj?$|s1X?HZ&XC|%=V-tG z%Yd)B;de*BUmj|V`NzAZ9Bn#~PpmJQrw6o}?C(EJtER)aSM+utDG#lMbm!uZw=D0N9zWRYv({{bkV@_gH=nPrM^0DnFyr4g%ya2i zl`b`B4s|^{7ZS-UysqpiEkvWIXwp5kXQ-#Cx>RlIS*i|I&sB@c)`u55q@D&R{a995huVN8oY} z&RXC{Y>@