diff --git a/CRAN-RELEASE b/CRAN-RELEASE deleted file mode 100644 index 7f90958..0000000 --- a/CRAN-RELEASE +++ /dev/null @@ -1,2 +0,0 @@ -This package was submitted to CRAN on 2021-11-12. -Once it is accepted, delete this file and tag the release (commit 3c40def). diff --git a/CRAN-SUBMISSION b/CRAN-SUBMISSION index ad75faf..ba8118b 100644 --- a/CRAN-SUBMISSION +++ b/CRAN-SUBMISSION @@ -1,3 +1,3 @@ -Version: 1.1.0 -Date: 2022-02-09 16:29:23 UTC -SHA: b21d92a5ae04e0190728231bf84d3529c758df21 +Version: 1.2.0 +Date: 2023-09-11 22:08:20 UTC +SHA: ba10aa7dd9c9ea57675350652ba4e2bf4df1b7cf diff --git a/DESCRIPTION b/DESCRIPTION index c96a641..e6ff805 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: zoid Title: Bayesian Zero-and-One Inflated Dirichlet Regression Modelling -Version: 1.1.0 +Version: 1.2.0 Authors@R: c(person(given = "Eric J.", family = "Ward", @@ -37,16 +37,16 @@ License: GPL (>=3) Encoding: UTF-8 LazyData: true Roxygen: list(markdown = TRUE) -RoxygenNote: 7.1.2 +RoxygenNote: 7.2.3 Biarch: true URL: https://nwfsc-cb.github.io/zoid/ BugReports: https://github.com/nwfsc-cb/zoid/issues Depends: R (>= 3.4.0) Imports: - methods, - gtools, compositions, + gtools, + methods, Rcpp (>= 0.12.0), RcppParallel (>= 5.0.1), rstan (>= 2.26.0), diff --git a/NEWS.md b/NEWS.md index bb28255..c9e7a6a 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,7 @@ +# zoid 1.2.0 + +* Updated Stan code to reflect changes needed to be compatible with Stan 2.26 (arrays) + # zoid 1.1.0 * Changed treatment of zeros in Stan model. Package consistent with Jensen et al. (in review) diff --git a/R/fit_prior.R b/R/fit_prior.R index 1ad6da6..a9d941e 100644 --- a/R/fit_prior.R +++ b/R/fit_prior.R @@ -30,7 +30,7 @@ fit_prior <- function(n_bins, n_draws = 10000, target = 1 / n_bins, iterations = n_draws = n_draws, method = "BFGS" ), silent = TRUE) - if (class(o) != "try-error") { + if(!identical(class(o), "try-error")) { if (o$value < best) { best <- o$value best_value <- list( diff --git a/src/Makevars b/src/Makevars index 1e5cb61..e2ddefe 100644 --- a/src/Makevars +++ b/src/Makevars @@ -3,8 +3,8 @@ STANHEADERS_SRC = $(shell "$(R_HOME)/bin$(R_ARCH_BIN)/Rscript" -e "message()" -e "cat(system.file('include', 'src', package = 'StanHeaders', mustWork = TRUE))" -e "message()" | grep "StanHeaders") STANC_FLAGS = $(shell "$(R_HOME)/bin$(R_ARCH_BIN)/Rscript" -e "cat(ifelse(utils::packageVersion('rstan') >= 2.26, '-DUSE_STANC3',''))") -PKG_CPPFLAGS = -I"../inst/include" -I"$(STANHEADERS_SRC)" -DBOOST_DISABLE_ASSERTS -DEIGEN_NO_DEBUG -DBOOST_MATH_OVERFLOW_ERROR_POLICY=errno_on_error $(STANC_FLAGS) +PKG_CPPFLAGS = -I"../inst/include" -I"$(STANHEADERS_SRC)" -DBOOST_DISABLE_ASSERTS -DEIGEN_NO_DEBUG -DBOOST_MATH_OVERFLOW_ERROR_POLICY=errno_on_error $(STANC_FLAGS) -D_HAS_AUTO_PTR_ETC=0 PKG_CXXFLAGS = $(shell "$(R_HOME)/bin$(R_ARCH_BIN)/Rscript" -e "RcppParallel::CxxFlags()") $(shell "$(R_HOME)/bin$(R_ARCH_BIN)/Rscript" -e "StanHeaders:::CxxFlags()") PKG_LIBS = $(shell "$(R_HOME)/bin$(R_ARCH_BIN)/Rscript" -e "RcppParallel::RcppParallelLibs()") $(shell "$(R_HOME)/bin$(R_ARCH_BIN)/Rscript" -e "StanHeaders:::LdFlags()") -CXX_STD = CXX14 +CXX_STD = CXX17 diff --git a/src/RcppExports.o b/src/RcppExports.o index cad082a..e6f2f7d 100644 Binary files a/src/RcppExports.o and b/src/RcppExports.o differ diff --git a/src/stanExports_dirichregmod.h b/src/stanExports_dirichregmod.h index bb22bca..df0dcf6 100644 --- a/src/stanExports_dirichregmod.h +++ b/src/stanExports_dirichregmod.h @@ -17,1240 +17,1766 @@ #ifndef MODELS_HPP #define MODELS_HPP #define STAN__SERVICES__COMMAND_HPP +#ifndef USE_STANC3 +#define USE_STANC3 +#endif #include -// Code generated by Stan version 2.21.0 +// Code generated by stanc v2.26.1-4-gd72b68b7-dirty #include namespace model_dirichregmod_namespace { +inline void validate_positive_index(const char* var_name, const char* expr, + int val) { + if (val < 1) { + std::stringstream msg; + msg << "Found dimension size less than one in simplex declaration" + << "; variable=" << var_name << "; dimension size expression=" << expr + << "; expression value=" << val; + std::string msg_str(msg.str()); + throw std::invalid_argument(msg_str.c_str()); + } +} +inline void validate_unit_vector_index(const char* var_name, const char* expr, + int val) { + if (val <= 1) { + std::stringstream msg; + if (val == 1) { + msg << "Found dimension size one in unit vector declaration." + << " One-dimensional unit vector is discrete" + << " but the target distribution must be continuous." + << " variable=" << var_name << "; dimension size expression=" << expr; + } else { + msg << "Found dimension size less than one in unit vector declaration" + << "; variable=" << var_name << "; dimension size expression=" << expr + << "; expression value=" << val; + } + std::string msg_str(msg.str()); + throw std::invalid_argument(msg_str.c_str()); + } +} using std::istream; using std::string; using std::stringstream; using std::vector; +using std::pow; using stan::io::dump; using stan::math::lgamma; -using stan::model::prob_grad; +using stan::model::model_base_crtp; +using stan::model::rvalue; +using stan::model::cons_list; +using stan::model::index_uni; +using stan::model::index_max; +using stan::model::index_min; +using stan::model::index_min_max; +using stan::model::index_multi; +using stan::model::index_omni; +using stan::model::nil_index_list; using namespace stan::math; -static int current_statement_begin__; -stan::io::program_reader prog_reader__() { - stan::io::program_reader reader; - reader.add_event(0, 0, "start", "model_dirichregmod"); - reader.add_event(171, 169, "end", "model_dirichregmod"); - return reader; -} +using stan::math::pow; +stan::math::profile_map profiles__; +static int current_statement__= 0; +static const std::vector locations_array__ = {" (found before start of program)", + " (in 'dirichregmod', line 54, column 2 to column 36)", + " (in 'dirichregmod', line 55, column 2 to column 36)", + " (in 'dirichregmod', line 58, column 2 to column 20)", + " (in 'dirichregmod', line 59, column 2 to column 52)", + " (in 'dirichregmod', line 60, column 2 to column 51)", + " (in 'dirichregmod', line 61, column 2 to column 30)", + " (in 'dirichregmod', line 62, column 2 to column 47)", + " (in 'dirichregmod', line 64, column 2 to column 10)", + " (in 'dirichregmod', line 65, column 19 to column 38)", + " (in 'dirichregmod', line 65, column 18 to column 39)", + " (in 'dirichregmod', line 65, column 2 to column 39)", + " (in 'dirichregmod', line 67, column 4 to column 25)", + " (in 'dirichregmod', line 66, column 23 to line 68, column 3)", + " (in 'dirichregmod', line 66, column 2 to line 68, column 3)", + " (in 'dirichregmod', line 71, column 6 to column 32)", + " (in 'dirichregmod', line 70, column 25 to line 72, column 5)", + " (in 'dirichregmod', line 70, column 4 to line 72, column 5)", + " (in 'dirichregmod', line 69, column 26 to line 73, column 3)", + " (in 'dirichregmod', line 69, column 2 to line 73, column 3)", + " (in 'dirichregmod', line 76, column 11 to column 17)", + " (in 'dirichregmod', line 76, column 4 to column 26)", + " (in 'dirichregmod', line 78, column 6 to column 53)", + " (in 'dirichregmod', line 77, column 23 to line 79, column 5)", + " (in 'dirichregmod', line 77, column 4 to line 79, column 5)", + " (in 'dirichregmod', line 80, column 4 to column 29)", + " (in 'dirichregmod', line 82, column 6 to column 26)", + " (in 'dirichregmod', line 81, column 23 to line 83, column 5)", + " (in 'dirichregmod', line 81, column 4 to line 83, column 5)", + " (in 'dirichregmod', line 75, column 25 to line 84, column 3)", + " (in 'dirichregmod', line 75, column 2 to line 84, column 3)", + " (in 'dirichregmod', line 87, column 6 to column 45)", + " (in 'dirichregmod', line 86, column 23 to line 90, column 5)", + " (in 'dirichregmod', line 86, column 4 to line 90, column 5)", + " (in 'dirichregmod', line 93, column 6 to column 70)", + " (in 'dirichregmod', line 92, column 23 to line 94, column 5)", + " (in 'dirichregmod', line 92, column 4 to line 94, column 5)", + " (in 'dirichregmod', line 85, column 24 to line 95, column 3)", + " (in 'dirichregmod', line 85, column 2 to line 95, column 3)", + " (in 'dirichregmod', line 123, column 2 to column 18)", + " (in 'dirichregmod', line 124, column 2 to column 17)", + " (in 'dirichregmod', line 125, column 2 to column 42)", + " (in 'dirichregmod', line 126, column 2 to column 57)", + " (in 'dirichregmod', line 127, column 2 to column 61)", + " (in 'dirichregmod', line 128, column 2 to column 60)", + " (in 'dirichregmod', line 129, column 2 to column 25)", + " (in 'dirichregmod', line 132, column 6 to column 23)", + " (in 'dirichregmod', line 133, column 6 to column 64)", + " (in 'dirichregmod', line 135, column 8 to column 40)", + " (in 'dirichregmod', line 136, column 8 to column 43)", + " (in 'dirichregmod', line 137, column 8 to column 181)", + " (in 'dirichregmod', line 134, column 32 to line 138, column 7)", + " (in 'dirichregmod', line 134, column 6 to line 138, column 7)", + " (in 'dirichregmod', line 141, column 8 to column 55)", + " (in 'dirichregmod', line 142, column 33 to column 49)", + " (in 'dirichregmod', line 142, column 8 to column 49)", + " (in 'dirichregmod', line 143, column 8 to column 53)", + " (in 'dirichregmod', line 144, column 32 to column 48)", + " (in 'dirichregmod', line 144, column 8 to column 48)", + " (in 'dirichregmod', line 145, column 8 to column 66)", + " (in 'dirichregmod', line 147, column 10 to column 38)", + " (in 'dirichregmod', line 148, column 10 to column 41)", + " (in 'dirichregmod', line 149, column 10 to column 54)", + " (in 'dirichregmod', line 150, column 10 to column 39)", + " (in 'dirichregmod', line 146, column 34 to line 151, column 9)", + " (in 'dirichregmod', line 146, column 8 to line 151, column 9)", + " (in 'dirichregmod', line 140, column 22 to line 152, column 7)", + " (in 'dirichregmod', line 140, column 6 to line 152, column 7)", + " (in 'dirichregmod', line 131, column 23 to line 153, column 5)", + " (in 'dirichregmod', line 131, column 4 to line 153, column 5)", + " (in 'dirichregmod', line 130, column 24 to line 154, column 3)", + " (in 'dirichregmod', line 130, column 2 to line 154, column 3)", + " (in 'dirichregmod', line 98, column 2 to column 18)", + " (in 'dirichregmod', line 99, column 2 to column 17)", + " (in 'dirichregmod', line 101, column 4 to column 26)", + " (in 'dirichregmod', line 100, column 18 to line 102, column 3)", + " (in 'dirichregmod', line 100, column 2 to line 102, column 3)", + " (in 'dirichregmod', line 106, column 6 to column 41)", + " (in 'dirichregmod', line 105, column 27 to line 107, column 5)", + " (in 'dirichregmod', line 105, column 4 to line 107, column 5)", + " (in 'dirichregmod', line 104, column 22 to line 108, column 3)", + " (in 'dirichregmod', line 104, column 2 to line 108, column 3)", + " (in 'dirichregmod', line 112, column 6 to column 58)", + " (in 'dirichregmod', line 114, column 8 to column 40)", + " (in 'dirichregmod', line 115, column 8 to column 43)", + " (in 'dirichregmod', line 117, column 8 to column 175)", + " (in 'dirichregmod', line 113, column 32 to line 118, column 7)", + " (in 'dirichregmod', line 113, column 6 to line 118, column 7)", + " (in 'dirichregmod', line 110, column 23 to line 119, column 5)", + " (in 'dirichregmod', line 110, column 4 to line 119, column 5)", + " (in 'dirichregmod', line 109, column 24 to line 120, column 3)", + " (in 'dirichregmod', line 109, column 2 to line 120, column 3)", + " (in 'dirichregmod', line 2, column 2 to column 16)", + " (in 'dirichregmod', line 3, column 2 to column 13)", + " (in 'dirichregmod', line 4, column 9 to column 18)", + " (in 'dirichregmod', line 4, column 20 to column 26)", + " (in 'dirichregmod', line 4, column 2 to column 30)", + " (in 'dirichregmod', line 5, column 2 to column 14)", + " (in 'dirichregmod', line 6, column 9 to column 18)", + " (in 'dirichregmod', line 6, column 20 to column 27)", + " (in 'dirichregmod', line 6, column 2 to column 38)", + " (in 'dirichregmod', line 7, column 8 to column 14)", + " (in 'dirichregmod', line 7, column 15 to column 23)", + " (in 'dirichregmod', line 7, column 2 to column 38)", + " (in 'dirichregmod', line 8, column 2 to column 15)", + " (in 'dirichregmod', line 9, column 2 to column 15)", + " (in 'dirichregmod', line 10, column 2 to column 16)", + " (in 'dirichregmod', line 13, column 8 to column 17)", + " (in 'dirichregmod', line 13, column 18 to column 24)", + " (in 'dirichregmod', line 13, column 2 to column 38)", + " (in 'dirichregmod', line 15, column 8 to column 17)", + " (in 'dirichregmod', line 15, column 18 to column 24)", + " (in 'dirichregmod', line 15, column 2 to column 44)", + " (in 'dirichregmod', line 16, column 9 to column 18)", + " (in 'dirichregmod', line 16, column 20 to column 26)", + " (in 'dirichregmod', line 16, column 2 to column 33)", + " (in 'dirichregmod', line 17, column 9 to column 18)", + " (in 'dirichregmod', line 17, column 20 to column 26)", + " (in 'dirichregmod', line 17, column 2 to column 34)", + " (in 'dirichregmod', line 18, column 9 to column 18)", + " (in 'dirichregmod', line 18, column 2 to column 24)", + " (in 'dirichregmod', line 19, column 9 to column 15)", + " (in 'dirichregmod', line 19, column 2 to column 22)", + " (in 'dirichregmod', line 22, column 4 to column 16)", + " (in 'dirichregmod', line 21, column 21 to line 23, column 3)", + " (in 'dirichregmod', line 21, column 2 to line 23, column 3)", + " (in 'dirichregmod', line 26, column 4 to column 15)", + " (in 'dirichregmod', line 28, column 6 to column 31)", + " (in 'dirichregmod', line 27, column 23 to line 29, column 5)", + " (in 'dirichregmod', line 27, column 4 to line 29, column 5)", + " (in 'dirichregmod', line 25, column 24 to line 30, column 3)", + " (in 'dirichregmod', line 25, column 2 to line 30, column 3)", + " (in 'dirichregmod', line 36, column 8 to column 25)", + " (in 'dirichregmod', line 35, column 13 to line 37, column 7)", + " (in 'dirichregmod', line 34, column 8 to column 25)", + " (in 'dirichregmod', line 33, column 20 to line 35, column 7)", + " (in 'dirichregmod', line 33, column 6 to line 37, column 7)", + " (in 'dirichregmod', line 39, column 8 to column 31)", + " (in 'dirichregmod', line 38, column 40 to line 40, column 7)", + " (in 'dirichregmod', line 38, column 6 to line 40, column 7)", + " (in 'dirichregmod', line 32, column 23 to line 41, column 5)", + " (in 'dirichregmod', line 32, column 4 to line 41, column 5)", + " (in 'dirichregmod', line 31, column 24 to line 42, column 3)", + " (in 'dirichregmod', line 31, column 2 to line 42, column 3)", + " (in 'dirichregmod', line 47, column 10 to column 34)", + " (in 'dirichregmod', line 48, column 10 to column 44)", + " (in 'dirichregmod', line 46, column 34 to line 49, column 9)", + " (in 'dirichregmod', line 46, column 8 to line 49, column 9)", + " (in 'dirichregmod', line 45, column 25 to line 50, column 7)", + " (in 'dirichregmod', line 45, column 6 to line 50, column 7)", + " (in 'dirichregmod', line 44, column 26 to line 51, column 5)", + " (in 'dirichregmod', line 44, column 4 to line 51, column 5)", + " (in 'dirichregmod', line 54, column 18 to column 26)", + " (in 'dirichregmod', line 55, column 9 to column 17)", + " (in 'dirichregmod', line 55, column 18 to column 25)", + " (in 'dirichregmod', line 59, column 26 to column 35)", + " (in 'dirichregmod', line 59, column 37 to column 43)", + " (in 'dirichregmod', line 60, column 26 to column 35)", + " (in 'dirichregmod', line 60, column 37 to column 43)", + " (in 'dirichregmod', line 61, column 9 to column 15)", + " (in 'dirichregmod', line 61, column 16 to column 23)", + " (in 'dirichregmod', line 62, column 26 to column 35)", + " (in 'dirichregmod', line 62, column 36 to column 42)", + " (in 'dirichregmod', line 125, column 8 to column 17)", + " (in 'dirichregmod', line 125, column 26 to column 32)", + " (in 'dirichregmod', line 126, column 8 to column 26)", + " (in 'dirichregmod', line 126, column 35 to column 50)", + " (in 'dirichregmod', line 127, column 8 to column 26)", + " (in 'dirichregmod', line 127, column 27 to column 42)", + " (in 'dirichregmod', line 128, column 8 to column 26)", + " (in 'dirichregmod', line 128, column 27 to column 42)"}; #include -class model_dirichregmod - : public stan::model::model_base_crtp { +class model_dirichregmod final : public model_base_crtp { private: - int N_samples; - int N_bins; - matrix_d X; - int N_covar; - matrix_d design_X; - std::vector > prod_idx; - int overdisp; - int postpred; - double prior_sd; - std::vector > is_zero; - std::vector > is_proportion; - matrix_d logX; - matrix_d logNX; - vector_d ESS; - vector_d ones; + int N_samples; + int N_bins; + Eigen::Matrix X; + int N_covar; + Eigen::Matrix design_X; + std::vector> prod_idx; + int overdisp; + int postpred; + double prior_sd; + std::vector> is_zero; + std::vector> is_proportion; + Eigen::Matrix logX; + Eigen::Matrix logNX; + Eigen::Matrix ESS; + Eigen::Matrix ones; + int beta_raw_1dim__; + int ynew_1dim__; + int ynew_2dim__; + int newy_is_zero_1dim__; + int newy_is_zero_2dim__; + int newy_is_one_1dim__; + int newy_is_one_2dim__; + public: - model_dirichregmod(stan::io::var_context& context__, - std::ostream* pstream__ = 0) - : model_base_crtp(0) { - ctor_body(context__, 0, pstream__); + ~model_dirichregmod() { } + + inline std::string model_name() const final { return "model_dirichregmod"; } + inline std::vector model_compile_info() const noexcept { + return std::vector{"stanc_version = stanc3 v2.26.1-4-gd72b68b7-dirty", "stancflags = "}; + } + + + model_dirichregmod(stan::io::var_context& context__, + unsigned int random_seed__ = 0, + std::ostream* pstream__ = nullptr) : model_base_crtp(0) { + using local_scalar_t__ = double ; + boost::ecuyer1988 base_rng__ = + stan::services::util::create_rng(random_seed__, 0); + (void) base_rng__; // suppress unused var warning + static const char* function__ = "model_dirichregmod_namespace::model_dirichregmod"; + (void) function__; // suppress unused var warning + local_scalar_t__ DUMMY_VAR__(std::numeric_limits::quiet_NaN()); + (void) DUMMY_VAR__; // suppress unused var warning + + try { + int pos__; + pos__ = std::numeric_limits::min(); + + pos__ = 1; + current_statement__ = 92; + context__.validate_dims("data initialization","N_samples","int", + context__.to_vec()); + N_samples = std::numeric_limits::min(); + + current_statement__ = 92; + N_samples = context__.vals_i("N_samples")[(1 - 1)]; + current_statement__ = 93; + context__.validate_dims("data initialization","N_bins","int", + context__.to_vec()); + N_bins = std::numeric_limits::min(); + + current_statement__ = 93; + N_bins = context__.vals_i("N_bins")[(1 - 1)]; + current_statement__ = 94; + validate_non_negative_index("X", "N_samples", N_samples); + current_statement__ = 95; + validate_non_negative_index("X", "N_bins", N_bins); + current_statement__ = 96; + context__.validate_dims("data initialization","X","double", + context__.to_vec(N_samples, N_bins)); + X = Eigen::Matrix(N_samples, N_bins); + stan::math::fill(X, std::numeric_limits::quiet_NaN()); + + { + std::vector X_flat__; + current_statement__ = 96; + assign(X_flat__, nil_index_list(), context__.vals_r("X"), + "assigning variable X_flat__"); + current_statement__ = 96; + pos__ = 1; + current_statement__ = 96; + for (int sym1__ = 1; sym1__ <= N_bins; ++sym1__) { + current_statement__ = 96; + for (int sym2__ = 1; sym2__ <= N_samples; ++sym2__) { + current_statement__ = 96; + assign(X, + cons_list(index_uni(sym2__), + cons_list(index_uni(sym1__), nil_index_list())), + X_flat__[(pos__ - 1)], "assigning variable X"); + current_statement__ = 96; + pos__ = (pos__ + 1);}} + } + current_statement__ = 97; + context__.validate_dims("data initialization","N_covar","int", + context__.to_vec()); + N_covar = std::numeric_limits::min(); + + current_statement__ = 97; + N_covar = context__.vals_i("N_covar")[(1 - 1)]; + current_statement__ = 98; + validate_non_negative_index("design_X", "N_samples", N_samples); + current_statement__ = 99; + validate_non_negative_index("design_X", "N_covar", N_covar); + current_statement__ = 100; + context__.validate_dims("data initialization","design_X","double", + context__.to_vec(N_samples, N_covar)); + design_X = Eigen::Matrix(N_samples, N_covar); + stan::math::fill(design_X, std::numeric_limits::quiet_NaN()); + + { + std::vector design_X_flat__; + current_statement__ = 100; + assign(design_X_flat__, nil_index_list(), + context__.vals_r("design_X"), "assigning variable design_X_flat__"); + current_statement__ = 100; + pos__ = 1; + current_statement__ = 100; + for (int sym1__ = 1; sym1__ <= N_covar; ++sym1__) { + current_statement__ = 100; + for (int sym2__ = 1; sym2__ <= N_samples; ++sym2__) { + current_statement__ = 100; + assign(design_X, + cons_list(index_uni(sym2__), + cons_list(index_uni(sym1__), nil_index_list())), + design_X_flat__[(pos__ - 1)], "assigning variable design_X"); + current_statement__ = 100; + pos__ = (pos__ + 1);}} + } + current_statement__ = 101; + validate_non_negative_index("prod_idx", "N_bins", N_bins); + current_statement__ = 102; + validate_non_negative_index("prod_idx", "N_bins - 1", (N_bins - 1)); + current_statement__ = 103; + context__.validate_dims("data initialization","prod_idx","int", + context__.to_vec(N_bins, (N_bins - 1))); + prod_idx = std::vector>(N_bins, std::vector( + (N_bins - 1), std::numeric_limits::min())); + + { + std::vector prod_idx_flat__; + current_statement__ = 103; + assign(prod_idx_flat__, nil_index_list(), + context__.vals_i("prod_idx"), "assigning variable prod_idx_flat__"); + current_statement__ = 103; + pos__ = 1; + current_statement__ = 103; + for (int sym1__ = 1; sym1__ <= (N_bins - 1); ++sym1__) { + current_statement__ = 103; + for (int sym2__ = 1; sym2__ <= N_bins; ++sym2__) { + current_statement__ = 103; + assign(prod_idx, + cons_list(index_uni(sym2__), + cons_list(index_uni(sym1__), nil_index_list())), + prod_idx_flat__[(pos__ - 1)], "assigning variable prod_idx"); + current_statement__ = 103; + pos__ = (pos__ + 1);}} + } + current_statement__ = 104; + context__.validate_dims("data initialization","overdisp","int", + context__.to_vec()); + overdisp = std::numeric_limits::min(); + + current_statement__ = 104; + overdisp = context__.vals_i("overdisp")[(1 - 1)]; + current_statement__ = 105; + context__.validate_dims("data initialization","postpred","int", + context__.to_vec()); + postpred = std::numeric_limits::min(); + + current_statement__ = 105; + postpred = context__.vals_i("postpred")[(1 - 1)]; + current_statement__ = 106; + context__.validate_dims("data initialization","prior_sd","double", + context__.to_vec()); + prior_sd = std::numeric_limits::quiet_NaN(); + + current_statement__ = 106; + prior_sd = context__.vals_r("prior_sd")[(1 - 1)]; + current_statement__ = 107; + validate_non_negative_index("is_zero", "N_samples", N_samples); + current_statement__ = 108; + validate_non_negative_index("is_zero", "N_bins", N_bins); + current_statement__ = 109; + is_zero = std::vector>(N_samples, std::vector(N_bins, std::numeric_limits::min())); + + current_statement__ = 110; + validate_non_negative_index("is_proportion", "N_samples", N_samples); + current_statement__ = 111; + validate_non_negative_index("is_proportion", "N_bins", N_bins); + current_statement__ = 112; + is_proportion = std::vector>(N_samples, std::vector(N_bins, std::numeric_limits::min())); + + current_statement__ = 113; + validate_non_negative_index("logX", "N_samples", N_samples); + current_statement__ = 114; + validate_non_negative_index("logX", "N_bins", N_bins); + current_statement__ = 115; + logX = Eigen::Matrix(N_samples, N_bins); + stan::math::fill(logX, std::numeric_limits::quiet_NaN()); + + current_statement__ = 116; + validate_non_negative_index("logNX", "N_samples", N_samples); + current_statement__ = 117; + validate_non_negative_index("logNX", "N_bins", N_bins); + current_statement__ = 118; + logNX = Eigen::Matrix(N_samples, N_bins); + stan::math::fill(logNX, std::numeric_limits::quiet_NaN()); + + current_statement__ = 119; + validate_non_negative_index("ESS", "N_samples", N_samples); + current_statement__ = 120; + ESS = Eigen::Matrix(N_samples); + stan::math::fill(ESS, std::numeric_limits::quiet_NaN()); + + current_statement__ = 121; + validate_non_negative_index("ones", "N_bins", N_bins); + current_statement__ = 122; + ones = Eigen::Matrix(N_bins); + stan::math::fill(ones, std::numeric_limits::quiet_NaN()); + + current_statement__ = 125; + for (int i = 1; i <= N_bins; ++i) { + current_statement__ = 123; + assign(ones, cons_list(index_uni(i), nil_index_list()), 1, + "assigning variable ones");} + current_statement__ = 131; + for (int i = 1; i <= N_samples; ++i) { + current_statement__ = 126; + assign(ESS, cons_list(index_uni(i), nil_index_list()), 0, + "assigning variable ESS"); + current_statement__ = 129; + for (int j = 1; j <= N_bins; ++j) { + current_statement__ = 127; + assign(ESS, cons_list(index_uni(i), nil_index_list()), + (ESS[(i - 1)] + + rvalue(X, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), "X")), + "assigning variable ESS");}} + current_statement__ = 143; + for (int i = 1; i <= N_samples; ++i) { + current_statement__ = 141; + for (int j = 1; j <= N_bins; ++j) { + current_statement__ = 136; + if (logical_eq( + rvalue(X, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), "X"), 0)) { + current_statement__ = 134; + assign(is_zero, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), 1, + "assigning variable is_zero"); + } else { + current_statement__ = 132; + assign(is_zero, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), 0, + "assigning variable is_zero"); + } + current_statement__ = 139; + if ((primitive_value( + logical_lt( + rvalue(X, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), "X"), + ESS[(i - 1)])) && primitive_value( + logical_gt( + rvalue(X, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), "X"), 0)))) { + current_statement__ = 137; + assign(is_proportion, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), 1, + "assigning variable is_proportion"); + } }} + current_statement__ = 151; + for (int i = 1; i <= N_samples; ++i) { + current_statement__ = 149; + for (int j = 1; j <= N_bins; ++j) { + current_statement__ = 147; + if (logical_eq(is_proportion[(i - 1)][(j - 1)], 1)) { + current_statement__ = 144; + assign(logX, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), + stan::math::log( + rvalue(X, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), "X")), + "assigning variable logX"); + current_statement__ = 145; + assign(logNX, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), + stan::math::log( + (ESS[(i - 1)] - + rvalue(X, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), "X"))), + "assigning variable logNX"); + } }} + current_statement__ = 152; + validate_non_negative_index("phi_inv", "overdisp", overdisp); + current_statement__ = 153; + beta_raw_1dim__ = std::numeric_limits::min(); + + current_statement__ = 153; + beta_raw_1dim__ = (N_bins - 1); + current_statement__ = 153; + validate_non_negative_index("beta_raw", "N_bins - 1", beta_raw_1dim__); + current_statement__ = 154; + validate_non_negative_index("beta_raw", "N_covar", N_covar); + current_statement__ = 155; + validate_non_negative_index("p_zero", "N_samples", N_samples); + current_statement__ = 156; + validate_non_negative_index("p_zero", "N_bins", N_bins); + current_statement__ = 157; + validate_non_negative_index("p_one", "N_samples", N_samples); + current_statement__ = 158; + validate_non_negative_index("p_one", "N_bins", N_bins); + current_statement__ = 159; + validate_non_negative_index("beta", "N_bins", N_bins); + current_statement__ = 160; + validate_non_negative_index("beta", "N_covar", N_covar); + current_statement__ = 161; + validate_non_negative_index("mu", "N_samples", N_samples); + current_statement__ = 162; + validate_non_negative_index("mu", "N_bins", N_bins); + current_statement__ = 163; + validate_non_negative_index("log_lik", "N_samples", N_samples); + current_statement__ = 164; + validate_non_negative_index("log_lik", "N_bins", N_bins); + current_statement__ = 165; + ynew_1dim__ = std::numeric_limits::min(); + + current_statement__ = 165; + ynew_1dim__ = (N_samples * postpred); + current_statement__ = 165; + validate_non_negative_index("ynew", "N_samples * postpred", ynew_1dim__); + current_statement__ = 166; + ynew_2dim__ = std::numeric_limits::min(); + + current_statement__ = 166; + ynew_2dim__ = (N_bins * postpred); + current_statement__ = 166; + validate_non_negative_index("ynew", "N_bins * postpred", ynew_2dim__); + current_statement__ = 167; + newy_is_zero_1dim__ = std::numeric_limits::min(); + + current_statement__ = 167; + newy_is_zero_1dim__ = (N_samples * postpred); + current_statement__ = 167; + validate_non_negative_index("newy_is_zero", "N_samples * postpred", + newy_is_zero_1dim__); + current_statement__ = 168; + newy_is_zero_2dim__ = std::numeric_limits::min(); + + current_statement__ = 168; + newy_is_zero_2dim__ = (N_bins * postpred); + current_statement__ = 168; + validate_non_negative_index("newy_is_zero", "N_bins * postpred", + newy_is_zero_2dim__); + current_statement__ = 169; + newy_is_one_1dim__ = std::numeric_limits::min(); + + current_statement__ = 169; + newy_is_one_1dim__ = (N_samples * postpred); + current_statement__ = 169; + validate_non_negative_index("newy_is_one", "N_samples * postpred", + newy_is_one_1dim__); + current_statement__ = 170; + newy_is_one_2dim__ = std::numeric_limits::min(); + + current_statement__ = 170; + newy_is_one_2dim__ = (N_bins * postpred); + current_statement__ = 170; + validate_non_negative_index("newy_is_one", "N_bins * postpred", + newy_is_one_2dim__); + } catch (const std::exception& e) { + stan::lang::rethrow_located(e, locations_array__[current_statement__]); + // Next line prevents compiler griping about no return + throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***"); } - model_dirichregmod(stan::io::var_context& context__, - unsigned int random_seed__, - std::ostream* pstream__ = 0) - : model_base_crtp(0) { - ctor_body(context__, random_seed__, pstream__); + num_params_r__ = 0U; + + try { + num_params_r__ += overdisp; + num_params_r__ += beta_raw_1dim__ * N_covar; + } catch (const std::exception& e) { + stan::lang::rethrow_located(e, locations_array__[current_statement__]); + // Next line prevents compiler griping about no return + throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***"); } - void ctor_body(stan::io::var_context& context__, - unsigned int random_seed__, - std::ostream* pstream__) { - typedef double local_scalar_t__; - boost::ecuyer1988 base_rng__ = - stan::services::util::create_rng(random_seed__, 0); - (void) base_rng__; // suppress unused var warning - current_statement_begin__ = -1; - static const char* function__ = "model_dirichregmod_namespace::model_dirichregmod"; - (void) function__; // dummy to suppress unused var warning - size_t pos__; - (void) pos__; // dummy to suppress unused var warning - std::vector vals_i__; - std::vector vals_r__; - local_scalar_t__ DUMMY_VAR__(std::numeric_limits::quiet_NaN()); - (void) DUMMY_VAR__; // suppress unused var warning - try { - // initialize data block variables from context__ - current_statement_begin__ = 2; - context__.validate_dims("data initialization", "N_samples", "int", context__.to_vec()); - N_samples = int(0); - vals_i__ = context__.vals_i("N_samples"); - pos__ = 0; - N_samples = vals_i__[pos__++]; - current_statement_begin__ = 3; - context__.validate_dims("data initialization", "N_bins", "int", context__.to_vec()); - N_bins = int(0); - vals_i__ = context__.vals_i("N_bins"); - pos__ = 0; - N_bins = vals_i__[pos__++]; - current_statement_begin__ = 4; - validate_non_negative_index("X", "N_samples", N_samples); - validate_non_negative_index("X", "N_bins", N_bins); - context__.validate_dims("data initialization", "X", "matrix_d", context__.to_vec(N_samples,N_bins)); - X = Eigen::Matrix(N_samples, N_bins); - vals_r__ = context__.vals_r("X"); - pos__ = 0; - size_t X_j_2_max__ = N_bins; - size_t X_j_1_max__ = N_samples; - for (size_t j_2__ = 0; j_2__ < X_j_2_max__; ++j_2__) { - for (size_t j_1__ = 0; j_1__ < X_j_1_max__; ++j_1__) { - X(j_1__, j_2__) = vals_r__[pos__++]; - } - } - current_statement_begin__ = 5; - context__.validate_dims("data initialization", "N_covar", "int", context__.to_vec()); - N_covar = int(0); - vals_i__ = context__.vals_i("N_covar"); - pos__ = 0; - N_covar = vals_i__[pos__++]; - current_statement_begin__ = 6; - validate_non_negative_index("design_X", "N_samples", N_samples); - validate_non_negative_index("design_X", "N_covar", N_covar); - context__.validate_dims("data initialization", "design_X", "matrix_d", context__.to_vec(N_samples,N_covar)); - design_X = Eigen::Matrix(N_samples, N_covar); - vals_r__ = context__.vals_r("design_X"); - pos__ = 0; - size_t design_X_j_2_max__ = N_covar; - size_t design_X_j_1_max__ = N_samples; - for (size_t j_2__ = 0; j_2__ < design_X_j_2_max__; ++j_2__) { - for (size_t j_1__ = 0; j_1__ < design_X_j_1_max__; ++j_1__) { - design_X(j_1__, j_2__) = vals_r__[pos__++]; - } - } - current_statement_begin__ = 7; - validate_non_negative_index("prod_idx", "N_bins", N_bins); - validate_non_negative_index("prod_idx", "(N_bins - 1)", (N_bins - 1)); - context__.validate_dims("data initialization", "prod_idx", "int", context__.to_vec(N_bins,(N_bins - 1))); - prod_idx = std::vector >(N_bins, std::vector((N_bins - 1), int(0))); - vals_i__ = context__.vals_i("prod_idx"); - pos__ = 0; - size_t prod_idx_k_0_max__ = N_bins; - size_t prod_idx_k_1_max__ = (N_bins - 1); - for (size_t k_1__ = 0; k_1__ < prod_idx_k_1_max__; ++k_1__) { - for (size_t k_0__ = 0; k_0__ < prod_idx_k_0_max__; ++k_0__) { - prod_idx[k_0__][k_1__] = vals_i__[pos__++]; - } - } - current_statement_begin__ = 8; - context__.validate_dims("data initialization", "overdisp", "int", context__.to_vec()); - overdisp = int(0); - vals_i__ = context__.vals_i("overdisp"); - pos__ = 0; - overdisp = vals_i__[pos__++]; - current_statement_begin__ = 9; - context__.validate_dims("data initialization", "postpred", "int", context__.to_vec()); - postpred = int(0); - vals_i__ = context__.vals_i("postpred"); - pos__ = 0; - postpred = vals_i__[pos__++]; - current_statement_begin__ = 10; - context__.validate_dims("data initialization", "prior_sd", "double", context__.to_vec()); - prior_sd = double(0); - vals_r__ = context__.vals_r("prior_sd"); - pos__ = 0; - prior_sd = vals_r__[pos__++]; - // initialize transformed data variables - current_statement_begin__ = 13; - validate_non_negative_index("is_zero", "N_samples", N_samples); - validate_non_negative_index("is_zero", "N_bins", N_bins); - is_zero = std::vector >(N_samples, std::vector(N_bins, int(0))); - stan::math::fill(is_zero, std::numeric_limits::min()); - current_statement_begin__ = 15; - validate_non_negative_index("is_proportion", "N_samples", N_samples); - validate_non_negative_index("is_proportion", "N_bins", N_bins); - is_proportion = std::vector >(N_samples, std::vector(N_bins, int(0))); - stan::math::fill(is_proportion, std::numeric_limits::min()); - current_statement_begin__ = 16; - validate_non_negative_index("logX", "N_samples", N_samples); - validate_non_negative_index("logX", "N_bins", N_bins); - logX = Eigen::Matrix(N_samples, N_bins); - stan::math::fill(logX, DUMMY_VAR__); - current_statement_begin__ = 17; - validate_non_negative_index("logNX", "N_samples", N_samples); - validate_non_negative_index("logNX", "N_bins", N_bins); - logNX = Eigen::Matrix(N_samples, N_bins); - stan::math::fill(logNX, DUMMY_VAR__); - current_statement_begin__ = 18; - validate_non_negative_index("ESS", "N_samples", N_samples); - ESS = Eigen::Matrix(N_samples); - stan::math::fill(ESS, DUMMY_VAR__); - current_statement_begin__ = 19; - validate_non_negative_index("ones", "N_bins", N_bins); - ones = Eigen::Matrix(N_bins); - stan::math::fill(ones, DUMMY_VAR__); - // execute transformed data statements - current_statement_begin__ = 22; - for (int i = 1; i <= N_bins; ++i) { - current_statement_begin__ = 23; - stan::model::assign(ones, - stan::model::cons_list(stan::model::index_uni(i), stan::model::nil_index_list()), - 1, - "assigning variable ones"); - } - current_statement_begin__ = 26; - for (int i = 1; i <= N_samples; ++i) { - current_statement_begin__ = 27; - stan::model::assign(ESS, - stan::model::cons_list(stan::model::index_uni(i), stan::model::nil_index_list()), - 0, - "assigning variable ESS"); - current_statement_begin__ = 28; - for (int j = 1; j <= N_bins; ++j) { - current_statement_begin__ = 29; - stan::model::assign(ESS, - stan::model::cons_list(stan::model::index_uni(i), stan::model::nil_index_list()), - (get_base1(ESS, i, "ESS", 1) + get_base1(X, i, j, "X", 1)), - "assigning variable ESS"); - } - } - current_statement_begin__ = 32; - for (int i = 1; i <= N_samples; ++i) { - current_statement_begin__ = 33; - for (int j = 1; j <= N_bins; ++j) { - current_statement_begin__ = 34; - if (as_bool(logical_eq(get_base1(X, i, j, "X", 1), 0))) { - current_statement_begin__ = 35; - stan::model::assign(is_zero, - stan::model::cons_list(stan::model::index_uni(i), stan::model::cons_list(stan::model::index_uni(j), stan::model::nil_index_list())), - 1, - "assigning variable is_zero"); - } else { - current_statement_begin__ = 37; - stan::model::assign(is_zero, - stan::model::cons_list(stan::model::index_uni(i), stan::model::cons_list(stan::model::index_uni(j), stan::model::nil_index_list())), - 0, - "assigning variable is_zero"); - } - current_statement_begin__ = 40; - if (as_bool((primitive_value(logical_lt(get_base1(X, i, j, "X", 1), get_base1(ESS, i, "ESS", 1))) && primitive_value(logical_gt(get_base1(X, i, j, "X", 1), 0))))) { - current_statement_begin__ = 41; - stan::model::assign(is_proportion, - stan::model::cons_list(stan::model::index_uni(i), stan::model::cons_list(stan::model::index_uni(j), stan::model::nil_index_list())), - 1, - "assigning variable is_proportion"); - } - } - } - current_statement_begin__ = 48; - for (int i = 1; i <= N_samples; ++i) { - current_statement_begin__ = 49; - for (int j = 1; j <= N_bins; ++j) { - current_statement_begin__ = 50; - if (as_bool(logical_eq(get_base1(get_base1(is_proportion, i, "is_proportion", 1), j, "is_proportion", 2), 1))) { - current_statement_begin__ = 51; - stan::model::assign(logX, - stan::model::cons_list(stan::model::index_uni(i), stan::model::cons_list(stan::model::index_uni(j), stan::model::nil_index_list())), - stan::math::log(get_base1(X, i, j, "X", 1)), - "assigning variable logX"); - current_statement_begin__ = 52; - stan::model::assign(logNX, - stan::model::cons_list(stan::model::index_uni(i), stan::model::cons_list(stan::model::index_uni(j), stan::model::nil_index_list())), - stan::math::log((get_base1(ESS, i, "ESS", 1) - get_base1(X, i, j, "X", 1))), - "assigning variable logNX"); - } - } - } - // validate transformed data - // validate, set parameter ranges - num_params_r__ = 0U; - param_ranges_i__.clear(); - current_statement_begin__ = 59; - validate_non_negative_index("phi_inv", "overdisp", overdisp); - num_params_r__ += overdisp; - current_statement_begin__ = 60; - validate_non_negative_index("beta_raw", "(N_bins - 1)", (N_bins - 1)); - validate_non_negative_index("beta_raw", "N_covar", N_covar); - num_params_r__ += ((N_bins - 1) * N_covar); - } catch (const std::exception& e) { - stan::lang::rethrow_located(e, current_statement_begin__, prog_reader__()); - // Next line prevents compiler griping about no return - throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***"); - } + } + template * = nullptr, stan::require_vector_like_vt* = nullptr> + inline stan::scalar_type_t log_prob_impl(VecR& params_r__, + VecI& params_i__, + std::ostream* pstream__ = nullptr) const { + using T__ = stan::scalar_type_t; + using local_scalar_t__ = T__; + T__ lp__(0.0); + stan::math::accumulator lp_accum__; + static const char* function__ = "model_dirichregmod_namespace::log_prob"; +(void) function__; // suppress unused var warning + stan::io::reader in__(params_r__, params_i__); + local_scalar_t__ DUMMY_VAR__(std::numeric_limits::quiet_NaN()); + (void) DUMMY_VAR__; // suppress unused var warning + + try { + Eigen::Matrix phi_inv; + phi_inv = Eigen::Matrix(overdisp); + stan::math::fill(phi_inv, DUMMY_VAR__); + + current_statement__ = 1; + phi_inv = in__.vector(overdisp); + current_statement__ = 1; + for (int sym1__ = 1; sym1__ <= overdisp; ++sym1__) { + current_statement__ = 1; + if (jacobian__) { + current_statement__ = 1; + assign(phi_inv, cons_list(index_uni(sym1__), nil_index_list()), + stan::math::lb_constrain(phi_inv[(sym1__ - 1)], 0, lp__), + "assigning variable phi_inv"); + } else { + current_statement__ = 1; + assign(phi_inv, cons_list(index_uni(sym1__), nil_index_list()), + stan::math::lb_constrain(phi_inv[(sym1__ - 1)], 0), + "assigning variable phi_inv"); + }} + Eigen::Matrix beta_raw; + beta_raw = Eigen::Matrix(beta_raw_1dim__, N_covar); + stan::math::fill(beta_raw, DUMMY_VAR__); + + current_statement__ = 2; + beta_raw = in__.matrix(beta_raw_1dim__, N_covar); + local_scalar_t__ phi; + phi = DUMMY_VAR__; + + Eigen::Matrix p_zero; + p_zero = Eigen::Matrix(N_samples, N_bins); + stan::math::fill(p_zero, DUMMY_VAR__); + + Eigen::Matrix p_one; + p_one = Eigen::Matrix(N_samples, N_bins); + stan::math::fill(p_one, DUMMY_VAR__); + + Eigen::Matrix beta; + beta = Eigen::Matrix(N_bins, N_covar); + stan::math::fill(beta, DUMMY_VAR__); + + Eigen::Matrix mu; + mu = Eigen::Matrix(N_samples, N_bins); + stan::math::fill(mu, DUMMY_VAR__); + + current_statement__ = 8; + phi = 1; + current_statement__ = 11; + if (logical_eq(overdisp, 1)) { + current_statement__ = 9; + phi = (1 / phi_inv[(1 - 1)]); + } + current_statement__ = 14; + for (int l = 1; l <= N_covar; ++l) { + current_statement__ = 12; + assign(beta, + cons_list(index_uni(N_bins), + cons_list(index_uni(l), nil_index_list())), 0.0, + "assigning variable beta");} + current_statement__ = 19; + for (int k = 1; k <= (N_bins - 1); ++k) { + current_statement__ = 17; + for (int l = 1; l <= N_covar; ++l) { + current_statement__ = 15; + assign(beta, + cons_list(index_uni(k), + cons_list(index_uni(l), nil_index_list())), + rvalue(beta_raw, + cons_list(index_uni(k), + cons_list(index_uni(l), nil_index_list())), "beta_raw"), + "assigning variable beta");}} + current_statement__ = 30; + for (int n = 1; n <= N_samples; ++n) { + current_statement__ = 20; + validate_non_negative_index("logits", "N_bins", N_bins); + Eigen::Matrix logits; + logits = Eigen::Matrix(N_bins); + stan::math::fill(logits, DUMMY_VAR__); + + current_statement__ = 24; + for (int m = 1; m <= N_bins; ++m) { + current_statement__ = 22; + assign(logits, cons_list(index_uni(m), nil_index_list()), + multiply( + rvalue(design_X, + cons_list(index_uni(n), + cons_list(index_omni(), nil_index_list())), "design_X"), + transpose( + rvalue(beta, + cons_list(index_uni(m), + cons_list(index_omni(), nil_index_list())), "beta"))), + "assigning variable logits");} + current_statement__ = 25; + assign(logits, nil_index_list(), + softmax(stan::model::deep_copy(logits)), + "assigning variable logits"); + current_statement__ = 28; + for (int m = 1; m <= N_bins; ++m) { + current_statement__ = 26; + assign(mu, + cons_list(index_uni(n), + cons_list(index_uni(m), nil_index_list())), logits[(m - 1)], + "assigning variable mu");}} + current_statement__ = 38; + for (int i = 1; i <= N_samples; ++i) { + current_statement__ = 33; + for (int j = 1; j <= N_bins; ++j) { + current_statement__ = 31; + assign(p_zero, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), + pow( + (1 - + rvalue(mu, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), "mu")), + (ESS[(i - 1)] * phi)), "assigning variable p_zero");} + current_statement__ = 36; + for (int j = 1; j <= N_bins; ++j) { + current_statement__ = 34; + assign(p_one, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), + ((1 - + rvalue(p_zero, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), "p_zero")) * + prod( + rvalue(p_zero, + cons_list(index_uni(i), + cons_list( + index_multi(rvalue(prod_idx, + cons_list(index_uni(j), + cons_list(index_omni(), + nil_index_list())), "prod_idx")), + nil_index_list())), "p_zero"))), + "assigning variable p_one");}} + current_statement__ = 3; + current_statement__ = 3; + check_greater_or_equal(function__, "phi", phi, 0); + current_statement__ = 4; + for (int sym1__ = 1; sym1__ <= N_samples; ++sym1__) { + current_statement__ = 4; + for (int sym2__ = 1; sym2__ <= N_bins; ++sym2__) { + current_statement__ = 4; + current_statement__ = 4; + check_greater_or_equal(function__, "p_zero[sym1__, sym2__]", + rvalue(p_zero, + cons_list(index_uni(sym1__), + cons_list(index_uni(sym2__), + nil_index_list())), "p_zero"), 0);}} + current_statement__ = 4; + for (int sym1__ = 1; sym1__ <= N_samples; ++sym1__) { + current_statement__ = 4; + for (int sym2__ = 1; sym2__ <= N_bins; ++sym2__) { + current_statement__ = 4; + current_statement__ = 4; + check_less_or_equal(function__, "p_zero[sym1__, sym2__]", + rvalue(p_zero, + cons_list(index_uni(sym1__), + cons_list(index_uni(sym2__), + nil_index_list())), "p_zero"), 1);}} + current_statement__ = 5; + for (int sym1__ = 1; sym1__ <= N_samples; ++sym1__) { + current_statement__ = 5; + for (int sym2__ = 1; sym2__ <= N_bins; ++sym2__) { + current_statement__ = 5; + current_statement__ = 5; + check_greater_or_equal(function__, "p_one[sym1__, sym2__]", + rvalue(p_one, + cons_list(index_uni(sym1__), + cons_list(index_uni(sym2__), + nil_index_list())), "p_one"), 0);}} + current_statement__ = 5; + for (int sym1__ = 1; sym1__ <= N_samples; ++sym1__) { + current_statement__ = 5; + for (int sym2__ = 1; sym2__ <= N_bins; ++sym2__) { + current_statement__ = 5; + current_statement__ = 5; + check_less_or_equal(function__, "p_one[sym1__, sym2__]", + rvalue(p_one, + cons_list(index_uni(sym1__), + cons_list(index_uni(sym2__), + nil_index_list())), "p_one"), 1);}} + current_statement__ = 7; + for (int sym1__ = 1; sym1__ <= N_samples; ++sym1__) { + current_statement__ = 7; + for (int sym2__ = 1; sym2__ <= N_bins; ++sym2__) { + current_statement__ = 7; + current_statement__ = 7; + check_greater_or_equal(function__, "mu[sym1__, sym2__]", + rvalue(mu, + cons_list(index_uni(sym1__), + cons_list(index_uni(sym2__), + nil_index_list())), "mu"), 0);}} + current_statement__ = 7; + for (int sym1__ = 1; sym1__ <= N_samples; ++sym1__) { + current_statement__ = 7; + for (int sym2__ = 1; sym2__ <= N_bins; ++sym2__) { + current_statement__ = 7; + current_statement__ = 7; + check_less_or_equal(function__, "mu[sym1__, sym2__]", + rvalue(mu, + cons_list(index_uni(sym1__), + cons_list(index_uni(sym2__), + nil_index_list())), "mu"), 1);}} + { + local_scalar_t__ alpha_temp; + alpha_temp = DUMMY_VAR__; + + local_scalar_t__ beta_temp; + beta_temp = DUMMY_VAR__; + + current_statement__ = 76; + if (logical_eq(overdisp, 1)) { + current_statement__ = 74; + lp_accum__.add(cauchy_lpdf(phi_inv, 0, 5)); + } + current_statement__ = 81; + for (int i = 1; i <= N_covar; ++i) { + current_statement__ = 79; + for (int j = 1; j <= (N_bins - 1); ++j) { + current_statement__ = 77; + lp_accum__.add( + normal_lpdf( + rvalue(beta_raw, + cons_list(index_uni(j), + cons_list(index_uni(i), nil_index_list())), "beta_raw"), + 0, prior_sd));}} + current_statement__ = 91; + for (int i = 1; i <= N_samples; ++i) { + current_statement__ = 89; + for (int j = 1; j <= N_bins; ++j) { + current_statement__ = 82; + lp_accum__.add( + bernoulli_lpmf(is_zero[(i - 1)][(j - 1)], + rvalue(p_zero, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), "p_zero"))); + current_statement__ = 87; + if (logical_eq(is_proportion[(i - 1)][(j - 1)], 1)) { + current_statement__ = 83; + alpha_temp = ((rvalue(mu, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), + "mu") * ESS[(i - 1)]) * phi); + current_statement__ = 84; + beta_temp = (((1 - + rvalue(mu, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), + "mu")) * ESS[(i - 1)]) * phi); + current_statement__ = 85; + lp_accum__.add( + ((((stan::math::log( + ((1 - + rvalue(p_zero, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), + "p_zero")) - + rvalue(p_one, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), + "p_one"))) + + ((alpha_temp - 1) * + rvalue(logX, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), + "logX"))) + + ((beta_temp - 1) * + rvalue(logNX, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), + "logNX"))) - + (((alpha_temp + beta_temp) - 1) * + stan::math::log(ESS[(i - 1)]))) - + lbeta(alpha_temp, beta_temp))); + } }} + } + } catch (const std::exception& e) { + stan::lang::rethrow_located(e, locations_array__[current_statement__]); + // Next line prevents compiler griping about no return + throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***"); } - ~model_dirichregmod() { } - void transform_inits(const stan::io::var_context& context__, - std::vector& params_i__, - std::vector& params_r__, - std::ostream* pstream__) const { - typedef double local_scalar_t__; - stan::io::writer writer__(params_r__, params_i__); - size_t pos__; - (void) pos__; // dummy call to supress warning - std::vector vals_r__; - std::vector vals_i__; - current_statement_begin__ = 59; - if (!(context__.contains_r("phi_inv"))) - stan::lang::rethrow_located(std::runtime_error(std::string("Variable phi_inv missing")), current_statement_begin__, prog_reader__()); - vals_r__ = context__.vals_r("phi_inv"); - pos__ = 0U; - validate_non_negative_index("phi_inv", "overdisp", overdisp); - context__.validate_dims("parameter initialization", "phi_inv", "vector_d", context__.to_vec(overdisp)); - Eigen::Matrix phi_inv(overdisp); - size_t phi_inv_j_1_max__ = overdisp; - for (size_t j_1__ = 0; j_1__ < phi_inv_j_1_max__; ++j_1__) { - phi_inv(j_1__) = vals_r__[pos__++]; - } - try { - writer__.vector_lb_unconstrain(0, phi_inv); - } catch (const std::exception& e) { - stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable phi_inv: ") + e.what()), current_statement_begin__, prog_reader__()); + lp_accum__.add(lp__); + return lp_accum__.sum(); + } // log_prob_impl() + + template * = nullptr, stan::require_vector_like_vt* = nullptr, stan::require_std_vector_vt* = nullptr> + inline void write_array_impl(RNG& base_rng__, VecR& params_r__, + VecI& params_i__, VecVar& vars__, + const bool emit_transformed_parameters__ = true, + const bool emit_generated_quantities__ = true, + std::ostream* pstream__ = nullptr) const { + using local_scalar_t__ = double; + vars__.resize(0); + stan::io::reader in__(params_r__, params_i__); + static const char* function__ = "model_dirichregmod_namespace::write_array"; +(void) function__; // suppress unused var warning + (void) function__; // suppress unused var warning + double lp__ = 0.0; + (void) lp__; // dummy to suppress unused var warning + stan::math::accumulator lp_accum__; + local_scalar_t__ DUMMY_VAR__(std::numeric_limits::quiet_NaN()); + (void) DUMMY_VAR__; // suppress unused var warning + + try { + Eigen::Matrix phi_inv; + phi_inv = Eigen::Matrix(overdisp); + stan::math::fill(phi_inv, std::numeric_limits::quiet_NaN()); + + current_statement__ = 1; + phi_inv = in__.vector(overdisp); + current_statement__ = 1; + for (int sym1__ = 1; sym1__ <= overdisp; ++sym1__) { + current_statement__ = 1; + assign(phi_inv, cons_list(index_uni(sym1__), nil_index_list()), + stan::math::lb_constrain(phi_inv[(sym1__ - 1)], 0), + "assigning variable phi_inv");} + Eigen::Matrix beta_raw; + beta_raw = Eigen::Matrix(beta_raw_1dim__, N_covar); + stan::math::fill(beta_raw, std::numeric_limits::quiet_NaN()); + + current_statement__ = 2; + beta_raw = in__.matrix(beta_raw_1dim__, N_covar); + double phi; + phi = std::numeric_limits::quiet_NaN(); + + Eigen::Matrix p_zero; + p_zero = Eigen::Matrix(N_samples, N_bins); + stan::math::fill(p_zero, std::numeric_limits::quiet_NaN()); + + Eigen::Matrix p_one; + p_one = Eigen::Matrix(N_samples, N_bins); + stan::math::fill(p_one, std::numeric_limits::quiet_NaN()); + + Eigen::Matrix beta; + beta = Eigen::Matrix(N_bins, N_covar); + stan::math::fill(beta, std::numeric_limits::quiet_NaN()); + + Eigen::Matrix mu; + mu = Eigen::Matrix(N_samples, N_bins); + stan::math::fill(mu, std::numeric_limits::quiet_NaN()); + + for (int sym1__ = 1; sym1__ <= overdisp; ++sym1__) { + vars__.emplace_back(phi_inv[(sym1__ - 1)]);} + for (int sym1__ = 1; sym1__ <= N_covar; ++sym1__) { + for (int sym2__ = 1; sym2__ <= beta_raw_1dim__; ++sym2__) { + vars__.emplace_back( + rvalue(beta_raw, + cons_list(index_uni(sym2__), + cons_list(index_uni(sym1__), nil_index_list())), "beta_raw")); + }} + if (logical_negation((primitive_value(emit_transformed_parameters__) || + primitive_value(emit_generated_quantities__)))) { + return ; + } + current_statement__ = 8; + phi = 1; + current_statement__ = 11; + if (logical_eq(overdisp, 1)) { + current_statement__ = 9; + phi = (1 / phi_inv[(1 - 1)]); + } + current_statement__ = 14; + for (int l = 1; l <= N_covar; ++l) { + current_statement__ = 12; + assign(beta, + cons_list(index_uni(N_bins), + cons_list(index_uni(l), nil_index_list())), 0.0, + "assigning variable beta");} + current_statement__ = 19; + for (int k = 1; k <= (N_bins - 1); ++k) { + current_statement__ = 17; + for (int l = 1; l <= N_covar; ++l) { + current_statement__ = 15; + assign(beta, + cons_list(index_uni(k), + cons_list(index_uni(l), nil_index_list())), + rvalue(beta_raw, + cons_list(index_uni(k), + cons_list(index_uni(l), nil_index_list())), "beta_raw"), + "assigning variable beta");}} + current_statement__ = 30; + for (int n = 1; n <= N_samples; ++n) { + current_statement__ = 20; + validate_non_negative_index("logits", "N_bins", N_bins); + Eigen::Matrix logits; + logits = Eigen::Matrix(N_bins); + stan::math::fill(logits, std::numeric_limits::quiet_NaN()); + + current_statement__ = 24; + for (int m = 1; m <= N_bins; ++m) { + current_statement__ = 22; + assign(logits, cons_list(index_uni(m), nil_index_list()), + multiply( + rvalue(design_X, + cons_list(index_uni(n), + cons_list(index_omni(), nil_index_list())), "design_X"), + transpose( + rvalue(beta, + cons_list(index_uni(m), + cons_list(index_omni(), nil_index_list())), "beta"))), + "assigning variable logits");} + current_statement__ = 25; + assign(logits, nil_index_list(), + softmax(stan::model::deep_copy(logits)), + "assigning variable logits"); + current_statement__ = 28; + for (int m = 1; m <= N_bins; ++m) { + current_statement__ = 26; + assign(mu, + cons_list(index_uni(n), + cons_list(index_uni(m), nil_index_list())), logits[(m - 1)], + "assigning variable mu");}} + current_statement__ = 38; + for (int i = 1; i <= N_samples; ++i) { + current_statement__ = 33; + for (int j = 1; j <= N_bins; ++j) { + current_statement__ = 31; + assign(p_zero, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), + pow( + (1 - + rvalue(mu, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), "mu")), + (ESS[(i - 1)] * phi)), "assigning variable p_zero");} + current_statement__ = 36; + for (int j = 1; j <= N_bins; ++j) { + current_statement__ = 34; + assign(p_one, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), + ((1 - + rvalue(p_zero, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), "p_zero")) * + prod( + rvalue(p_zero, + cons_list(index_uni(i), + cons_list( + index_multi(rvalue(prod_idx, + cons_list(index_uni(j), + cons_list(index_omni(), + nil_index_list())), "prod_idx")), + nil_index_list())), "p_zero"))), + "assigning variable p_one");}} + current_statement__ = 3; + current_statement__ = 3; + check_greater_or_equal(function__, "phi", phi, 0); + current_statement__ = 4; + for (int sym1__ = 1; sym1__ <= N_samples; ++sym1__) { + current_statement__ = 4; + for (int sym2__ = 1; sym2__ <= N_bins; ++sym2__) { + current_statement__ = 4; + current_statement__ = 4; + check_greater_or_equal(function__, "p_zero[sym1__, sym2__]", + rvalue(p_zero, + cons_list(index_uni(sym1__), + cons_list(index_uni(sym2__), + nil_index_list())), "p_zero"), 0);}} + current_statement__ = 4; + for (int sym1__ = 1; sym1__ <= N_samples; ++sym1__) { + current_statement__ = 4; + for (int sym2__ = 1; sym2__ <= N_bins; ++sym2__) { + current_statement__ = 4; + current_statement__ = 4; + check_less_or_equal(function__, "p_zero[sym1__, sym2__]", + rvalue(p_zero, + cons_list(index_uni(sym1__), + cons_list(index_uni(sym2__), + nil_index_list())), "p_zero"), 1);}} + current_statement__ = 5; + for (int sym1__ = 1; sym1__ <= N_samples; ++sym1__) { + current_statement__ = 5; + for (int sym2__ = 1; sym2__ <= N_bins; ++sym2__) { + current_statement__ = 5; + current_statement__ = 5; + check_greater_or_equal(function__, "p_one[sym1__, sym2__]", + rvalue(p_one, + cons_list(index_uni(sym1__), + cons_list(index_uni(sym2__), + nil_index_list())), "p_one"), 0);}} + current_statement__ = 5; + for (int sym1__ = 1; sym1__ <= N_samples; ++sym1__) { + current_statement__ = 5; + for (int sym2__ = 1; sym2__ <= N_bins; ++sym2__) { + current_statement__ = 5; + current_statement__ = 5; + check_less_or_equal(function__, "p_one[sym1__, sym2__]", + rvalue(p_one, + cons_list(index_uni(sym1__), + cons_list(index_uni(sym2__), + nil_index_list())), "p_one"), 1);}} + current_statement__ = 7; + for (int sym1__ = 1; sym1__ <= N_samples; ++sym1__) { + current_statement__ = 7; + for (int sym2__ = 1; sym2__ <= N_bins; ++sym2__) { + current_statement__ = 7; + current_statement__ = 7; + check_greater_or_equal(function__, "mu[sym1__, sym2__]", + rvalue(mu, + cons_list(index_uni(sym1__), + cons_list(index_uni(sym2__), + nil_index_list())), "mu"), 0);}} + current_statement__ = 7; + for (int sym1__ = 1; sym1__ <= N_samples; ++sym1__) { + current_statement__ = 7; + for (int sym2__ = 1; sym2__ <= N_bins; ++sym2__) { + current_statement__ = 7; + current_statement__ = 7; + check_less_or_equal(function__, "mu[sym1__, sym2__]", + rvalue(mu, + cons_list(index_uni(sym1__), + cons_list(index_uni(sym2__), + nil_index_list())), "mu"), 1);}} + if (emit_transformed_parameters__) { + vars__.emplace_back(phi); + for (int sym1__ = 1; sym1__ <= N_bins; ++sym1__) { + for (int sym2__ = 1; sym2__ <= N_samples; ++sym2__) { + vars__.emplace_back( + rvalue(p_zero, + cons_list(index_uni(sym2__), + cons_list(index_uni(sym1__), nil_index_list())), "p_zero")); + }} + for (int sym1__ = 1; sym1__ <= N_bins; ++sym1__) { + for (int sym2__ = 1; sym2__ <= N_samples; ++sym2__) { + vars__.emplace_back( + rvalue(p_one, + cons_list(index_uni(sym2__), + cons_list(index_uni(sym1__), nil_index_list())), "p_one")); + }} + for (int sym1__ = 1; sym1__ <= N_covar; ++sym1__) { + for (int sym2__ = 1; sym2__ <= N_bins; ++sym2__) { + vars__.emplace_back( + rvalue(beta, + cons_list(index_uni(sym2__), + cons_list(index_uni(sym1__), nil_index_list())), "beta"));} } - current_statement_begin__ = 60; - if (!(context__.contains_r("beta_raw"))) - stan::lang::rethrow_located(std::runtime_error(std::string("Variable beta_raw missing")), current_statement_begin__, prog_reader__()); - vals_r__ = context__.vals_r("beta_raw"); - pos__ = 0U; - validate_non_negative_index("beta_raw", "(N_bins - 1)", (N_bins - 1)); - validate_non_negative_index("beta_raw", "N_covar", N_covar); - context__.validate_dims("parameter initialization", "beta_raw", "matrix_d", context__.to_vec((N_bins - 1),N_covar)); - Eigen::Matrix beta_raw((N_bins - 1), N_covar); - size_t beta_raw_j_2_max__ = N_covar; - size_t beta_raw_j_1_max__ = (N_bins - 1); - for (size_t j_2__ = 0; j_2__ < beta_raw_j_2_max__; ++j_2__) { - for (size_t j_1__ = 0; j_1__ < beta_raw_j_1_max__; ++j_1__) { - beta_raw(j_1__, j_2__) = vals_r__[pos__++]; - } - } - try { - writer__.matrix_unconstrain(beta_raw); - } catch (const std::exception& e) { - stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable beta_raw: ") + e.what()), current_statement_begin__, prog_reader__()); - } - params_r__ = writer__.data_r(); - params_i__ = writer__.data_i(); + for (int sym1__ = 1; sym1__ <= N_bins; ++sym1__) { + for (int sym2__ = 1; sym2__ <= N_samples; ++sym2__) { + vars__.emplace_back( + rvalue(mu, + cons_list(index_uni(sym2__), + cons_list(index_uni(sym1__), nil_index_list())), "mu"));}} + } + if (logical_negation(emit_generated_quantities__)) { + return ; + } + double alpha_temp; + alpha_temp = std::numeric_limits::quiet_NaN(); + + double beta_temp; + beta_temp = std::numeric_limits::quiet_NaN(); + + std::vector> log_lik; + log_lik = std::vector>(N_samples, Eigen::Matrix(N_bins)); + stan::math::fill(log_lik, std::numeric_limits::quiet_NaN()); + + std::vector> ynew; + ynew = std::vector>(ynew_1dim__, Eigen::Matrix(ynew_2dim__)); + stan::math::fill(ynew, std::numeric_limits::quiet_NaN()); + + std::vector> newy_is_zero; + newy_is_zero = std::vector>(newy_is_zero_1dim__, std::vector(newy_is_zero_2dim__, std::numeric_limits::min())); + + std::vector> newy_is_one; + newy_is_one = std::vector>(newy_is_one_1dim__, std::vector(newy_is_one_2dim__, std::numeric_limits::min())); + + int newy_is_proportion; + newy_is_proportion = std::numeric_limits::min(); + + current_statement__ = 71; + for (int i = 1; i <= N_samples; ++i) { + current_statement__ = 69; + for (int j = 1; j <= N_bins; ++j) { + current_statement__ = 46; + assign(log_lik, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), 0, + "assigning variable log_lik"); + current_statement__ = 47; + assign(log_lik, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), + (log_lik[(i - 1)][(j - 1)] + + bernoulli_lpmf(is_zero[(i - 1)][(j - 1)], + rvalue(p_zero, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), "p_zero"))), + "assigning variable log_lik"); + current_statement__ = 52; + if (logical_eq(is_proportion[(i - 1)][(j - 1)], 1)) { + current_statement__ = 48; + alpha_temp = ((rvalue(mu, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), + "mu") * ESS[(i - 1)]) * phi); + current_statement__ = 49; + beta_temp = (((1 - + rvalue(mu, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), + "mu")) * ESS[(i - 1)]) * phi); + current_statement__ = 50; + assign(log_lik, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), + (log_lik[(i - 1)][(j - 1)] + + ((((stan::math::log( + ((1 - + rvalue(p_zero, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), + "p_zero")) - + rvalue(p_one, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), + "p_one"))) + + ((alpha_temp - 1) * + rvalue(logX, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), + "logX"))) + + ((beta_temp - 1) * + rvalue(logNX, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), + "logNX"))) - + (((alpha_temp + beta_temp) - 1) * + stan::math::log(ESS[(i - 1)]))) - + lbeta(alpha_temp, beta_temp))), + "assigning variable log_lik"); + } + current_statement__ = 67; + if (logical_eq(postpred, 1)) { + current_statement__ = 53; + assign(newy_is_zero, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), + bernoulli_rng( + rvalue(p_zero, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), "p_zero"), + base_rng__), "assigning variable newy_is_zero"); + current_statement__ = 55; + if (logical_eq(newy_is_zero[(i - 1)][(j - 1)], 1)) { + current_statement__ = 54; + assign(ynew, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), 0.0, + "assigning variable ynew"); + } + current_statement__ = 56; + assign(newy_is_one, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), + bernoulli_rng( + rvalue(p_one, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), "p_one"), + base_rng__), "assigning variable newy_is_one"); + current_statement__ = 58; + if (logical_eq(newy_is_one[(i - 1)][(j - 1)], 1)) { + current_statement__ = 57; + assign(ynew, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), 1.0, + "assigning variable ynew"); + } + current_statement__ = 59; + newy_is_proportion = (newy_is_zero[(i - 1)][(j - 1)] + + newy_is_one[(i - 1)][(j - 1)]); + current_statement__ = 65; + if (logical_eq(newy_is_proportion, 0)) { + current_statement__ = 60; + alpha_temp = (rvalue(mu, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), + "mu") * ESS[(i - 1)]); + current_statement__ = 61; + beta_temp = ((1 - + rvalue(mu, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), + "mu")) * ESS[(i - 1)]); + current_statement__ = 62; + assign(ynew, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), + beta_rng(alpha_temp, beta_temp, base_rng__), + "assigning variable ynew"); + current_statement__ = 63; + assign(ynew, + cons_list(index_uni(i), + cons_list(index_uni(j), nil_index_list())), + (ynew[(i - 1)][(j - 1)] * ESS[(i - 1)]), + "assigning variable ynew"); + } + } }} + vars__.emplace_back(alpha_temp); + vars__.emplace_back(beta_temp); + for (int sym1__ = 1; sym1__ <= N_bins; ++sym1__) { + for (int sym2__ = 1; sym2__ <= N_samples; ++sym2__) { + vars__.emplace_back(log_lik[(sym2__ - 1)][(sym1__ - 1)]);}} + for (int sym1__ = 1; sym1__ <= ynew_2dim__; ++sym1__) { + for (int sym2__ = 1; sym2__ <= ynew_1dim__; ++sym2__) { + vars__.emplace_back(ynew[(sym2__ - 1)][(sym1__ - 1)]);}} + for (int sym1__ = 1; sym1__ <= newy_is_zero_2dim__; ++sym1__) { + for (int sym2__ = 1; sym2__ <= newy_is_zero_1dim__; ++sym2__) { + vars__.emplace_back(newy_is_zero[(sym2__ - 1)][(sym1__ - 1)]);}} + for (int sym1__ = 1; sym1__ <= newy_is_one_2dim__; ++sym1__) { + for (int sym2__ = 1; sym2__ <= newy_is_one_1dim__; ++sym2__) { + vars__.emplace_back(newy_is_one[(sym2__ - 1)][(sym1__ - 1)]);}} + vars__.emplace_back(newy_is_proportion); + } catch (const std::exception& e) { + stan::lang::rethrow_located(e, locations_array__[current_statement__]); + // Next line prevents compiler griping about no return + throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***"); } - void transform_inits(const stan::io::var_context& context, - Eigen::Matrix& params_r, - std::ostream* pstream__) const { - std::vector params_r_vec; - std::vector params_i_vec; - transform_inits(context, params_i_vec, params_r_vec, pstream__); - params_r.resize(params_r_vec.size()); - for (int i = 0; i < params_r.size(); ++i) - params_r(i) = params_r_vec[i]; + } // write_array_impl() + + template * = nullptr, stan::require_vector_like_vt* = nullptr> + inline void transform_inits_impl(const stan::io::var_context& context__, + VecI& params_i__, VecVar& vars__, + std::ostream* pstream__ = nullptr) const { + using local_scalar_t__ = double; + vars__.clear(); + vars__.reserve(num_params_r__); + + try { + int pos__; + pos__ = std::numeric_limits::min(); + + pos__ = 1; + Eigen::Matrix phi_inv; + phi_inv = Eigen::Matrix(overdisp); + stan::math::fill(phi_inv, std::numeric_limits::quiet_NaN()); + + { + std::vector phi_inv_flat__; + current_statement__ = 1; + assign(phi_inv_flat__, nil_index_list(), context__.vals_r("phi_inv"), + "assigning variable phi_inv_flat__"); + current_statement__ = 1; + pos__ = 1; + current_statement__ = 1; + for (int sym1__ = 1; sym1__ <= overdisp; ++sym1__) { + current_statement__ = 1; + assign(phi_inv, cons_list(index_uni(sym1__), nil_index_list()), + phi_inv_flat__[(pos__ - 1)], "assigning variable phi_inv"); + current_statement__ = 1; + pos__ = (pos__ + 1);} + } + Eigen::Matrix phi_inv_free__; + phi_inv_free__ = Eigen::Matrix(overdisp); + stan::math::fill(phi_inv_free__, std::numeric_limits::quiet_NaN()); + + current_statement__ = 1; + for (int sym1__ = 1; sym1__ <= overdisp; ++sym1__) { + current_statement__ = 1; + assign(phi_inv_free__, + cons_list(index_uni(sym1__), nil_index_list()), + stan::math::lb_free(phi_inv[(sym1__ - 1)], 0), + "assigning variable phi_inv_free__");} + Eigen::Matrix beta_raw; + beta_raw = Eigen::Matrix(beta_raw_1dim__, N_covar); + stan::math::fill(beta_raw, std::numeric_limits::quiet_NaN()); + + { + std::vector beta_raw_flat__; + current_statement__ = 2; + assign(beta_raw_flat__, nil_index_list(), + context__.vals_r("beta_raw"), "assigning variable beta_raw_flat__"); + current_statement__ = 2; + pos__ = 1; + current_statement__ = 2; + for (int sym1__ = 1; sym1__ <= N_covar; ++sym1__) { + current_statement__ = 2; + for (int sym2__ = 1; sym2__ <= beta_raw_1dim__; ++sym2__) { + current_statement__ = 2; + assign(beta_raw, + cons_list(index_uni(sym2__), + cons_list(index_uni(sym1__), nil_index_list())), + beta_raw_flat__[(pos__ - 1)], "assigning variable beta_raw"); + current_statement__ = 2; + pos__ = (pos__ + 1);}} + } + for (int sym1__ = 1; sym1__ <= overdisp; ++sym1__) { + vars__.emplace_back(phi_inv_free__[(sym1__ - 1)]);} + for (int sym1__ = 1; sym1__ <= N_covar; ++sym1__) { + for (int sym2__ = 1; sym2__ <= beta_raw_1dim__; ++sym2__) { + vars__.emplace_back( + rvalue(beta_raw, + cons_list(index_uni(sym2__), + cons_list(index_uni(sym1__), nil_index_list())), "beta_raw")); + }} + } catch (const std::exception& e) { + stan::lang::rethrow_located(e, locations_array__[current_statement__]); + // Next line prevents compiler griping about no return + throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***"); } - template - T__ log_prob(std::vector& params_r__, - std::vector& params_i__, - std::ostream* pstream__ = 0) const { - typedef T__ local_scalar_t__; - local_scalar_t__ DUMMY_VAR__(std::numeric_limits::quiet_NaN()); - (void) DUMMY_VAR__; // dummy to suppress unused var warning - T__ lp__(0.0); - stan::math::accumulator lp_accum__; - try { - stan::io::reader in__(params_r__, params_i__); - // model parameters - current_statement_begin__ = 59; - Eigen::Matrix phi_inv; - (void) phi_inv; // dummy to suppress unused var warning - if (jacobian__) - phi_inv = in__.vector_lb_constrain(0, overdisp, lp__); - else - phi_inv = in__.vector_lb_constrain(0, overdisp); - current_statement_begin__ = 60; - Eigen::Matrix beta_raw; - (void) beta_raw; // dummy to suppress unused var warning - if (jacobian__) - beta_raw = in__.matrix_constrain((N_bins - 1), N_covar, lp__); - else - beta_raw = in__.matrix_constrain((N_bins - 1), N_covar); - // transformed parameters - current_statement_begin__ = 63; - local_scalar_t__ phi; - (void) phi; // dummy to suppress unused var warning - stan::math::initialize(phi, DUMMY_VAR__); - stan::math::fill(phi, DUMMY_VAR__); - current_statement_begin__ = 64; - validate_non_negative_index("p_zero", "N_samples", N_samples); - validate_non_negative_index("p_zero", "N_bins", N_bins); - Eigen::Matrix p_zero(N_samples, N_bins); - stan::math::initialize(p_zero, DUMMY_VAR__); - stan::math::fill(p_zero, DUMMY_VAR__); - current_statement_begin__ = 65; - validate_non_negative_index("p_one", "N_samples", N_samples); - validate_non_negative_index("p_one", "N_bins", N_bins); - Eigen::Matrix p_one(N_samples, N_bins); - stan::math::initialize(p_one, DUMMY_VAR__); - stan::math::fill(p_one, DUMMY_VAR__); - current_statement_begin__ = 66; - validate_non_negative_index("beta", "N_bins", N_bins); - validate_non_negative_index("beta", "N_covar", N_covar); - Eigen::Matrix beta(N_bins, N_covar); - stan::math::initialize(beta, DUMMY_VAR__); - stan::math::fill(beta, DUMMY_VAR__); - current_statement_begin__ = 67; - validate_non_negative_index("mu", "N_samples", N_samples); - validate_non_negative_index("mu", "N_bins", N_bins); - Eigen::Matrix mu(N_samples, N_bins); - stan::math::initialize(mu, DUMMY_VAR__); - stan::math::fill(mu, DUMMY_VAR__); - // transformed parameters block statements - current_statement_begin__ = 69; - stan::math::assign(phi, 1); - current_statement_begin__ = 70; - if (as_bool(logical_eq(overdisp, 1))) { - current_statement_begin__ = 70; - stan::math::assign(phi, (1 / get_base1(phi_inv, 1, "phi_inv", 1))); - } - current_statement_begin__ = 71; - for (int l = 1; l <= N_covar; ++l) { - current_statement_begin__ = 72; - stan::model::assign(beta, - stan::model::cons_list(stan::model::index_uni(N_bins), stan::model::cons_list(stan::model::index_uni(l), stan::model::nil_index_list())), - 0.0, - "assigning variable beta"); - } - current_statement_begin__ = 74; - for (int k = 1; k <= (N_bins - 1); ++k) { - current_statement_begin__ = 75; - for (int l = 1; l <= N_covar; ++l) { - current_statement_begin__ = 76; - stan::model::assign(beta, - stan::model::cons_list(stan::model::index_uni(k), stan::model::cons_list(stan::model::index_uni(l), stan::model::nil_index_list())), - get_base1(beta_raw, k, l, "beta_raw", 1), - "assigning variable beta"); - } - } - current_statement_begin__ = 81; - for (int n = 1; n <= N_samples; ++n) { - { - current_statement_begin__ = 82; - validate_non_negative_index("logits", "N_bins", N_bins); - Eigen::Matrix logits(N_bins); - stan::math::initialize(logits, DUMMY_VAR__); - stan::math::fill(logits, DUMMY_VAR__); - current_statement_begin__ = 83; - for (int m = 1; m <= N_bins; ++m) { - current_statement_begin__ = 84; - stan::model::assign(logits, - stan::model::cons_list(stan::model::index_uni(m), stan::model::nil_index_list()), - multiply(stan::model::rvalue(design_X, stan::model::cons_list(stan::model::index_uni(n), stan::model::cons_list(stan::model::index_omni(), stan::model::nil_index_list())), "design_X"), transpose(stan::model::rvalue(beta, stan::model::cons_list(stan::model::index_uni(m), stan::model::cons_list(stan::model::index_omni(), stan::model::nil_index_list())), "beta"))), - "assigning variable logits"); - } - current_statement_begin__ = 86; - stan::math::assign(logits, softmax(logits)); - current_statement_begin__ = 87; - for (int m = 1; m <= N_bins; ++m) { - current_statement_begin__ = 88; - stan::model::assign(mu, - stan::model::cons_list(stan::model::index_uni(n), stan::model::cons_list(stan::model::index_uni(m), stan::model::nil_index_list())), - get_base1(logits, m, "logits", 1), - "assigning variable mu"); - } - } - } - current_statement_begin__ = 92; - for (int i = 1; i <= N_samples; ++i) { - current_statement_begin__ = 93; - for (int j = 1; j <= N_bins; ++j) { - current_statement_begin__ = 94; - stan::model::assign(p_zero, - stan::model::cons_list(stan::model::index_uni(i), stan::model::cons_list(stan::model::index_uni(j), stan::model::nil_index_list())), - pow((1 - get_base1(mu, i, j, "mu", 1)), (get_base1(ESS, i, "ESS", 1) * phi)), - "assigning variable p_zero"); - } - current_statement_begin__ = 99; - for (int j = 1; j <= N_bins; ++j) { - current_statement_begin__ = 100; - stan::model::assign(p_one, - stan::model::cons_list(stan::model::index_uni(i), stan::model::cons_list(stan::model::index_uni(j), stan::model::nil_index_list())), - ((1 - get_base1(p_zero, i, j, "p_zero", 1)) * prod(stan::model::rvalue(p_zero, stan::model::cons_list(stan::model::index_uni(i), stan::model::cons_list(stan::model::index_multi(stan::model::rvalue(prod_idx, stan::model::cons_list(stan::model::index_uni(j), stan::model::cons_list(stan::model::index_omni(), stan::model::nil_index_list())), "prod_idx")), stan::model::nil_index_list())), "p_zero"))), - "assigning variable p_one"); - } - } - // validate transformed parameters - const char* function__ = "validate transformed params"; - (void) function__; // dummy to suppress unused var warning - current_statement_begin__ = 63; - if (stan::math::is_uninitialized(phi)) { - std::stringstream msg__; - msg__ << "Undefined transformed parameter: phi"; - stan::lang::rethrow_located(std::runtime_error(std::string("Error initializing variable phi: ") + msg__.str()), current_statement_begin__, prog_reader__()); - } - check_greater_or_equal(function__, "phi", phi, 0); - current_statement_begin__ = 64; - size_t p_zero_j_1_max__ = N_samples; - size_t p_zero_j_2_max__ = N_bins; - for (size_t j_1__ = 0; j_1__ < p_zero_j_1_max__; ++j_1__) { - for (size_t j_2__ = 0; j_2__ < p_zero_j_2_max__; ++j_2__) { - if (stan::math::is_uninitialized(p_zero(j_1__, j_2__))) { - std::stringstream msg__; - msg__ << "Undefined transformed parameter: p_zero" << "(" << j_1__ << ", " << j_2__ << ")"; - stan::lang::rethrow_located(std::runtime_error(std::string("Error initializing variable p_zero: ") + msg__.str()), current_statement_begin__, prog_reader__()); - } - } - } - check_greater_or_equal(function__, "p_zero", p_zero, 0); - check_less_or_equal(function__, "p_zero", p_zero, 1); - current_statement_begin__ = 65; - size_t p_one_j_1_max__ = N_samples; - size_t p_one_j_2_max__ = N_bins; - for (size_t j_1__ = 0; j_1__ < p_one_j_1_max__; ++j_1__) { - for (size_t j_2__ = 0; j_2__ < p_one_j_2_max__; ++j_2__) { - if (stan::math::is_uninitialized(p_one(j_1__, j_2__))) { - std::stringstream msg__; - msg__ << "Undefined transformed parameter: p_one" << "(" << j_1__ << ", " << j_2__ << ")"; - stan::lang::rethrow_located(std::runtime_error(std::string("Error initializing variable p_one: ") + msg__.str()), current_statement_begin__, prog_reader__()); - } - } - } - check_greater_or_equal(function__, "p_one", p_one, 0); - check_less_or_equal(function__, "p_one", p_one, 1); - current_statement_begin__ = 66; - size_t beta_j_1_max__ = N_bins; - size_t beta_j_2_max__ = N_covar; - for (size_t j_1__ = 0; j_1__ < beta_j_1_max__; ++j_1__) { - for (size_t j_2__ = 0; j_2__ < beta_j_2_max__; ++j_2__) { - if (stan::math::is_uninitialized(beta(j_1__, j_2__))) { - std::stringstream msg__; - msg__ << "Undefined transformed parameter: beta" << "(" << j_1__ << ", " << j_2__ << ")"; - stan::lang::rethrow_located(std::runtime_error(std::string("Error initializing variable beta: ") + msg__.str()), current_statement_begin__, prog_reader__()); - } - } - } - current_statement_begin__ = 67; - size_t mu_j_1_max__ = N_samples; - size_t mu_j_2_max__ = N_bins; - for (size_t j_1__ = 0; j_1__ < mu_j_1_max__; ++j_1__) { - for (size_t j_2__ = 0; j_2__ < mu_j_2_max__; ++j_2__) { - if (stan::math::is_uninitialized(mu(j_1__, j_2__))) { - std::stringstream msg__; - msg__ << "Undefined transformed parameter: mu" << "(" << j_1__ << ", " << j_2__ << ")"; - stan::lang::rethrow_located(std::runtime_error(std::string("Error initializing variable mu: ") + msg__.str()), current_statement_begin__, prog_reader__()); - } - } - } - check_greater_or_equal(function__, "mu", mu, 0); - check_less_or_equal(function__, "mu", mu, 1); - // model body + } // transform_inits_impl() + + inline void get_param_names(std::vector& names__) const { + + names__.clear(); + names__.emplace_back("phi_inv"); + names__.emplace_back("beta_raw"); + names__.emplace_back("phi"); + names__.emplace_back("p_zero"); + names__.emplace_back("p_one"); + names__.emplace_back("beta"); + names__.emplace_back("mu"); + names__.emplace_back("alpha_temp"); + names__.emplace_back("beta_temp"); + names__.emplace_back("log_lik"); + names__.emplace_back("ynew"); + names__.emplace_back("newy_is_zero"); + names__.emplace_back("newy_is_one"); + names__.emplace_back("newy_is_proportion"); + } // get_param_names() + + inline void get_dims(std::vector>& dimss__) const { + dimss__.clear(); + dimss__.emplace_back(std::vector{static_cast(overdisp)}); + + dimss__.emplace_back(std::vector{ + static_cast(beta_raw_1dim__) + , static_cast(N_covar)}); + + dimss__.emplace_back(std::vector{}); + + dimss__.emplace_back(std::vector{static_cast(N_samples), + static_cast(N_bins)}); + + dimss__.emplace_back(std::vector{static_cast(N_samples), + static_cast(N_bins)}); + + dimss__.emplace_back(std::vector{static_cast(N_bins), + static_cast(N_covar)}); + + dimss__.emplace_back(std::vector{static_cast(N_samples), + static_cast(N_bins)}); + + dimss__.emplace_back(std::vector{}); + + dimss__.emplace_back(std::vector{}); + + dimss__.emplace_back(std::vector{static_cast(N_samples), + static_cast(N_bins)}); + + dimss__.emplace_back(std::vector{static_cast(ynew_1dim__) + , + static_cast(ynew_2dim__)}); + + dimss__.emplace_back(std::vector{ + static_cast(newy_is_zero_1dim__) + , + static_cast(newy_is_zero_2dim__) + }); + + dimss__.emplace_back(std::vector{ + static_cast(newy_is_one_1dim__) + , + static_cast(newy_is_one_2dim__) + }); + + dimss__.emplace_back(std::vector{}); + + } // get_dims() + + inline void constrained_param_names( + std::vector& param_names__, + bool emit_transformed_parameters__ = true, + bool emit_generated_quantities__ = true) const + final { + + for (int sym1__ = 1; sym1__ <= overdisp; ++sym1__) { + { + param_names__.emplace_back(std::string() + "phi_inv" + '.' + std::to_string(sym1__)); + }} + for (int sym1__ = 1; sym1__ <= N_covar; ++sym1__) { + { + for (int sym2__ = 1; sym2__ <= beta_raw_1dim__; ++sym2__) { + { + param_names__.emplace_back(std::string() + "beta_raw" + '.' + std::to_string(sym2__) + '.' + std::to_string(sym1__)); + }} + }} + if (emit_transformed_parameters__) { + param_names__.emplace_back(std::string() + "phi"); + for (int sym1__ = 1; sym1__ <= N_bins; ++sym1__) { + { + for (int sym2__ = 1; sym2__ <= N_samples; ++sym2__) { { - current_statement_begin__ = 106; - local_scalar_t__ alpha_temp(DUMMY_VAR__); - (void) alpha_temp; // dummy to suppress unused var warning - stan::math::initialize(alpha_temp, DUMMY_VAR__); - stan::math::fill(alpha_temp, DUMMY_VAR__); - current_statement_begin__ = 107; - local_scalar_t__ beta_temp(DUMMY_VAR__); - (void) beta_temp; // dummy to suppress unused var warning - stan::math::initialize(beta_temp, DUMMY_VAR__); - stan::math::fill(beta_temp, DUMMY_VAR__); - current_statement_begin__ = 109; - if (as_bool(logical_eq(overdisp, 1))) { - current_statement_begin__ = 110; - lp_accum__.add(cauchy_log(phi_inv, 0, 5)); - } - current_statement_begin__ = 113; - for (int i = 1; i <= N_covar; ++i) { - current_statement_begin__ = 114; - for (int j = 1; j <= (N_bins - 1); ++j) { - current_statement_begin__ = 115; - lp_accum__.add(normal_log(get_base1(beta_raw, j, i, "beta_raw", 1), 0, prior_sd)); - } - } - current_statement_begin__ = 119; - for (int i = 1; i <= N_samples; ++i) { - current_statement_begin__ = 120; - for (int j = 1; j <= N_bins; ++j) { - current_statement_begin__ = 122; - lp_accum__.add(bernoulli_log(get_base1(get_base1(is_zero, i, "is_zero", 1), j, "is_zero", 2), get_base1(p_zero, i, j, "p_zero", 1))); - current_statement_begin__ = 124; - if (as_bool(logical_eq(get_base1(get_base1(is_proportion, i, "is_proportion", 1), j, "is_proportion", 2), 1))) { - current_statement_begin__ = 125; - stan::math::assign(alpha_temp, ((get_base1(mu, i, j, "mu", 1) * get_base1(ESS, i, "ESS", 1)) * phi)); - current_statement_begin__ = 126; - stan::math::assign(beta_temp, (((1 - get_base1(mu, i, j, "mu", 1)) * get_base1(ESS, i, "ESS", 1)) * phi)); - current_statement_begin__ = 128; - lp_accum__.add(((((stan::math::log(((1 - get_base1(p_zero, i, j, "p_zero", 1)) - get_base1(p_one, i, j, "p_one", 1))) + ((alpha_temp - 1) * get_base1(logX, i, j, "logX", 1))) + ((beta_temp - 1) * get_base1(logNX, i, j, "logNX", 1))) - (((alpha_temp + beta_temp) - 1) * stan::math::log(get_base1(ESS, i, "ESS", 1)))) - lbeta(alpha_temp, beta_temp))); - } - } - } - } - } catch (const std::exception& e) { - stan::lang::rethrow_located(e, current_statement_begin__, prog_reader__()); - // Next line prevents compiler griping about no return - throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***"); - } - lp_accum__.add(lp__); - return lp_accum__.sum(); - } // log_prob() - template - T_ log_prob(Eigen::Matrix& params_r, - std::ostream* pstream = 0) const { - std::vector vec_params_r; - vec_params_r.reserve(params_r.size()); - for (int i = 0; i < params_r.size(); ++i) - vec_params_r.push_back(params_r(i)); - std::vector vec_params_i; - return log_prob(vec_params_r, vec_params_i, pstream); + param_names__.emplace_back(std::string() + "p_zero" + '.' + std::to_string(sym2__) + '.' + std::to_string(sym1__)); + }} + }} + for (int sym1__ = 1; sym1__ <= N_bins; ++sym1__) { + { + for (int sym2__ = 1; sym2__ <= N_samples; ++sym2__) { + { + param_names__.emplace_back(std::string() + "p_one" + '.' + std::to_string(sym2__) + '.' + std::to_string(sym1__)); + }} + }} + for (int sym1__ = 1; sym1__ <= N_covar; ++sym1__) { + { + for (int sym2__ = 1; sym2__ <= N_bins; ++sym2__) { + { + param_names__.emplace_back(std::string() + "beta" + '.' + std::to_string(sym2__) + '.' + std::to_string(sym1__)); + }} + }} + for (int sym1__ = 1; sym1__ <= N_bins; ++sym1__) { + { + for (int sym2__ = 1; sym2__ <= N_samples; ++sym2__) { + { + param_names__.emplace_back(std::string() + "mu" + '.' + std::to_string(sym2__) + '.' + std::to_string(sym1__)); + }} + }} } - void get_param_names(std::vector& names__) const { - names__.resize(0); - names__.push_back("phi_inv"); - names__.push_back("beta_raw"); - names__.push_back("phi"); - names__.push_back("p_zero"); - names__.push_back("p_one"); - names__.push_back("beta"); - names__.push_back("mu"); - names__.push_back("alpha_temp"); - names__.push_back("beta_temp"); - names__.push_back("log_lik"); - names__.push_back("ynew"); - names__.push_back("newy_is_zero"); - names__.push_back("newy_is_one"); - names__.push_back("newy_is_proportion"); + + if (emit_generated_quantities__) { + param_names__.emplace_back(std::string() + "alpha_temp"); + param_names__.emplace_back(std::string() + "beta_temp"); + for (int sym1__ = 1; sym1__ <= N_bins; ++sym1__) { + { + for (int sym2__ = 1; sym2__ <= N_samples; ++sym2__) { + { + param_names__.emplace_back(std::string() + "log_lik" + '.' + std::to_string(sym2__) + '.' + std::to_string(sym1__)); + }} + }} + for (int sym1__ = 1; sym1__ <= ynew_2dim__; ++sym1__) { + { + for (int sym2__ = 1; sym2__ <= ynew_1dim__; ++sym2__) { + { + param_names__.emplace_back(std::string() + "ynew" + '.' + std::to_string(sym2__) + '.' + std::to_string(sym1__)); + }} + }} + for (int sym1__ = 1; sym1__ <= newy_is_zero_2dim__; ++sym1__) { + { + for (int sym2__ = 1; sym2__ <= newy_is_zero_1dim__; ++sym2__) { + { + param_names__.emplace_back(std::string() + "newy_is_zero" + '.' + std::to_string(sym2__) + '.' + std::to_string(sym1__)); + }} + }} + for (int sym1__ = 1; sym1__ <= newy_is_one_2dim__; ++sym1__) { + { + for (int sym2__ = 1; sym2__ <= newy_is_one_1dim__; ++sym2__) { + { + param_names__.emplace_back(std::string() + "newy_is_one" + '.' + std::to_string(sym2__) + '.' + std::to_string(sym1__)); + }} + }} + param_names__.emplace_back(std::string() + "newy_is_proportion"); } - void get_dims(std::vector >& dimss__) const { - dimss__.resize(0); - std::vector dims__; - dims__.resize(0); - dims__.push_back(overdisp); - dimss__.push_back(dims__); - dims__.resize(0); - dims__.push_back((N_bins - 1)); - dims__.push_back(N_covar); - dimss__.push_back(dims__); - dims__.resize(0); - dimss__.push_back(dims__); - dims__.resize(0); - dims__.push_back(N_samples); - dims__.push_back(N_bins); - dimss__.push_back(dims__); - dims__.resize(0); - dims__.push_back(N_samples); - dims__.push_back(N_bins); - dimss__.push_back(dims__); - dims__.resize(0); - dims__.push_back(N_bins); - dims__.push_back(N_covar); - dimss__.push_back(dims__); - dims__.resize(0); - dims__.push_back(N_samples); - dims__.push_back(N_bins); - dimss__.push_back(dims__); - dims__.resize(0); - dimss__.push_back(dims__); - dims__.resize(0); - dimss__.push_back(dims__); - dims__.resize(0); - dims__.push_back(N_samples); - dims__.push_back(N_bins); - dimss__.push_back(dims__); - dims__.resize(0); - dims__.push_back((N_samples * postpred)); - dims__.push_back((N_bins * postpred)); - dimss__.push_back(dims__); - dims__.resize(0); - dims__.push_back((N_samples * postpred)); - dims__.push_back((N_bins * postpred)); - dimss__.push_back(dims__); - dims__.resize(0); - dims__.push_back((N_samples * postpred)); - dims__.push_back((N_bins * postpred)); - dimss__.push_back(dims__); - dims__.resize(0); - dimss__.push_back(dims__); + + } // constrained_param_names() + + inline void unconstrained_param_names( + std::vector& param_names__, + bool emit_transformed_parameters__ = true, + bool emit_generated_quantities__ = true) const + final { + + for (int sym1__ = 1; sym1__ <= overdisp; ++sym1__) { + { + param_names__.emplace_back(std::string() + "phi_inv" + '.' + std::to_string(sym1__)); + }} + for (int sym1__ = 1; sym1__ <= N_covar; ++sym1__) { + { + for (int sym2__ = 1; sym2__ <= beta_raw_1dim__; ++sym2__) { + { + param_names__.emplace_back(std::string() + "beta_raw" + '.' + std::to_string(sym2__) + '.' + std::to_string(sym1__)); + }} + }} + if (emit_transformed_parameters__) { + param_names__.emplace_back(std::string() + "phi"); + for (int sym1__ = 1; sym1__ <= N_bins; ++sym1__) { + { + for (int sym2__ = 1; sym2__ <= N_samples; ++sym2__) { + { + param_names__.emplace_back(std::string() + "p_zero" + '.' + std::to_string(sym2__) + '.' + std::to_string(sym1__)); + }} + }} + for (int sym1__ = 1; sym1__ <= N_bins; ++sym1__) { + { + for (int sym2__ = 1; sym2__ <= N_samples; ++sym2__) { + { + param_names__.emplace_back(std::string() + "p_one" + '.' + std::to_string(sym2__) + '.' + std::to_string(sym1__)); + }} + }} + for (int sym1__ = 1; sym1__ <= N_covar; ++sym1__) { + { + for (int sym2__ = 1; sym2__ <= N_bins; ++sym2__) { + { + param_names__.emplace_back(std::string() + "beta" + '.' + std::to_string(sym2__) + '.' + std::to_string(sym1__)); + }} + }} + for (int sym1__ = 1; sym1__ <= N_bins; ++sym1__) { + { + for (int sym2__ = 1; sym2__ <= N_samples; ++sym2__) { + { + param_names__.emplace_back(std::string() + "mu" + '.' + std::to_string(sym2__) + '.' + std::to_string(sym1__)); + }} + }} } - template - void write_array(RNG& base_rng__, - std::vector& params_r__, - std::vector& params_i__, - std::vector& vars__, - bool include_tparams__ = true, - bool include_gqs__ = true, - std::ostream* pstream__ = 0) const { - typedef double local_scalar_t__; - vars__.resize(0); - stan::io::reader in__(params_r__, params_i__); - static const char* function__ = "model_dirichregmod_namespace::write_array"; - (void) function__; // dummy to suppress unused var warning - // read-transform, write parameters - Eigen::Matrix phi_inv = in__.vector_lb_constrain(0, overdisp); - size_t phi_inv_j_1_max__ = overdisp; - for (size_t j_1__ = 0; j_1__ < phi_inv_j_1_max__; ++j_1__) { - vars__.push_back(phi_inv(j_1__)); - } - Eigen::Matrix beta_raw = in__.matrix_constrain((N_bins - 1), N_covar); - size_t beta_raw_j_2_max__ = N_covar; - size_t beta_raw_j_1_max__ = (N_bins - 1); - for (size_t j_2__ = 0; j_2__ < beta_raw_j_2_max__; ++j_2__) { - for (size_t j_1__ = 0; j_1__ < beta_raw_j_1_max__; ++j_1__) { - vars__.push_back(beta_raw(j_1__, j_2__)); - } - } - double lp__ = 0.0; - (void) lp__; // dummy to suppress unused var warning - stan::math::accumulator lp_accum__; - local_scalar_t__ DUMMY_VAR__(std::numeric_limits::quiet_NaN()); - (void) DUMMY_VAR__; // suppress unused var warning - if (!include_tparams__ && !include_gqs__) return; - try { - // declare and define transformed parameters - current_statement_begin__ = 63; - double phi; - (void) phi; // dummy to suppress unused var warning - stan::math::initialize(phi, DUMMY_VAR__); - stan::math::fill(phi, DUMMY_VAR__); - current_statement_begin__ = 64; - validate_non_negative_index("p_zero", "N_samples", N_samples); - validate_non_negative_index("p_zero", "N_bins", N_bins); - Eigen::Matrix p_zero(N_samples, N_bins); - stan::math::initialize(p_zero, DUMMY_VAR__); - stan::math::fill(p_zero, DUMMY_VAR__); - current_statement_begin__ = 65; - validate_non_negative_index("p_one", "N_samples", N_samples); - validate_non_negative_index("p_one", "N_bins", N_bins); - Eigen::Matrix p_one(N_samples, N_bins); - stan::math::initialize(p_one, DUMMY_VAR__); - stan::math::fill(p_one, DUMMY_VAR__); - current_statement_begin__ = 66; - validate_non_negative_index("beta", "N_bins", N_bins); - validate_non_negative_index("beta", "N_covar", N_covar); - Eigen::Matrix beta(N_bins, N_covar); - stan::math::initialize(beta, DUMMY_VAR__); - stan::math::fill(beta, DUMMY_VAR__); - current_statement_begin__ = 67; - validate_non_negative_index("mu", "N_samples", N_samples); - validate_non_negative_index("mu", "N_bins", N_bins); - Eigen::Matrix mu(N_samples, N_bins); - stan::math::initialize(mu, DUMMY_VAR__); - stan::math::fill(mu, DUMMY_VAR__); - // do transformed parameters statements - current_statement_begin__ = 69; - stan::math::assign(phi, 1); - current_statement_begin__ = 70; - if (as_bool(logical_eq(overdisp, 1))) { - current_statement_begin__ = 70; - stan::math::assign(phi, (1 / get_base1(phi_inv, 1, "phi_inv", 1))); - } - current_statement_begin__ = 71; - for (int l = 1; l <= N_covar; ++l) { - current_statement_begin__ = 72; - stan::model::assign(beta, - stan::model::cons_list(stan::model::index_uni(N_bins), stan::model::cons_list(stan::model::index_uni(l), stan::model::nil_index_list())), - 0.0, - "assigning variable beta"); - } - current_statement_begin__ = 74; - for (int k = 1; k <= (N_bins - 1); ++k) { - current_statement_begin__ = 75; - for (int l = 1; l <= N_covar; ++l) { - current_statement_begin__ = 76; - stan::model::assign(beta, - stan::model::cons_list(stan::model::index_uni(k), stan::model::cons_list(stan::model::index_uni(l), stan::model::nil_index_list())), - get_base1(beta_raw, k, l, "beta_raw", 1), - "assigning variable beta"); - } - } - current_statement_begin__ = 81; - for (int n = 1; n <= N_samples; ++n) { - { - current_statement_begin__ = 82; - validate_non_negative_index("logits", "N_bins", N_bins); - Eigen::Matrix logits(N_bins); - stan::math::initialize(logits, DUMMY_VAR__); - stan::math::fill(logits, DUMMY_VAR__); - current_statement_begin__ = 83; - for (int m = 1; m <= N_bins; ++m) { - current_statement_begin__ = 84; - stan::model::assign(logits, - stan::model::cons_list(stan::model::index_uni(m), stan::model::nil_index_list()), - multiply(stan::model::rvalue(design_X, stan::model::cons_list(stan::model::index_uni(n), stan::model::cons_list(stan::model::index_omni(), stan::model::nil_index_list())), "design_X"), transpose(stan::model::rvalue(beta, stan::model::cons_list(stan::model::index_uni(m), stan::model::cons_list(stan::model::index_omni(), stan::model::nil_index_list())), "beta"))), - "assigning variable logits"); - } - current_statement_begin__ = 86; - stan::math::assign(logits, softmax(logits)); - current_statement_begin__ = 87; - for (int m = 1; m <= N_bins; ++m) { - current_statement_begin__ = 88; - stan::model::assign(mu, - stan::model::cons_list(stan::model::index_uni(n), stan::model::cons_list(stan::model::index_uni(m), stan::model::nil_index_list())), - get_base1(logits, m, "logits", 1), - "assigning variable mu"); - } - } - } - current_statement_begin__ = 92; - for (int i = 1; i <= N_samples; ++i) { - current_statement_begin__ = 93; - for (int j = 1; j <= N_bins; ++j) { - current_statement_begin__ = 94; - stan::model::assign(p_zero, - stan::model::cons_list(stan::model::index_uni(i), stan::model::cons_list(stan::model::index_uni(j), stan::model::nil_index_list())), - pow((1 - get_base1(mu, i, j, "mu", 1)), (get_base1(ESS, i, "ESS", 1) * phi)), - "assigning variable p_zero"); - } - current_statement_begin__ = 99; - for (int j = 1; j <= N_bins; ++j) { - current_statement_begin__ = 100; - stan::model::assign(p_one, - stan::model::cons_list(stan::model::index_uni(i), stan::model::cons_list(stan::model::index_uni(j), stan::model::nil_index_list())), - ((1 - get_base1(p_zero, i, j, "p_zero", 1)) * prod(stan::model::rvalue(p_zero, stan::model::cons_list(stan::model::index_uni(i), stan::model::cons_list(stan::model::index_multi(stan::model::rvalue(prod_idx, stan::model::cons_list(stan::model::index_uni(j), stan::model::cons_list(stan::model::index_omni(), stan::model::nil_index_list())), "prod_idx")), stan::model::nil_index_list())), "p_zero"))), - "assigning variable p_one"); - } - } - if (!include_gqs__ && !include_tparams__) return; - // validate transformed parameters - const char* function__ = "validate transformed params"; - (void) function__; // dummy to suppress unused var warning - current_statement_begin__ = 63; - check_greater_or_equal(function__, "phi", phi, 0); - current_statement_begin__ = 64; - check_greater_or_equal(function__, "p_zero", p_zero, 0); - check_less_or_equal(function__, "p_zero", p_zero, 1); - current_statement_begin__ = 65; - check_greater_or_equal(function__, "p_one", p_one, 0); - check_less_or_equal(function__, "p_one", p_one, 1); - current_statement_begin__ = 67; - check_greater_or_equal(function__, "mu", mu, 0); - check_less_or_equal(function__, "mu", mu, 1); - // write transformed parameters - if (include_tparams__) { - vars__.push_back(phi); - size_t p_zero_j_2_max__ = N_bins; - size_t p_zero_j_1_max__ = N_samples; - for (size_t j_2__ = 0; j_2__ < p_zero_j_2_max__; ++j_2__) { - for (size_t j_1__ = 0; j_1__ < p_zero_j_1_max__; ++j_1__) { - vars__.push_back(p_zero(j_1__, j_2__)); - } - } - size_t p_one_j_2_max__ = N_bins; - size_t p_one_j_1_max__ = N_samples; - for (size_t j_2__ = 0; j_2__ < p_one_j_2_max__; ++j_2__) { - for (size_t j_1__ = 0; j_1__ < p_one_j_1_max__; ++j_1__) { - vars__.push_back(p_one(j_1__, j_2__)); - } - } - size_t beta_j_2_max__ = N_covar; - size_t beta_j_1_max__ = N_bins; - for (size_t j_2__ = 0; j_2__ < beta_j_2_max__; ++j_2__) { - for (size_t j_1__ = 0; j_1__ < beta_j_1_max__; ++j_1__) { - vars__.push_back(beta(j_1__, j_2__)); - } - } - size_t mu_j_2_max__ = N_bins; - size_t mu_j_1_max__ = N_samples; - for (size_t j_2__ = 0; j_2__ < mu_j_2_max__; ++j_2__) { - for (size_t j_1__ = 0; j_1__ < mu_j_1_max__; ++j_1__) { - vars__.push_back(mu(j_1__, j_2__)); - } - } - } - if (!include_gqs__) return; - // declare and define generated quantities - current_statement_begin__ = 135; - double alpha_temp; - (void) alpha_temp; // dummy to suppress unused var warning - stan::math::initialize(alpha_temp, DUMMY_VAR__); - stan::math::fill(alpha_temp, DUMMY_VAR__); - current_statement_begin__ = 136; - double beta_temp; - (void) beta_temp; // dummy to suppress unused var warning - stan::math::initialize(beta_temp, DUMMY_VAR__); - stan::math::fill(beta_temp, DUMMY_VAR__); - current_statement_begin__ = 137; - validate_non_negative_index("log_lik", "N_bins", N_bins); - validate_non_negative_index("log_lik", "N_samples", N_samples); - std::vector > log_lik(N_samples, Eigen::Matrix(N_bins)); - stan::math::initialize(log_lik, DUMMY_VAR__); - stan::math::fill(log_lik, DUMMY_VAR__); - current_statement_begin__ = 138; - validate_non_negative_index("ynew", "(N_bins * postpred)", (N_bins * postpred)); - validate_non_negative_index("ynew", "(N_samples * postpred)", (N_samples * postpred)); - std::vector > ynew((N_samples * postpred), Eigen::Matrix((N_bins * postpred))); - stan::math::initialize(ynew, DUMMY_VAR__); - stan::math::fill(ynew, DUMMY_VAR__); - current_statement_begin__ = 139; - validate_non_negative_index("newy_is_zero", "(N_samples * postpred)", (N_samples * postpred)); - validate_non_negative_index("newy_is_zero", "(N_bins * postpred)", (N_bins * postpred)); - std::vector > newy_is_zero((N_samples * postpred), std::vector((N_bins * postpred), int(0))); - stan::math::fill(newy_is_zero, std::numeric_limits::min()); - current_statement_begin__ = 140; - validate_non_negative_index("newy_is_one", "(N_samples * postpred)", (N_samples * postpred)); - validate_non_negative_index("newy_is_one", "(N_bins * postpred)", (N_bins * postpred)); - std::vector > newy_is_one((N_samples * postpred), std::vector((N_bins * postpred), int(0))); - stan::math::fill(newy_is_one, std::numeric_limits::min()); - current_statement_begin__ = 141; - int newy_is_proportion; - (void) newy_is_proportion; // dummy to suppress unused var warning - stan::math::fill(newy_is_proportion, std::numeric_limits::min()); - // generated quantities statements - current_statement_begin__ = 143; - for (int i = 1; i <= N_samples; ++i) { - current_statement_begin__ = 144; - for (int j = 1; j <= N_bins; ++j) { - current_statement_begin__ = 145; - stan::model::assign(log_lik, - stan::model::cons_list(stan::model::index_uni(i), stan::model::cons_list(stan::model::index_uni(j), stan::model::nil_index_list())), - 0, - "assigning variable log_lik"); - current_statement_begin__ = 146; - stan::model::assign(log_lik, - stan::model::cons_list(stan::model::index_uni(i), stan::model::cons_list(stan::model::index_uni(j), stan::model::nil_index_list())), - (stan::model::rvalue(log_lik, stan::model::cons_list(stan::model::index_uni(i), stan::model::cons_list(stan::model::index_uni(j), stan::model::nil_index_list())), "log_lik") + bernoulli_log(get_base1(get_base1(is_zero, i, "is_zero", 1), j, "is_zero", 2), get_base1(p_zero, i, j, "p_zero", 1))), - "assigning variable log_lik"); - current_statement_begin__ = 147; - if (as_bool(logical_eq(get_base1(get_base1(is_proportion, i, "is_proportion", 1), j, "is_proportion", 2), 1))) { - current_statement_begin__ = 148; - stan::math::assign(alpha_temp, ((get_base1(mu, i, j, "mu", 1) * get_base1(ESS, i, "ESS", 1)) * phi)); - current_statement_begin__ = 149; - stan::math::assign(beta_temp, (((1 - get_base1(mu, i, j, "mu", 1)) * get_base1(ESS, i, "ESS", 1)) * phi)); - current_statement_begin__ = 150; - stan::model::assign(log_lik, - stan::model::cons_list(stan::model::index_uni(i), stan::model::cons_list(stan::model::index_uni(j), stan::model::nil_index_list())), - (stan::model::rvalue(log_lik, stan::model::cons_list(stan::model::index_uni(i), stan::model::cons_list(stan::model::index_uni(j), stan::model::nil_index_list())), "log_lik") + ((((stan::math::log(((1 - get_base1(p_zero, i, j, "p_zero", 1)) - get_base1(p_one, i, j, "p_one", 1))) + ((alpha_temp - 1) * get_base1(logX, i, j, "logX", 1))) + ((beta_temp - 1) * get_base1(logNX, i, j, "logNX", 1))) - (((alpha_temp + beta_temp) - 1) * stan::math::log(get_base1(ESS, i, "ESS", 1)))) - lbeta(alpha_temp, beta_temp))), - "assigning variable log_lik"); - } - current_statement_begin__ = 154; - if (as_bool(logical_eq(postpred, 1))) { - current_statement_begin__ = 155; - stan::model::assign(newy_is_zero, - stan::model::cons_list(stan::model::index_uni(i), stan::model::cons_list(stan::model::index_uni(j), stan::model::nil_index_list())), - bernoulli_rng(get_base1(p_zero, i, j, "p_zero", 1), base_rng__), - "assigning variable newy_is_zero"); - current_statement_begin__ = 156; - if (as_bool(logical_eq(get_base1(get_base1(newy_is_zero, i, "newy_is_zero", 1), j, "newy_is_zero", 2), 1))) { - current_statement_begin__ = 156; - stan::model::assign(ynew, - stan::model::cons_list(stan::model::index_uni(i), stan::model::cons_list(stan::model::index_uni(j), stan::model::nil_index_list())), - 0.0, - "assigning variable ynew"); - } - current_statement_begin__ = 157; - stan::model::assign(newy_is_one, - stan::model::cons_list(stan::model::index_uni(i), stan::model::cons_list(stan::model::index_uni(j), stan::model::nil_index_list())), - bernoulli_rng(get_base1(p_one, i, j, "p_one", 1), base_rng__), - "assigning variable newy_is_one"); - current_statement_begin__ = 158; - if (as_bool(logical_eq(get_base1(get_base1(newy_is_one, i, "newy_is_one", 1), j, "newy_is_one", 2), 1))) { - current_statement_begin__ = 158; - stan::model::assign(ynew, - stan::model::cons_list(stan::model::index_uni(i), stan::model::cons_list(stan::model::index_uni(j), stan::model::nil_index_list())), - 1.0, - "assigning variable ynew"); - } - current_statement_begin__ = 159; - stan::math::assign(newy_is_proportion, (get_base1(get_base1(newy_is_zero, i, "newy_is_zero", 1), j, "newy_is_zero", 2) + get_base1(get_base1(newy_is_one, i, "newy_is_one", 1), j, "newy_is_one", 2))); - current_statement_begin__ = 160; - if (as_bool(logical_eq(newy_is_proportion, 0))) { - current_statement_begin__ = 161; - stan::math::assign(alpha_temp, (get_base1(mu, i, j, "mu", 1) * get_base1(ESS, i, "ESS", 1))); - current_statement_begin__ = 162; - stan::math::assign(beta_temp, ((1 - get_base1(mu, i, j, "mu", 1)) * get_base1(ESS, i, "ESS", 1))); - current_statement_begin__ = 163; - stan::model::assign(ynew, - stan::model::cons_list(stan::model::index_uni(i), stan::model::cons_list(stan::model::index_uni(j), stan::model::nil_index_list())), - beta_rng(alpha_temp, beta_temp, base_rng__), - "assigning variable ynew"); - current_statement_begin__ = 164; - stan::model::assign(ynew, - stan::model::cons_list(stan::model::index_uni(i), stan::model::cons_list(stan::model::index_uni(j), stan::model::nil_index_list())), - (get_base1(get_base1(ynew, i, "ynew", 1), j, "ynew", 2) * get_base1(ESS, i, "ESS", 1)), - "assigning variable ynew"); - } - } - } - } - // validate, write generated quantities - current_statement_begin__ = 135; - vars__.push_back(alpha_temp); - current_statement_begin__ = 136; - vars__.push_back(beta_temp); - current_statement_begin__ = 137; - size_t log_lik_j_1_max__ = N_bins; - size_t log_lik_k_0_max__ = N_samples; - for (size_t j_1__ = 0; j_1__ < log_lik_j_1_max__; ++j_1__) { - for (size_t k_0__ = 0; k_0__ < log_lik_k_0_max__; ++k_0__) { - vars__.push_back(log_lik[k_0__](j_1__)); - } - } - current_statement_begin__ = 138; - size_t ynew_j_1_max__ = (N_bins * postpred); - size_t ynew_k_0_max__ = (N_samples * postpred); - for (size_t j_1__ = 0; j_1__ < ynew_j_1_max__; ++j_1__) { - for (size_t k_0__ = 0; k_0__ < ynew_k_0_max__; ++k_0__) { - vars__.push_back(ynew[k_0__](j_1__)); - } - } - current_statement_begin__ = 139; - size_t newy_is_zero_k_0_max__ = (N_samples * postpred); - size_t newy_is_zero_k_1_max__ = (N_bins * postpred); - for (size_t k_1__ = 0; k_1__ < newy_is_zero_k_1_max__; ++k_1__) { - for (size_t k_0__ = 0; k_0__ < newy_is_zero_k_0_max__; ++k_0__) { - vars__.push_back(newy_is_zero[k_0__][k_1__]); - } - } - current_statement_begin__ = 140; - size_t newy_is_one_k_0_max__ = (N_samples * postpred); - size_t newy_is_one_k_1_max__ = (N_bins * postpred); - for (size_t k_1__ = 0; k_1__ < newy_is_one_k_1_max__; ++k_1__) { - for (size_t k_0__ = 0; k_0__ < newy_is_one_k_0_max__; ++k_0__) { - vars__.push_back(newy_is_one[k_0__][k_1__]); - } - } - current_statement_begin__ = 141; - vars__.push_back(newy_is_proportion); - } catch (const std::exception& e) { - stan::lang::rethrow_located(e, current_statement_begin__, prog_reader__()); - // Next line prevents compiler griping about no return - throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***"); - } + + if (emit_generated_quantities__) { + param_names__.emplace_back(std::string() + "alpha_temp"); + param_names__.emplace_back(std::string() + "beta_temp"); + for (int sym1__ = 1; sym1__ <= N_bins; ++sym1__) { + { + for (int sym2__ = 1; sym2__ <= N_samples; ++sym2__) { + { + param_names__.emplace_back(std::string() + "log_lik" + '.' + std::to_string(sym2__) + '.' + std::to_string(sym1__)); + }} + }} + for (int sym1__ = 1; sym1__ <= ynew_2dim__; ++sym1__) { + { + for (int sym2__ = 1; sym2__ <= ynew_1dim__; ++sym2__) { + { + param_names__.emplace_back(std::string() + "ynew" + '.' + std::to_string(sym2__) + '.' + std::to_string(sym1__)); + }} + }} + for (int sym1__ = 1; sym1__ <= newy_is_zero_2dim__; ++sym1__) { + { + for (int sym2__ = 1; sym2__ <= newy_is_zero_1dim__; ++sym2__) { + { + param_names__.emplace_back(std::string() + "newy_is_zero" + '.' + std::to_string(sym2__) + '.' + std::to_string(sym1__)); + }} + }} + for (int sym1__ = 1; sym1__ <= newy_is_one_2dim__; ++sym1__) { + { + for (int sym2__ = 1; sym2__ <= newy_is_one_1dim__; ++sym2__) { + { + param_names__.emplace_back(std::string() + "newy_is_one" + '.' + std::to_string(sym2__) + '.' + std::to_string(sym1__)); + }} + }} + param_names__.emplace_back(std::string() + "newy_is_proportion"); } + + } // unconstrained_param_names() + + inline std::string get_constrained_sizedtypes() const { + stringstream s__; + s__ << "[{\"name\":\"phi_inv\",\"type\":{\"name\":\"vector\",\"length\":" << overdisp << "},\"block\":\"parameters\"},{\"name\":\"beta_raw\",\"type\":{\"name\":\"matrix\",\"rows\":" << beta_raw_1dim__ << ",\"cols\":" << N_covar << "},\"block\":\"parameters\"},{\"name\":\"phi\",\"type\":{\"name\":\"real\"},\"block\":\"transformed_parameters\"},{\"name\":\"p_zero\",\"type\":{\"name\":\"matrix\",\"rows\":" << N_samples << ",\"cols\":" << N_bins << "},\"block\":\"transformed_parameters\"},{\"name\":\"p_one\",\"type\":{\"name\":\"matrix\",\"rows\":" << N_samples << ",\"cols\":" << N_bins << "},\"block\":\"transformed_parameters\"},{\"name\":\"beta\",\"type\":{\"name\":\"matrix\",\"rows\":" << N_bins << ",\"cols\":" << N_covar << "},\"block\":\"transformed_parameters\"},{\"name\":\"mu\",\"type\":{\"name\":\"matrix\",\"rows\":" << N_samples << ",\"cols\":" << N_bins << "},\"block\":\"transformed_parameters\"},{\"name\":\"alpha_temp\",\"type\":{\"name\":\"real\"},\"block\":\"generated_quantities\"},{\"name\":\"beta_temp\",\"type\":{\"name\":\"real\"},\"block\":\"generated_quantities\"},{\"name\":\"log_lik\",\"type\":{\"name\":\"array\",\"length\":" << N_samples << ",\"element_type\":{\"name\":\"vector\",\"length\":" << N_bins << "}},\"block\":\"generated_quantities\"},{\"name\":\"ynew\",\"type\":{\"name\":\"array\",\"length\":" << ynew_1dim__ << ",\"element_type\":{\"name\":\"vector\",\"length\":" << ynew_2dim__ << "}},\"block\":\"generated_quantities\"},{\"name\":\"newy_is_zero\",\"type\":{\"name\":\"array\",\"length\":" << newy_is_zero_1dim__ << ",\"element_type\":{\"name\":\"array\",\"length\":" << newy_is_zero_2dim__ << ",\"element_type\":{\"name\":\"int\"}}},\"block\":\"generated_quantities\"},{\"name\":\"newy_is_one\",\"type\":{\"name\":\"array\",\"length\":" << newy_is_one_1dim__ << ",\"element_type\":{\"name\":\"array\",\"length\":" << newy_is_one_2dim__ << ",\"element_type\":{\"name\":\"int\"}}},\"block\":\"generated_quantities\"},{\"name\":\"newy_is_proportion\",\"type\":{\"name\":\"int\"},\"block\":\"generated_quantities\"}]"; + return s__.str(); + } // get_constrained_sizedtypes() + + inline std::string get_unconstrained_sizedtypes() const { + stringstream s__; + s__ << "[{\"name\":\"phi_inv\",\"type\":{\"name\":\"vector\",\"length\":" << overdisp << "},\"block\":\"parameters\"},{\"name\":\"beta_raw\",\"type\":{\"name\":\"matrix\",\"rows\":" << beta_raw_1dim__ << ",\"cols\":" << N_covar << "},\"block\":\"parameters\"},{\"name\":\"phi\",\"type\":{\"name\":\"real\"},\"block\":\"transformed_parameters\"},{\"name\":\"p_zero\",\"type\":{\"name\":\"matrix\",\"rows\":" << N_samples << ",\"cols\":" << N_bins << "},\"block\":\"transformed_parameters\"},{\"name\":\"p_one\",\"type\":{\"name\":\"matrix\",\"rows\":" << N_samples << ",\"cols\":" << N_bins << "},\"block\":\"transformed_parameters\"},{\"name\":\"beta\",\"type\":{\"name\":\"matrix\",\"rows\":" << N_bins << ",\"cols\":" << N_covar << "},\"block\":\"transformed_parameters\"},{\"name\":\"mu\",\"type\":{\"name\":\"matrix\",\"rows\":" << N_samples << ",\"cols\":" << N_bins << "},\"block\":\"transformed_parameters\"},{\"name\":\"alpha_temp\",\"type\":{\"name\":\"real\"},\"block\":\"generated_quantities\"},{\"name\":\"beta_temp\",\"type\":{\"name\":\"real\"},\"block\":\"generated_quantities\"},{\"name\":\"log_lik\",\"type\":{\"name\":\"array\",\"length\":" << N_samples << ",\"element_type\":{\"name\":\"vector\",\"length\":" << N_bins << "}},\"block\":\"generated_quantities\"},{\"name\":\"ynew\",\"type\":{\"name\":\"array\",\"length\":" << ynew_1dim__ << ",\"element_type\":{\"name\":\"vector\",\"length\":" << ynew_2dim__ << "}},\"block\":\"generated_quantities\"},{\"name\":\"newy_is_zero\",\"type\":{\"name\":\"array\",\"length\":" << newy_is_zero_1dim__ << ",\"element_type\":{\"name\":\"array\",\"length\":" << newy_is_zero_2dim__ << ",\"element_type\":{\"name\":\"int\"}}},\"block\":\"generated_quantities\"},{\"name\":\"newy_is_one\",\"type\":{\"name\":\"array\",\"length\":" << newy_is_one_1dim__ << ",\"element_type\":{\"name\":\"array\",\"length\":" << newy_is_one_2dim__ << ",\"element_type\":{\"name\":\"int\"}}},\"block\":\"generated_quantities\"},{\"name\":\"newy_is_proportion\",\"type\":{\"name\":\"int\"},\"block\":\"generated_quantities\"}]"; + return s__.str(); + } // get_unconstrained_sizedtypes() + + + // Begin method overload boilerplate template - void write_array(RNG& base_rng, - Eigen::Matrix& params_r, - Eigen::Matrix& vars, - bool include_tparams = true, - bool include_gqs = true, - std::ostream* pstream = 0) const { - std::vector params_r_vec(params_r.size()); - for (int i = 0; i < params_r.size(); ++i) - params_r_vec[i] = params_r(i); - std::vector vars_vec; - std::vector params_i_vec; - write_array(base_rng, params_r_vec, params_i_vec, vars_vec, include_tparams, include_gqs, pstream); + inline void write_array(RNG& base_rng, + Eigen::Matrix& params_r, + Eigen::Matrix& vars, + const bool emit_transformed_parameters = true, + const bool emit_generated_quantities = true, + std::ostream* pstream = nullptr) const { + std::vector vars_vec(vars.size()); + std::vector params_i; + write_array_impl(base_rng, params_r, params_i, vars_vec, + emit_transformed_parameters, emit_generated_quantities, pstream); vars.resize(vars_vec.size()); - for (int i = 0; i < vars.size(); ++i) - vars(i) = vars_vec[i]; + for (int i = 0; i < vars.size(); ++i) { + vars.coeffRef(i) = vars_vec[i]; + } } - std::string model_name() const { - return "model_dirichregmod"; + template + inline void write_array(RNG& base_rng, std::vector& params_r, + std::vector& params_i, + std::vector& vars, + bool emit_transformed_parameters = true, + bool emit_generated_quantities = true, + std::ostream* pstream = nullptr) const { + write_array_impl(base_rng, params_r, params_i, vars, emit_transformed_parameters, emit_generated_quantities, pstream); } - void constrained_param_names(std::vector& param_names__, - bool include_tparams__ = true, - bool include_gqs__ = true) const { - std::stringstream param_name_stream__; - size_t phi_inv_j_1_max__ = overdisp; - for (size_t j_1__ = 0; j_1__ < phi_inv_j_1_max__; ++j_1__) { - param_name_stream__.str(std::string()); - param_name_stream__ << "phi_inv" << '.' << j_1__ + 1; - param_names__.push_back(param_name_stream__.str()); - } - size_t beta_raw_j_2_max__ = N_covar; - size_t beta_raw_j_1_max__ = (N_bins - 1); - for (size_t j_2__ = 0; j_2__ < beta_raw_j_2_max__; ++j_2__) { - for (size_t j_1__ = 0; j_1__ < beta_raw_j_1_max__; ++j_1__) { - param_name_stream__.str(std::string()); - param_name_stream__ << "beta_raw" << '.' << j_1__ + 1 << '.' << j_2__ + 1; - param_names__.push_back(param_name_stream__.str()); - } - } - if (!include_gqs__ && !include_tparams__) return; - if (include_tparams__) { - param_name_stream__.str(std::string()); - param_name_stream__ << "phi"; - param_names__.push_back(param_name_stream__.str()); - size_t p_zero_j_2_max__ = N_bins; - size_t p_zero_j_1_max__ = N_samples; - for (size_t j_2__ = 0; j_2__ < p_zero_j_2_max__; ++j_2__) { - for (size_t j_1__ = 0; j_1__ < p_zero_j_1_max__; ++j_1__) { - param_name_stream__.str(std::string()); - param_name_stream__ << "p_zero" << '.' << j_1__ + 1 << '.' << j_2__ + 1; - param_names__.push_back(param_name_stream__.str()); - } - } - size_t p_one_j_2_max__ = N_bins; - size_t p_one_j_1_max__ = N_samples; - for (size_t j_2__ = 0; j_2__ < p_one_j_2_max__; ++j_2__) { - for (size_t j_1__ = 0; j_1__ < p_one_j_1_max__; ++j_1__) { - param_name_stream__.str(std::string()); - param_name_stream__ << "p_one" << '.' << j_1__ + 1 << '.' << j_2__ + 1; - param_names__.push_back(param_name_stream__.str()); - } - } - size_t beta_j_2_max__ = N_covar; - size_t beta_j_1_max__ = N_bins; - for (size_t j_2__ = 0; j_2__ < beta_j_2_max__; ++j_2__) { - for (size_t j_1__ = 0; j_1__ < beta_j_1_max__; ++j_1__) { - param_name_stream__.str(std::string()); - param_name_stream__ << "beta" << '.' << j_1__ + 1 << '.' << j_2__ + 1; - param_names__.push_back(param_name_stream__.str()); - } - } - size_t mu_j_2_max__ = N_bins; - size_t mu_j_1_max__ = N_samples; - for (size_t j_2__ = 0; j_2__ < mu_j_2_max__; ++j_2__) { - for (size_t j_1__ = 0; j_1__ < mu_j_1_max__; ++j_1__) { - param_name_stream__.str(std::string()); - param_name_stream__ << "mu" << '.' << j_1__ + 1 << '.' << j_2__ + 1; - param_names__.push_back(param_name_stream__.str()); - } - } - } - if (!include_gqs__) return; - param_name_stream__.str(std::string()); - param_name_stream__ << "alpha_temp"; - param_names__.push_back(param_name_stream__.str()); - param_name_stream__.str(std::string()); - param_name_stream__ << "beta_temp"; - param_names__.push_back(param_name_stream__.str()); - size_t log_lik_j_1_max__ = N_bins; - size_t log_lik_k_0_max__ = N_samples; - for (size_t j_1__ = 0; j_1__ < log_lik_j_1_max__; ++j_1__) { - for (size_t k_0__ = 0; k_0__ < log_lik_k_0_max__; ++k_0__) { - param_name_stream__.str(std::string()); - param_name_stream__ << "log_lik" << '.' << k_0__ + 1 << '.' << j_1__ + 1; - param_names__.push_back(param_name_stream__.str()); - } - } - size_t ynew_j_1_max__ = (N_bins * postpred); - size_t ynew_k_0_max__ = (N_samples * postpred); - for (size_t j_1__ = 0; j_1__ < ynew_j_1_max__; ++j_1__) { - for (size_t k_0__ = 0; k_0__ < ynew_k_0_max__; ++k_0__) { - param_name_stream__.str(std::string()); - param_name_stream__ << "ynew" << '.' << k_0__ + 1 << '.' << j_1__ + 1; - param_names__.push_back(param_name_stream__.str()); - } - } - size_t newy_is_zero_k_0_max__ = (N_samples * postpred); - size_t newy_is_zero_k_1_max__ = (N_bins * postpred); - for (size_t k_1__ = 0; k_1__ < newy_is_zero_k_1_max__; ++k_1__) { - for (size_t k_0__ = 0; k_0__ < newy_is_zero_k_0_max__; ++k_0__) { - param_name_stream__.str(std::string()); - param_name_stream__ << "newy_is_zero" << '.' << k_0__ + 1 << '.' << k_1__ + 1; - param_names__.push_back(param_name_stream__.str()); - } - } - size_t newy_is_one_k_0_max__ = (N_samples * postpred); - size_t newy_is_one_k_1_max__ = (N_bins * postpred); - for (size_t k_1__ = 0; k_1__ < newy_is_one_k_1_max__; ++k_1__) { - for (size_t k_0__ = 0; k_0__ < newy_is_one_k_0_max__; ++k_0__) { - param_name_stream__.str(std::string()); - param_name_stream__ << "newy_is_one" << '.' << k_0__ + 1 << '.' << k_1__ + 1; - param_names__.push_back(param_name_stream__.str()); - } - } - param_name_stream__.str(std::string()); - param_name_stream__ << "newy_is_proportion"; - param_names__.push_back(param_name_stream__.str()); + template + inline T_ log_prob(Eigen::Matrix& params_r, + std::ostream* pstream = nullptr) const { + Eigen::Matrix params_i; + return log_prob_impl(params_r, params_i, pstream); } - void unconstrained_param_names(std::vector& param_names__, - bool include_tparams__ = true, - bool include_gqs__ = true) const { - std::stringstream param_name_stream__; - size_t phi_inv_j_1_max__ = overdisp; - for (size_t j_1__ = 0; j_1__ < phi_inv_j_1_max__; ++j_1__) { - param_name_stream__.str(std::string()); - param_name_stream__ << "phi_inv" << '.' << j_1__ + 1; - param_names__.push_back(param_name_stream__.str()); - } - size_t beta_raw_j_2_max__ = N_covar; - size_t beta_raw_j_1_max__ = (N_bins - 1); - for (size_t j_2__ = 0; j_2__ < beta_raw_j_2_max__; ++j_2__) { - for (size_t j_1__ = 0; j_1__ < beta_raw_j_1_max__; ++j_1__) { - param_name_stream__.str(std::string()); - param_name_stream__ << "beta_raw" << '.' << j_1__ + 1 << '.' << j_2__ + 1; - param_names__.push_back(param_name_stream__.str()); - } - } - if (!include_gqs__ && !include_tparams__) return; - if (include_tparams__) { - param_name_stream__.str(std::string()); - param_name_stream__ << "phi"; - param_names__.push_back(param_name_stream__.str()); - size_t p_zero_j_2_max__ = N_bins; - size_t p_zero_j_1_max__ = N_samples; - for (size_t j_2__ = 0; j_2__ < p_zero_j_2_max__; ++j_2__) { - for (size_t j_1__ = 0; j_1__ < p_zero_j_1_max__; ++j_1__) { - param_name_stream__.str(std::string()); - param_name_stream__ << "p_zero" << '.' << j_1__ + 1 << '.' << j_2__ + 1; - param_names__.push_back(param_name_stream__.str()); - } - } - size_t p_one_j_2_max__ = N_bins; - size_t p_one_j_1_max__ = N_samples; - for (size_t j_2__ = 0; j_2__ < p_one_j_2_max__; ++j_2__) { - for (size_t j_1__ = 0; j_1__ < p_one_j_1_max__; ++j_1__) { - param_name_stream__.str(std::string()); - param_name_stream__ << "p_one" << '.' << j_1__ + 1 << '.' << j_2__ + 1; - param_names__.push_back(param_name_stream__.str()); - } - } - size_t beta_j_2_max__ = N_covar; - size_t beta_j_1_max__ = N_bins; - for (size_t j_2__ = 0; j_2__ < beta_j_2_max__; ++j_2__) { - for (size_t j_1__ = 0; j_1__ < beta_j_1_max__; ++j_1__) { - param_name_stream__.str(std::string()); - param_name_stream__ << "beta" << '.' << j_1__ + 1 << '.' << j_2__ + 1; - param_names__.push_back(param_name_stream__.str()); - } - } - size_t mu_j_2_max__ = N_bins; - size_t mu_j_1_max__ = N_samples; - for (size_t j_2__ = 0; j_2__ < mu_j_2_max__; ++j_2__) { - for (size_t j_1__ = 0; j_1__ < mu_j_1_max__; ++j_1__) { - param_name_stream__.str(std::string()); - param_name_stream__ << "mu" << '.' << j_1__ + 1 << '.' << j_2__ + 1; - param_names__.push_back(param_name_stream__.str()); - } - } - } - if (!include_gqs__) return; - param_name_stream__.str(std::string()); - param_name_stream__ << "alpha_temp"; - param_names__.push_back(param_name_stream__.str()); - param_name_stream__.str(std::string()); - param_name_stream__ << "beta_temp"; - param_names__.push_back(param_name_stream__.str()); - size_t log_lik_j_1_max__ = N_bins; - size_t log_lik_k_0_max__ = N_samples; - for (size_t j_1__ = 0; j_1__ < log_lik_j_1_max__; ++j_1__) { - for (size_t k_0__ = 0; k_0__ < log_lik_k_0_max__; ++k_0__) { - param_name_stream__.str(std::string()); - param_name_stream__ << "log_lik" << '.' << k_0__ + 1 << '.' << j_1__ + 1; - param_names__.push_back(param_name_stream__.str()); - } - } - size_t ynew_j_1_max__ = (N_bins * postpred); - size_t ynew_k_0_max__ = (N_samples * postpred); - for (size_t j_1__ = 0; j_1__ < ynew_j_1_max__; ++j_1__) { - for (size_t k_0__ = 0; k_0__ < ynew_k_0_max__; ++k_0__) { - param_name_stream__.str(std::string()); - param_name_stream__ << "ynew" << '.' << k_0__ + 1 << '.' << j_1__ + 1; - param_names__.push_back(param_name_stream__.str()); - } - } - size_t newy_is_zero_k_0_max__ = (N_samples * postpred); - size_t newy_is_zero_k_1_max__ = (N_bins * postpred); - for (size_t k_1__ = 0; k_1__ < newy_is_zero_k_1_max__; ++k_1__) { - for (size_t k_0__ = 0; k_0__ < newy_is_zero_k_0_max__; ++k_0__) { - param_name_stream__.str(std::string()); - param_name_stream__ << "newy_is_zero" << '.' << k_0__ + 1 << '.' << k_1__ + 1; - param_names__.push_back(param_name_stream__.str()); - } - } - size_t newy_is_one_k_0_max__ = (N_samples * postpred); - size_t newy_is_one_k_1_max__ = (N_bins * postpred); - for (size_t k_1__ = 0; k_1__ < newy_is_one_k_1_max__; ++k_1__) { - for (size_t k_0__ = 0; k_0__ < newy_is_one_k_0_max__; ++k_0__) { - param_name_stream__.str(std::string()); - param_name_stream__ << "newy_is_one" << '.' << k_0__ + 1 << '.' << k_1__ + 1; - param_names__.push_back(param_name_stream__.str()); - } - } - param_name_stream__.str(std::string()); - param_name_stream__ << "newy_is_proportion"; - param_names__.push_back(param_name_stream__.str()); + template + inline T__ log_prob(std::vector& params_r, + std::vector& params_i, + std::ostream* pstream = nullptr) const { + return log_prob_impl(params_r, params_i, pstream); } -}; // model -} // namespace -typedef model_dirichregmod_namespace::model_dirichregmod stan_model; + + inline void transform_inits(const stan::io::var_context& context, + Eigen::Matrix& params_r, + std::ostream* pstream = nullptr) const final { + std::vector params_r_vec(params_r.size()); + std::vector params_i; + transform_inits_impl(context, params_i, params_r_vec, pstream); + params_r.resize(params_r_vec.size()); + for (int i = 0; i < params_r.size(); ++i) { + params_r.coeffRef(i) = params_r_vec[i]; + } + } + inline void transform_inits(const stan::io::var_context& context, + std::vector& params_i, + std::vector& vars, + std::ostream* pstream = nullptr) const final { + transform_inits_impl(context, params_i, vars, pstream); + } +}; +} +using stan_model = model_dirichregmod_namespace::model_dirichregmod; #ifndef USING_R +// Boilerplate stan::model::model_base& new_model( stan::io::var_context& data_context, unsigned int seed, @@ -1258,5 +1784,8 @@ stan::model::model_base& new_model( stan_model* m = new stan_model(data_context, seed, msg_stream); return *m; } +stan::math::profile_map& get_stan_profile_data() { + return model_dirichregmod_namespace::profiles__; +} #endif #endif diff --git a/src/stanExports_dirichregmod.o b/src/stanExports_dirichregmod.o index 6a3ab75..4bfc2bb 100644 Binary files a/src/stanExports_dirichregmod.o and b/src/stanExports_dirichregmod.o differ diff --git a/src/zoid.so b/src/zoid.so index 31a381d..3384d00 100755 Binary files a/src/zoid.so and b/src/zoid.so differ