From b887a1d872d292fd1593512b7491545222127edc Mon Sep 17 00:00:00 2001 From: Carlos Parada Date: Thu, 26 Aug 2021 19:04:43 -0700 Subject: [PATCH 1/3] Performance improvement Only sort the tail values of the array, rather than the whole thing. --- R/psis.R | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/R/psis.R b/R/psis.R index e82a1ef9..0eec1c0a 100644 --- a/R/psis.R +++ b/R/psis.R @@ -205,8 +205,9 @@ do_psis_i <- function(log_ratios_i, tail_len_i, ...) { khat <- Inf if (enough_tail_samples(tail_len_i)) { - ord <- sort.int(lw_i, index.return = TRUE) - tail_ids <- seq(S - tail_len_i + 1, S) + tail_start <- S - tail_len_i + 1 + ord <- sort.int(lw_i, index.return = TRUE, partial=tail_start:S) + tail_ids <- seq(tail_start, S) lw_tail <- ord$x[tail_ids] if (abs(max(lw_tail) - min(lw_tail)) < .Machine$double.eps/100) { warning( From 6e070e1d8b29b2308081a2a48f22af304a23bc39 Mon Sep 17 00:00:00 2001 From: Closed-Limelike-Curves Date: Mon, 20 Sep 2021 11:57:07 -0500 Subject: [PATCH 2/3] off-by-one --- .lh/R/psis.R.json | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) create mode 100644 .lh/R/psis.R.json diff --git a/.lh/R/psis.R.json b/.lh/R/psis.R.json new file mode 100644 index 00000000..0978361e --- /dev/null +++ b/.lh/R/psis.R.json @@ -0,0 +1,18 @@ +{ + "sourceFile": "R/psis.R", + "activeCommit": 0, + "commits": [ + { + "activePatchIndex": 0, + "patches": [ + { + "date": 1632157006627, + "content": "Index: \n===================================================================\n--- \n+++ \n" + } + ], + "date": 1632157006627, + "name": "Commit-0", + "content": "#' Pareto smoothed importance sampling (PSIS)\n#'\n#' Implementation of Pareto smoothed importance sampling (PSIS), a method for\n#' stabilizing importance ratios. The version of PSIS implemented here\n#' corresponds to the algorithm presented in Vehtari, Simpson, Gelman, Yao,\n#' and Gabry (2019).\n#' For PSIS diagnostics see the [pareto-k-diagnostic] page.\n#'\n#' @export\n#' @param log_ratios An array, matrix, or vector of importance ratios on the log\n#' scale (for PSIS-LOO these are *negative* log-likelihood values). See the\n#' **Methods (by class)** section below for a detailed description of how\n#' to specify the inputs for each method.\n#' @param ... Arguments passed on to the various methods.\n#' @template cores\n#' @param r_eff Vector of relative effective sample size estimates containing\n#' one element per observation. The values provided should be the relative\n#' effective sample sizes of `1/exp(log_ratios)` (i.e., `1/ratios`).\n#' This is related to the relative efficiency of estimating the normalizing\n#' term in self-normalizing importance sampling. If `r_eff` is not\n#' provided then the reported PSIS effective sample sizes and Monte Carlo\n#' error estimates will be over-optimistic. See the [relative_eff()]\n#' helper function for computing `r_eff`. If using `psis` with\n#' draws of the `log_ratios` not obtained from MCMC then the warning\n#' message thrown when not specifying `r_eff` can be disabled by\n#' setting `r_eff` to `NA`.\n#'\n#' @return The `psis()` methods return an object of class `\"psis\"`,\n#' which is a named list with the following components:\n#'\n#' \\describe{\n#' \\item{`log_weights`}{\n#' Vector or matrix of smoothed (and truncated) but *unnormalized* log\n#' weights. To get normalized weights use the\n#' [`weights()`][weights.importance_sampling] method provided for objects of\n#' class `\"psis\"`.\n#' }\n#' \\item{`diagnostics`}{\n#' A named list containing two vectors:\n#' * `pareto_k`: Estimates of the shape parameter \\eqn{k} of the\n#' generalized Pareto distribution. See the [pareto-k-diagnostic]\n#' page for details.\n#' * `n_eff`: PSIS effective sample size estimates.\n#' }\n#' }\n#'\n#' Objects of class `\"psis\"` also have the following [attributes][attributes()]:\n#' \\describe{\n#' \\item{`norm_const_log`}{\n#' Vector of precomputed values of `colLogSumExps(log_weights)` that are\n#' used internally by the `weights` method to normalize the log weights.\n#' }\n#' \\item{`tail_len`}{\n#' Vector of tail lengths used for fitting the generalized Pareto distribution.\n#' }\n#' \\item{`r_eff`}{\n#' If specified, the user's `r_eff` argument.\n#' }\n#' \\item{`dims`}{\n#' Integer vector of length 2 containing `S` (posterior sample size)\n#' and `N` (number of observations).\n#' }\n#' \\item{`method`}{\n#' Method used for importance sampling, here `psis`.\n#' }\n#' }\n#'\n#' @seealso\n#' * [loo()] for approximate LOO-CV using PSIS.\n#' * [pareto-k-diagnostic] for PSIS diagnostics.\n#'\n#' @template loo-and-psis-references\n#'\n#' @examples\n#' log_ratios <- -1 * example_loglik_array()\n#' r_eff <- relative_eff(exp(-log_ratios))\n#' psis_result <- psis(log_ratios, r_eff = r_eff)\n#' str(psis_result)\n#' plot(psis_result)\n#'\n#' # extract smoothed weights\n#' lw <- weights(psis_result) # default args are log=TRUE, normalize=TRUE\n#' ulw <- weights(psis_result, normalize=FALSE) # unnormalized log-weights\n#'\n#' w <- weights(psis_result, log=FALSE) # normalized weights (not log-weights)\n#' uw <- weights(psis_result, log=FALSE, normalize = FALSE) # unnormalized weights\n#'\n#'\n#'\npsis <- function(log_ratios, ...) UseMethod(\"psis\")\n\n#' @export\n#' @templateVar fn psis\n#' @template array\n#'\npsis.array <-\n function(log_ratios, ...,\n r_eff = NULL,\n cores = getOption(\"mc.cores\", 1)) {\n importance_sampling.array(log_ratios = log_ratios, ...,\n r_eff = r_eff,\n cores = cores,\n method = \"psis\")\n }\n\n\n#' @export\n#' @templateVar fn psis\n#' @template matrix\n#'\npsis.matrix <-\n function(log_ratios,\n ...,\n r_eff = NULL,\n cores = getOption(\"mc.cores\", 1)) {\n importance_sampling.matrix(log_ratios,\n ...,\n r_eff = r_eff,\n cores = cores,\n method = \"psis\")\n }\n\n#' @export\n#' @templateVar fn psis\n#' @template vector\n#'\npsis.default <-\n function(log_ratios, ..., r_eff = NULL) {\n importance_sampling.default(log_ratios = log_ratios, ...,\n r_eff = r_eff,\n method = \"psis\")\n }\n\n\n#' @rdname psis\n#' @export\n#' @param x For `is.psis()`, an object to check.\nis.psis <- function(x) {\n inherits(x, \"psis\") && is.list(x)\n}\n\n\n# internal ----------------------------------------------------------------\n\n#' @noRd\n#' @seealso importance_sampling_object\npsis_object <-\n function(unnormalized_log_weights,\n pareto_k,\n tail_len,\n r_eff) {\n importance_sampling_object(unnormalized_log_weights = unnormalized_log_weights,\n pareto_k = pareto_k,\n tail_len = tail_len,\n r_eff = r_eff,\n method = \"psis\")\n }\n\n\n#' @noRd\n#' @seealso do_importance_sampling\ndo_psis <- function(log_ratios, r_eff, cores, method){\n do_importance_sampling(log_ratios = log_ratios,\n r_eff = r_eff,\n cores = cores,\n method = \"psis\")\n}\n\n#' Extract named components from each list in the list of lists obtained by\n#' parallelizing `do_psis_i()`\n#'\n#' @noRd\n#' @param x List of lists.\n#' @param item String naming the component or attribute to pull out of each list\n#' (or list-like object).\n#' @param fun,fun.val passed to `vapply()`'s `FUN` and `FUN.VALUE` arguments.\n#' @return Numeric vector or matrix.\n#'\npsis_apply <- function(x, item, fun = c(\"[[\", \"attr\"), fun_val = numeric(1)) {\n if (!is.list(x)) stop(\"Internal error ('x' must be a list for psis_apply)\")\n vapply(x, FUN = match.arg(fun), FUN.VALUE = fun_val, item)\n}\n\n#' PSIS on a single vector\n#'\n#' @noRd\n#' @param log_ratios_i A vector of log importance ratios (for `loo()`, negative\n#' log likelihoods).\n#' @param tail_len_i An integer tail length.\n#' @param ... Not used. Included to conform to API for differen IS methods.\n#'\n#' @details\n#' * If there are enough tail samples then the tail is smoothed with PSIS\n#' * The log weights (or log ratios if no smoothing) larger than the largest raw\n#' ratio are set to the largest raw ratio\n#'\n#' @return A named list containing:\n#' * `lw`: vector of unnormalized log weights\n#' * `pareto_k`: scalar Pareto k estimate.\n#'\ndo_psis_i <- function(log_ratios_i, tail_len_i, ...) {\n S <- length(log_ratios_i)\n # shift log ratios for safer exponentation\n lw_i <- log_ratios_i - max(log_ratios_i)\n khat <- Inf\n\n if (enough_tail_samples(tail_len_i)) {\n tail_start <- S - tail_len_i + 1\n ord <- sort.int(lw_i, index.return = TRUE, partial = (tail_start - 1):S)\n tail_ids <- seq(tail_start, S)\n lw_tail <- ord$x[tail_ids]\n if (abs(max(lw_tail) - min(lw_tail)) < .Machine$double.eps/100) {\n warning(\n \"Can't fit generalized Pareto distribution \",\n \"because all tail values are the same.\",\n call. = FALSE\n )\n } else {\n cutoff <- ord$x[tail_start - 1] # largest value smaller than tail values\n smoothed <- psis_smooth_tail(lw_tail, cutoff)\n khat <- smoothed$k\n lw_i[ord$ix[tail_ids]] <- smoothed$tail\n }\n }\n\n # truncate at max of raw wts (i.e., 0 since max has been subtracted)\n lw_i[lw_i > 0] <- 0\n # shift log weights back so that the smallest log weights remain unchanged\n lw_i <- lw_i + max(log_ratios_i)\n\n list(log_weights = lw_i, pareto_k = khat)\n}\n\n#' PSIS tail smoothing for a single vector\n#'\n#' @noRd\n#' @param x Vector of tail elements already sorted in ascending order.\n#' @return A named list containing:\n#' * `tail`: vector same size as `x` containing the logs of the\n#' order statistics of the generalized pareto distribution.\n#' * `k`: scalar shape parameter estimate.\n#'\npsis_smooth_tail <- function(x, cutoff) {\n len <- length(x)\n exp_cutoff <- exp(cutoff)\n\n # save time not sorting since x already sorted\n fit <- gpdfit(exp(x) - exp_cutoff, sort_x = FALSE)\n k <- fit$k\n sigma <- fit$sigma\n if (is.finite(k)) {\n p <- (seq_len(len) - 0.5) / len\n qq <- qgpd(p, k, sigma) + exp_cutoff\n tail <- log(qq)\n } else {\n tail <- x\n }\n list(tail = tail, k = k)\n}\n\n\n#' Calculate tail lengths to use for fitting the GPD\n#'\n#' The number of weights (i.e., tail length) used to fit the generalized Pareto\n#' distribution is now decreasing with the number of posterior draws S, and is\n#' also adjusted based on the relative MCMC neff for `exp(log_lik)`. This will\n#' answer the questions about the asymptotic properties, works better for thick\n#' tailed proposal distributions, and is adjusted based on dependent Markov chain\n#' samples. Specifically, the tail length is now `3*sqrt(S)/r_eff` but capped at\n#' 20% of the total number of weights.\n#'\n#' @noRd\n#' @param r_eff A N-vector of relative MCMC effective sample sizes of `exp(log-lik matrix)`.\n#' @param S The (integer) size of posterior sample.\n#' @return An N-vector of tail lengths.\n#'\nn_pareto <- function(r_eff, S) {\n ceiling(pmin(0.2 * S, 3 * sqrt(S / r_eff)))\n}\n\n#' Check for enough tail samples to fit GPD\n#'\n#' @noRd\n#' @param tail_len Integer tail length.\n#' @param min_len The minimum allowed tail length.\n#' @return `TRUE` or `FALSE`\n#'\nenough_tail_samples <- function(tail_len, min_len = 5) {\n tail_len >= min_len\n}\n\n\n#' Throw warnings about pareto k estimates\n#'\n#' @noRd\n#' @param k A vector of pareto k estimates.\n#' @param high The value at which to warn about slighly high estimates.\n#' @param too_high The value at which to warn about very high estimates.\n#' @return Nothing, just possibly throws warnings.\n#'\nthrow_pareto_warnings <- function(k, high = 0.5, too_high = 0.7) {\n if (any(k > too_high)) {\n .warn(\"Some Pareto k diagnostic values are too high. \", .k_help())\n } else if (any(k > high)) {\n .warn(\"Some Pareto k diagnostic values are slightly high. \", .k_help())\n }\n}\n\n#' Warn if not enough tail samples to fit GPD\n#'\n#' @noRd\n#' @param tail_lengths Vector of tail lengths.\n#' @return `tail_lengths`, invisibly.\n#'\nthrow_tail_length_warnings <- function(tail_lengths) {\n tail_len_bad <- !sapply(tail_lengths, enough_tail_samples)\n if (any(tail_len_bad)) {\n if (length(tail_lengths) == 1) {\n warning(\n \"Not enough tail samples to fit the generalized Pareto distribution.\",\n call. = FALSE, immediate. = TRUE\n )\n } else {\n bad <- which(tail_len_bad)\n Nbad <- length(bad)\n warning(\n \"Not enough tail samples to fit the generalized Pareto distribution \",\n \"in some or all columns of matrix of log importance ratios. \",\n \"Skipping the following columns: \",\n paste(if (Nbad <= 10) bad else bad[1:10], collapse = \", \"),\n if (Nbad > 10) paste0(\", ... [\", Nbad - 10, \" more not printed].\\n\") else \"\\n\",\n call. = FALSE,\n immediate. = TRUE\n )\n }\n }\n invisible(tail_lengths)\n}\n\n#' Prepare `r_eff` to pass to `psis()` and throw warnings/errors if necessary\n#'\n#' @noRd\n#' @param r_eff User's `r_eff` argument.\n#' @param len The length `r_eff` should have if not `NULL` or `NA`.\n#' @return\n#' * If `r_eff` has length `len` then `r_eff` is returned.\n#' * If `r_eff` is `NULL` then a warning is thrown and `rep(1, len)` is returned.\n#' * If `r_eff` is `NA` then the warning is skipped and\n#' `rep(1, len)` is returned.\n#' * If `r_eff` has length `len` but has `NA`s then an error is thrown.\n#'\nprepare_psis_r_eff <- function(r_eff, len) {\n if (isTRUE(is.null(r_eff) || all(is.na(r_eff)))) {\n if (!called_from_loo() && is.null(r_eff)) {\n throw_psis_r_eff_warning()\n }\n r_eff <- rep(1, len)\n } else if (length(r_eff) != len) {\n stop(\"'r_eff' must have one value per observation.\", call. = FALSE)\n } else if (anyNA(r_eff)) {\n stop(\"Can't mix NA and not NA values in 'r_eff'.\", call. = FALSE)\n }\n return(r_eff)\n}\n\n#' Check if `psis()` was called from one of the loo methods\n#'\n#' @noRd\n#' @return `TRUE` if the `loo()` array, matrix, or function method is found in\n#' the active call list, `FALSE` otherwise.\n#'\ncalled_from_loo <- function() {\n calls <- sys.calls()\n txt <- unlist(lapply(calls, deparse))\n patts <- \"loo.array\\\\(|loo.matrix\\\\(|loo.function\\\\(\"\n check <- sapply(txt, function(x) grepl(patts, x))\n isTRUE(any(check))\n}\n\n#' Warning message about missing `r_eff` argument\n#' @noRd\nthrow_psis_r_eff_warning <- function() {\n warning(\n \"Relative effective sample sizes ('r_eff' argument) not specified. \",\n \"PSIS n_eff will not be adjusted based on MCMC n_eff.\",\n call. = FALSE\n )\n}\n\n" + } + ] +} \ No newline at end of file From 788c88981b08d5e6468026bb5c6735c1930d522e Mon Sep 17 00:00:00 2001 From: Closed-Limelike-Curves Date: Mon, 20 Sep 2021 16:07:51 -0500 Subject: [PATCH 3/3] Whoops --- .lh/R/psis.R.json | 18 ------------------ R/psis.R | 4 ++-- 2 files changed, 2 insertions(+), 20 deletions(-) delete mode 100644 .lh/R/psis.R.json diff --git a/.lh/R/psis.R.json b/.lh/R/psis.R.json deleted file mode 100644 index 0978361e..00000000 --- a/.lh/R/psis.R.json +++ /dev/null @@ -1,18 +0,0 @@ -{ - "sourceFile": "R/psis.R", - "activeCommit": 0, - "commits": [ - { - "activePatchIndex": 0, - "patches": [ - { - "date": 1632157006627, - "content": "Index: \n===================================================================\n--- \n+++ \n" - } - ], - "date": 1632157006627, - "name": "Commit-0", - "content": "#' Pareto smoothed importance sampling (PSIS)\n#'\n#' Implementation of Pareto smoothed importance sampling (PSIS), a method for\n#' stabilizing importance ratios. The version of PSIS implemented here\n#' corresponds to the algorithm presented in Vehtari, Simpson, Gelman, Yao,\n#' and Gabry (2019).\n#' For PSIS diagnostics see the [pareto-k-diagnostic] page.\n#'\n#' @export\n#' @param log_ratios An array, matrix, or vector of importance ratios on the log\n#' scale (for PSIS-LOO these are *negative* log-likelihood values). See the\n#' **Methods (by class)** section below for a detailed description of how\n#' to specify the inputs for each method.\n#' @param ... Arguments passed on to the various methods.\n#' @template cores\n#' @param r_eff Vector of relative effective sample size estimates containing\n#' one element per observation. The values provided should be the relative\n#' effective sample sizes of `1/exp(log_ratios)` (i.e., `1/ratios`).\n#' This is related to the relative efficiency of estimating the normalizing\n#' term in self-normalizing importance sampling. If `r_eff` is not\n#' provided then the reported PSIS effective sample sizes and Monte Carlo\n#' error estimates will be over-optimistic. See the [relative_eff()]\n#' helper function for computing `r_eff`. If using `psis` with\n#' draws of the `log_ratios` not obtained from MCMC then the warning\n#' message thrown when not specifying `r_eff` can be disabled by\n#' setting `r_eff` to `NA`.\n#'\n#' @return The `psis()` methods return an object of class `\"psis\"`,\n#' which is a named list with the following components:\n#'\n#' \\describe{\n#' \\item{`log_weights`}{\n#' Vector or matrix of smoothed (and truncated) but *unnormalized* log\n#' weights. To get normalized weights use the\n#' [`weights()`][weights.importance_sampling] method provided for objects of\n#' class `\"psis\"`.\n#' }\n#' \\item{`diagnostics`}{\n#' A named list containing two vectors:\n#' * `pareto_k`: Estimates of the shape parameter \\eqn{k} of the\n#' generalized Pareto distribution. See the [pareto-k-diagnostic]\n#' page for details.\n#' * `n_eff`: PSIS effective sample size estimates.\n#' }\n#' }\n#'\n#' Objects of class `\"psis\"` also have the following [attributes][attributes()]:\n#' \\describe{\n#' \\item{`norm_const_log`}{\n#' Vector of precomputed values of `colLogSumExps(log_weights)` that are\n#' used internally by the `weights` method to normalize the log weights.\n#' }\n#' \\item{`tail_len`}{\n#' Vector of tail lengths used for fitting the generalized Pareto distribution.\n#' }\n#' \\item{`r_eff`}{\n#' If specified, the user's `r_eff` argument.\n#' }\n#' \\item{`dims`}{\n#' Integer vector of length 2 containing `S` (posterior sample size)\n#' and `N` (number of observations).\n#' }\n#' \\item{`method`}{\n#' Method used for importance sampling, here `psis`.\n#' }\n#' }\n#'\n#' @seealso\n#' * [loo()] for approximate LOO-CV using PSIS.\n#' * [pareto-k-diagnostic] for PSIS diagnostics.\n#'\n#' @template loo-and-psis-references\n#'\n#' @examples\n#' log_ratios <- -1 * example_loglik_array()\n#' r_eff <- relative_eff(exp(-log_ratios))\n#' psis_result <- psis(log_ratios, r_eff = r_eff)\n#' str(psis_result)\n#' plot(psis_result)\n#'\n#' # extract smoothed weights\n#' lw <- weights(psis_result) # default args are log=TRUE, normalize=TRUE\n#' ulw <- weights(psis_result, normalize=FALSE) # unnormalized log-weights\n#'\n#' w <- weights(psis_result, log=FALSE) # normalized weights (not log-weights)\n#' uw <- weights(psis_result, log=FALSE, normalize = FALSE) # unnormalized weights\n#'\n#'\n#'\npsis <- function(log_ratios, ...) UseMethod(\"psis\")\n\n#' @export\n#' @templateVar fn psis\n#' @template array\n#'\npsis.array <-\n function(log_ratios, ...,\n r_eff = NULL,\n cores = getOption(\"mc.cores\", 1)) {\n importance_sampling.array(log_ratios = log_ratios, ...,\n r_eff = r_eff,\n cores = cores,\n method = \"psis\")\n }\n\n\n#' @export\n#' @templateVar fn psis\n#' @template matrix\n#'\npsis.matrix <-\n function(log_ratios,\n ...,\n r_eff = NULL,\n cores = getOption(\"mc.cores\", 1)) {\n importance_sampling.matrix(log_ratios,\n ...,\n r_eff = r_eff,\n cores = cores,\n method = \"psis\")\n }\n\n#' @export\n#' @templateVar fn psis\n#' @template vector\n#'\npsis.default <-\n function(log_ratios, ..., r_eff = NULL) {\n importance_sampling.default(log_ratios = log_ratios, ...,\n r_eff = r_eff,\n method = \"psis\")\n }\n\n\n#' @rdname psis\n#' @export\n#' @param x For `is.psis()`, an object to check.\nis.psis <- function(x) {\n inherits(x, \"psis\") && is.list(x)\n}\n\n\n# internal ----------------------------------------------------------------\n\n#' @noRd\n#' @seealso importance_sampling_object\npsis_object <-\n function(unnormalized_log_weights,\n pareto_k,\n tail_len,\n r_eff) {\n importance_sampling_object(unnormalized_log_weights = unnormalized_log_weights,\n pareto_k = pareto_k,\n tail_len = tail_len,\n r_eff = r_eff,\n method = \"psis\")\n }\n\n\n#' @noRd\n#' @seealso do_importance_sampling\ndo_psis <- function(log_ratios, r_eff, cores, method){\n do_importance_sampling(log_ratios = log_ratios,\n r_eff = r_eff,\n cores = cores,\n method = \"psis\")\n}\n\n#' Extract named components from each list in the list of lists obtained by\n#' parallelizing `do_psis_i()`\n#'\n#' @noRd\n#' @param x List of lists.\n#' @param item String naming the component or attribute to pull out of each list\n#' (or list-like object).\n#' @param fun,fun.val passed to `vapply()`'s `FUN` and `FUN.VALUE` arguments.\n#' @return Numeric vector or matrix.\n#'\npsis_apply <- function(x, item, fun = c(\"[[\", \"attr\"), fun_val = numeric(1)) {\n if (!is.list(x)) stop(\"Internal error ('x' must be a list for psis_apply)\")\n vapply(x, FUN = match.arg(fun), FUN.VALUE = fun_val, item)\n}\n\n#' PSIS on a single vector\n#'\n#' @noRd\n#' @param log_ratios_i A vector of log importance ratios (for `loo()`, negative\n#' log likelihoods).\n#' @param tail_len_i An integer tail length.\n#' @param ... Not used. Included to conform to API for differen IS methods.\n#'\n#' @details\n#' * If there are enough tail samples then the tail is smoothed with PSIS\n#' * The log weights (or log ratios if no smoothing) larger than the largest raw\n#' ratio are set to the largest raw ratio\n#'\n#' @return A named list containing:\n#' * `lw`: vector of unnormalized log weights\n#' * `pareto_k`: scalar Pareto k estimate.\n#'\ndo_psis_i <- function(log_ratios_i, tail_len_i, ...) {\n S <- length(log_ratios_i)\n # shift log ratios for safer exponentation\n lw_i <- log_ratios_i - max(log_ratios_i)\n khat <- Inf\n\n if (enough_tail_samples(tail_len_i)) {\n tail_start <- S - tail_len_i + 1\n ord <- sort.int(lw_i, index.return = TRUE, partial = (tail_start - 1):S)\n tail_ids <- seq(tail_start, S)\n lw_tail <- ord$x[tail_ids]\n if (abs(max(lw_tail) - min(lw_tail)) < .Machine$double.eps/100) {\n warning(\n \"Can't fit generalized Pareto distribution \",\n \"because all tail values are the same.\",\n call. = FALSE\n )\n } else {\n cutoff <- ord$x[tail_start - 1] # largest value smaller than tail values\n smoothed <- psis_smooth_tail(lw_tail, cutoff)\n khat <- smoothed$k\n lw_i[ord$ix[tail_ids]] <- smoothed$tail\n }\n }\n\n # truncate at max of raw wts (i.e., 0 since max has been subtracted)\n lw_i[lw_i > 0] <- 0\n # shift log weights back so that the smallest log weights remain unchanged\n lw_i <- lw_i + max(log_ratios_i)\n\n list(log_weights = lw_i, pareto_k = khat)\n}\n\n#' PSIS tail smoothing for a single vector\n#'\n#' @noRd\n#' @param x Vector of tail elements already sorted in ascending order.\n#' @return A named list containing:\n#' * `tail`: vector same size as `x` containing the logs of the\n#' order statistics of the generalized pareto distribution.\n#' * `k`: scalar shape parameter estimate.\n#'\npsis_smooth_tail <- function(x, cutoff) {\n len <- length(x)\n exp_cutoff <- exp(cutoff)\n\n # save time not sorting since x already sorted\n fit <- gpdfit(exp(x) - exp_cutoff, sort_x = FALSE)\n k <- fit$k\n sigma <- fit$sigma\n if (is.finite(k)) {\n p <- (seq_len(len) - 0.5) / len\n qq <- qgpd(p, k, sigma) + exp_cutoff\n tail <- log(qq)\n } else {\n tail <- x\n }\n list(tail = tail, k = k)\n}\n\n\n#' Calculate tail lengths to use for fitting the GPD\n#'\n#' The number of weights (i.e., tail length) used to fit the generalized Pareto\n#' distribution is now decreasing with the number of posterior draws S, and is\n#' also adjusted based on the relative MCMC neff for `exp(log_lik)`. This will\n#' answer the questions about the asymptotic properties, works better for thick\n#' tailed proposal distributions, and is adjusted based on dependent Markov chain\n#' samples. Specifically, the tail length is now `3*sqrt(S)/r_eff` but capped at\n#' 20% of the total number of weights.\n#'\n#' @noRd\n#' @param r_eff A N-vector of relative MCMC effective sample sizes of `exp(log-lik matrix)`.\n#' @param S The (integer) size of posterior sample.\n#' @return An N-vector of tail lengths.\n#'\nn_pareto <- function(r_eff, S) {\n ceiling(pmin(0.2 * S, 3 * sqrt(S / r_eff)))\n}\n\n#' Check for enough tail samples to fit GPD\n#'\n#' @noRd\n#' @param tail_len Integer tail length.\n#' @param min_len The minimum allowed tail length.\n#' @return `TRUE` or `FALSE`\n#'\nenough_tail_samples <- function(tail_len, min_len = 5) {\n tail_len >= min_len\n}\n\n\n#' Throw warnings about pareto k estimates\n#'\n#' @noRd\n#' @param k A vector of pareto k estimates.\n#' @param high The value at which to warn about slighly high estimates.\n#' @param too_high The value at which to warn about very high estimates.\n#' @return Nothing, just possibly throws warnings.\n#'\nthrow_pareto_warnings <- function(k, high = 0.5, too_high = 0.7) {\n if (any(k > too_high)) {\n .warn(\"Some Pareto k diagnostic values are too high. \", .k_help())\n } else if (any(k > high)) {\n .warn(\"Some Pareto k diagnostic values are slightly high. \", .k_help())\n }\n}\n\n#' Warn if not enough tail samples to fit GPD\n#'\n#' @noRd\n#' @param tail_lengths Vector of tail lengths.\n#' @return `tail_lengths`, invisibly.\n#'\nthrow_tail_length_warnings <- function(tail_lengths) {\n tail_len_bad <- !sapply(tail_lengths, enough_tail_samples)\n if (any(tail_len_bad)) {\n if (length(tail_lengths) == 1) {\n warning(\n \"Not enough tail samples to fit the generalized Pareto distribution.\",\n call. = FALSE, immediate. = TRUE\n )\n } else {\n bad <- which(tail_len_bad)\n Nbad <- length(bad)\n warning(\n \"Not enough tail samples to fit the generalized Pareto distribution \",\n \"in some or all columns of matrix of log importance ratios. \",\n \"Skipping the following columns: \",\n paste(if (Nbad <= 10) bad else bad[1:10], collapse = \", \"),\n if (Nbad > 10) paste0(\", ... [\", Nbad - 10, \" more not printed].\\n\") else \"\\n\",\n call. = FALSE,\n immediate. = TRUE\n )\n }\n }\n invisible(tail_lengths)\n}\n\n#' Prepare `r_eff` to pass to `psis()` and throw warnings/errors if necessary\n#'\n#' @noRd\n#' @param r_eff User's `r_eff` argument.\n#' @param len The length `r_eff` should have if not `NULL` or `NA`.\n#' @return\n#' * If `r_eff` has length `len` then `r_eff` is returned.\n#' * If `r_eff` is `NULL` then a warning is thrown and `rep(1, len)` is returned.\n#' * If `r_eff` is `NA` then the warning is skipped and\n#' `rep(1, len)` is returned.\n#' * If `r_eff` has length `len` but has `NA`s then an error is thrown.\n#'\nprepare_psis_r_eff <- function(r_eff, len) {\n if (isTRUE(is.null(r_eff) || all(is.na(r_eff)))) {\n if (!called_from_loo() && is.null(r_eff)) {\n throw_psis_r_eff_warning()\n }\n r_eff <- rep(1, len)\n } else if (length(r_eff) != len) {\n stop(\"'r_eff' must have one value per observation.\", call. = FALSE)\n } else if (anyNA(r_eff)) {\n stop(\"Can't mix NA and not NA values in 'r_eff'.\", call. = FALSE)\n }\n return(r_eff)\n}\n\n#' Check if `psis()` was called from one of the loo methods\n#'\n#' @noRd\n#' @return `TRUE` if the `loo()` array, matrix, or function method is found in\n#' the active call list, `FALSE` otherwise.\n#'\ncalled_from_loo <- function() {\n calls <- sys.calls()\n txt <- unlist(lapply(calls, deparse))\n patts <- \"loo.array\\\\(|loo.matrix\\\\(|loo.function\\\\(\"\n check <- sapply(txt, function(x) grepl(patts, x))\n isTRUE(any(check))\n}\n\n#' Warning message about missing `r_eff` argument\n#' @noRd\nthrow_psis_r_eff_warning <- function() {\n warning(\n \"Relative effective sample sizes ('r_eff' argument) not specified. \",\n \"PSIS n_eff will not be adjusted based on MCMC n_eff.\",\n call. = FALSE\n )\n}\n\n" - } - ] -} \ No newline at end of file diff --git a/R/psis.R b/R/psis.R index 0eec1c0a..c59f50a2 100644 --- a/R/psis.R +++ b/R/psis.R @@ -206,7 +206,7 @@ do_psis_i <- function(log_ratios_i, tail_len_i, ...) { if (enough_tail_samples(tail_len_i)) { tail_start <- S - tail_len_i + 1 - ord <- sort.int(lw_i, index.return = TRUE, partial=tail_start:S) + ord <- sort.int(lw_i, index.return = TRUE, partial=(tail_start - 1):S) tail_ids <- seq(tail_start, S) lw_tail <- ord$x[tail_ids] if (abs(max(lw_tail) - min(lw_tail)) < .Machine$double.eps/100) { @@ -216,7 +216,7 @@ do_psis_i <- function(log_ratios_i, tail_len_i, ...) { call. = FALSE ) } else { - cutoff <- ord$x[min(tail_ids) - 1] # largest value smaller than tail values + cutoff <- ord$x[tail_start - 1] # largest value smaller than tail values smoothed <- psis_smooth_tail(lw_tail, cutoff) khat <- smoothed$k lw_i[ord$ix[tail_ids]] <- smoothed$tail