-
-
Notifications
You must be signed in to change notification settings - Fork 187
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.
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
Feature/issue 962 bivar norm #1009
Closed
Closed
Changes from 1 commit
Commits
Show all changes
41 commits
Select commit
Hold shift + click to select a range
f1964a3
Binormal integral using Owens T functions and test
rtrangucci 8976a72
Binormal tests passing
rtrangucci 26c5a9c
More accurate numeric evaluation of bivariate normal integral
rtrangucci 21f2333
Rev tests
rtrangucci c6f9d34
gradient tests
rtrangucci 5b7f6f8
Varis on stack test
rtrangucci d85f493
Fwd tests added
rtrangucci b427d4f
More efficient gradient calculations
rtrangucci 611636b
Accuracy tested vs. MM acceptable
rtrangucci e459c7c
Generalized function specification, added testing to test higher-orde…
rtrangucci ad8dbbc
Removed old test
rtrangucci f078373
Doc for lcdf
rtrangucci d530927
Binormal lpdf added
rtrangucci a14dcd4
Fixed typo in lpdf referring to sigma
rtrangucci e01221e
Rewrite of lcdf and tests to take 2-vector y args
rtrangucci 53f31ff
Added all permutations of tests for vectorized calls
rtrangucci 36e85c9
Added finite diff testing to reverse mode AD
rtrangucci 257b1fe
Higher order tests and mixed mode tests
rtrangucci 2b8a660
Removed old files
rtrangucci 2a36726
Removed old test file
rtrangucci abf12bc
Reverting cdf log changes
rtrangucci 82aaaff
Adding lcdf to prim mat header file
rtrangucci 315b830
Doc added
rtrangucci d4c9ad4
Added includes for header tests
rtrangucci 1a88b14
Lint changes
rtrangucci b358949
Lint fixes and fixes test errors
rtrangucci ff53f15
Merge commit '14b8c56e7a6d4f21b5346e255e63e3df09517e76' into HEAD
yashikno 4199f64
[Jenkins] auto-formatting by clang-format version 6.0.0 (tags/google/…
stan-buildbot a2047c3
Lint fixes
rtrangucci 464a9b7
Inlined double std_binormal_integral
rtrangucci 786e75f
Fixes sign-comparison warnings
rtrangucci 0ea95ec
Updated gradients to handle boundary conditions
rtrangucci a6dd957
Merge commit '09a2aa82a06c4f702dd8583e6683d461194dcdc1' into HEAD
yashikno 132dba6
[Jenkins] auto-formatting by clang-format version 5.0.2-svn328729-1~e…
stan-buildbot 73a2c3e
Updates to tests for gradients at boundaries
rtrangucci 1963d81
Merge branch 'feature/issue-962-bivar-norm' of ssh://github.com/stan-…
rtrangucci 5c99911
[Jenkins] auto-formatting by clang-format version 5.0.2-svn328729-1~e…
stan-buildbot e8ad567
Fixes errant cout, adds square
rtrangucci c5cb5a2
[Jenkins] auto-formatting by clang-format version 5.0.2-svn328729-1~e…
stan-buildbot 36f16f1
Fix rho initialization problem
rtrangucci 0c0ea11
Merge branch 'develop' into feature/issue-962-bivar-norm
syclik File filter
Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Updated gradients to handle boundary conditions
- Loading branch information
commit 0ea95ec23220c93234ab8b866a05be4b399b7866
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -8,6 +8,7 @@ | |
#include <stan/math/prim/mat/meta/operands_and_partials.hpp> | ||
#include <stan/math/prim/mat/meta/vector_seq_view.hpp> | ||
#include <stan/math/prim/scal/err/check_consistent_sizes.hpp> | ||
#include <stan/math/prim/scal/err/check_finite.hpp> | ||
#include <stan/math/prim/scal/err/check_not_nan.hpp> | ||
#include <stan/math/prim/scal/err/check_bounded.hpp> | ||
#include <stan/math/prim/scal/fun/size_zero.hpp> | ||
|
@@ -22,6 +23,7 @@ | |
#include <stan/math/prim/scal/prob/std_normal_lpdf.hpp> | ||
#include <stan/math/prim/scal/fun/std_binormal_integral.hpp> | ||
#include <stan/math/prim/scal/fun/inv.hpp> | ||
#include <stan/math/prim/scal/fun/value_of_rec.hpp> | ||
#include <stan/math/prim/scal/meta/max_size_mvt.hpp> | ||
#include <stan/math/prim/mat/meta/value_type.hpp> | ||
|
||
|
@@ -65,6 +67,7 @@ typename return_type<T_y, T_rho>::type std_binormal_lcdf(const T_y& y, | |
using std::log; | ||
using std::max; | ||
using std::sqrt; | ||
using stan::math::value_of_rec; | ||
|
||
check_bounded(function, "Correlation parameter", rho, -1.0, 1.0); | ||
if (stan::is_vector_like<T_y_child_type>::value | ||
|
@@ -109,32 +112,44 @@ typename return_type<T_y, T_rho>::type std_binormal_lcdf(const T_y& y, | |
cdf_log += log(cdf_); | ||
|
||
if (contains_nonconstant_struct<T_y, T_rho>::value) { | ||
const T_partials_return inv_cdf_ | ||
= cdf_ > 0 ? inv(cdf_) : std::numeric_limits<double>::infinity(); | ||
const T_partials_return inv_cdf_ = 1.0 / cdf_; | ||
const T_partials_return one_minus_rho_sq = (1 + rho_dbl) * (1 - rho_dbl); | ||
const T_partials_return sqrt_one_minus_rho_sq = sqrt(one_minus_rho_sq); | ||
const T_partials_return rho_times_y2 = rho_dbl * y2_dbl; | ||
const T_partials_return y1_minus_rho_times_y2 = y1_dbl - rho_times_y2; | ||
const bool y1_isfinite = std::isfinite(value_of_rec(y1_dbl)); | ||
const bool y2_isfinite = std::isfinite(value_of_rec(y2_dbl)); | ||
const bool y1_y2_arefinite = y1_isfinite && y2_isfinite; | ||
const bool rho_lt_one = fabs(rho_dbl) < 1; | ||
if (!is_constant_struct<T_y>::value) { | ||
ops_partials.edge1_.partials_vec_[n](0) | ||
+= cdf_ > 0 && cdf_ < 1 ? inv_cdf_ * exp(std_normal_lpdf(y1_dbl)) | ||
* Phi((y2_dbl - rho_dbl * y1_dbl) | ||
/ sqrt_one_minus_rho_sq) | ||
: cdf_ > 0 ? 1 : 0; | ||
+= cdf_ > 0 && rho_lt_one && y1_y2_arefinite | ||
? inv_cdf_ * exp(std_normal_lpdf(y1_dbl)) | ||
* Phi((y2_dbl - rho_dbl * y1_dbl) | ||
/ sqrt_one_minus_rho_sq) | ||
: (!y2_isfinite && y1_isfinite && cdf_ > 0) | ||
|| (rho_dbl == 1 && y1_dbl < y2_dbl && cdf_ > 0) | ||
|| (rho_dbl == -1 && y2_dbl > -y1_dbl && cdf_ > 0) | ||
? inv_cdf_ * exp(std_normal_lpdf(y1_dbl)) : | ||
cdf_ > 0 ? 0 : inv_cdf_; | ||
ops_partials.edge1_.partials_vec_[n](1) | ||
+= cdf_ > 0 && cdf_ < 1 | ||
? inv_cdf_ * exp(std_normal_lpdf(y2_dbl)) | ||
+= cdf_ > 0 && rho_lt_one && y1_y2_arefinite | ||
? inv_cdf_ * exp(std_normal_lpdf(y2_dbl)) | ||
* Phi(y1_minus_rho_times_y2 / sqrt_one_minus_rho_sq) | ||
: cdf_ > 0 ? 1 : 0; | ||
: (!y1_isfinite && y2_isfinite && cdf_ > 0) | ||
|| (rho_dbl == 1 && y2_dbl < y1_dbl && cdf_ > 0) | ||
|| (rho_dbl == -1 && y2_dbl > -y1_dbl && cdf_ > 0) | ||
? inv_cdf_ * exp(std_normal_lpdf(y2_dbl)) : | ||
cdf_ > 0 ? 0 : inv_cdf_; | ||
} | ||
if (!is_constant_struct<T_rho>::value) | ||
ops_partials.edge2_.partials_[n] | ||
+= cdf_ > 0 && cdf_ < 1 | ||
+= cdf_ > 0 && y1_y2_arefinite && rho_lt_one | ||
? inv_cdf_ * 0.5 / (stan::math::pi() * sqrt_one_minus_rho_sq) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Arithmetic is left associative, so all the constants should be pushed together, here |
||
* exp(-0.5 / one_minus_rho_sq * y1_minus_rho_times_y2 | ||
* y1_minus_rho_times_y2 | ||
- 0.5 * y2_dbl * y2_dbl) | ||
: cdf_ > 0 ? 1 : 0; | ||
- 0.5 * y2_dbl * y2_dbl) : | ||
cdf_ > 0 ? 0 : inv_cdf_; | ||
} | ||
} | ||
return ops_partials.build(cdf_log); | ||
|
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Stan's
is_finite
will work directly ony1
, I think---I'm not sure why there's a recursive value_of here.There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I don't think that Stan has an is_finite implementation that works for arbitrary types. Here's an implementation for
vars
: https://github.com/stan-dev/math/blob/502487c511594ccb93eb979df5b8fe163becb417/stan/math/rev/scal/fun/boost_isfinite.hpp but as far as I can tell this is it. I have thevalue_of_rec
to handle the nested AD types cases.There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thanks for looking. That's a separate shortcoming of the math lib you don't need to address here.