/
sis.R
147 lines (139 loc) · 4.72 KB
/
sis.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
#' Standard importance sampling (SIS)
#'
#' Implementation of standard importance sampling (SIS).
#'
#' @param log_ratios An array, matrix, or vector of importance ratios on the log
#' scale (for Importance sampling LOO, these are *negative* log-likelihood
#' values). See the **Methods (by class)** section below for a detailed
#' description of how to specify the inputs for each method.
#' @template cores
#' @param ... Arguments passed on to the various methods.
#' @param r_eff Vector of relative effective sample size estimates containing
#' one element per observation. The values provided should be the relative
#' effective sample sizes of `1/exp(log_ratios)` (i.e., `1/ratios`).
#' This is related to the relative efficiency of estimating the normalizing
#' term in self-normalizing importance sampling. See the [relative_eff()]
#' helper function for computing `r_eff`. If using `psis` with
#' draws of the `log_ratios` not obtained from MCMC then the warning
#' message thrown when not specifying `r_eff` can be disabled by
#' setting `r_eff` to `NA`.
#'
#' @return The `sis()` methods return an object of class `"sis"`,
#' which is a named list with the following components:
#'
#' \describe{
#' \item{`log_weights`}{
#' Vector or matrix of smoothed but *unnormalized* log
#' weights. To get normalized weights use the
#' [`weights()`][weights.importance_sampling] method provided for objects of
#' class `sis`.
#' }
#' \item{`diagnostics`}{
#' A named list containing one vector:
#' * `pareto_k`: Not used in `sis`, all set to 0.
#' * `n_eff`: effective sample size estimates.
#' }
#' }
#'
#' Objects of class `"sis"` also have the following [attributes][attributes()]:
#' \describe{
#' \item{`norm_const_log`}{
#' Vector of precomputed values of `colLogSumExps(log_weights)` that are
#' used internally by the `weights` method to normalize the log weights.
#' }
#' \item{`r_eff`}{
#' If specified, the user's `r_eff` argument.
#' }
#' \item{`tail_len`}{
#' Not used for `sis`.
#' }
#' \item{`dims`}{
#' Integer vector of length 2 containing `S` (posterior sample size)
#' and `N` (number of observations).
#' }
#' \item{`method`}{
#' Method used for importance sampling, here `sis`.
#' }
#' }
#'
#' @seealso
#' * [psis()] for approximate LOO-CV using PSIS.
#' * [loo()] for approximate LOO-CV.
#' * [pareto-k-diagnostic] for PSIS diagnostics.
#'
#' @template loo-and-psis-references
#'
#' @examples
#' log_ratios <- -1 * example_loglik_array()
#' r_eff <- relative_eff(exp(-log_ratios))
#' sis_result <- sis(log_ratios, r_eff = r_eff)
#' str(sis_result)
#'
#' # extract smoothed weights
#' lw <- weights(sis_result) # default args are log=TRUE, normalize=TRUE
#' ulw <- weights(sis_result, normalize=FALSE) # unnormalized log-weights
#'
#' w <- weights(sis_result, log=FALSE) # normalized weights (not log-weights)
#' uw <- weights(sis_result, log=FALSE, normalize = FALSE) # unnormalized weights
#'
#' @export
sis <- function(log_ratios, ...) UseMethod("sis")
#' @export
#' @templateVar fn sis
#' @template array
#'
sis.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 = "sis")
}
#' @export
#' @templateVar fn sis
#' @template matrix
#'
sis.matrix <-
function(log_ratios,
...,
r_eff = NULL,
cores = getOption("mc.cores", 1)) {
importance_sampling.matrix(log_ratios,
...,
r_eff = r_eff,
cores = cores,
method = "sis")
}
#' @export
#' @templateVar fn sis
#' @template vector
#'
sis.default <-
function(log_ratios, ..., r_eff = NULL) {
importance_sampling.default(log_ratios = log_ratios, ...,
r_eff = r_eff,
method = "sis")
}
#' @rdname psis
#' @export
is.sis <- function(x) {
inherits(x, "sis") && is.list(x)
}
# internal ----------------------------------------------------------------
#' Standard IS on a single vector
#'
#' @noRd
#' @param log_ratios_i A vector of log importance ratios (for `loo()`, negative
#' log likelihoods).
#' @param ... Not used. Included to conform to PSIS API.
#'
#' @details Implementation standard importance sampling.
#' @return A named list containing:
#' * `lw`: vector of unnormalized log weights
#' * `pareto_k`: scalar Pareto k estimate. For IS, this defaults to 0.
do_sis_i <- function(log_ratios_i, ...) {
S <- length(log_ratios_i)
list(log_weights = log_ratios_i, pareto_k = 0)
}