Feature/log bessel i #640

bgoodri · 2017-10-07T17:55:16Z

Submission Checklist

Run unit tests: ./runTests.py test/unit
Run cpplint: make cpplint
Declare copyright holder and open-source license: see below

Summary:

Introduce log_modified_bessel_first_kind(T1 v, T2 z) function that uses log_sum_exp to avoid overflow and can be evaluated for non-integer v.

Intended Effect:

Not overflowing while expanding the number of situations where the Bessel I function can be used.

How to Verify:

Run the unit tests

Side Effects:

Inconsistency with other Bessel functions, which currently are restricted to integer orders and do not have log_ versions.

Fixes #639

Documentation:

I am not sure log_modified_bessel_first_kind should necessarily be exposed in the Stan language since the only place it is used in Stan is von_mises_lpdf, which only calls it with v = 0 and double-precision z. Thus, I have not implemented rev or fwd specializations because the derivative with respect to v is a lot of work for no known use case. It probably wouldn't hurt anything to expose it in the Stan language, but any derivatives would be calculated by autodiffing a truncated infinite sum.

Copyright and Licensing

Please list the copyright holder for the work you are submitting (this will be you or your assignee, such as a university or company): Trustees of Columbia University

By submitting this pull request, the copyright holder is agreeing to 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/)

by using log_modified_bessel_first_kind

bob-carpenter

Overall this looks great. The changes requested are some minor efficiency issues, some needless includes, and abit of doc on template parameters.

bob-carpenter · 2017-10-07T21:41:31Z

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

+      check_greater_or_equal("log_modified_bessel_first_kind", "v", v, -1);
+
+      using std::log;
+      using stan::math::log_sum_exp;


Symbols in the same namespace get brought in automatically---so you don't need the stan::math:: using statements. They won't hurt anything being there.

bob-carpenter · 2017-10-07T21:44:22Z

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

+      }
+
+      T log_half_z = log(0.5 * z);
+      T lgam = lgamma(v + 1.0);


Shouldn't lgam be set inside the v > -1 branch? The other branch sets it to 0.

Boost's lgamma function stupidly throws an exception if v = -1 and thus lgamma is evaluated at 0. It should be returning positive infinity, in which case the first term in the denominator of the Bessel I function is just zero. So, we avoid that, handle the m = 1 case manually, and start the infinity sum with the m = 2 case.

bob-carpenter · 2017-10-07T21:46:26Z

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

