Feature/issue 962 bivar norm #1009

rtrangucci · 2018-08-27T22:07:18Z

Summary

Adds std_binormal_lcdf and std_binormal_integral to compute the log-CDF and the CDF of the bivariate standard normal density.

Tests

Added tests for prim, rev, fwd, mix for std_binormal_integral and tests for prim, rev, mix for std_binormal_lcdf. Gradient and higher-order tests are in rev and mix for both functions.

./runTests.py test/unit/math/prim/scal/fun/std_binormal_integral_test.cpp 
./runTests.py test/unit/math/rev/scal/fun/std_binormal_integral_test.cpp 
./runTests.py test/unit/math/fwd/scal/fun/std_binormal_integral_test.cpp 
./runTests.py test/unit/math/mix/mat/fun/std_binormal_integral_test.cpp 

./runTests.py test/unit/math/prim/mat/fun/std_binormal_lcdf_test.cpp 
./runTests.py test/unit/math/rev/mat/fun/std_binormal_lcdf_test.cpp 
./runTests.py test/unit/math/mix/mat/fun/std_binormal_lcdf_test.cpp

Side Effects

None

Checklist

Math issue Bivariate normal distribution #962
Copyright holder: Rob Trangucci

The copyright holder is typically you or your assignee, such as a university or company. By submitting this pull request, the copyright holder is agreeing to the license the submitted work under the following licenses:
- Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
- Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)
the basic tests are passing
- unit tests pass (to run, use: ./runTests.py test/unit)
- header checks pass, (make test-headers)
- docs build, (make doxygen)
- code passes the built in C++ standards checks (make cpplint)
the code is written in idiomatic C++ and changes are documented in the doxygen
the new changes are tested

…r derivatives, fixed segfaults for lcdf

…stable/2017-11-14)

…xp1~20180509124008.99 (branches/release_50)

…dev/math into feature/issue-962-bivar-norm

…xp1~20180509124008.99 (branches/release_50)

stan/math/prim/scal/fun/std_binormal_integral.hpp

…xp1~20180509124008.99 (branches/release_50)

rtrangucci · 2018-09-13T20:54:13Z

@bob-carpenter this is ready for a review when you've got some time

bob-carpenter

Thanks. This is basically OK. I'd be a lot happier if all of the cut-and-paste code in the tests could be cosolidated.

stan/math/prim/mat/prob/std_binormal_lcdf.hpp

bob-carpenter · 2018-09-14T14:00:23Z

stan/math/prim/mat/prob/std_binormal_lcdf.hpp

+      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));


Stan's is_finite will work directly on y1, I think---I'm not sure why there's a recursive value_of here.

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 the value_of_rec to handle the nested AD types cases.

Thanks for looking. That's a separate shortcoming of the math lib you don't need to address here.

bob-carpenter · 2018-09-14T14:02:16Z

stan/math/prim/mat/prob/std_binormal_lcdf.hpp

+      if (!is_constant_struct<T_rho>::value)
+        ops_partials.edge2_.partials_[n]
+            += cdf_ > 0 && y1_y2_arefinite && rho_lt_one
+                   ? inv_cdf_ * 0.5 / (stan::math::pi() * sqrt_one_minus_rho_sq)


Arithmetic is left associative, so all the constants should be pushed together, here (0.5 / ...).

bob-carpenter · 2018-09-14T14:03:56Z

stan/math/prim/mat/prob/std_binormal_lcdf.hpp

+                   ? inv_cdf_ * 0.5 / (stan::math::pi() * sqrt_one_minus_rho_sq)
+                         * exp(-0.5 / one_minus_rho_sq * y1_minus_rho_times_y2
+                                   * y1_minus_rho_times_y2
+                               - 0.5 * y2_dbl * y2_dbl)


The -0.5 should go on the outside to save a multiplication,

exp(-0.5 * (square(y1_minus_rho_times_y2) / one_minus_rho_sq + square(y2_dbl)));

Is square better to use than direct multiplication? i.e. x * x vs. square(x)?

Only if the variable can be an autodiff variable.

bob-carpenter · 2018-09-14T14:25:04Z

test/unit/math/mix/mat/prob/std_binormal_lcdf_test.cpp

+  }
+};
+
+// vector<Row Vector>, real


These should also be templatable. No big deal if you don't combine, but it's easier to maintain with fewer cut-and-paste functions.

Yeah, I'm OK templating these tests more. I did so for the copula PR, so redoing these should be easier.

bob-carpenter · 2018-09-14T14:25:51Z

test/unit/math/mix/mat/prob/std_binormal_lcdf_test.cpp

+  template <typename T>
+  inline T operator()(
+      const Eigen::Matrix<T, Eigen::Dynamic, 1>& inp_vec) const {
+    vector<Matrix<T, Dynamic, 1>> y;


These four lines of code are repeated other than for type, so they should be pullable into a function.

bob-carpenter · 2018-09-14T14:26:13Z

test/unit/math/mix/mat/prob/std_binormal_lcdf_test.cpp

+};
+
+// Row Vector, double
+struct RV_D_std_binorm_lcdf {


Consolidate vetor and row vector versions using templating.

bob-carpenter · 2018-09-14T14:27:24Z

test/unit/math/mix/mat/prob/std_binormal_lcdf_test.cpp

+  }
+};
+
+// double vector<Row Vector>, real


Basically, way too many of these. Figure out how to consolidate using templates if at all possible.

If not, document the scheme up top.

bob-carpenter · 2018-09-14T14:45:27Z

