From a0602adcc8cff061244dfa2eba3a935dea6114ee Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Paul-Christian=20B=C3=BCrkner?= Date: Thu, 5 Dec 2024 17:54:21 +0100 Subject: [PATCH] fix issue #1716 --- DESCRIPTION | 4 ++-- R/stan-likelihood.R | 32 +++++++++----------------------- tests/testthat/tests.stancode.R | 25 ++++++++++++++----------- 3 files changed, 25 insertions(+), 36 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 5e87587d..5e7b2bc7 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,8 +2,8 @@ Package: brms Encoding: UTF-8 Type: Package Title: Bayesian Regression Models using 'Stan' -Version: 2.22.7 -Date: 2024-11-26 +Version: 2.22.8 +Date: 2024-12-05 Authors@R: c(person("Paul-Christian", "Bürkner", email = "paul.buerkner@gmail.com", role = c("aut", "cre")), diff --git a/R/stan-likelihood.R b/R/stan-likelihood.R index 7884b27a..f4555b64 100644 --- a/R/stan-likelihood.R +++ b/R/stan-likelihood.R @@ -354,13 +354,8 @@ stan_log_lik_add_se <- function(sigma, bterms, reqn = stan_log_lik_adj(bterms), # multiply 'dpar' by the 'rate' denominator within the Stan likelihood # @param log add the rate denominator on the log scale if sensible? -# @param req_dot_multiply Censoring may turn non-vectorized into vectorized -# statements later on (see stan_log_lik_cens) which then makes the * operator -# invalid and requires .* instead. Accordingly, req_dot_multiply should be -# FALSE if [n] is required only because of censoring. stan_log_lik_multiply_rate_denom <- function( dpar, bterms, reqn = stan_log_lik_adj(bterms), - req_dot_multiply = stan_log_lik_adj(bterms, c("trunc", "weights")), log = FALSE, transform = NULL, threads = NULL, ...) { dpar_transform <- dpar @@ -381,7 +376,7 @@ stan_log_lik_multiply_rate_denom <- function( # dpar without resp name or index dpar_clean <- sub("(_|\\[).*", "", dpar) is_pred <- dpar_clean %in% c("mu", names(bterms$dpars)) - operator <- str_if(req_dot_multiply || !is_pred, "*", ".*") + operator <- ".*" } glue("{dpar_transform} {operator} {denom}") } @@ -631,13 +626,11 @@ stan_log_lik_binomial <- function(bterms, ...) { stan_log_lik_beta_binomial <- function(bterms, ...) { p <- stan_log_lik_dpars(bterms) p$trials <- stan_log_lik_advars(bterms, "trials", ...)$trials - req_dot_multiply <- !stan_log_lik_adj(bterms) && is_pred_dpar(bterms, "phi") - multiply <- str_if(req_dot_multiply, " .* ", " * ") sdist( "beta_binomial", p$trials, - paste0(p$mu, multiply, p$phi), - paste0("(1 - ", p$mu, ")", multiply, p$phi) + paste0(p$mu, " .* ", p$phi), + paste0("(1 - ", p$mu, ") .* ", p$phi) ) } @@ -665,11 +658,10 @@ stan_log_lik_com_poisson <- function(bterms, ...) { } stan_log_lik_gamma <- function(bterms, ...) { - reqn <- stan_log_lik_adj(bterms) || is_pred_dpar(bterms, "shape") + reqn <- stan_log_lik_adj(bterms) p <- stan_log_lik_dpars(bterms, reqn = reqn) # Stan uses shape-rate parameterization with rate = shape / mean - div_op <- str_if(reqn, " / ", " ./ ") - sdist("gamma", p$shape, paste0(p$shape, div_op, p$mu)) + sdist("gamma", p$shape, paste0(p$shape, " ./ ", p$mu)) } stan_log_lik_exponential <- function(bterms, ...) { @@ -681,18 +673,14 @@ stan_log_lik_exponential <- function(bterms, ...) { stan_log_lik_weibull <- function(bterms, ...) { p <- stan_log_lik_dpars(bterms) # Stan uses shape-scale parameterization for weibull - need_dot_div <- !stan_log_lik_adj(bterms) && is_pred_dpar(bterms, "shape") - div_op <- str_if(need_dot_div, " ./ ", " / ") - p$scale <- paste0(p$mu, div_op, "tgamma(1 + 1", div_op, p$shape, ")") + p$scale <- paste0(p$mu, " ./ tgamma(1 + 1 ./ ", p$shape, ")") sdist("weibull", p$shape, p$scale) } stan_log_lik_frechet <- function(bterms, ...) { p <- stan_log_lik_dpars(bterms) # Stan uses shape-scale parameterization for frechet - need_dot_div <- !stan_log_lik_adj(bterms) && is_pred_dpar(bterms, "nu") - div_op <- str_if(need_dot_div, " ./ ", " / ") - p$scale <- paste0(p$mu, div_op, "tgamma(1 - 1", div_op, p$nu, ")") + p$scale <- paste0(p$mu, " ./ tgamma(1 - 1 ./ ", p$nu, ")") sdist("frechet", p$nu, p$scale) } @@ -724,11 +712,9 @@ stan_log_lik_wiener <- function(bterms, ...) { stan_log_lik_beta <- function(bterms, ...) { p <- stan_log_lik_dpars(bterms) - req_dot_multiply <- !stan_log_lik_adj(bterms) && is_pred_dpar(bterms, "phi") - multiply <- str_if(req_dot_multiply, " .* ", " * ") sdist("beta", - paste0(p$mu, multiply, p$phi), - paste0("(1 - ", p$mu, ")", multiply, p$phi) + paste0(p$mu, " .* ", p$phi), + paste0("(1 - ", p$mu, ") .* ", p$phi) ) } diff --git a/tests/testthat/tests.stancode.R b/tests/testthat/tests.stancode.R index 5b120a07..88b221c4 100644 --- a/tests/testthat/tests.stancode.R +++ b/tests/testthat/tests.stancode.R @@ -410,8 +410,8 @@ test_that("customized covariances appear in the Stan code", { test_that("truncation appears in the Stan code", { scode <- stancode(time | trunc(0) ~ age + sex + disease, data = kidney, family = "gamma") - expect_match2(scode, "target += gamma_lpdf(Y[n] | shape, shape / mu[n]) -") - expect_match2(scode, "gamma_lccdf(lb[n] | shape, shape / mu[n]);") + expect_match2(scode, "target += gamma_lpdf(Y[n] | shape, shape ./ mu[n]) -") + expect_match2(scode, "gamma_lccdf(lb[n] | shape, shape ./ mu[n]);") scode <- stancode(time | trunc(ub = 100) ~ age + sex + disease, data = kidney, family = student("log")) @@ -620,7 +620,7 @@ test_that("Stan code for multivariate models is correct", { # multivariate weibull models bform <- bform + weibull() scode <- stancode(bform, dat) - expect_match2(scode, "weibull_lpdf(Y_g | shape_g, mu_g / tgamma(1 + 1 / shape_g));") + expect_match2(scode, "weibull_lpdf(Y_g | shape_g, mu_g ./ tgamma(1 + 1 ./ shape_g));") }) test_that("Stan code for categorical models is correct", { @@ -1419,6 +1419,9 @@ test_that("weighted, censored, and truncated likelihoods are correct", { scode <- stancode(y | cens(x) ~ 1, dat, family = asym_laplace()) expect_match2(scode, "target += asym_laplace_lccdf(Y[n] | mu[n], sigma, quantile);") + scode <- stancode(bf(y | cens(x) ~ 1, shape ~ 1), dat, family = Gamma()) + expect_match2(scode, "target += gamma_lpdf(Y[Jevent[1:Nevent]] | shape[Jevent[1:Nevent]], shape[Jevent[1:Nevent]] ./ mu[Jevent[1:Nevent]]);") + dat$x[1] <- 2 scode <- stancode(y | cens(x, y2) ~ 1, dat, family = asym_laplace()) expect_match2(scode, "target += log_diff_exp(\n") @@ -1442,15 +1445,15 @@ test_that("weighted, censored, and truncated likelihoods are correct", { expect_match2( stancode(y | trials(y2) + weights(y2) ~ 1, dat, beta_binomial()), - "target += weights[n] * (beta_binomial_lpmf(Y[n] | trials[n], mu[n] * phi," + "target += weights[n] * (beta_binomial_lpmf(Y[n] | trials[n], mu[n] .* phi," ) expect_match2( stancode(y | trials(y2) + trunc(0, 30) ~ 1, dat, beta_binomial()), - "log_diff_exp(beta_binomial_lcdf(ub[n] | trials[n], mu[n] * phi," + "log_diff_exp(beta_binomial_lcdf(ub[n] | trials[n], mu[n] .* phi," ) expect_match2( stancode(y | trials(y2) + cens(x, y2) ~ 1, dat, beta_binomial()), - "beta_binomial_lcdf(rcens[n] | trials[n], mu[n] * phi," + "beta_binomial_lcdf(rcens[n] | trials[n], mu[n] .* phi," ) }) @@ -1684,7 +1687,7 @@ test_that("Stan code of addition term 'rate' is correct", { expect_match2(scode, "target += poisson_lpmf(Y | mu .* denom);") scode <- stancode(y | rate(time) ~ x, data, negbinomial()) - expect_match2(scode, "target += neg_binomial_2_log_lpmf(Y | mu + log_denom, shape * denom);") + expect_match2(scode, "target += neg_binomial_2_log_lpmf(Y | mu + log_denom, shape .* denom);") bform <- bf(y | rate(time) ~ mi(x), shape ~ mi(x), family = negbinomial()) + bf(x | mi() ~ 1, family = gaussian()) @@ -1692,11 +1695,11 @@ test_that("Stan code of addition term 'rate' is correct", { expect_match2(scode, "target += neg_binomial_2_log_lpmf(Y_y | mu_y + log_denom_y, shape_y .* denom_y);") scode <- stancode(y | rate(time) ~ x, data, brmsfamily("negbinomial2")) - expect_match2(scode, "target += neg_binomial_2_log_lpmf(Y | mu + log_denom, inv(sigma) * denom);") + expect_match2(scode, "target += neg_binomial_2_log_lpmf(Y | mu + log_denom, inv(sigma) .* denom);") scode <- stancode(y | rate(time) + cens(1) ~ x, data, geometric()) expect_match2(scode, - "target += neg_binomial_2_lpmf(Y[Jevent[1:Nevent]] | mu[Jevent[1:Nevent]] .* denom[Jevent[1:Nevent]], 1 * denom[Jevent[1:Nevent]]);" + "target += neg_binomial_2_lpmf(Y[Jevent[1:Nevent]] | mu[Jevent[1:Nevent]] .* denom[Jevent[1:Nevent]], 1 .* denom[Jevent[1:Nevent]]);" ) }) @@ -1783,7 +1786,7 @@ test_that("Stan code of mixture model is correct", { data = data, mixture(Gamma("log"), weibull)) expect_match(scode, "data \\{[^\\}]*real theta1;") expect_match(scode, "data \\{[^\\}]*real theta2;") - expect_match2(scode, "ps[1] = log(theta1) + gamma_lpdf(Y[n] | shape1[n], shape1[n] / mu1[n]);") + expect_match2(scode, "ps[1] = log(theta1) + gamma_lpdf(Y[n] | shape1[n], shape1[n] ./ mu1[n]);") expect_match2(scode, "target += weights[n] * log_sum_exp(ps);") scode <- stancode(bf(abs(y) | se(c) ~ x), data = data, @@ -2146,7 +2149,7 @@ test_that("Stan code for missing value terms works correctly", { scode <- stancode(bform, dat) expect_match2(scode, "vector[Nmi_x] Ymi_x;") expect_match2(scode, - "target += beta_lpdf(Y_x[Jevent_x[1:Nevent_x]] | mu_x[Jevent_x[1:Nevent_x]] * phi_x, (1 - mu_x[Jevent_x[1:Nevent_x]]) * phi_x);" + "target += beta_lpdf(Y_x[Jevent_x[1:Nevent_x]] | mu_x[Jevent_x[1:Nevent_x]] .* phi_x, (1 - mu_x[Jevent_x[1:Nevent_x]]) .* phi_x);" ) # tests #1608