+      double m = 2;
+      T lfac = 0;
+      T old_out;
+      do {


Why isn't the test just out != old_out.

I always find the do-while convention takes me a while to get my head around, but I could get used to it.

You are right. I forgot that NaN != NaN is false.

Upon further review, I was right the first time: NaN != NaN is true because NaN == NaN is false. So, if out and old_out somehow both go NaN then the loop would never terminate.

Thanks for reviewing this. I somehow thought any comparisons involving NaN were automatically false, but you're right that != is an exception to this general rule and always returns true.

bob-carpenter · 2017-10-07T21:49:13Z

stan/math/prim/scal/prob/von_mises_lpdf.hpp

@@ -67,7 +68,7 @@ namespace stan {
        kappa_dbl[i] = value_of(kappa_vec[i]);
        if (include_summand<propto, T_scale>::value)
          log_bessel0[i]
-            = log(modified_bessel_first_kind(0, value_of(kappa_vec[i])));
+            = log_modified_bessel_first_kind(0, value_of(kappa_vec[i]));


bob-carpenter · 2017-10-07T21:52:55Z

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

+    /* Log of the modified Bessel function of the first kind,
+     * which is better known as the Bessel I function. See
+     * modified_bessel_first_kind.hpp for the function definition.
+     *


You want to add doc for the template parameters,

@tparam T1 type of the order (v) @tparam T2 type of argument (z)

bob-carpenter · 2017-10-07T21:53:59Z

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

+    template <typename T1, typename T2>
+    inline typename boost::math::tools::promote_args<T1, T2, double>::type
+    log_modified_bessel_first_kind(const T1 v, const T2 z) {
+      check_not_nan("log_modified_bessel_first_kind", "v", v);


These comments can't really use v since that's not what the user is necessarily going to have typed into their program, so it gets confusing.

It's better if you can replace v with first argument and z with second argument (order) in the error messages.

bob-carpenter · 2017-10-07T21:56:28Z

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

+      if (stan::math::is_inf(z)) return z;
+      if (z > 100) {
+        // Boost does something like this in asymptotic_bessel_i_large_x
+        T lim = (4 * square(v) + 10) / (8 * z);


This would be more efficient as (square(v) + 2.5) / (2 * z).

bob-carpenter · 2017-10-07T21:57:55Z

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

+        }
+      }
+
+      T log_half_z = log(0.5 * z);


The type can be tighter here as T2. That way, if T2 == double, this won't be overpromoted.

Same thing goes for the next one (lgam) being just T2. You only need the promoted T when v and z are involved (directly or indirectly).

I think you have to metaprogram here. If T1 or T2 is int then we want the left-hand side to be a double.

Good point about integers; it'll need to be

typename boost::math::tools::promote_args<T1>::type log_half_z

bob-carpenter · 2017-10-07T21:58:08Z

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

+          num *= mu - 25;
+          denom *= ex * 3;
+          s -= num / denom;
+          s = z - log(std::sqrt(2 * z * M_PI)) + log(s);


This needs to be using std::sqrt and then just sqrt here if z is going to be anything other than double or int.

bob-carpenter · 2017-10-07T21:59:16Z

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

+      T out;
+      if (v > -1) {
+        out = log_sum_exp(v * log_half_z - lgam, lcons - lgamma(v + 2));
+        lgam += log(v + 1.0);


The righthand side can be log1p(v) --- this may need an include and a using std::log1p if there's not a primitive stan::math::log1p.

works faster for v = 0 or v = 1

bgoodri · 2017-10-07T23:50:19Z

I think I responded to all of your comments, but in the latest push I also copied special cases from Boost when v = 0 or v = 1.

bob-carpenter · 2017-10-08T22:47:25Z

The failing test is a failed instantiation of the log-Bessel for fvar<var>, I think.

bob-carpenter · 2017-10-08T22:50:13Z

Thanks. Just the one issue with instantiation to sort out and this should be good to go. I also responded on the metaprogramming --- promote_args does what you need.

…

On Oct 7, 2017, at 7:50 PM, bgoodri ***@***.***> wrote: I think I responded to all of your comments, but in the latest push I also copied special cases from Boost when v = 0 or v = 1. — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or mute the thread.

bgoodri · 2017-10-08T23:47:58Z

I don't understand the error message. It just says over and over again that it can't do the template instantiation. But it has T1 and T2, so I don't know what is stopping it from working in fwd mode.

bob-carpenter · 2017-10-09T00:18:55Z

I'll check it out and have a look. Probably just a missing include somewhere.

…re/log_besselI

bob-carpenter · 2017-10-09T00:43:11Z

I fixed the issue with this and will push as soon as I generalize the testing a bit.

…

On Oct 8, 2017, at 7:47 PM, bgoodri ***@***.***> wrote: I don't understand the error message. It just says over and over again that it can't do the template instantiation. But it has T1 and T2, so I don't know what is stopping it from working in fwd mode. — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or mute the thread.

bob-carpenter · 2017-10-09T01:32:12Z

I pushed a fix for the compilation problem. I then plugged in the general testing framework with a test in: test/unit/math/mix/scal/fun/log_modified_bessel_first_kind_test.cpp I had to update the testing framework to deal with exceptions, but now there's an issue with the autodiff test failing. It's failing at all orders vs. finite diffs. Is this something you might expect to be way off? I didn't see any other tests for gradients---should've caught that the first time!

bgoodri · 2017-10-09T02:27:39Z

It won't autodiff correctly at v = 0 or v = 1 but other than that it should be fine.

bob-carpenter · 2017-10-09T03:43:48Z

OK, let me write a finer test for this case, then, and see if that works.

bob-carpenter · 2017-10-09T04:11:58Z

I'm pushing a version that passes tests locally for me.

There are still some known issues and then a bunch of cases where the double version throws an exception but the autodiff version does not.

  // TODO(carpenter):  need to fix boundary conditions in function
  //   50 : tolerance too tight for autodiff, but basically OK
  //   0 : known that won't work for finite diffs test
  //   -0.5  : autodiff doesn't throw
  //   stan::math::not_a_number() : autodiff doesn't throw
  //   stan::math::negative_infinity() : autodiff doesn't throw

bob-carpenter · 2017-10-09T16:07:28Z

Jenkins, retest this please.

bob-carpenter · 2017-10-09T16:08:34Z

@seantalts --- any idea what's up? My request to Jenkins to retest this didn't kick off a new testing round. Looks like the failure's a Jenkins issue of some kind.

seantalts · 2017-10-09T16:10:58Z

It did kick off a new testing round. It's just queued. I don't think the link will work until the job actually begins. You can see it in the job queue on the home page of Jenkins and you can see its description by hovering over it.

…re/log_besselI

bob-carpenter · 2017-10-09T16:59:13Z

Thanks --- they will show that mustard-colored ball when it gets registered even before it starts. Just took a while to show up, which meant given that it's a computer, I figured it was broken.

seantalts · 2017-10-09T19:42:05Z

Jenkins, retest this please.

seantalts · 2017-10-13T00:22:15Z

My bad. Jenkins, retest this please.

syclik · 2017-10-22T03:16:06Z

@bob-carpenter, can you change your review to approved? @bgoodri went through and updated the pull request.

seantalts · 2017-10-23T15:17:28Z

Hey @bgoodri, looks like there's a compiler warning in develop on Stan during test-headers that's breaking the build: http://d1m1s1b1.stat.columbia.edu:8080/job/Stan/job/develop/11/execution/node/83/log/

This wasn't caught because our current Math build doesn't test the full Stan test suite when it tests Stan Upstream. The new pipelines fix this but aren't in yet.

Can you make a PR to fix this and we can roll forward?

aaronjg · 2018-03-28T01:23:57Z

It looks like this is currently not exposed. I just wanted to mention that the log_bessel_k is used in other distributions beyond the von Misses, such as the Generalized Inverse Gaussian Distribution, and having log_bessel_i exposed could be helpful in implementing these.

bgoodri · 2018-03-28T01:48:08Z

Yeah, we need to expose it on the Stan side.

…

On Tue, Mar 27, 2018 at 9:23 PM, aaronjg ***@***.***> wrote: It looks like this is currently not exposed. I just wanted to mention that the log_bessel_k is used in other distributions beyond the von Misses, such as the Generalized Inverse Gaussian Distribution, and having log_bessel_i exposed could be helpful in implementing these. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#640 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ADOrqkdiveMv2LaAnbw8WcJzwExTw8jTks5tiuYugaJpZM4Pxb6y> .

bgoodri added 5 commits June 7, 2017 19:18

add and test log_modified_bessel_first_kind

ad17474

change to templates

f607c1e

improve tests

0960589

make von_mises_lpdf numerically stable

8aea60e

by using log_modified_bessel_first_kind

Merge branch 'develop' into feature/log_besselI

2d0cdf8

bob-carpenter requested changes Oct 7, 2017

View reviewed changes

bgoodri added 2 commits October 7, 2017 18:58

copy and test Boost's specializations

89d78e8

works faster for v = 0 or v = 1

respond to Bob's comments

e302b73

Merge branch 'develop' of https://github.com/stan-dev/math into featu…

466d2c1

…re/log_besselI

log-mod-bessel-1st compile fix; ad test exceptions

42e089f

upgrade autodiff tests, doc for log mod bessel I

b5d19d2

Merge branch 'develop' of https://github.com/stan-dev/math into featu…

e713101

…re/log_besselI

bob-carpenter approved these changes Oct 22, 2017

View reviewed changes

bob-carpenter merged commit 62a8c51 into develop Oct 22, 2017

bob-carpenter deleted the feature/log_besselI branch October 22, 2017 03:41

a-hurst mentioned this pull request Dec 11, 2018

"numeric overflow" error for models with von Mises distributions and large kappa values #1076

Open

Feature/log bessel i #640

Feature/log bessel i #640

Conversation

bgoodri commented Oct 7, 2017 • edited by bob-carpenter Loading

Submission Checklist

Summary:

Intended Effect:

How to Verify:

Side Effects:

Documentation:

Copyright and Licensing

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

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

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

bgoodri commented Oct 7, 2017

bob-carpenter commented Oct 8, 2017

bob-carpenter commented Oct 8, 2017 via email

bgoodri commented Oct 8, 2017

bob-carpenter commented Oct 9, 2017

bob-carpenter commented Oct 9, 2017 via email

bob-carpenter commented Oct 9, 2017 via email

bgoodri commented Oct 9, 2017

bob-carpenter commented Oct 9, 2017

bob-carpenter commented Oct 9, 2017

bob-carpenter commented Oct 9, 2017

bob-carpenter commented Oct 9, 2017

seantalts commented Oct 9, 2017

bob-carpenter commented Oct 9, 2017 via email

seantalts commented Oct 9, 2017

seantalts commented Oct 13, 2017

syclik commented Oct 22, 2017

seantalts commented Oct 23, 2017

aaronjg commented Mar 28, 2018

bgoodri commented Mar 28, 2018 via email

bgoodri commented Oct 7, 2017 •

edited by bob-carpenter

Loading