test/unit/math/prim/mat/prob/std_binormal_lcdf_test.cpp

+                       vector<Eigen::Matrix<T, Eigen::Dynamic, 1>>& y, T& rho) {
+  int cntr = 0;
+  for (int i = 0; i < N_y; ++i) {
+    Eigen::Matrix<T, Eigen::Dynamic, 1> el(2);


At the very least, these repeated pieces need to be turned into functions.

Does the y need to be passed in here as an argument since it's just being constructed? How about returning it rather than passing it in?

syclik · 2019-02-03T01:51:52Z

Closing due to inactivity.

rtrangucci added 24 commits August 13, 2018 17:49

Binormal integral using Owens T functions and test

f1964a3

Binormal tests passing

8976a72

More accurate numeric evaluation of bivariate normal integral

26c5a9c

Rev tests

21f2333

gradient tests

c6f9d34

Varis on stack test

5b7f6f8

Fwd tests added

d85f493

More efficient gradient calculations

b427d4f

Accuracy tested vs. MM acceptable

611636b

Generalized function specification, added testing to test higher-orde…

e459c7c

…r derivatives, fixed segfaults for lcdf

Removed old test

ad8dbbc

Doc for lcdf

f078373

Binormal lpdf added

d530927

Fixed typo in lpdf referring to sigma

a14dcd4

Rewrite of lcdf and tests to take 2-vector y args

e01221e

Added all permutations of tests for vectorized calls

53f31ff

Added finite diff testing to reverse mode AD

36e85c9

Higher order tests and mixed mode tests

257b1fe

Removed old files

2b8a660

Removed old test file

2a36726

Reverting cdf log changes

abf12bc

Adding lcdf to prim mat header file

82aaaff

Doc added

315b830

Added includes for header tests

d4c9ad4

rtrangucci requested review from bgoodri and bob-carpenter August 27, 2018 22:09

rtrangucci added the feature label Aug 27, 2018

rtrangucci and others added 3 commits August 27, 2018 18:29

Lint changes

1a88b14

Lint fixes and fixes test errors

b358949

Merge commit '14b8c56e7a6d4f21b5346e255e63e3df09517e76' into HEAD

ff53f15

stan-buildbot and others added 10 commits August 27, 2018 18:51

[Jenkins] auto-formatting by clang-format version 6.0.0 (tags/google/…

4199f64

…stable/2017-11-14)

Lint fixes

a2047c3

Inlined double std_binormal_integral

464a9b7

Fixes sign-comparison warnings

786e75f

Updated gradients to handle boundary conditions

0ea95ec

Merge commit '09a2aa82a06c4f702dd8583e6683d461194dcdc1' into HEAD

a6dd957

[Jenkins] auto-formatting by clang-format version 5.0.2-svn328729-1~e…

132dba6

…xp1~20180509124008.99 (branches/release_50)

Updates to tests for gradients at boundaries

73a2c3e

Merge branch 'feature/issue-962-bivar-norm' of ssh://github.com/stan-…

1963d81

…dev/math into feature/issue-962-bivar-norm

[Jenkins] auto-formatting by clang-format version 5.0.2-svn328729-1~e…

5c99911

…xp1~20180509124008.99 (branches/release_50)

bgoodri reviewed Sep 6, 2018

View reviewed changes

stan/math/prim/scal/fun/std_binormal_integral.hpp Outdated Show resolved Hide resolved

bgoodri reviewed Sep 6, 2018

View reviewed changes

stan/math/prim/scal/fun/std_binormal_integral.hpp Outdated Show resolved Hide resolved

rtrangucci and others added 3 commits September 6, 2018 21:09

Fixes errant cout, adds square

e8ad567

[Jenkins] auto-formatting by clang-format version 5.0.2-svn328729-1~e…

c5cb5a2

…xp1~20180509124008.99 (branches/release_50)

Fix rho initialization problem

36f16f1

bob-carpenter requested changes Sep 14, 2018

View reviewed changes

Merge branch 'develop' into feature/issue-962-bivar-norm

0c0ea11

syclik closed this Feb 3, 2019

syclik deleted the feature/issue-962-bivar-norm branch September 15, 2022 14:10

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Feature/issue 962 bivar norm #1009

Feature/issue 962 bivar norm #1009

rtrangucci commented Aug 27, 2018 •

edited

Loading

rtrangucci commented Sep 13, 2018

bob-carpenter left a comment

bob-carpenter Sep 14, 2018

rtrangucci Oct 2, 2018

bob-carpenter Oct 2, 2018

bob-carpenter Sep 14, 2018

bob-carpenter Sep 14, 2018

rtrangucci Oct 2, 2018 •

edited

Loading

bob-carpenter Oct 2, 2018

bob-carpenter Sep 14, 2018

rtrangucci Sep 15, 2018

bob-carpenter Sep 14, 2018

bob-carpenter Sep 14, 2018

bob-carpenter Sep 14, 2018

bob-carpenter Sep 14, 2018

syclik commented Feb 3, 2019

Feature/issue 962 bivar norm #1009

Feature/issue 962 bivar norm #1009

Conversation

rtrangucci commented Aug 27, 2018 • edited Loading

Summary

Tests

Side Effects

Checklist

rtrangucci commented Sep 13, 2018

bob-carpenter left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

rtrangucci Oct 2, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

syclik commented Feb 3, 2019

rtrangucci commented Aug 27, 2018 •

edited

Loading

rtrangucci Oct 2, 2018 •

edited

Loading