Fix accidental precision loss (#814)

[danluu]

Submission Checklist

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

Summary:

There's accidental precision loss due to the use of integer instead of floating point divide. The offending function also does redundant loads and possibly a redundant divide, even with optimizations enabled. See #814 for more details and discussion.

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):

Dan Luu

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/)

[bob-carpenter]

Thanks.

We can come back and clean up the potential catastrophic arithmetic and raw calls to alloc in another pass.

[bob-carpenter]

This should get a test that fails in the previous release and passes now. We want to try to stop regressions on these things. This function was presumably missing good tests or we would've caught that case.

[danluu]

@bob-carpenter Where should the test live in the test hierarchy? I don't think I understand the directory structure.

[danluu]

Oh wait, I see that there's unit/math/rev/mat/fun/variance_test.cpp, so perhaps there.

[bob-carpenter]

Yes, that's where the gradient test for vector and row vector inputs should go for now.

@seantalts is going to collapse the prim / arr / mat distinction soon---everyone's on board with that and the current division of directories is just too painful.

[bbbales2]

@danluu I can finish up this pull request if you don't have time. This looks good and we should get it in.

This is test code that fails with the current implementation, but works in the new one (I basically wrote up the example in the issue #814).

It should go in test/unit/math/rev/mat/fun/variance_test.cpp

TEST(AgradRevMatrix, variance_avoid_precision_loss) {
  using stan::math::variance;
  using stan::math::vector_d;
  using stan::math::vector_v;

  vector_v v(4);
  v << 1.0, 2.0, 4.0, 5.0;
  stan::math::var output = variance(v);
  output.grad();

  EXPECT_FLOAT_EQ(10.0 / 3.0, output.val());

  for (int i = 0; i < v.size(); ++i) {
    EXPECT_FLOAT_EQ(2.0 * (v(i).val() - 3.0) / 3.0, v(i).adj());
  }
}

[danluu]

Sorry, got really busy recently. I certainly don't mind finishing it, but I'm not sure when I'd get around to it. If you want to finish it up, please go ahead :-).

[bbbales2]

@danluu Hahaha, it worked! I saw @seantalts do this, but apparently since I have permissions on this repo, I also get free permissions on your fork! So I merged in the new develop and added a test to your branch.

That probably violates all forms of repo etiquette. I probably should have asked first haha.

[bbbales2]

@bob-carpenter Test added, and I merged in the latest develop. This should be good to go