From 1374757553c87814cd368f6f04fa3761f6b02292 Mon Sep 17 00:00:00 2001 From: MatthijsV Date: Thu, 19 Oct 2017 12:04:30 +0100 Subject: [PATCH 01/14] Normal Linear GLM working --- stan/math/prim/mat.hpp | 1 + .../prim/mat/prob/normal_linear_glm_lpdf.hpp | 159 ++++++++++++++++++ 2 files changed, 160 insertions(+) create mode 100644 stan/math/prim/mat/prob/normal_linear_glm_lpdf.hpp diff --git a/stan/math/prim/mat.hpp b/stan/math/prim/mat.hpp index 0d80084eb8c..0e515d91179 100644 --- a/stan/math/prim/mat.hpp +++ b/stan/math/prim/mat.hpp @@ -288,6 +288,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/prim/mat/prob/normal_linear_glm_lpdf.hpp b/stan/math/prim/mat/prob/normal_linear_glm_lpdf.hpp new file mode 100644 index 00000000000..6ad86678b2a --- /dev/null +++ b/stan/math/prim/mat/prob/normal_linear_glm_lpdf.hpp @@ -0,0 +1,159 @@ +#ifndef STAN_MATH_PRIM_MAT_PROB_NORMAL_LINEAR_GLM_LPDF_HPP +#define STAN_MATH_PRIM_MAT_PROB_NORMAL_LINEAR_GLM_LPDF_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { + namespace math { + + /** + * Returns the log PDF of the Generalized Linear Model (GLM) + * with Normal distribution and linear link function. + * If containers are supplied, returns the log sum of the probabilities. + * @tparam T_n type of real vector of dependent variables (labels); + * this can also be a single real value; + * @tparam T_x type of the matrix of independent variables (features); this + * should be an Eigen::Matrix type whose number of rows should match the + * length of n and whose number of columns should match the length of beta + * @tparam T_beta type of the weight vector; + * this can also be a single double value; + * @tparam T_alpha type of the intercept; + * this can either be a vector of doubles of a single double value; + * @tparam T_scale type of the scale vector. + * @param n real vector parameter + * @param x design matrix + * @param beta weight vector + * @param alpha intercept (in log odds) + * @param sigma (Sequence of) scale parameters for the normal + * distribution. + * @return log probability or log sum of probabilities + * @throw std::domain_error if x, beta or alpha is infinite. + * @throw std::invalid_argument if container sizes mismatch. + * @throw std::domain_error if the scale is not positive. + */ + template + typename return_type::type + normal_linear_glm_lpdf(const T_n &n, const T_x &x, const T_beta &beta, + const T_alpha &alpha, const T_scale& sigma) { + static const std::string function = "normal_linear_glm_lpdf"; + typedef typename stan::partials_return_type::type + T_partials_return; + + using std::exp; + using Eigen::Dynamic; + using Eigen::Matrix; + using Eigen::Array; + + if (!(stan::length(n) && stan::length(x) && stan::length(beta) + && stan::length(sigma))) + return 0.0; + + T_partials_return logp(0.0); + + check_not_nan(function, "Vector of dependent variables", n); + check_not_nan(function, "Matrix of independent variables", x); + check_not_nan(function, "Weight vector", beta); + check_not_nan(function, "Intercept", alpha); + check_positive(function, "Scale vector", sigma); + check_consistent_sizes(function, + "Rows in matrix of independent variables", + x.col(0), "Vector of dependent variables", n); + check_consistent_sizes(function, + "Columns in matrix of independent variables", + x.row(0), "Weight vector", beta); + + if (!include_summand::value) + return 0.0; + + const size_t N = x.col(0).size(); + const size_t M = x.row(0).size(); + + Array inv_sigma(N, 1); + Array log_sigma(N, 1); + Array n_dbl(N, 1); + { + scalar_seq_view n_vec(n); + scalar_seq_view sigma_vec(sigma); + for (size_t n = 0; n < N; ++n) { + inv_sigma[n] = 1/value_of(sigma_vec[n]); + log_sigma[n] = log(value_of(sigma_vec[n])); + n_dbl[n] = n_vec[n]; + } + } + Matrix beta_dbl(M, 1); + { + scalar_seq_view beta_vec(beta); + for (size_t m = 0; m < M; ++m) { + beta_dbl[m] = value_of(beta_vec[m]); + } + } + Matrix x_dbl = value_of(x); + + Array mu_dbl = (x_dbl * beta_dbl + Matrix::Ones(N, 1) + * value_of(alpha)).array(); + Array n_minus_mu_over_sigma = + (n_dbl - mu_dbl) * inv_sigma; + Array n_minus_mu_over_sigma_squared + = n_minus_mu_over_sigma.square(); + + if (include_summand::value) + logp += NEG_LOG_SQRT_TWO_PI * N; + if (include_summand::value) + logp -= log_sigma.sum(); + if (include_summand::value) + logp -= 0.5 * n_minus_mu_over_sigma_squared.sum(); + + // Compute the necessary derivatives. + operands_and_partials ops_partials(x, beta, alpha, sigma); + if (!(is_constant_struct::value && is_constant_struct::value + && is_constant_struct::value)) { + Matrix mu_derivative = (inv_sigma + * n_minus_mu_over_sigma).matrix(); + if (!is_constant_struct::value) { + ops_partials.edge2_.partials_ = x_dbl.transpose() * mu_derivative; + } + if (!is_constant_struct::value) { + ops_partials.edge1_.partials_ = mu_derivative + * beta_dbl.transpose(); + } + if (!is_constant_struct::value) { + ops_partials.edge3_.partials_[0] = mu_derivative.trace(); + } + if (!is_constant_struct::value) { + ops_partials.edge4_.partials_ = ((inv_sigma - Array::Ones(N, 1)) * n_minus_mu_over_sigma_squared).matrix(); + } + } + + return ops_partials.build(logp); + } + + template + inline + typename return_type::type + normal_linear_glm_lpdf(const T_n &n, const T_x &x, const T_beta &beta, + const T_alpha &alpha, const T_scale& sigma) { + return normal_linear_glm_lpdf(n, x, beta, alpha, sigma); + } + } // namespace math +} // namespace stan +#endif From fd6518201d4d3a5b893328ed8055006c4a457eb0 Mon Sep 17 00:00:00 2001 From: MatthijsV Date: Thu, 19 Oct 2017 12:06:17 +0100 Subject: [PATCH 02/14] Added Test, which passes --- .../mat/prob/normal_linear_glm_lpdf_test.cpp | 159 ++++++++++++++++++ 1 file changed, 159 insertions(+) create mode 100644 test/unit/math/rev/mat/prob/normal_linear_glm_lpdf_test.cpp diff --git a/test/unit/math/rev/mat/prob/normal_linear_glm_lpdf_test.cpp b/test/unit/math/rev/mat/prob/normal_linear_glm_lpdf_test.cpp new file mode 100644 index 00000000000..311beff6d4b --- /dev/null +++ b/test/unit/math/rev/mat/prob/normal_linear_glm_lpdf_test.cpp @@ -0,0 +1,159 @@ +#include +#include +#include +#include + +using stan::math::var; +using Eigen::Dynamic; +using Eigen::Matrix; + +typedef std::chrono::high_resolution_clock::time_point TimeVar; +#define duration(a) std::chrono::duration_cast(a).count() +#define timeNow() std::chrono::high_resolution_clock::now() + +// We check that the values of the new regression match those of one built +// from existing primitives. +TEST(ProbDistributionsNormalLinearGLM, glm_matches_normal_linear_doubles) +{ + Matrix n(3, 1); + n << 51, 32, 12; + Matrix x(3, 2); + x << -12, 46, -42, + 24, 25, 27; + Matrix beta(2, 1); + beta << 0.3, 2; + double alpha = 0.3; + Matrix alphavec = alpha * Matrix::Ones(); + Matrix theta(3, 1); + theta = x * beta + alphavec; + Matrix sigma(3, 1); + sigma << 10, 4, 6; + + EXPECT_FLOAT_EQ((stan::math::normal_lpdf(n, theta, sigma)), + (stan::math::normal_linear_glm_lpdf(n, x, beta, alpha, + sigma))); + EXPECT_FLOAT_EQ((stan::math::normal_lpdf(n, theta, sigma)), + (stan::math::normal_linear_glm_lpdf(n, x, beta, alpha, + sigma))); + EXPECT_FLOAT_EQ((stan::math::normal_lpdf(n, theta, sigma)), + (stan::math::normal_linear_glm_lpdf(n, x, beta, alpha, + sigma))); + EXPECT_FLOAT_EQ((stan::math::normal_lpdf>(n, + theta, sigma)), + (stan::math::normal_linear_glm_lpdf>(n, x, beta, alpha, sigma))); + EXPECT_FLOAT_EQ((stan::math::normal_lpdf>(n, theta, sigma)), + (stan::math::normal_linear_glm_lpdf>(n, x, beta, alpha, sigma))); + EXPECT_FLOAT_EQ((stan::math::normal_lpdf>(n, theta, sigma)), + (stan::math::normal_linear_glm_lpdf>(n, x, beta, alpha, sigma))); +} + +// We check that the gradients of the new regression match those of one built +// from existing primitives. +TEST(ProbDistributionsNormalLinearGLM, glm_matches_normal_linear_vars) +{ + + Matrix n(3, 1); + n << 14, 32, 21; + Matrix x(3, 2); + x << -12, 46, -42, + 24, 25, 27; + Matrix beta(2, 1); + beta << 0.3, 2; + var alpha = 0.3; + Matrix sigma(3, 1); + sigma << 10, 4, 6; + Matrix alphavec = alpha * Matrix::Ones(); + Matrix theta(3, 1); + theta = x * beta + alphavec; + + var lp = stan::math::normal_lpdf(n, theta, sigma); + lp.grad(); + + stan::math::recover_memory(); + + Matrix n2(3, 1); + n2 << 14, 32, 21; + Matrix x2(3, 2); + x2 << -12, 46, -42, + 24, 25, 27; + Matrix beta2(2, 1); + beta2 << 0.3, 2; + var alpha2 = 0.3; + Matrix sigma2(3, 1); + sigma2 << 10, 4, 6; + + var lp2 = stan::math::normal_linear_glm_lpdf(n2, x2, beta2, alpha2, sigma2); + lp2.grad(); + + EXPECT_FLOAT_EQ(lp.val(), + lp2.val()); + for (size_t i = 0; i < 2; i++) + { + EXPECT_FLOAT_EQ(beta[i].adj(), beta2[i].adj()); + } + EXPECT_FLOAT_EQ(alpha.adj(), alpha2.adj()); + for (size_t j = 0; j < 3; j++) + { + EXPECT_FLOAT_EQ(sigma[j].adj(), sigma2[j].adj()); + for (size_t i = 0; i < 2; i++) + { + EXPECT_FLOAT_EQ(x(j, i).adj(), x2(j, i).adj()); + } + } +} + +// Here, we compare the speed of the new regression to that of one built from +// existing primitives. + +/* +TEST(ProbDistributionsNormalLinearGLM, glm_matches_normal_linear_speed) { + const int R = 30000; + const int C = 1000; + + Matrix n(R, 1); + for (size_t i = 0; i < R; i++) { + n[i] = rand()%2000; + } + + int T1 = 0; + int T2 = 0; + + for (size_t testnumber = 0; testnumber < 30; testnumber++){ + Matrix xreal = Matrix::Random(R, C); + Matrix betareal = Matrix::Random(C, 1); + Matrix sigmareal = Matrix::Random(R, 1) + + Matrix::Ones(); // Random Matrix has entries between -1 and 1. We add 1 to it to get positive entries. + Matrix alphareal = Matrix::Random(1, 1); + Matrix alpharealvec = Matrix::Ones() * alphareal; + + Matrix beta = betareal; + Matrix sigma = sigmareal; + Matrix theta(R, 1); + + + TimeVar t1 = timeNow(); + theta = (xreal * beta) + alpharealvec; + var lp = stan::math::normal_lpdf(n, theta, sigma); + + lp.grad(); + TimeVar t2 = timeNow(); + + stan::math::recover_memory(); + + Matrix beta2 = betareal; + Matrix sigma2 = sigmareal; + + TimeVar t3 = timeNow(); + var lp2 = stan::math::normal_linear_glm_lpdf(n, xreal, beta2, alphareal[0], + sigma2); + lp2.grad(); + TimeVar t4 = timeNow(); + stan::math::recover_memory(); + T1 += duration(t2 - t1); + T2 += duration(t4 - t3); + + } + + std::cout << "Existing Primitives:" << std::endl << T1 << std::endl << "New Primitives:" << std::endl << T2 << std::endl; +} +*/ From de35b3bbe6b47ba3673a3176573590938fcf9d9d Mon Sep 17 00:00:00 2001 From: MatthijsV Date: Thu, 19 Oct 2017 12:08:59 +0100 Subject: [PATCH 03/14] Removed poissonglm from mat.hpp, to help merging --- stan/math/prim/mat.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/stan/math/prim/mat.hpp b/stan/math/prim/mat.hpp index 0e515d91179..84b20e1f8a1 100644 --- a/stan/math/prim/mat.hpp +++ b/stan/math/prim/mat.hpp @@ -292,7 +292,6 @@ #include #include #include -#include #include #include #include From 19b0b417538291cea39d467efd23b8123b47e288 Mon Sep 17 00:00:00 2001 From: MatthijsV Date: Thu, 19 Oct 2017 12:11:15 +0100 Subject: [PATCH 04/14] Removed poisson regression files to help merging --- .../prim/mat/prob/poisson_log_glm_lpmf.hpp | 159 ------------------ .../mat/prob/poisson_log_glm_lpmf_test.cpp | 144 ---------------- 2 files changed, 303 deletions(-) delete mode 100644 stan/math/prim/mat/prob/poisson_log_glm_lpmf.hpp delete mode 100644 test/unit/math/rev/mat/prob/poisson_log_glm_lpmf_test.cpp diff --git a/stan/math/prim/mat/prob/poisson_log_glm_lpmf.hpp b/stan/math/prim/mat/prob/poisson_log_glm_lpmf.hpp deleted file mode 100644 index e10837e7bc9..00000000000 --- a/stan/math/prim/mat/prob/poisson_log_glm_lpmf.hpp +++ /dev/null @@ -1,159 +0,0 @@ -#ifndef STAN_MATH_PRIM_MAT_PROB_POISSON_LOG_GLM_LPMF_HPP -#define STAN_MATH_PRIM_MAT_PROB_POISSON_LOG_GLM_LPMF_HPP - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -namespace stan { - namespace math { - - /** - * Returns the log PMF of the Generalized Linear Model (GLM) - * with Poisson distribution and log link function. - * If containers are supplied, returns the log sum of the probabilities. - * @tparam T_n type of vector of dependent variables (labels), integers >=0; - * this can also be a single positive integer; - * @tparam T_x type of the matrix of independent variables (features); this - * should be an Eigen::Matrix type whose number of rows should match the - * length of n and whose number of columns should match the length of beta - * @tparam T_beta type of the weight vector; - * this can also be a single double value; - * @tparam T_alpha type of the intercept; - * this can either be a vector of doubles of a single double value; - * @param n positive integer vector parameter - * @param x design matrix - * @param beta weight vector - * @param alpha intercept (in log odds) - * @return log probability or log sum of probabilities - * @throw std::domain_error if theta is infinite. - * @throw std::invalid_argument if container sizes mismatch. - */ - template - typename return_type::type - poisson_log_glm_lpmf(const T_n &n, const T_x &x, const T_beta &beta, - const T_alpha &alpha) { - static const std::string function = "poisson_log_glm_lpmf"; - typedef typename stan::partials_return_type::type - T_partials_return; - - using std::exp; - using Eigen::Dynamic; - using Eigen::Matrix; - - if (!(stan::length(n) && stan::length(x) && stan::length(beta))) - return 0.0; - - T_partials_return logp(0.0); - - check_nonnegative(function, "Vector of dependent variables", n); - check_not_nan(function, "Matrix of independent variables", x); - check_not_nan(function, "Weight vector", beta); - check_not_nan(function, "Intercept", alpha); - check_consistent_sizes(function, - "Rows in matrix of independent variables", - x.col(0), "Vector of dependent variables", n); - check_consistent_sizes(function, - "Columns in matrix of independent variables", - x.row(0), "Weight vector", beta); - - if (!include_summand::value) - return 0.0; - - const size_t N = x.col(0).size(); - const size_t M = x.row(0).size(); - - /* - for (size_t i = 0; i < N; i++) - if ((std::numeric_limits::infinity() == theta_vec[i]) || - (-std::numeric_limits::infinity() == theta_vec[i] - && n_vec[i] != 0)) - return LOG_ZERO; - */ - // Should we do something like this? - - Matrix n_vec(N, 1); - { - scalar_seq_view n_seq_view(n); - for (size_t n = 0; n < N; ++n) { - n_vec[n] = n_seq_view[n]; - } - } - - Matrix beta_dbl(M, 1); - { - scalar_seq_view beta_vec(beta); - for (size_t m = 0; m < M; ++m) { - beta_dbl[m] = value_of(beta_vec[m]); - } - } - Matrix x_dbl = value_of(x); - - Matrix theta_dbl = (x_dbl * beta_dbl - + Matrix::Ones(N, 1) - * value_of(alpha)); - Matrix exp_theta = - theta_dbl.array().exp().matrix(); - - for (size_t i = 0; i < N; i++) { - // Compute the log-density and handle extreme values gracefully - // using Taylor approximations. - if (!(theta_dbl[i] == -std::numeric_limits::infinity() - && n_vec[i] == 0)) { - if (include_summand::value) - logp -= lgamma(n_vec[i] + 1.0); - if (include_summand::value) - logp += n_vec[i] * theta_dbl[i] - exp_theta[i]; - } - } - - // Compute the necessary derivatives. - operands_and_partials ops_partials(x, beta, alpha); - if (!(is_constant_struct::value && is_constant_struct::value - && is_constant_struct::value)) { - Matrix theta_derivative = - n_vec - exp_theta; - if (!is_constant_struct::value) { - ops_partials.edge2_.partials_ = x_dbl.transpose() * theta_derivative; - } - if (!is_constant_struct::value) { - ops_partials.edge1_.partials_ = theta_derivative - * beta_dbl.transpose(); - } - if (!is_constant_struct::value) { - ops_partials.edge3_.partials_[0] = theta_derivative.trace(); - } - } - return ops_partials.build(logp); - } - - template - inline - typename return_type::type - poisson_log_glm_lpmf(const T_n &n, const T_x &x, const T_beta &beta, - const T_alpha &alpha) { - return poisson_log_glm_lpmf(n, x, beta, alpha); - } - } // namespace math -} // namespace stan -#endif diff --git a/test/unit/math/rev/mat/prob/poisson_log_glm_lpmf_test.cpp b/test/unit/math/rev/mat/prob/poisson_log_glm_lpmf_test.cpp deleted file mode 100644 index b9afa49a1be..00000000000 --- a/test/unit/math/rev/mat/prob/poisson_log_glm_lpmf_test.cpp +++ /dev/null @@ -1,144 +0,0 @@ -#include -#include -#include -#include - -using stan::math::var; -using Eigen::Dynamic; -using Eigen::Matrix; - -typedef std::chrono::high_resolution_clock::time_point TimeVar; -#define duration(a) std::chrono::duration_cast(a).count() -#define timeNow() std::chrono::high_resolution_clock::now() - -// We check that the values of the new regression match those of one built -// from existing primitives. -TEST(ProbDistributionsPoissonLogGLM, glm_matches_poisson_log_doubles) -{ - Matrix n(3, 1); - n << 15, 3, 5; - Matrix x(3, 2); - x << -12, 46, -42, - 24, 25, 27; - Matrix beta(2, 1); - beta << 0.3, 2; - double alpha = 0.3; - Matrix alphavec = alpha * Matrix::Ones(); - Matrix theta(3, 1); - theta = x * beta + alphavec; - - EXPECT_FLOAT_EQ((stan::math::poisson_log_lpmf(n, theta)), - (stan::math::poisson_log_glm_lpmf(n, x, beta, alpha))); - EXPECT_FLOAT_EQ((stan::math::poisson_log_lpmf(n, theta)), - (stan::math::poisson_log_glm_lpmf(n, x, beta, alpha))); - EXPECT_FLOAT_EQ((stan::math::poisson_log_lpmf(n, theta)), - (stan::math::poisson_log_glm_lpmf(n, x, beta, alpha))); - EXPECT_FLOAT_EQ((stan::math::poisson_log_lpmf>(n, theta)), - (stan::math::poisson_log_glm_lpmf>(n, x, beta, alpha))); - EXPECT_FLOAT_EQ((stan::math::poisson_log_lpmf>(n, theta)), - (stan::math::poisson_log_glm_lpmf>(n, x, beta, alpha))); - EXPECT_FLOAT_EQ((stan::math::poisson_log_lpmf>(n, theta)), - (stan::math::poisson_log_glm_lpmf>(n, x, beta, alpha))); -} - -// We check that the gradients of the new regression match those of one built -// from existing primitives. - -TEST(ProbDistributionsPoissonLogGLM, glm_matches_poisson_log_vars) -{ - - Matrix n(3, 1); - n << 14, 2, 5; - Matrix x(3, 2); - x << -12, 46, -42, - 24, 25, 27; - Matrix beta(2, 1); - beta << 0.3, 2; - var alpha = 0.3; - Matrix alphavec = alpha * Matrix::Ones(); - Matrix theta(3, 1); - theta = x * beta + alphavec; - - var lp = stan::math::poisson_log_lpmf(n, theta); - lp.grad(); - - stan::math::recover_memory(); - - Matrix n2(3, 1); - n2 << 14, 2, 5; - Matrix x2(3, 2); - x2 << -12, 46, -42, - 24, 25, 27; - Matrix beta2(2, 1); - beta2 << 0.3, 2; - var alpha2 = 0.3; - - var lp2 = stan::math::poisson_log_glm_lpmf(n2, x2, beta2, alpha2); - lp2.grad(); - - EXPECT_FLOAT_EQ(lp.val(), - lp2.val()); - for (size_t i = 0; i < 2; i++) - { - EXPECT_FLOAT_EQ(beta[i].adj(), beta2[i].adj()); - } - EXPECT_FLOAT_EQ(alpha.adj(), alpha2.adj()); - for (size_t j = 0; j < 3; j++) - { - for (size_t i = 0; i < 2; i++) - { - EXPECT_FLOAT_EQ(x(j, i).adj(), x2(j, i).adj()); - } - } -} - - -// Here, we compare the speed of the new regression to that of one built from -// existing primitives. -/* -TEST(ProbDistributionsPoissonLogGLM, glm_matches_poisson_log_speed) { - const int R = 30000; - const int C = 1000; - - Matrix n(R, 1); - for (size_t i = 0; i < R; i++) { - n[i] = rand()%200; - } - - int T1 = 0; - int T2 = 0; - - for (size_t testnumber = 0; testnumber < 30; testnumber++){ - Matrix xreal = Matrix::Random(R, C); - Matrix betareal = Matrix::Random(C, 1); - Matrix alphareal = Matrix::Random(1, 1); - Matrix alpharealvec = Matrix::Ones() * alphareal; - - Matrix beta = betareal; - Matrix theta(R, 1); - - - TimeVar t1 = timeNow(); - theta = (xreal * beta) + alpharealvec; - var lp = stan::math::poisson_log_lpmf(n, theta); - - lp.grad(); - TimeVar t2 = timeNow(); - - stan::math::recover_memory(); - - Matrix beta2 = betareal; - - TimeVar t3 = timeNow(); - var lp2 = stan::math::poisson_log_glm_lpmf(n, xreal, beta2, alphareal[0]); - lp2.grad(); - TimeVar t4 = timeNow(); - stan::math::recover_memory(); - T1 += duration(t2 - t1); - T2 += duration(t4 - t3); - - } - - std::cout << "Existing Primitives:" << std::endl << T1 << std::endl << "New Primitives:" << std::endl << T2 << std::endl; -} -*/ From d9acb80acd7dcb698bde0d999cdecc1c7b3c94c7 Mon Sep 17 00:00:00 2001 From: MatthijsV Date: Thu, 19 Oct 2017 13:29:40 +0100 Subject: [PATCH 05/14] Satisfied cpplint + did performance tests - 4.5 performance gain over hand-built linear regression --- .../prim/mat/prob/normal_linear_glm_lpdf.hpp | 19 ++++++++++--------- .../mat/prob/normal_linear_glm_lpdf_test.cpp | 4 ++-- 2 files changed, 12 insertions(+), 11 deletions(-) diff --git a/stan/math/prim/mat/prob/normal_linear_glm_lpdf.hpp b/stan/math/prim/mat/prob/normal_linear_glm_lpdf.hpp index 6ad86678b2a..e55bd6b7bbc 100644 --- a/stan/math/prim/mat/prob/normal_linear_glm_lpdf.hpp +++ b/stan/math/prim/mat/prob/normal_linear_glm_lpdf.hpp @@ -86,18 +86,17 @@ namespace stan { const size_t N = x.col(0).size(); const size_t M = x.row(0).size(); - Array inv_sigma(N, 1); - Array log_sigma(N, 1); - Array n_dbl(N, 1); + Array sigma_dbl(N, 1); + Array n_dbl(N, 1); { scalar_seq_view n_vec(n); scalar_seq_view sigma_vec(sigma); for (size_t n = 0; n < N; ++n) { - inv_sigma[n] = 1/value_of(sigma_vec[n]); - log_sigma[n] = log(value_of(sigma_vec[n])); + sigma_dbl[n] = value_of(sigma_vec[n]); n_dbl[n] = n_vec[n]; } } + Array inv_sigma = 1/sigma_dbl; Matrix beta_dbl(M, 1); { scalar_seq_view beta_vec(beta); @@ -107,7 +106,8 @@ namespace stan { } Matrix x_dbl = value_of(x); - Array mu_dbl = (x_dbl * beta_dbl + Matrix::Ones(N, 1) + Array mu_dbl = (x_dbl * beta_dbl + + Matrix::Ones(N, 1) * value_of(alpha)).array(); Array n_minus_mu_over_sigma = (n_dbl - mu_dbl) * inv_sigma; @@ -117,12 +117,13 @@ namespace stan { if (include_summand::value) logp += NEG_LOG_SQRT_TWO_PI * N; if (include_summand::value) - logp -= log_sigma.sum(); + logp -= sigma_dbl.log().sum(); if (include_summand::value) logp -= 0.5 * n_minus_mu_over_sigma_squared.sum(); - + // Compute the necessary derivatives. - operands_and_partials ops_partials(x, beta, alpha, sigma); + operands_and_partials ops_partials(x, beta, + alpha, sigma); if (!(is_constant_struct::value && is_constant_struct::value && is_constant_struct::value)) { Matrix mu_derivative = (inv_sigma diff --git a/test/unit/math/rev/mat/prob/normal_linear_glm_lpdf_test.cpp b/test/unit/math/rev/mat/prob/normal_linear_glm_lpdf_test.cpp index 311beff6d4b..95bf2fdabb1 100644 --- a/test/unit/math/rev/mat/prob/normal_linear_glm_lpdf_test.cpp +++ b/test/unit/math/rev/mat/prob/normal_linear_glm_lpdf_test.cpp @@ -105,7 +105,7 @@ TEST(ProbDistributionsNormalLinearGLM, glm_matches_normal_linear_vars) // Here, we compare the speed of the new regression to that of one built from // existing primitives. -/* + TEST(ProbDistributionsNormalLinearGLM, glm_matches_normal_linear_speed) { const int R = 30000; const int C = 1000; @@ -156,4 +156,4 @@ TEST(ProbDistributionsNormalLinearGLM, glm_matches_normal_linear_speed) { std::cout << "Existing Primitives:" << std::endl << T1 << std::endl << "New Primitives:" << std::endl << T2 << std::endl; } -*/ + From 5d155f26029ba4c577a32ee6283085c1475ca829 Mon Sep 17 00:00:00 2001 From: MatthijsV Date: Thu, 19 Oct 2017 13:31:07 +0100 Subject: [PATCH 06/14] commented out performance test --- test/unit/math/rev/mat/prob/normal_linear_glm_lpdf_test.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/unit/math/rev/mat/prob/normal_linear_glm_lpdf_test.cpp b/test/unit/math/rev/mat/prob/normal_linear_glm_lpdf_test.cpp index 95bf2fdabb1..311beff6d4b 100644 --- a/test/unit/math/rev/mat/prob/normal_linear_glm_lpdf_test.cpp +++ b/test/unit/math/rev/mat/prob/normal_linear_glm_lpdf_test.cpp @@ -105,7 +105,7 @@ TEST(ProbDistributionsNormalLinearGLM, glm_matches_normal_linear_vars) // Here, we compare the speed of the new regression to that of one built from // existing primitives. - +/* TEST(ProbDistributionsNormalLinearGLM, glm_matches_normal_linear_speed) { const int R = 30000; const int C = 1000; @@ -156,4 +156,4 @@ TEST(ProbDistributionsNormalLinearGLM, glm_matches_normal_linear_speed) { std::cout << "Existing Primitives:" << std::endl << T1 << std::endl << "New Primitives:" << std::endl << T2 << std::endl; } - +*/ From 7e4e5fd8db39174566fcc898f004ad2b7ae28443 Mon Sep 17 00:00:00 2001 From: MatthijsV Date: Tue, 24 Oct 2017 19:34:03 +0100 Subject: [PATCH 07/14] initial stab at neg_binomial_2_log_glm --- stan/math/prim/mat.hpp | 2 +- .../mat/prob/bernoulli_logit_glm_lpmf.hpp | 2 +- .../mat/prob/neg_binomial_2_log_glm_lpmf.hpp | 180 ++++++++++++++++++ .../prim/mat/prob/normal_linear_glm_lpdf.hpp | 160 ---------------- ...p => neg_binomial_2_log_glm_lpmf_test.cpp} | 77 ++++---- 5 files changed, 218 insertions(+), 203 deletions(-) create mode 100644 stan/math/prim/mat/prob/neg_binomial_2_log_glm_lpmf.hpp delete mode 100644 stan/math/prim/mat/prob/normal_linear_glm_lpdf.hpp rename test/unit/math/rev/mat/prob/{normal_linear_glm_lpdf_test.cpp => neg_binomial_2_log_glm_lpmf_test.cpp} (57%) diff --git a/stan/math/prim/mat.hpp b/stan/math/prim/mat.hpp index 84b20e1f8a1..6549a3fa434 100644 --- a/stan/math/prim/mat.hpp +++ b/stan/math/prim/mat.hpp @@ -288,7 +288,7 @@ #include #include #include -#include +#include #include #include #include diff --git a/stan/math/prim/mat/prob/bernoulli_logit_glm_lpmf.hpp b/stan/math/prim/mat/prob/bernoulli_logit_glm_lpmf.hpp index 57be9c29e89..55acb382c8f 100644 --- a/stan/math/prim/mat/prob/bernoulli_logit_glm_lpmf.hpp +++ b/stan/math/prim/mat/prob/bernoulli_logit_glm_lpmf.hpp @@ -79,7 +79,7 @@ namespace stan { const size_t N = x.col(0).size(); const size_t M = x.row(0).size(); - Matrix signs(N, 1); + Matrix signs(N, 1); { scalar_seq_view n_vec(n); for (size_t n = 0; n < N; ++n) { diff --git a/stan/math/prim/mat/prob/neg_binomial_2_log_glm_lpmf.hpp b/stan/math/prim/mat/prob/neg_binomial_2_log_glm_lpmf.hpp new file mode 100644 index 00000000000..6adc07d4e95 --- /dev/null +++ b/stan/math/prim/mat/prob/neg_binomial_2_log_glm_lpmf.hpp @@ -0,0 +1,180 @@ +#ifndef STAN_MATH_PRIM_MAT_PROB_NEG_BINOMIAL_2_LOG_GLM_LPMF_HPP +#define STAN_MATH_PRIM_MAT_PROB_NEG_BINOMIAL_2_LOG_GLM_LPMF_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { + namespace math { + + /** + * Returns the log PMF of the Generalized Linear Model (GLM) + * with Negative-Binomial-2 distribution and log link function. + * If containers are supplied, returns the log sum of the probabilities. + * @tparam T_n type of positive int vector of dependent variables (labels); + * this can also be a single positive integer value; + * @tparam T_x type of the matrix of independent variables (features); this + * should be an Eigen::Matrix type whose number of rows should match the + * length of n and whose number of columns should match the length of beta + * @tparam T_beta type of the weight vector; + * this can also be a single double value; + * @tparam T_alpha type of the intercept; + * this can either be a vector of doubles of a single double value; + * @tparam T_precision type of the precision vector; + * this can also be a single double value; + * @param n failures count vector parameter + * @param x design matrix + * @param beta weight vector + * @param alpha intercept (in log odds) + * @param phi (vector of) precision parameters + * @return log probability or log sum of probabilities + * @throw std::invalid_argument if container sizes mismatch. + */ + template + typename return_type::type + neg_binomial_2_log_glm_lpmf(const T_n &n, const T_x &x, const T_beta &beta, + const T_alpha &alpha, const T_precision &phi) { + static const std::string function = "neg_binomial_2_log_glm_lpmf"; + typedef typename stan::partials_return_type::type + T_partials_return; + + using std::exp; + using Eigen::Dynamic; + using Eigen::Matrix; + using Eigen::Array; + + if (!(stan::length(n) && stan::length(x) && stan::length(beta) + && stan::length(phi))) + return 0.0; + + T_partials_return logp(0.0); + + check_nonnegative(function, "Failures variables", n); + check_not_nan(function, "Matrix of independent variables", x); + check_not_nan(function, "Weight vector", beta); + check_not_nan(function, "Intercept", alpha); + check_positive_finite(function, "Precision parameter", phi); + check_consistent_sizes(function, + "Rows in matrix of independent variables", + x.col(0), "Vector of dependent variables", n); + check_consistent_sizes(function, + "Rows in matrix of independent variables", + x.col(0), "Vector of failure counts", phi); + check_consistent_sizes(function, + "Columns in matrix of independent variables", + x.row(0), "Weight vector", beta); + + if (!include_summand::value) + return 0.0; + + const size_t N = x.col(0).size(); + const size_t M = x.row(0).size(); + + Array n_arr(N, 1); + Array phi_arr(N, 1); + { + scalar_seq_view n_vec(n); + scalar_seq_view phi_vec(n); + for (size_t n = 0; n < N; ++n) { + n_arr[n] = n_vec[n]; + phi_arr[n] = value_of(phi_vec[n]); + } + } + Matrix beta_dbl(M, 1); + { + scalar_seq_view beta_vec(beta); + for (size_t m = 0; m < M; ++m) { + beta_dbl[m] = value_of(beta_vec[m]); + } + } + Matrix x_dbl = value_of(x); + + Array theta_dbl = (x_dbl * beta_dbl + + Matrix::Ones(N, 1) * value_of(alpha)).array(); + Array log_phi = phi_arr.log(); + Array logsumexp_eta_logphi = + theta_dbl.binaryExpr(log_phi, std::ref(log_sum_exp)); + Array n_plus_phi = n_arr + phi_arr; + + // Compute the log-density. + if (include_summand::value) + logp -= (n_vec + + Array::Ones(N, 1)) + .unaryExpr(std::ref(lgamma)).sum(); + if (include_summand::value) + logp += phi_arr.binaryExpr(phi_arr, std::ref(multiply_log)) + - phi_arr.unaryExpr(std::ref(lgamma)).sum(); + if (include_summand::value) + logp -= (n_plus_phi * logsumexp_eta_logphi).sum(); + if (include_summand::value) + logp += (n_arr * theta_dbl).sum(); + if (include_summand::value) + logp += n_plus_phi.unaryExpr(std::ref(lgamma)).sum(); + + // Compute the necessary derivatives. + operands_and_partials ops_partials(x, beta, alpha, phi); + if (!(is_constant_struct::value && is_constant_struct::value + && is_constant_struct::value + && is_constant_struct)) { + Matrix theta_derivative(N, 1); + theta_derivative = (n_arr - (n_plus_phi / (phi_arr / (theta_dbl.exp()) + + Array::Ones(N, 1)))).matrix(); + if (!is_constant_struct::value) { + ops_partials.edge2_.partials_ = x_dbl.transpose() * theta_derivative; + } + if (!is_constant_struct::value) { + ops_partials.edge1_.partials_ = theta_derivative + * beta_dbl.transpose(); + } + if (!is_constant_struct::value) { + ops_partials.edge3_.partials_[0] = theta_derivative.trace(); + } + if (!is_constant_struct::value) { + ops_partials.edge4_.partials_ = (Array::Ones(N, 1) + - n_plus_phi / (theta_dbl.exp() + phi_arr) + log_phi + - logsumexp_eta_logphi - phi.unaryExpr(std::ref(digamma)) + + n_plus_phi.unaryExpr(std::ref(digamma))).matrix(); + } + } + return ops_partials.build(logp); + } + + template + inline + typename return_type::type + neg_binomial_2_log_glm_lpmf(const T_n &n, const T_x &x, + const T_beta &beta, const T_alpha &alpha, + const T_precision& phi) { + return neg_binomial_2_log_glm_lpmf(n, x, beta, alpha, phi); + } + } // namespace math +} // namespace stan +#endif diff --git a/stan/math/prim/mat/prob/normal_linear_glm_lpdf.hpp b/stan/math/prim/mat/prob/normal_linear_glm_lpdf.hpp deleted file mode 100644 index e55bd6b7bbc..00000000000 --- a/stan/math/prim/mat/prob/normal_linear_glm_lpdf.hpp +++ /dev/null @@ -1,160 +0,0 @@ -#ifndef STAN_MATH_PRIM_MAT_PROB_NORMAL_LINEAR_GLM_LPDF_HPP -#define STAN_MATH_PRIM_MAT_PROB_NORMAL_LINEAR_GLM_LPDF_HPP - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -namespace stan { - namespace math { - - /** - * Returns the log PDF of the Generalized Linear Model (GLM) - * with Normal distribution and linear link function. - * If containers are supplied, returns the log sum of the probabilities. - * @tparam T_n type of real vector of dependent variables (labels); - * this can also be a single real value; - * @tparam T_x type of the matrix of independent variables (features); this - * should be an Eigen::Matrix type whose number of rows should match the - * length of n and whose number of columns should match the length of beta - * @tparam T_beta type of the weight vector; - * this can also be a single double value; - * @tparam T_alpha type of the intercept; - * this can either be a vector of doubles of a single double value; - * @tparam T_scale type of the scale vector. - * @param n real vector parameter - * @param x design matrix - * @param beta weight vector - * @param alpha intercept (in log odds) - * @param sigma (Sequence of) scale parameters for the normal - * distribution. - * @return log probability or log sum of probabilities - * @throw std::domain_error if x, beta or alpha is infinite. - * @throw std::invalid_argument if container sizes mismatch. - * @throw std::domain_error if the scale is not positive. - */ - template - typename return_type::type - normal_linear_glm_lpdf(const T_n &n, const T_x &x, const T_beta &beta, - const T_alpha &alpha, const T_scale& sigma) { - static const std::string function = "normal_linear_glm_lpdf"; - typedef typename stan::partials_return_type::type - T_partials_return; - - using std::exp; - using Eigen::Dynamic; - using Eigen::Matrix; - using Eigen::Array; - - if (!(stan::length(n) && stan::length(x) && stan::length(beta) - && stan::length(sigma))) - return 0.0; - - T_partials_return logp(0.0); - - check_not_nan(function, "Vector of dependent variables", n); - check_not_nan(function, "Matrix of independent variables", x); - check_not_nan(function, "Weight vector", beta); - check_not_nan(function, "Intercept", alpha); - check_positive(function, "Scale vector", sigma); - check_consistent_sizes(function, - "Rows in matrix of independent variables", - x.col(0), "Vector of dependent variables", n); - check_consistent_sizes(function, - "Columns in matrix of independent variables", - x.row(0), "Weight vector", beta); - - if (!include_summand::value) - return 0.0; - - const size_t N = x.col(0).size(); - const size_t M = x.row(0).size(); - - Array sigma_dbl(N, 1); - Array n_dbl(N, 1); - { - scalar_seq_view n_vec(n); - scalar_seq_view sigma_vec(sigma); - for (size_t n = 0; n < N; ++n) { - sigma_dbl[n] = value_of(sigma_vec[n]); - n_dbl[n] = n_vec[n]; - } - } - Array inv_sigma = 1/sigma_dbl; - Matrix beta_dbl(M, 1); - { - scalar_seq_view beta_vec(beta); - for (size_t m = 0; m < M; ++m) { - beta_dbl[m] = value_of(beta_vec[m]); - } - } - Matrix x_dbl = value_of(x); - - Array mu_dbl = (x_dbl * beta_dbl - + Matrix::Ones(N, 1) - * value_of(alpha)).array(); - Array n_minus_mu_over_sigma = - (n_dbl - mu_dbl) * inv_sigma; - Array n_minus_mu_over_sigma_squared - = n_minus_mu_over_sigma.square(); - - if (include_summand::value) - logp += NEG_LOG_SQRT_TWO_PI * N; - if (include_summand::value) - logp -= sigma_dbl.log().sum(); - if (include_summand::value) - logp -= 0.5 * n_minus_mu_over_sigma_squared.sum(); - - // Compute the necessary derivatives. - operands_and_partials ops_partials(x, beta, - alpha, sigma); - if (!(is_constant_struct::value && is_constant_struct::value - && is_constant_struct::value)) { - Matrix mu_derivative = (inv_sigma - * n_minus_mu_over_sigma).matrix(); - if (!is_constant_struct::value) { - ops_partials.edge2_.partials_ = x_dbl.transpose() * mu_derivative; - } - if (!is_constant_struct::value) { - ops_partials.edge1_.partials_ = mu_derivative - * beta_dbl.transpose(); - } - if (!is_constant_struct::value) { - ops_partials.edge3_.partials_[0] = mu_derivative.trace(); - } - if (!is_constant_struct::value) { - ops_partials.edge4_.partials_ = ((inv_sigma - Array::Ones(N, 1)) * n_minus_mu_over_sigma_squared).matrix(); - } - } - - return ops_partials.build(logp); - } - - template - inline - typename return_type::type - normal_linear_glm_lpdf(const T_n &n, const T_x &x, const T_beta &beta, - const T_alpha &alpha, const T_scale& sigma) { - return normal_linear_glm_lpdf(n, x, beta, alpha, sigma); - } - } // namespace math -} // namespace stan -#endif diff --git a/test/unit/math/rev/mat/prob/normal_linear_glm_lpdf_test.cpp b/test/unit/math/rev/mat/prob/neg_binomial_2_log_glm_lpmf_test.cpp similarity index 57% rename from test/unit/math/rev/mat/prob/normal_linear_glm_lpdf_test.cpp rename to test/unit/math/rev/mat/prob/neg_binomial_2_log_glm_lpmf_test.cpp index 311beff6d4b..b195d9508e7 100644 --- a/test/unit/math/rev/mat/prob/normal_linear_glm_lpdf_test.cpp +++ b/test/unit/math/rev/mat/prob/neg_binomial_2_log_glm_lpmf_test.cpp @@ -13,76 +13,72 @@ typedef std::chrono::high_resolution_clock::time_point TimeVar; // We check that the values of the new regression match those of one built // from existing primitives. -TEST(ProbDistributionsNormalLinearGLM, glm_matches_normal_linear_doubles) +TEST(ProbDistributionsNegBinomial2LogGLM, glm_matches_neg_binomial_2_log_doubles) { Matrix n(3, 1); - n << 51, 32, 12; + n << 1, 0, 1; Matrix x(3, 2); x << -12, 46, -42, 24, 25, 27; Matrix beta(2, 1); beta << 0.3, 2; double alpha = 0.3; - Matrix alphavec = alpha * Matrix::Ones(); + Matrix alphavec = alpha * Matrix::Ones(); Matrix theta(3, 1); theta = x * beta + alphavec; - Matrix sigma(3, 1); - sigma << 10, 4, 6; - - EXPECT_FLOAT_EQ((stan::math::normal_lpdf(n, theta, sigma)), - (stan::math::normal_linear_glm_lpdf(n, x, beta, alpha, - sigma))); - EXPECT_FLOAT_EQ((stan::math::normal_lpdf(n, theta, sigma)), - (stan::math::normal_linear_glm_lpdf(n, x, beta, alpha, - sigma))); - EXPECT_FLOAT_EQ((stan::math::normal_lpdf(n, theta, sigma)), - (stan::math::normal_linear_glm_lpdf(n, x, beta, alpha, - sigma))); - EXPECT_FLOAT_EQ((stan::math::normal_lpdf>(n, - theta, sigma)), - (stan::math::normal_linear_glm_lpdf>(n, x, beta, alpha, sigma))); - EXPECT_FLOAT_EQ((stan::math::normal_lpdf>(n, theta, sigma)), - (stan::math::normal_linear_glm_lpdf>(n, x, beta, alpha, sigma))); - EXPECT_FLOAT_EQ((stan::math::normal_lpdf>(n, theta, sigma)), - (stan::math::normal_linear_glm_lpdf>(n, x, beta, alpha, sigma))); + Matrix phi(3, 1); + phi << 2, 1, 0.2; + + EXPECT_FLOAT_EQ((stan::math::neg_binomial_2_log_lpmf(n, theta, phi)), + (stan::math::neg_binomial_2_log_glm_lpmf(n, x, beta, alpha, phi))); + EXPECT_FLOAT_EQ((stan::math::neg_binomial_2_log_lpmf(n, theta, phi)), + (stan::math::neg_binomial_2_log_glm_lpmf(n, x, beta, alpha, phi))); + EXPECT_FLOAT_EQ((stan::math::neg_binomial_2_log_lpmf(n, theta, phi)), + (stan::math::neg_binomial_2_log_glm_lpmf(n, x, beta, alpha, phi))); + EXPECT_FLOAT_EQ((stan::math::neg_binomial_2_log_lpmf>(n, theta, phi)), + (stan::math::neg_binomial_2_log_glm_lpmf>(n, x, beta, alpha, phi))); + EXPECT_FLOAT_EQ((stan::math::neg_binomial_2_log_lpmf>(n, theta, phi)), + (stan::math::neg_binomial_2_log_glm_lpmf>(n, x, beta, alpha, phi))); + EXPECT_FLOAT_EQ((stan::math::neg_binomial_2_log_lpmf>(n, theta, phi)), + (stan::math::neg_binomial_2_log_glm_lpmf>(n, x, beta, alpha, phi))); } // We check that the gradients of the new regression match those of one built // from existing primitives. -TEST(ProbDistributionsNormalLinearGLM, glm_matches_normal_linear_vars) +TEST(ProbDistributionsNegBinomial2LogGLM, glm_matches_neg_binomial_2_log_vars) { Matrix n(3, 1); - n << 14, 32, 21; + n << 1, 0, 1; Matrix x(3, 2); x << -12, 46, -42, 24, 25, 27; Matrix beta(2, 1); beta << 0.3, 2; var alpha = 0.3; - Matrix sigma(3, 1); - sigma << 10, 4, 6; Matrix alphavec = alpha * Matrix::Ones(); Matrix theta(3, 1); theta = x * beta + alphavec; + Matrix phi(3, 1); + phi << 2, 1, 0.2; - var lp = stan::math::normal_lpdf(n, theta, sigma); + var lp = stan::math::neg_binomial_2_log_lpmf(n, theta, phi); lp.grad(); stan::math::recover_memory(); Matrix n2(3, 1); - n2 << 14, 32, 21; + n2 << 1, 0, 1; Matrix x2(3, 2); x2 << -12, 46, -42, 24, 25, 27; Matrix beta2(2, 1); beta2 << 0.3, 2; var alpha2 = 0.3; - Matrix sigma2(3, 1); - sigma2 << 10, 4, 6; + Matrix phi2(3, 1); + phi2 << 2, 1, 0.2; - var lp2 = stan::math::normal_linear_glm_lpdf(n2, x2, beta2, alpha2, sigma2); + var lp2 = stan::math::neg_binomial_2_log_glm_lpmf(n2, x2, beta2, alpha2, phi2); lp2.grad(); EXPECT_FLOAT_EQ(lp.val(), @@ -94,7 +90,7 @@ TEST(ProbDistributionsNormalLinearGLM, glm_matches_normal_linear_vars) EXPECT_FLOAT_EQ(alpha.adj(), alpha2.adj()); for (size_t j = 0; j < 3; j++) { - EXPECT_FLOAT_EQ(sigma[j].adj(), sigma2[j].adj()); + EXPECT_FLOAT_EQ(phi[i].adj(), phi2[i].adj()); for (size_t i = 0; i < 2; i++) { EXPECT_FLOAT_EQ(x(j, i).adj(), x2(j, i).adj()); @@ -106,13 +102,13 @@ TEST(ProbDistributionsNormalLinearGLM, glm_matches_normal_linear_vars) // existing primitives. /* -TEST(ProbDistributionsNormalLinearGLM, glm_matches_normal_linear_speed) { +TEST(ProbDistributionsNegBinomial2LogGLM, glm_matches_neg_binomial_2_log_speed) { const int R = 30000; const int C = 1000; Matrix n(R, 1); for (size_t i = 0; i < R; i++) { - n[i] = rand()%2000; + n[i] = rand()%2; } int T1 = 0; @@ -121,19 +117,18 @@ TEST(ProbDistributionsNormalLinearGLM, glm_matches_normal_linear_speed) { for (size_t testnumber = 0; testnumber < 30; testnumber++){ Matrix xreal = Matrix::Random(R, C); Matrix betareal = Matrix::Random(C, 1); - Matrix sigmareal = Matrix::Random(R, 1) - + Matrix::Ones(); // Random Matrix has entries between -1 and 1. We add 1 to it to get positive entries. Matrix alphareal = Matrix::Random(1, 1); Matrix alpharealvec = Matrix::Ones() * alphareal; + Matrix phireal = Matrix::Random(R, 1); Matrix beta = betareal; - Matrix sigma = sigmareal; + Matrix phi = phireal; Matrix theta(R, 1); TimeVar t1 = timeNow(); theta = (xreal * beta) + alpharealvec; - var lp = stan::math::normal_lpdf(n, theta, sigma); + var lp = stan::math::neg_binomial_2_log_lpmf(n, theta, phi); lp.grad(); TimeVar t2 = timeNow(); @@ -141,11 +136,10 @@ TEST(ProbDistributionsNormalLinearGLM, glm_matches_normal_linear_speed) { stan::math::recover_memory(); Matrix beta2 = betareal; - Matrix sigma2 = sigmareal; + Matrix phi2 = phireal; TimeVar t3 = timeNow(); - var lp2 = stan::math::normal_linear_glm_lpdf(n, xreal, beta2, alphareal[0], - sigma2); + var lp2 = stan::math::neg_binomial_2_log_glm_lpmf(n, xreal, beta2, alphareal[0], phi2); lp2.grad(); TimeVar t4 = timeNow(); stan::math::recover_memory(); @@ -157,3 +151,4 @@ TEST(ProbDistributionsNormalLinearGLM, glm_matches_normal_linear_speed) { std::cout << "Existing Primitives:" << std::endl << T1 << std::endl << "New Primitives:" << std::endl << T2 << std::endl; } */ + From 95f1a74475ab941ee21961a46441523e42e1181c Mon Sep 17 00:00:00 2001 From: MatthijsV Date: Tue, 24 Oct 2017 20:59:27 +0100 Subject: [PATCH 08/14] neg_binomial_2_log_glm finished and working --- .../mat/prob/neg_binomial_2_log_glm_lpmf.hpp | 47 ++++++++++++++----- .../prob/neg_binomial_2_log_glm_lpmf_test.cpp | 9 ++-- 2 files changed, 39 insertions(+), 17 deletions(-) diff --git a/stan/math/prim/mat/prob/neg_binomial_2_log_glm_lpmf.hpp b/stan/math/prim/mat/prob/neg_binomial_2_log_glm_lpmf.hpp index 6adc07d4e95..78395978066 100644 --- a/stan/math/prim/mat/prob/neg_binomial_2_log_glm_lpmf.hpp +++ b/stan/math/prim/mat/prob/neg_binomial_2_log_glm_lpmf.hpp @@ -90,9 +90,9 @@ namespace stan { "Columns in matrix of independent variables", x.row(0), "Weight vector", beta); - if (!include_summand::value) + if (!include_summand::value) return 0.0; - + const size_t N = x.col(0).size(); const size_t M = x.row(0).size(); @@ -100,12 +100,13 @@ namespace stan { Array phi_arr(N, 1); { scalar_seq_view n_vec(n); - scalar_seq_view phi_vec(n); + scalar_seq_view phi_vec(phi); for (size_t n = 0; n < N; ++n) { n_arr[n] = n_vec[n]; phi_arr[n] = value_of(phi_vec[n]); } } + Matrix beta_dbl(M, 1); { scalar_seq_view beta_vec(beta); @@ -114,35 +115,50 @@ namespace stan { } } Matrix x_dbl = value_of(x); - + Array theta_dbl = (x_dbl * beta_dbl + Matrix::Ones(N, 1) * value_of(alpha)).array(); Array log_phi = phi_arr.log(); Array logsumexp_eta_logphi = - theta_dbl.binaryExpr(log_phi, std::ref(log_sum_exp)); + theta_dbl.binaryExpr(log_phi, + [](T_partials_return xx, T_partials_return yy) { + return log_sum_exp(xx, yy); + }); Array n_plus_phi = n_arr + phi_arr; // Compute the log-density. + if (include_summand::value) - logp -= (n_vec + logp -= (n_arr + Array::Ones(N, 1)) - .unaryExpr(std::ref(lgamma)).sum(); + .unaryExpr([](T_partials_return xx) { + return lgamma(xx); + }).sum(); if (include_summand::value) - logp += phi_arr.binaryExpr(phi_arr, std::ref(multiply_log)) - - phi_arr.unaryExpr(std::ref(lgamma)).sum(); + logp += (phi_arr.binaryExpr(phi_arr, [](T_partials_return xx, + T_partials_return yy) { + return multiply_log(xx, yy); + }) + - phi_arr.unaryExpr([](T_partials_return xx) { + return lgamma(xx); + })).sum(); if (include_summand::value) logp -= (n_plus_phi * logsumexp_eta_logphi).sum(); if (include_summand::value) logp += (n_arr * theta_dbl).sum(); if (include_summand::value) - logp += n_plus_phi.unaryExpr(std::ref(lgamma)).sum(); + logp += n_plus_phi.unaryExpr([](T_partials_return xx) { + return lgamma(xx); + }).sum(); + // Compute the necessary derivatives. operands_and_partials ops_partials(x, beta, alpha, phi); + if (!(is_constant_struct::value && is_constant_struct::value && is_constant_struct::value - && is_constant_struct)) { + && is_constant_struct::value)) { Matrix theta_derivative(N, 1); theta_derivative = (n_arr - (n_plus_phi / (phi_arr / (theta_dbl.exp()) + Array::Ones(N, 1)))).matrix(); @@ -159,10 +175,15 @@ namespace stan { if (!is_constant_struct::value) { ops_partials.edge4_.partials_ = (Array::Ones(N, 1) - n_plus_phi / (theta_dbl.exp() + phi_arr) + log_phi - - logsumexp_eta_logphi - phi.unaryExpr(std::ref(digamma)) - + n_plus_phi.unaryExpr(std::ref(digamma))).matrix(); + - logsumexp_eta_logphi - phi_arr.unaryExpr([](T_partials_return xx) { + return digamma(xx); + }) + + n_plus_phi.unaryExpr([](T_partials_return xx) { + return digamma(xx); + })).matrix(); } } + return ops_partials.build(logp); } diff --git a/test/unit/math/rev/mat/prob/neg_binomial_2_log_glm_lpmf_test.cpp b/test/unit/math/rev/mat/prob/neg_binomial_2_log_glm_lpmf_test.cpp index b195d9508e7..4ee27539597 100644 --- a/test/unit/math/rev/mat/prob/neg_binomial_2_log_glm_lpmf_test.cpp +++ b/test/unit/math/rev/mat/prob/neg_binomial_2_log_glm_lpmf_test.cpp @@ -31,6 +31,7 @@ TEST(ProbDistributionsNegBinomial2LogGLM, glm_matches_neg_binomial_2_log_doubles EXPECT_FLOAT_EQ((stan::math::neg_binomial_2_log_lpmf(n, theta, phi)), (stan::math::neg_binomial_2_log_glm_lpmf(n, x, beta, alpha, phi))); + EXPECT_FLOAT_EQ((stan::math::neg_binomial_2_log_lpmf(n, theta, phi)), (stan::math::neg_binomial_2_log_glm_lpmf(n, x, beta, alpha, phi))); EXPECT_FLOAT_EQ((stan::math::neg_binomial_2_log_lpmf(n, theta, phi)), @@ -90,7 +91,7 @@ TEST(ProbDistributionsNegBinomial2LogGLM, glm_matches_neg_binomial_2_log_vars) EXPECT_FLOAT_EQ(alpha.adj(), alpha2.adj()); for (size_t j = 0; j < 3; j++) { - EXPECT_FLOAT_EQ(phi[i].adj(), phi2[i].adj()); + EXPECT_FLOAT_EQ(phi[j].adj(), phi2[j].adj()); for (size_t i = 0; i < 2; i++) { EXPECT_FLOAT_EQ(x(j, i).adj(), x2(j, i).adj()); @@ -101,7 +102,7 @@ TEST(ProbDistributionsNegBinomial2LogGLM, glm_matches_neg_binomial_2_log_vars) // Here, we compare the speed of the new regression to that of one built from // existing primitives. -/* + TEST(ProbDistributionsNegBinomial2LogGLM, glm_matches_neg_binomial_2_log_speed) { const int R = 30000; const int C = 1000; @@ -119,7 +120,7 @@ TEST(ProbDistributionsNegBinomial2LogGLM, glm_matches_neg_binomial_2_log_speed) Matrix betareal = Matrix::Random(C, 1); Matrix alphareal = Matrix::Random(1, 1); Matrix alpharealvec = Matrix::Ones() * alphareal; - Matrix phireal = Matrix::Random(R, 1); + Matrix phireal = Matrix::Random(R, 1) + Matrix::Ones(R, 1); // We want phireal to be positive. Matrix beta = betareal; Matrix phi = phireal; @@ -150,5 +151,5 @@ TEST(ProbDistributionsNegBinomial2LogGLM, glm_matches_neg_binomial_2_log_speed) std::cout << "Existing Primitives:" << std::endl << T1 << std::endl << "New Primitives:" << std::endl << T2 << std::endl; } -*/ + From 6f71ca9cc0acc8cbb5e9ed657090a3aef8e7fd92 Mon Sep 17 00:00:00 2001 From: MatthijsV Date: Wed, 25 Oct 2017 09:40:38 +0100 Subject: [PATCH 09/14] minor speedup --- .../mat/prob/neg_binomial_2_log_glm_lpmf.hpp | 24 +++++++++---------- .../prob/neg_binomial_2_log_glm_lpmf_test.cpp | 4 ++-- 2 files changed, 13 insertions(+), 15 deletions(-) diff --git a/stan/math/prim/mat/prob/neg_binomial_2_log_glm_lpmf.hpp b/stan/math/prim/mat/prob/neg_binomial_2_log_glm_lpmf.hpp index 78395978066..f5d291b7383 100644 --- a/stan/math/prim/mat/prob/neg_binomial_2_log_glm_lpmf.hpp +++ b/stan/math/prim/mat/prob/neg_binomial_2_log_glm_lpmf.hpp @@ -157,8 +157,7 @@ namespace stan { T_precision> ops_partials(x, beta, alpha, phi); if (!(is_constant_struct::value && is_constant_struct::value - && is_constant_struct::value - && is_constant_struct::value)) { + && is_constant_struct::value)) { Matrix theta_derivative(N, 1); theta_derivative = (n_arr - (n_plus_phi / (phi_arr / (theta_dbl.exp()) + Array::Ones(N, 1)))).matrix(); @@ -172,18 +171,17 @@ namespace stan { if (!is_constant_struct::value) { ops_partials.edge3_.partials_[0] = theta_derivative.trace(); } - if (!is_constant_struct::value) { - ops_partials.edge4_.partials_ = (Array::Ones(N, 1) - - n_plus_phi / (theta_dbl.exp() + phi_arr) + log_phi - - logsumexp_eta_logphi - phi_arr.unaryExpr([](T_partials_return xx) { - return digamma(xx); - }) - + n_plus_phi.unaryExpr([](T_partials_return xx) { - return digamma(xx); - })).matrix(); - } } - + if (!is_constant_struct::value) { + ops_partials.edge4_.partials_ = (Array::Ones(N, 1) + - n_plus_phi / (theta_dbl.exp() + phi_arr) + log_phi + - logsumexp_eta_logphi - phi_arr.unaryExpr([](T_partials_return xx) { + return digamma(xx); + }) + + n_plus_phi.unaryExpr([](T_partials_return xx) { + return digamma(xx); + })).matrix(); + } return ops_partials.build(logp); } diff --git a/test/unit/math/rev/mat/prob/neg_binomial_2_log_glm_lpmf_test.cpp b/test/unit/math/rev/mat/prob/neg_binomial_2_log_glm_lpmf_test.cpp index 4ee27539597..a3bc203891b 100644 --- a/test/unit/math/rev/mat/prob/neg_binomial_2_log_glm_lpmf_test.cpp +++ b/test/unit/math/rev/mat/prob/neg_binomial_2_log_glm_lpmf_test.cpp @@ -102,7 +102,7 @@ TEST(ProbDistributionsNegBinomial2LogGLM, glm_matches_neg_binomial_2_log_vars) // Here, we compare the speed of the new regression to that of one built from // existing primitives. - +/* TEST(ProbDistributionsNegBinomial2LogGLM, glm_matches_neg_binomial_2_log_speed) { const int R = 30000; const int C = 1000; @@ -151,5 +151,5 @@ TEST(ProbDistributionsNegBinomial2LogGLM, glm_matches_neg_binomial_2_log_speed) std::cout << "Existing Primitives:" << std::endl << T1 << std::endl << "New Primitives:" << std::endl << T2 << std::endl; } - +*/ From 39065624a808c581f4fa21a27b92ae5e54d6ad1d Mon Sep 17 00:00:00 2001 From: MatthijsV Date: Wed, 1 Nov 2017 11:14:24 +0000 Subject: [PATCH 10/14] done for review --- .../mat/prob/bernoulli_logit_glm_lpmf.hpp | 17 +++--- .../mat/prob/neg_binomial_2_log_glm_lpmf.hpp | 60 +++++++++---------- 2 files changed, 35 insertions(+), 42 deletions(-) diff --git a/stan/math/prim/mat/prob/bernoulli_logit_glm_lpmf.hpp b/stan/math/prim/mat/prob/bernoulli_logit_glm_lpmf.hpp index 55acb382c8f..0f8afcd5b49 100644 --- a/stan/math/prim/mat/prob/bernoulli_logit_glm_lpmf.hpp +++ b/stan/math/prim/mat/prob/bernoulli_logit_glm_lpmf.hpp @@ -7,14 +7,10 @@ #include #include #include -#include #include -#include #include #include #include -#include -#include #include #include #include @@ -32,15 +28,16 @@ namespace stan { * should be an Eigen::Matrix type whose number of rows should match the * length of n and whose number of columns should match the length of beta * @tparam T_beta type of the weight vector; - * this can also be a single double value; + * this can also be a single value; * @tparam T_alpha type of the intercept; - * this can either be a vector of doubles of a single double value; + * this has to be single value; * @param n binary vector parameter * @param x design matrix * @param beta weight vector * @param alpha intercept (in log odds) * @return log probability or log sum of probabilities - * @throw std::domain_error if theta is infinite. + * @throw std::domain_error if x, beta or alpha is infinite. + * @throw std::domain_error if n is not binary. * @throw std::invalid_argument if container sizes mismatch. */ template #include #include -#include -#include #include -#include #include -#include #include #include #include @@ -21,10 +17,6 @@ #include #include #include -#include -#include -#include -#include #include #include @@ -41,11 +33,11 @@ namespace stan { * should be an Eigen::Matrix type whose number of rows should match the * length of n and whose number of columns should match the length of beta * @tparam T_beta type of the weight vector; - * this can also be a single double value; + * this can also be a single value; * @tparam T_alpha type of the intercept; - * this can either be a vector of doubles of a single double value; - * @tparam T_precision type of the precision vector; - * this can also be a single double value; + * this should be a single value; + * @tparam T_precision type of the (positive) precision vector phi; + * this can also be a single value; * @param n failures count vector parameter * @param x design matrix * @param beta weight vector @@ -53,6 +45,9 @@ namespace stan { * @param phi (vector of) precision parameters * @return log probability or log sum of probabilities * @throw std::invalid_argument if container sizes mismatch. + * @throw std::domain_error if x, beta or alpha is infinite. + * @throw std::domain_error if phi is infinite or non-positive. + * @throw std::domain_error if n is negative. */ template @@ -76,9 +71,9 @@ namespace stan { T_partials_return logp(0.0); check_nonnegative(function, "Failures variables", n); - check_not_nan(function, "Matrix of independent variables", x); - check_not_nan(function, "Weight vector", beta); - check_not_nan(function, "Intercept", alpha); + check_finite(function, "Matrix of independent variables", x); + check_finite(function, "Weight vector", beta); + check_finite(function, "Intercept", alpha); check_positive_finite(function, "Precision parameter", phi); check_consistent_sizes(function, "Rows in matrix of independent variables", @@ -92,7 +87,7 @@ namespace stan { if (!include_summand::value) return 0.0; - + const size_t N = x.col(0).size(); const size_t M = x.row(0).size(); @@ -106,7 +101,7 @@ namespace stan { phi_arr[n] = value_of(phi_vec[n]); } } - + Matrix beta_dbl(M, 1); { scalar_seq_view beta_vec(beta); @@ -115,7 +110,7 @@ namespace stan { } } Matrix x_dbl = value_of(x); - + Array theta_dbl = (x_dbl * beta_dbl + Matrix::Ones(N, 1) * value_of(alpha)).array(); Array log_phi = phi_arr.log(); @@ -127,35 +122,36 @@ namespace stan { Array n_plus_phi = n_arr + phi_arr; // Compute the log-density. - - if (include_summand::value) + if (include_summand::value) { logp -= (n_arr + Array::Ones(N, 1)) .unaryExpr([](T_partials_return xx) { - return lgamma(xx); - }).sum(); - if (include_summand::value) + return lgamma(xx); + }).sum(); + } + if (include_summand::value) { logp += (phi_arr.binaryExpr(phi_arr, [](T_partials_return xx, - T_partials_return yy) { - return multiply_log(xx, yy); - }) + T_partials_return yy) { + return multiply_log(xx, yy); + }) - phi_arr.unaryExpr([](T_partials_return xx) { return lgamma(xx); })).sum(); + } if (include_summand::value) logp -= (n_plus_phi * logsumexp_eta_logphi).sum(); if (include_summand::value) logp += (n_arr * theta_dbl).sum(); - if (include_summand::value) + if (include_summand::value) { logp += n_plus_phi.unaryExpr([](T_partials_return xx) { return lgamma(xx); }).sum(); - + } // Compute the necessary derivatives. - operands_and_partials ops_partials(x, beta, alpha, phi); - + if (!(is_constant_struct::value && is_constant_struct::value && is_constant_struct::value)) { Matrix theta_derivative(N, 1); @@ -176,8 +172,8 @@ namespace stan { ops_partials.edge4_.partials_ = (Array::Ones(N, 1) - n_plus_phi / (theta_dbl.exp() + phi_arr) + log_phi - logsumexp_eta_logphi - phi_arr.unaryExpr([](T_partials_return xx) { - return digamma(xx); - }) + return digamma(xx); + }) + n_plus_phi.unaryExpr([](T_partials_return xx) { return digamma(xx); })).matrix(); From 625ace15809f7cbd062cd0db55695dfceb75574b Mon Sep 17 00:00:00 2001 From: MatthijsV Date: Wed, 1 Nov 2017 15:41:52 +0000 Subject: [PATCH 11/14] satisfied cpplint for test --- .../prob/neg_binomial_2_log_glm_lpmf_test.cpp | 69 +++++++++++-------- 1 file changed, 39 insertions(+), 30 deletions(-) diff --git a/test/unit/math/rev/mat/prob/neg_binomial_2_log_glm_lpmf_test.cpp b/test/unit/math/rev/mat/prob/neg_binomial_2_log_glm_lpmf_test.cpp index a3bc203891b..d1edf76809d 100644 --- a/test/unit/math/rev/mat/prob/neg_binomial_2_log_glm_lpmf_test.cpp +++ b/test/unit/math/rev/mat/prob/neg_binomial_2_log_glm_lpmf_test.cpp @@ -1,20 +1,16 @@ #include #include #include -#include + using stan::math::var; using Eigen::Dynamic; using Eigen::Matrix; -typedef std::chrono::high_resolution_clock::time_point TimeVar; -#define duration(a) std::chrono::duration_cast(a).count() -#define timeNow() std::chrono::high_resolution_clock::now() - // We check that the values of the new regression match those of one built // from existing primitives. -TEST(ProbDistributionsNegBinomial2LogGLM, glm_matches_neg_binomial_2_log_doubles) -{ +TEST(ProbDistributionsNegBinomial2LogGLM, + glm_matches_neg_binomial_2_log_doubles) { Matrix n(3, 1); n << 1, 0, 1; Matrix x(3, 2); @@ -23,32 +19,42 @@ TEST(ProbDistributionsNegBinomial2LogGLM, glm_matches_neg_binomial_2_log_doubles Matrix beta(2, 1); beta << 0.3, 2; double alpha = 0.3; - Matrix alphavec = alpha * Matrix::Ones(); + Matrix alphavec = alpha * Matrix::Ones(); Matrix theta(3, 1); theta = x * beta + alphavec; Matrix phi(3, 1); phi << 2, 1, 0.2; EXPECT_FLOAT_EQ((stan::math::neg_binomial_2_log_lpmf(n, theta, phi)), - (stan::math::neg_binomial_2_log_glm_lpmf(n, x, beta, alpha, phi))); + (stan::math::neg_binomial_2_log_glm_lpmf(n, x, beta, alpha, + phi))); EXPECT_FLOAT_EQ((stan::math::neg_binomial_2_log_lpmf(n, theta, phi)), - (stan::math::neg_binomial_2_log_glm_lpmf(n, x, beta, alpha, phi))); + (stan::math::neg_binomial_2_log_glm_lpmf(n, x, beta, + alpha, phi))); EXPECT_FLOAT_EQ((stan::math::neg_binomial_2_log_lpmf(n, theta, phi)), - (stan::math::neg_binomial_2_log_glm_lpmf(n, x, beta, alpha, phi))); - EXPECT_FLOAT_EQ((stan::math::neg_binomial_2_log_lpmf>(n, theta, phi)), - (stan::math::neg_binomial_2_log_glm_lpmf>(n, x, beta, alpha, phi))); - EXPECT_FLOAT_EQ((stan::math::neg_binomial_2_log_lpmf>(n, theta, phi)), - (stan::math::neg_binomial_2_log_glm_lpmf>(n, x, beta, alpha, phi))); - EXPECT_FLOAT_EQ((stan::math::neg_binomial_2_log_lpmf>(n, theta, phi)), - (stan::math::neg_binomial_2_log_glm_lpmf>(n, x, beta, alpha, phi))); + (stan::math::neg_binomial_2_log_glm_lpmf(n, x, beta, + alpha, phi))); + EXPECT_FLOAT_EQ( + (stan::math::neg_binomial_2_log_lpmf + >(n, theta, phi)), + (stan::math::neg_binomial_2_log_glm_lpmf + >(n, x, beta, alpha, phi))); + EXPECT_FLOAT_EQ( + (stan::math::neg_binomial_2_log_lpmf + >(n, theta, phi)), + (stan::math::neg_binomial_2_log_glm_lpmf + >(n, x, beta, alpha, phi))); + EXPECT_FLOAT_EQ( + (stan::math::neg_binomial_2_log_lpmf + >(n, theta, phi)), + (stan::math::neg_binomial_2_log_glm_lpmf + >(n, x, beta, alpha, phi))); } // We check that the gradients of the new regression match those of one built // from existing primitives. -TEST(ProbDistributionsNegBinomial2LogGLM, glm_matches_neg_binomial_2_log_vars) -{ - +TEST(ProbDistributionsNegBinomial2LogGLM, glm_matches_neg_binomial_2_log_vars) { Matrix n(3, 1); n << 1, 0, 1; Matrix x(3, 2); @@ -57,7 +63,7 @@ TEST(ProbDistributionsNegBinomial2LogGLM, glm_matches_neg_binomial_2_log_vars) Matrix beta(2, 1); beta << 0.3, 2; var alpha = 0.3; - Matrix alphavec = alpha * Matrix::Ones(); + Matrix alphavec = alpha * Matrix::Ones(); Matrix theta(3, 1); theta = x * beta + alphavec; Matrix phi(3, 1); @@ -78,22 +84,20 @@ TEST(ProbDistributionsNegBinomial2LogGLM, glm_matches_neg_binomial_2_log_vars) var alpha2 = 0.3; Matrix phi2(3, 1); phi2 << 2, 1, 0.2; - - var lp2 = stan::math::neg_binomial_2_log_glm_lpmf(n2, x2, beta2, alpha2, phi2); + + var lp2 = stan::math::neg_binomial_2_log_glm_lpmf(n2, x2, beta2, alpha2, + phi2); lp2.grad(); EXPECT_FLOAT_EQ(lp.val(), lp2.val()); - for (size_t i = 0; i < 2; i++) - { + for (size_t i = 0; i < 2; i++) { EXPECT_FLOAT_EQ(beta[i].adj(), beta2[i].adj()); } EXPECT_FLOAT_EQ(alpha.adj(), alpha2.adj()); - for (size_t j = 0; j < 3; j++) - { + for (size_t j = 0; j < 3; j++) { EXPECT_FLOAT_EQ(phi[j].adj(), phi2[j].adj()); - for (size_t i = 0; i < 2; i++) - { + for (size_t i = 0; i < 2; i++) { EXPECT_FLOAT_EQ(x(j, i).adj(), x2(j, i).adj()); } } @@ -103,6 +107,12 @@ TEST(ProbDistributionsNegBinomial2LogGLM, glm_matches_neg_binomial_2_log_vars) // existing primitives. /* +#include +typedef std::chrono::high_resolution_clock::time_point TimeVar; +#define duration(a) \ + std::chrono::duration_cast(a).count() +#define timeNow() std::chrono::high_resolution_clock::now() + TEST(ProbDistributionsNegBinomial2LogGLM, glm_matches_neg_binomial_2_log_speed) { const int R = 30000; const int C = 1000; @@ -152,4 +162,3 @@ TEST(ProbDistributionsNegBinomial2LogGLM, glm_matches_neg_binomial_2_log_speed) std::cout << "Existing Primitives:" << std::endl << T1 << std::endl << "New Primitives:" << std::endl << T2 << std::endl; } */ - From dde66fddd7b761c9677354461ceacf665cc59be9 Mon Sep 17 00:00:00 2001 From: MatthijsV Date: Fri, 10 Nov 2017 16:44:48 +0000 Subject: [PATCH 12/14] minor --- .../prim/mat/prob/neg_binomial_2_log_glm_lpmf.hpp | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/stan/math/prim/mat/prob/neg_binomial_2_log_glm_lpmf.hpp b/stan/math/prim/mat/prob/neg_binomial_2_log_glm_lpmf.hpp index 53a0a6d9f16..b6938071e71 100644 --- a/stan/math/prim/mat/prob/neg_binomial_2_log_glm_lpmf.hpp +++ b/stan/math/prim/mat/prob/neg_binomial_2_log_glm_lpmf.hpp @@ -18,7 +18,6 @@ #include #include #include -#include namespace stan { namespace math { @@ -27,17 +26,17 @@ namespace stan { * Returns the log PMF of the Generalized Linear Model (GLM) * with Negative-Binomial-2 distribution and log link function. * If containers are supplied, returns the log sum of the probabilities. - * @tparam T_n type of positive int vector of dependent variables (labels); + * @tparam T_n type of positive int vector of variates (labels); * this can also be a single positive integer value; - * @tparam T_x type of the matrix of independent variables (features); this + * @tparam T_x type of the matrix of covariates (features); this * should be an Eigen::Matrix type whose number of rows should match the * length of n and whose number of columns should match the length of beta * @tparam T_beta type of the weight vector; - * this can also be a single value; + * this can also be a scalar; * @tparam T_alpha type of the intercept; - * this should be a single value; + * this should be a scalar; * @tparam T_precision type of the (positive) precision vector phi; - * this can also be a single value; + * this can also be a scalar; * @param n failures count vector parameter * @param x design matrix * @param beta weight vector @@ -54,7 +53,7 @@ namespace stan { typename return_type::type neg_binomial_2_log_glm_lpmf(const T_n &n, const T_x &x, const T_beta &beta, const T_alpha &alpha, const T_precision &phi) { - static const std::string function = "neg_binomial_2_log_glm_lpmf"; + static const char* function = "neg_binomial_2_log_glm_lpmf"; typedef typename stan::partials_return_type::type T_partials_return; From 178adb815dbf41fd915c8601cb40073c289911cd Mon Sep 17 00:00:00 2001 From: MatthijsV Date: Mon, 13 Nov 2017 16:19:42 -0500 Subject: [PATCH 13/14] fixed passing by reference in lambdas --- .../mat/prob/neg_binomial_2_log_glm_lpmf.hpp | 22 ++++++++++--------- .../prob/neg_binomial_2_log_glm_lpmf_test.cpp | 6 +++-- 2 files changed, 16 insertions(+), 12 deletions(-) diff --git a/stan/math/prim/mat/prob/neg_binomial_2_log_glm_lpmf.hpp b/stan/math/prim/mat/prob/neg_binomial_2_log_glm_lpmf.hpp index b6938071e71..7d8ec443bbb 100644 --- a/stan/math/prim/mat/prob/neg_binomial_2_log_glm_lpmf.hpp +++ b/stan/math/prim/mat/prob/neg_binomial_2_log_glm_lpmf.hpp @@ -115,7 +115,8 @@ namespace stan { Array log_phi = phi_arr.log(); Array logsumexp_eta_logphi = theta_dbl.binaryExpr(log_phi, - [](T_partials_return xx, T_partials_return yy) { + [](const T_partials_return& xx, + const T_partials_return& yy) { return log_sum_exp(xx, yy); }); Array n_plus_phi = n_arr + phi_arr; @@ -124,16 +125,16 @@ namespace stan { if (include_summand::value) { logp -= (n_arr + Array::Ones(N, 1)) - .unaryExpr([](T_partials_return xx) { + .unaryExpr([](const T_partials_return& xx) { return lgamma(xx); }).sum(); } if (include_summand::value) { - logp += (phi_arr.binaryExpr(phi_arr, [](T_partials_return xx, - T_partials_return yy) { + logp += (phi_arr.binaryExpr(phi_arr, [](const T_partials_return& xx, + const T_partials_return& yy) { return multiply_log(xx, yy); }) - - phi_arr.unaryExpr([](T_partials_return xx) { + - phi_arr.unaryExpr([](const T_partials_return& xx) { return lgamma(xx); })).sum(); } @@ -142,7 +143,7 @@ namespace stan { if (include_summand::value) logp += (n_arr * theta_dbl).sum(); if (include_summand::value) { - logp += n_plus_phi.unaryExpr([](T_partials_return xx) { + logp += n_plus_phi.unaryExpr([](const T_partials_return& xx) { return lgamma(xx); }).sum(); } @@ -170,10 +171,11 @@ namespace stan { if (!is_constant_struct::value) { ops_partials.edge4_.partials_ = (Array::Ones(N, 1) - n_plus_phi / (theta_dbl.exp() + phi_arr) + log_phi - - logsumexp_eta_logphi - phi_arr.unaryExpr([](T_partials_return xx) { - return digamma(xx); - }) - + n_plus_phi.unaryExpr([](T_partials_return xx) { + - logsumexp_eta_logphi + - phi_arr.unaryExpr([](const T_partials_return& xx) { + return digamma(xx); + }) + + n_plus_phi.unaryExpr([](const T_partials_return& xx) { return digamma(xx); })).matrix(); } diff --git a/test/unit/math/rev/mat/prob/neg_binomial_2_log_glm_lpmf_test.cpp b/test/unit/math/rev/mat/prob/neg_binomial_2_log_glm_lpmf_test.cpp index d1edf76809d..e5bcf224811 100644 --- a/test/unit/math/rev/mat/prob/neg_binomial_2_log_glm_lpmf_test.cpp +++ b/test/unit/math/rev/mat/prob/neg_binomial_2_log_glm_lpmf_test.cpp @@ -114,8 +114,9 @@ typedef std::chrono::high_resolution_clock::time_point TimeVar; #define timeNow() std::chrono::high_resolution_clock::now() TEST(ProbDistributionsNegBinomial2LogGLM, glm_matches_neg_binomial_2_log_speed) { - const int R = 30000; - const int C = 1000; + for (size_t exponent = 7; exponent < 10; exponent++) { + const int R = 10000; + const int C = std::pow(3, exponent); Matrix n(R, 1); for (size_t i = 0; i < R; i++) { @@ -161,4 +162,5 @@ TEST(ProbDistributionsNegBinomial2LogGLM, glm_matches_neg_binomial_2_log_speed) std::cout << "Existing Primitives:" << std::endl << T1 << std::endl << "New Primitives:" << std::endl << T2 << std::endl; } +} */ From f168bfdf39e8baad8f4ef24e38a7527a0c54fc38 Mon Sep 17 00:00:00 2001 From: MatthijsV Date: Tue, 14 Nov 2017 14:20:58 -0500 Subject: [PATCH 14/14] minor --- stan/math/prim/mat/prob/bernoulli_logit_glm_lpmf.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/stan/math/prim/mat/prob/bernoulli_logit_glm_lpmf.hpp b/stan/math/prim/mat/prob/bernoulli_logit_glm_lpmf.hpp index cc8073e7a3a..adb0371b0a6 100644 --- a/stan/math/prim/mat/prob/bernoulli_logit_glm_lpmf.hpp +++ b/stan/math/prim/mat/prob/bernoulli_logit_glm_lpmf.hpp @@ -1,6 +1,7 @@ #ifndef STAN_MATH_PRIM_MAT_PROB_BERNOULLI_LOGIT_GLM_LPMF_HPP #define STAN_MATH_PRIM_MAT_PROB_BERNOULLI_LOGIT_GLM_LPMF_HPP + #include #include #include @@ -11,6 +12,8 @@ #include #include #include +#include +#include #include #include