Skip to content

Reuse intermediate computations in distributions part 2 #1752

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
mcol wants to merge 9 commits into from
Closed

Reuse intermediate computations in distributions part 2 #1752

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
41af3de
Reuse intermediate computations in chi_square and [scaled_]inv_chi_sq…
mcol Feb 28, 2020
0abd9f1
Reuse intermediate computations in gamma and inv_gamma
mcol Feb 28, 2020
e851f2c
Move size_zero checks after arguments checks and further cleanups
mcol Mar 2, 2020
cae9b46
Merge commit 'f0341173706be97c95f5b76af06f27a0fa597791' into HEAD
yashikno Mar 2, 2020
0680496
[Jenkins] auto-formatting by clang-format version 6.0.0 (tags/google/…
stan-buildbot Mar 2, 2020
b597f07
Merge branch 'develop' into cleanup/1230-better-intermediates-prob-p2
mcol Mar 6, 2020
6cd5d46
Remove duplicated using statements introduced by a bad merge
mcol Mar 6, 2020
356de66
Merge remote-tracking branch 'origin' into cleanup/1230-better-interm…
bbbales2 Jul 19, 2020
60a574d
Add using std::double to gamma cdf (Issue #1763)
bbbales2 Jul 19, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
38 changes: 19 additions & 19 deletions stan/math/prim/prob/chi_square_cdf.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -53,6 +53,7 @@ return_type_t<T_y, T_dof> chi_square_cdf(const T_y& y, const T_dof& nu) {

scalar_seq_view<T_y> y_vec(y);
scalar_seq_view<T_dof> nu_vec(nu);
size_t size_nu = stan::math::size(nu);
size_t N = max_size(y, nu);

// Explicit return for extreme values
Expand All @@ -63,16 +64,17 @@ return_type_t<T_y, T_dof> chi_square_cdf(const T_y& y, const T_dof& nu) {
}
}

VectorBuilder<!is_constant_all<T_y, T_dof>::value, T_partials_return, T_dof>
tgamma_vec(size_nu);
VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
gamma_vec(size(nu));
VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
digamma_vec(size(nu));

if (!is_constant_all<T_dof>::value) {
for (size_t i = 0; i < stan::math::size(nu); i++) {
const T_partials_return alpha_dbl = value_of(nu_vec[i]) * 0.5;
gamma_vec[i] = tgamma(alpha_dbl);
digamma_vec[i] = digamma(alpha_dbl);
digamma_vec(size_nu);
if (!is_constant_all<T_y, T_dof>::value) {
for (size_t i = 0; i < size_nu; i++) {
const T_partials_return half_nu_dbl = 0.5 * value_of(nu_vec[i]);
tgamma_vec[i] = tgamma(half_nu_dbl);
if (!is_constant_all<T_dof>::value) {
digamma_vec[i] = digamma(half_nu_dbl);
}
}
}

Expand All @@ -83,23 +85,21 @@ return_type_t<T_y, T_dof> chi_square_cdf(const T_y& y, const T_dof& nu) {
continue;
}

const T_partials_return y_dbl = value_of(y_vec[n]);
const T_partials_return alpha_dbl = value_of(nu_vec[n]) * 0.5;
const T_partials_return beta_dbl = 0.5;

const T_partials_return Pn = gamma_p(alpha_dbl, beta_dbl * y_dbl);
const T_partials_return half_nu_dbl = 0.5 * value_of(nu_vec[n]);
const T_partials_return half_y_dbl = 0.5 * value_of(y_vec[n]);
const T_partials_return Pn = gamma_p(half_nu_dbl, half_y_dbl);

cdf *= Pn;

if (!is_constant_all<T_y>::value) {
ops_partials.edge1_.partials_[n] += beta_dbl * exp(-beta_dbl * y_dbl)
* pow(beta_dbl * y_dbl, alpha_dbl - 1)
/ tgamma(alpha_dbl) / Pn;
ops_partials.edge1_.partials_[n] += 0.5 * exp(-half_y_dbl)
* pow(half_y_dbl, half_nu_dbl - 1)
/ (tgamma_vec[n] * Pn);
}
if (!is_constant_all<T_dof>::value) {
ops_partials.edge2_.partials_[n]
-= 0.5
* grad_reg_inc_gamma(alpha_dbl, beta_dbl * y_dbl, gamma_vec[n],
* grad_reg_inc_gamma(half_nu_dbl, half_y_dbl, tgamma_vec[n],
digamma_vec[n])
/ Pn;
}
Expand All @@ -111,7 +111,7 @@ return_type_t<T_y, T_dof> chi_square_cdf(const T_y& y, const T_dof& nu) {
}
}
if (!is_constant_all<T_dof>::value) {
for (size_t n = 0; n < stan::math::size(nu); ++n) {
for (size_t n = 0; n < size_nu; ++n) {
ops_partials.edge2_.partials_[n] *= cdf;
}
}
Expand Down
49 changes: 24 additions & 25 deletions stan/math/prim/prob/chi_square_lccdf.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,53 +55,52 @@ return_type_t<T_y, T_dof> chi_square_lccdf(const T_y& y, const T_dof& nu) {

scalar_seq_view<T_y> y_vec(y);
scalar_seq_view<T_dof> nu_vec(nu);
size_t size_nu = stan::math::size(nu);
size_t N = max_size(y, nu);

// Explicit return for extreme values
// The gradients are technically ill-defined, but treated as zero
for (size_t i = 0; i < stan::math::size(y); i++) {
if (value_of(y_vec[i]) == 0) {
const T_partials_return y_dbl = value_of(y_vec[i]);
if (y_dbl == 0) {
return ops_partials.build(0.0);
}
if (y_dbl == INFTY) {
return ops_partials.build(NEGATIVE_INFTY);
}
}

VectorBuilder<true, T_partials_return, T_dof> half_nu(size_nu);
VectorBuilder<!is_constant_all<T_y, T_dof>::value, T_partials_return, T_dof>
tgamma_vec(size_nu);
VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
gamma_vec(size(nu));
VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
digamma_vec(size(nu));

if (!is_constant_all<T_dof>::value) {
for (size_t i = 0; i < stan::math::size(nu); i++) {
const T_partials_return alpha_dbl = value_of(nu_vec[i]) * 0.5;
gamma_vec[i] = tgamma(alpha_dbl);
digamma_vec[i] = digamma(alpha_dbl);
digamma_vec(size_nu);
for (size_t i = 0; i < size_nu; i++) {
const T_partials_return half_nu_dbl = 0.5 * value_of(nu_vec[i]);
half_nu[i] = half_nu_dbl;
if (!is_constant_all<T_y, T_dof>::value) {
tgamma_vec[i] = tgamma(half_nu[i]);
}
if (!is_constant_all<T_dof>::value) {
digamma_vec[i] = digamma(half_nu_dbl);
}
}

for (size_t n = 0; n < N; n++) {
// Explicit results for extreme values
// The gradients are technically ill-defined, but treated as zero
if (value_of(y_vec[n]) == INFTY) {
return ops_partials.build(negative_infinity());
}

const T_partials_return y_dbl = value_of(y_vec[n]);
const T_partials_return alpha_dbl = value_of(nu_vec[n]) * 0.5;
const T_partials_return beta_dbl = 0.5;

const T_partials_return Pn = gamma_q(alpha_dbl, beta_dbl * y_dbl);
const T_partials_return half_y_dbl = 0.5 * value_of(y_vec[n]);
const T_partials_return Pn = gamma_q(half_nu[n], half_y_dbl);

ccdf_log += log(Pn);

if (!is_constant_all<T_y>::value) {
ops_partials.edge1_.partials_[n] -= beta_dbl * exp(-beta_dbl * y_dbl)
* pow(beta_dbl * y_dbl, alpha_dbl - 1)
/ tgamma(alpha_dbl) / Pn;
ops_partials.edge1_.partials_[n] -= 0.5 * exp(-half_y_dbl)
* pow(half_y_dbl, half_nu[n] - 1)
/ (tgamma_vec[n] * Pn);
}
if (!is_constant_all<T_dof>::value) {
ops_partials.edge2_.partials_[n]
+= 0.5
* grad_reg_inc_gamma(alpha_dbl, beta_dbl * y_dbl, gamma_vec[n],
* grad_reg_inc_gamma(half_nu[n], half_y_dbl, tgamma_vec[n],
digamma_vec[n])
/ Pn;
}
Expand Down
51 changes: 25 additions & 26 deletions stan/math/prim/prob/chi_square_lcdf.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,53 +55,52 @@ return_type_t<T_y, T_dof> chi_square_lcdf(const T_y& y, const T_dof& nu) {

scalar_seq_view<T_y> y_vec(y);
scalar_seq_view<T_dof> nu_vec(nu);
size_t size_nu = stan::math::size(nu);
size_t N = max_size(y, nu);

// Explicit return for extreme values
// The gradients are technically ill-defined, but treated as zero
for (size_t i = 0; i < stan::math::size(y); i++) {
if (value_of(y_vec[i]) == 0) {
return ops_partials.build(negative_infinity());
const T_partials_return y_dbl = value_of(y_vec[i]);
if (y_dbl == 0) {
return ops_partials.build(NEGATIVE_INFTY);
}
if (y_dbl == INFTY) {
return ops_partials.build(0.0);
}
}

VectorBuilder<true, T_partials_return, T_dof> half_nu(size_nu);
VectorBuilder<!is_constant_all<T_y, T_dof>::value, T_partials_return, T_dof>
tgamma_vec(size_nu);
VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
gamma_vec(size(nu));
VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
digamma_vec(size(nu));

if (!is_constant_all<T_dof>::value) {
for (size_t i = 0; i < stan::math::size(nu); i++) {
const T_partials_return alpha_dbl = value_of(nu_vec[i]) * 0.5;
gamma_vec[i] = tgamma(alpha_dbl);
digamma_vec[i] = digamma(alpha_dbl);
digamma_vec(size_nu);
for (size_t i = 0; i < size_nu; i++) {
const T_partials_return half_nu_dbl = 0.5 * value_of(nu_vec[i]);
half_nu[i] = half_nu_dbl;
if (!is_constant_all<T_y, T_dof>::value) {
tgamma_vec[i] = tgamma(half_nu[i]);
}
if (!is_constant_all<T_dof>::value) {
digamma_vec[i] = digamma(half_nu_dbl);
}
}

for (size_t n = 0; n < N; n++) {
// Explicit results for extreme values
// The gradients are technically ill-defined, but treated as zero
if (value_of(y_vec[n]) == INFTY) {
return ops_partials.build(0.0);
}

const T_partials_return y_dbl = value_of(y_vec[n]);
const T_partials_return alpha_dbl = value_of(nu_vec[n]) * 0.5;
const T_partials_return beta_dbl = 0.5;

const T_partials_return Pn = gamma_p(alpha_dbl, beta_dbl * y_dbl);
const T_partials_return half_y_dbl = 0.5 * value_of(y_vec[n]);
const T_partials_return Pn = gamma_p(half_nu[n], half_y_dbl);

cdf_log += log(Pn);

if (!is_constant_all<T_y>::value) {
ops_partials.edge1_.partials_[n] += beta_dbl * exp(-beta_dbl * y_dbl)
* pow(beta_dbl * y_dbl, alpha_dbl - 1)
/ tgamma(alpha_dbl) / Pn;
ops_partials.edge1_.partials_[n] += 0.5 * exp(-half_y_dbl)
* pow(half_y_dbl, half_nu[n] - 1)
/ (tgamma_vec[n] * Pn);
}
if (!is_constant_all<T_dof>::value) {
ops_partials.edge2_.partials_[n]
-= 0.5
* grad_reg_inc_gamma(alpha_dbl, beta_dbl * y_dbl, gamma_vec[n],
* grad_reg_inc_gamma(half_nu[n], half_y_dbl, tgamma_vec[n],
digamma_vec[n])
/ Pn;
}
Expand Down
18 changes: 0 additions & 18 deletions stan/math/prim/prob/chi_square_log.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,25 +8,7 @@ namespace stan {
namespace math {

/** \ingroup prob_dists
* The log of a chi-squared density for y with the specified
* degrees of freedom parameter.
* The degrees of freedom parameter must be greater than 0.
* y must be greater than or equal to 0.
*
\f{eqnarray*}{
y &\sim& \chi^2_\nu \\
\log (p (y \, |\, \nu)) &=& \log \left( \frac{2^{-\nu / 2}}{\Gamma (\nu / 2)}
y^{\nu / 2 - 1} \exp^{- y / 2} \right) \\
&=& - \frac{\nu}{2} \log(2) - \log (\Gamma (\nu / 2)) + (\frac{\nu}{2} - 1)
\log(y) - \frac{y}{2} \\ & & \mathrm{ where } \; y \ge 0 \f}
*
* @deprecated use <code>chi_square_lpdf</code>
* @param y A scalar variable.
* @param nu Degrees of freedom.
* @throw std::domain_error if nu is not greater than or equal to 0
* @throw std::domain_error if y is not greater than or equal to 0.
* @tparam T_y Type of scalar.
* @tparam T_dof Type of degrees of freedom.
*/
template <bool propto, typename T_y, typename T_dof>
return_type_t<T_y, T_dof> chi_square_log(const T_y& y, const T_dof& nu) {
Expand Down
52 changes: 24 additions & 28 deletions stan/math/prim/prob/chi_square_lpdf.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
#include <stan/math/prim/err.hpp>
#include <stan/math/prim/fun/constants.hpp>
#include <stan/math/prim/fun/digamma.hpp>
#include <stan/math/prim/fun/inv.hpp>
#include <stan/math/prim/fun/lgamma.hpp>
#include <stan/math/prim/fun/log.hpp>
#include <stan/math/prim/fun/max_size.hpp>
Expand Down Expand Up @@ -60,64 +61,59 @@ return_type_t<T_y, T_dof> chi_square_lpdf(const T_y& y, const T_dof& nu) {

scalar_seq_view<T_y> y_vec(y);
scalar_seq_view<T_dof> nu_vec(nu);
size_t size_y = stan::math::size(y);
size_t size_nu = stan::math::size(nu);
size_t N = max_size(y, nu);

for (size_t n = 0; n < stan::math::size(y); n++) {
for (size_t n = 0; n < size_y; n++) {
if (value_of(y_vec[n]) < 0) {
return LOG_ZERO;
}
}

VectorBuilder<include_summand<propto, T_y, T_dof>::value, T_partials_return,
T_y>
log_y(size(y));
for (size_t i = 0; i < stan::math::size(y); i++) {
log_y[i] = log(value_of(y_vec[i]));
}

VectorBuilder<true, T_partials_return, T_y> log_y(size_y);
VectorBuilder<include_summand<propto, T_y>::value, T_partials_return, T_y>
inv_y(size(y));
for (size_t i = 0; i < stan::math::size(y); i++) {
inv_y(size_y);
for (size_t i = 0; i < size_y; i++) {
const T_partials_return y_dbl = value_of(y_vec[i]);
log_y[i] = log(y_dbl);
if (include_summand<propto, T_y>::value) {
inv_y[i] = 1.0 / value_of(y_vec[i]);
inv_y[i] = inv(y_dbl);
}
}

VectorBuilder<include_summand<propto, T_dof>::value, T_partials_return, T_dof>
lgamma_half_nu(size(nu));
lgamma_half_nu(size_nu);
VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
digamma_half_nu_over_two(size(nu));

for (size_t i = 0; i < stan::math::size(nu); i++) {
T_partials_return half_nu = 0.5 * value_of(nu_vec[i]);
if (include_summand<propto, T_dof>::value) {
digamma_half_nu(size_nu);
if (include_summand<propto, T_dof>::value) {
for (size_t i = 0; i < size_nu; i++) {
T_partials_return half_nu = 0.5 * value_of(nu_vec[i]);
lgamma_half_nu[i] = lgamma(half_nu);
}
if (!is_constant_all<T_dof>::value) {
digamma_half_nu_over_two[i] = digamma(half_nu) * 0.5;
if (!is_constant_all<T_dof>::value) {
digamma_half_nu[i] = digamma(half_nu);
}
}
}

for (size_t n = 0; n < N; n++) {
const T_partials_return y_dbl = value_of(y_vec[n]);
const T_partials_return half_y = 0.5 * y_dbl;
const T_partials_return nu_dbl = value_of(nu_vec[n]);
const T_partials_return half_nu = 0.5 * nu_dbl;
const T_partials_return half_nu_m1 = 0.5 * nu_dbl - 1;

if (include_summand<propto, T_dof>::value) {
logp -= nu_dbl * HALF_LOG_TWO + lgamma_half_nu[n];
}
logp += (half_nu - 1.0) * log_y[n];

if (include_summand<propto, T_y>::value) {
logp -= half_y;
logp -= 0.5 * value_of(y_vec[n]);
}
logp += half_nu_m1 * log_y[n];

if (!is_constant_all<T_y>::value) {
ops_partials.edge1_.partials_[n] += (half_nu - 1.0) * inv_y[n] - 0.5;
ops_partials.edge1_.partials_[n] += half_nu_m1 * inv_y[n] - 0.5;
}
if (!is_constant_all<T_dof>::value) {
ops_partials.edge2_.partials_[n]
-= HALF_LOG_TWO + digamma_half_nu_over_two[n] - log_y[n] * 0.5;
-= 0.5 * (LOG_TWO + digamma_half_nu[n] - log_y[n]);
}
}
return ops_partials.build(logp);
Expand Down
Loading