From e2d7c7fd03afcde79d8c5f57c2c9cb77059113b0 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Sun, 4 Aug 2024 18:14:08 +0300 Subject: [PATCH 1/2] avoid under and overflows in stacking --- R/loo_model_weights.R | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/R/loo_model_weights.R b/R/loo_model_weights.R index 861eae64..3d9aba05 100644 --- a/R/loo_model_weights.R +++ b/R/loo_model_weights.R @@ -257,15 +257,12 @@ stacking_weights <- stop("At least two models are required for stacking weights.") } - exp_lpd_point <- exp(lpd_point) negative_log_score_loo <- function(w) { # objective function: log score stopifnot(length(w) == K - 1) - w_full <- c(w, 1 - sum(w)) - sum <- 0 - for (i in 1:N) { - sum <- sum + log(exp(lpd_point[i, ]) %*% w_full) - } + w_full <- (c(w, 1 - sum(w))) + # avoid over- and underflows using log weights and rowLogSumExps + sum <- sum(matrixStats::rowLogSumExps(sweep(lpd_point[1:N,], 2, log(w_full), '+'))) return(-as.numeric(sum)) } @@ -274,11 +271,11 @@ stacking_weights <- stopifnot(length(w) == K - 1) w_full <- c(w, 1 - sum(w)) grad <- rep(0, K - 1) + # avoid over- and underflows using log weights, rowLogSumExps, + # and by subtracting the row maximum of lpd_point + mlpd <- matrixStats::rowMaxs(lpd_point) for (k in 1:(K - 1)) { - for (i in 1:N) { - grad[k] <- grad[k] + - (exp_lpd_point[i, k] - exp_lpd_point[i, K]) / (exp_lpd_point[i,] %*% w_full) - } + grad[k] <- sum((exp(lpd_point[, k]-mlpd) - exp(lpd_point[, K]-mlpd)) / exp(matrixStats::rowLogSumExps(sweep(lpd_point, 2, log(w_full), '+'))-mlpd)) } return(-grad) } From 568f29b2ad856fe72c64966c0ad260a5c4c40a73 Mon Sep 17 00:00:00 2001 From: jgabry Date: Mon, 5 Aug 2024 13:17:22 -0600 Subject: [PATCH 2/2] remove extra parentheses --- R/loo_model_weights.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/loo_model_weights.R b/R/loo_model_weights.R index 3d9aba05..da899988 100644 --- a/R/loo_model_weights.R +++ b/R/loo_model_weights.R @@ -260,7 +260,7 @@ stacking_weights <- negative_log_score_loo <- function(w) { # objective function: log score stopifnot(length(w) == K - 1) - w_full <- (c(w, 1 - sum(w))) + w_full <- c(w, 1 - sum(w)) # avoid over- and underflows using log weights and rowLogSumExps sum <- sum(matrixStats::rowLogSumExps(sweep(lpd_point[1:N,], 2, log(w_full), '+'))) return(-as.numeric(sum)) @@ -275,7 +275,7 @@ stacking_weights <- # and by subtracting the row maximum of lpd_point mlpd <- matrixStats::rowMaxs(lpd_point) for (k in 1:(K - 1)) { - grad[k] <- sum((exp(lpd_point[, k]-mlpd) - exp(lpd_point[, K]-mlpd)) / exp(matrixStats::rowLogSumExps(sweep(lpd_point, 2, log(w_full), '+'))-mlpd)) + grad[k] <- sum((exp(lpd_point[, k] - mlpd) - exp(lpd_point[, K] - mlpd)) / exp(matrixStats::rowLogSumExps(sweep(lpd_point, 2, log(w_full), '+')) - mlpd)) } return(-grad) }