Skip to content

Commit

Permalink
Merge pull request #597 from stan-dev/array_syntax
Browse files Browse the repository at this point in the history 
Array syntax
  • Loading branch information
@jgabry
jgabry authored Sep 12, 2023
2 parents 2fade36 + 9a97ebd commit 80af4b6
Show file tree
Hide file tree
Showing 36 changed files with 1,099 additions and 910 deletions.
6 changes: 3 additions & 3 deletions DESCRIPTION
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
Package: rstanarm
Type: Package
Title: Bayesian Applied Regression Modeling via Stan
Version: 2.21.4
Date: 2023-04-03
Version: 2.26.1
Date: 2023-09-09
Encoding: UTF-8
Authors@R: c(person("Jonah", "Gabry", email = "[email protected]", role = "aut"),
person("Imad", "Ali", role = "ctb"),
Expand Down Expand Up @@ -42,7 +42,7 @@ Imports:
Matrix (>= 1.2-13),
nlme (>= 3.1-124),
posterior,
rstan (>= 2.21.1),
rstan (>= 2.26.1),
rstantools (>= 2.1.0),
shinystan (>= 2.3.0),
stats,
Expand Down
2 changes: 1 addition & 1 deletion R/doc-datasets.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -171,7 +171,7 @@
#'
#' @references
#' Carpenter, B. (2009) Bayesian estimators for the beta-binomial model of
#' batting ability. \url{https://lingpipe-blog.com/2009/09/23/}
#' batting ability. \url{https://web.archive.org/web/20220618114439/https://lingpipe-blog.com/2009/09/23/}
#'
#' Efron, B. and Morris, C. (1975) Data analysis using Stein's estimator and its
#' generalizations. \emph{Journal of the American Statistical Association}
Expand Down
214 changes: 119 additions & 95 deletions src/stan_files/bernoulli.stan
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,139 +3,157 @@

// GLM for a Bernoulli outcome
functions {
#include /functions/common_functions.stan
#include /functions/bernoulli_likelihoods.stan
#include /functions/common_functions.stan
#include /functions/bernoulli_likelihoods.stan
}
data {
// dimensions
int<lower=0> K; // number of predictors
int<lower=0> N[2]; // number of observations where y = 0 and y = 1 respectively
vector[K] xbar; // vector of column-means of rbind(X0, X1)
int<lower=0,upper=1> dense_X; // flag for dense vs. sparse
matrix[N[1],K] X0[dense_X]; // centered (by xbar) predictor matrix | y = 0
matrix[N[2],K] X1[dense_X]; // centered (by xbar) predictor matrix | y = 1

int<lower=0> K; // number of predictors
array[2] int<lower=0> N; // number of observations where y = 0 and y = 1 respectively
vector[K] xbar; // vector of column-means of rbind(X0, X1)
int<lower=0, upper=1> dense_X; // flag for dense vs. sparse
array[dense_X] matrix[N[1], K] X0; // centered (by xbar) predictor matrix | y = 0
array[dense_X] matrix[N[2], K] X1; // centered (by xbar) predictor matrix | y = 1
int<lower=0, upper=1> clogit; // 1 iff the number of successes is fixed in each stratum
int<lower=0> J; // number of strata (possibly zero)
int<lower=1,upper=J> strata[clogit == 1 ? N[1] + N[2] : 0];

array[clogit == 1 ? N[1] + N[2] : 0] int<lower=1, upper=J> strata;
// stuff for the sparse case
int<lower=0> nnz_X0; // number of non-zero elements in the implicit X0 matrix
vector[nnz_X0] w_X0; // non-zero elements in the implicit X0 matrix
int<lower=1, upper = K> v_X0[nnz_X0]; // column indices for w_X0
int<lower=0> nnz_X0; // number of non-zero elements in the implicit X0 matrix
vector[nnz_X0] w_X0; // non-zero elements in the implicit X0 matrix
array[nnz_X0] int<lower=1, upper=K> v_X0; // column indices for w_X0
// where the non-zeros start in each row of X0
int<lower=1, upper = rows(w_X0) + 1> u_X0[dense_X ? 0 : N[1] + 1];
int<lower=0> nnz_X1; // number of non-zero elements in the implicit X1 matrix
vector[nnz_X1] w_X1; // non-zero elements in the implicit X1 matrix
int<lower=1, upper = K> v_X1[nnz_X1]; // column indices for w_X1
array[dense_X ? 0 : N[1] + 1] int<lower=1, upper=rows(w_X0) + 1> u_X0;
int<lower=0> nnz_X1; // number of non-zero elements in the implicit X1 matrix
vector[nnz_X1] w_X1; // non-zero elements in the implicit X1 matrix
array[nnz_X1] int<lower=1, upper=K> v_X1; // column indices for w_X1
// where the non-zeros start in each row of X1
int<lower=1, upper = rows(w_X1) + 1> u_X1[dense_X ? 0 : N[2] + 1];
array[dense_X ? 0 : N[2] + 1] int<lower=1, upper=rows(w_X1) + 1> u_X1;
// declares prior_PD, has_intercept, link, prior_dist, prior_dist_for_intercept
#include /data/data_glm.stan

#include /data/data_glm.stan
int<lower=0> K_smooth;
matrix[N[1], K_smooth] S0;
matrix[N[2], K_smooth] S1;
int<lower=1> smooth_map[K_smooth];

int<lower=5,upper=5> family;

array[K_smooth] int<lower=1> smooth_map;
int<lower=5, upper=5> family;
// weights
int<lower=0,upper=1> has_weights; // 0 = No, 1 = Yes
int<lower=0, upper=1> has_weights; // 0 = No, 1 = Yes
vector[has_weights ? N[1] : 0] weights0;
vector[has_weights ? N[2] : 0] weights1;

// offset
int<lower=0,upper=1> has_offset; // 0 = No, 1 = Yes
int<lower=0, upper=1> has_offset; // 0 = No, 1 = Yes
vector[has_offset ? N[1] : 0] offset0;
vector[has_offset ? N[2] : 0] offset1;

// declares prior_{mean, scale, df}, prior_{mean, scale, df}_for_intercept, prior_{mean, scale, df}_for_aux
#include /data/hyperparameters.stan
#include /data/hyperparameters.stan

// declares t, p[t], l[t], q, len_theta_L, shape, scale, {len_}concentration, {len_}regularization
#include /data/glmer_stuff.stan

#include /data/glmer_stuff.stan
// more glmer stuff
int<lower=0> num_non_zero[2]; // number of non-zero elements in the Z matrices
vector[num_non_zero[1]] w0; // non-zero elements in the implicit Z0 matrix
vector[num_non_zero[2]] w1; // non-zero elements in the implicit Z1 matrix
int<lower=1, upper = q> v0[num_non_zero[1]]; // column indices for w0
int<lower=1, upper = q> v1[num_non_zero[2]]; // column indices for w1
array[2] int<lower=0> num_non_zero; // number of non-zero elements in the Z matrices
vector[num_non_zero[1]] w0; // non-zero elements in the implicit Z0 matrix
vector[num_non_zero[2]] w1; // non-zero elements in the implicit Z1 matrix
array[num_non_zero[1]] int<lower=1, upper=q> v0; // column indices for w0
array[num_non_zero[2]] int<lower=1, upper=q> v1; // column indices for w1
// where the non-zeros start in each row of Z0
int<lower=1, upper = rows(w0) + 1> u0[t > 0 ? N[1] + 1 : 0];
array[t > 0 ? N[1] + 1 : 0] int<lower=1, upper=rows(w0) + 1> u0;
// where the non-zeros start in each row of Z1
int<lower=1, upper = rows(w1) + 1> u1[t > 0 ? N[2] + 1 : 0];
int<lower=0, upper=1> special_case; // whether we only have to deal with (1|group)
array[t > 0 ? N[2] + 1 : 0] int<lower=1, upper=rows(w1) + 1> u1;
int<lower=0, upper=1> special_case; // whether we only have to deal with (1|group)
}
transformed data {
int NN = N[1] + N[2];
real aux = not_a_number();
int<lower=1> V0[special_case ? t : 0,N[1]] = make_V(N[1], special_case ? t : 0, v0);
int<lower=1> V1[special_case ? t : 0,N[2]] = make_V(N[2], special_case ? t : 0, v1);
int<lower=0> successes[clogit ? J : 0];
int<lower=0> failures[clogit ? J : 0];
int<lower=0> observations[clogit ? J : 0];

int can_do_bernoullilogitglm = K != 0 && // remove K!=0 after rstan includes this Stan bugfix: https://github.com/stan-dev/math/issues/1398
link == 1 && clogit == 0 && has_offset == 0 &&
prior_PD == 0 && dense_X == 1 && has_weights == 0 && t == 0;
matrix[can_do_bernoullilogitglm ? NN : 0, can_do_bernoullilogitglm ? K + K_smooth : 0] XS;
int y[can_do_bernoullilogitglm ? NN : 0];

array[special_case ? t : 0, N[1]] int<lower=1> V0 = make_V(N[1],
special_case ? t
: 0, v0);
array[special_case ? t : 0, N[2]] int<lower=1> V1 = make_V(N[2],
special_case ? t
: 0, v1);
array[clogit ? J : 0] int<lower=0> successes;
array[clogit ? J : 0] int<lower=0> failures;
array[clogit ? J : 0] int<lower=0> observations;

int can_do_bernoullilogitglm = K != 0
&& // remove K!=0 after rstan includes this Stan bugfix: https://github.com/stan-dev/math/issues/1398
link == 1 && clogit == 0 && has_offset == 0
&& prior_PD == 0 && dense_X == 1
&& has_weights == 0 && t == 0;
matrix[can_do_bernoullilogitglm ? NN : 0,
can_do_bernoullilogitglm ? K + K_smooth : 0] XS;
array[can_do_bernoullilogitglm ? NN : 0] int y;

// defines hs, len_z_T, len_var_group, delta, pos
#include /tdata/tdata_glm.stan
for (j in 1:J) {
#include /tdata/tdata_glm.stan

for (j in 1 : J) {
successes[j] = 0;
failures[j] = 0;
}
if (J > 0) for (i in 1:N[2]) successes[strata[i]] += 1;
if (J > 0) for (i in (N[2] + 1):NN) failures[strata[i]] += 1;
for (j in 1:J) observations[j] = failures[j] + successes[j];

if (J > 0)
for (i in 1 : N[2])
successes[strata[i]] += 1;
if (J > 0)
for (i in (N[2] + 1) : NN)
failures[strata[i]] += 1;
for (j in 1 : J)
observations[j] = failures[j] + successes[j];

if (can_do_bernoullilogitglm) {
XS = K_smooth > 0 ? append_col(append_row(X0[1], X1[1]), append_row(S0, S1)) : append_row(X0[1], X1[1]);
XS = K_smooth > 0
? append_col(append_row(X0[1], X1[1]), append_row(S0, S1))
: append_row(X0[1], X1[1]);
y = append_array(rep_array(0, N[1]), rep_array(1, N[2]));
}
}
parameters {
real<upper=(link == 4 ? 0.0 : positive_infinity())> gamma[has_intercept];
array[has_intercept] real<upper=(link == 4 ? 0.0 : positive_infinity())> gamma;
// declares z_beta, global, local, z_b, z_T, rho, zeta, tau
#include /parameters/parameters_glm.stan
#include /parameters/parameters_glm.stan
}
transformed parameters {
// defines beta, b, theta_L
#include /tparameters/tparameters_glm.stan
#include /tparameters/tparameters_glm.stan

if (t > 0) {
if (special_case) {
int start = 1;
theta_L = scale .* tau;
if (t == 1) b = theta_L[1] * z_b;
else for (i in 1:t) {
int end = start + l[i] - 1;
b[start:end] = theta_L[i] * z_b[start:end];
start = end + 1;
}
}
else {
theta_L = make_theta_L(len_theta_L, p,
1.0, tau, scale, zeta, rho, z_T);
if (t == 1)
b = theta_L[1] * z_b;
else
for (i in 1 : t) {
int end = start + l[i] - 1;
b[start : end] = theta_L[i] * z_b[start : end];
start = end + 1;
}
} else {
theta_L = make_theta_L(len_theta_L, p, 1.0, tau, scale, zeta, rho, z_T);
b = make_b(z_b, theta_L, p, l);
}
}
}
model {
if (can_do_bernoullilogitglm) {
vector[K + K_smooth] coeff = K_smooth > 0 ? append_row(beta, beta_smooth) : beta;
vector[K + K_smooth] coeff = K_smooth > 0 ? append_row(beta, beta_smooth)
: beta;
target += bernoulli_logit_glm_lpmf(y | XS, has_intercept ? gamma[1] : 0.0, coeff);
} else if (prior_PD == 0) {
// defines eta0, eta1
#include /model/make_eta_bern.stan
#include /model/make_eta_bern.stan

if (has_intercept == 1) {
if (link != 4) {
eta0 += gamma[1];
eta1 += gamma[1];
}
else {
} else {
real shift = fmax(max(eta0), max(eta1));
eta0 += gamma[1] - shift;
eta1 += gamma[1] - shift;
Expand All @@ -144,56 +162,62 @@ model {
// Log-likelihood
if (clogit) {
target += clogit_lpdf(eta0 | eta1, successes, failures, observations);
}
else if (has_weights == 0) {
} else if (has_weights == 0) {
target += bern_lpdf(eta0 | eta1, link, N);
}
else { // weighted log-likelihoods
} else {
// weighted log-likelihoods
target += dot_product(weights0, pw_bern(0, eta0, link));
target += dot_product(weights1, pw_bern(1, eta1, link));
}
}

#include /model/priors_glm.stan

#include /model/priors_glm.stan

if (t > 0) {
target += decov_lpdf(z_b | z_T, rho, zeta, tau,
regularization, delta, shape, t, p);
target += decov_lpdf(z_b | z_T, rho, zeta, tau, regularization, delta, shape, t, p);
}
}
generated quantities {
real mean_PPD = compute_mean_PPD ? 0 : negative_infinity();
real alpha[has_intercept];

array[has_intercept] real alpha;
if (has_intercept == 1) {
if (dense_X) alpha[1] = gamma[1] - dot_product(xbar, beta);
else alpha[1] = gamma[1];
if (dense_X)
alpha[1] = gamma[1] - dot_product(xbar, beta);
else
alpha[1] = gamma[1];
}

if (compute_mean_PPD) {
vector[N[1]] pi0;
vector[N[2]] pi1;
// defines eta0, eta1
#include /model/make_eta_bern.stan
#include /model/make_eta_bern.stan

if (has_intercept == 1) {
if (link != 4) {
eta0 += gamma[1];
eta1 += gamma[1];
}
else {
} else {
real shift;
shift = fmax(max(eta0), max(eta1));
eta0 += gamma[1] - shift;
eta1 += gamma[1] - shift;
alpha[1] -= shift;
}
}
if (clogit) for (j in 1:J) mean_PPD += successes[j]; // fixed by design
if (clogit)
for (j in 1 : J)
mean_PPD += successes[j]; // fixed by design
else {
pi0 = linkinv_bern(eta0, link);
pi1 = linkinv_bern(eta1, link);
for (n in 1:N[1]) mean_PPD += bernoulli_rng(pi0[n]);
for (n in 1:N[2]) mean_PPD += bernoulli_rng(pi1[n]);
for (n in 1 : N[1])
mean_PPD += bernoulli_rng(pi0[n]);
for (n in 1 : N[2])
mean_PPD += bernoulli_rng(pi1[n]);
}
mean_PPD /= NN;
}
}

Loading

0 comments on commit 80af4b6

Please sign in to comment.