From 1dd2dcac15e87f629c30bf460a648668debbf309 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sat, 20 Jan 2024 17:39:02 +0200 Subject: [PATCH 1/2] Use bessel identity to improve grad calc --- stan/math/prim/prob/von_mises_lpdf.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/prim/prob/von_mises_lpdf.hpp b/stan/math/prim/prob/von_mises_lpdf.hpp index 7e4bfb7d2a4..f8e5e1aa303 100644 --- a/stan/math/prim/prob/von_mises_lpdf.hpp +++ b/stan/math/prim/prob/von_mises_lpdf.hpp @@ -81,7 +81,7 @@ return_type_t von_mises_lpdf(T_y const& y, T_loc const& mu, if (!is_constant_all::value) { edge<2>(ops_partials).partials_ = cos_mu_minus_y - - modified_bessel_first_kind(-1, kappa_val) + - modified_bessel_first_kind(1, kappa_val) / modified_bessel_first_kind(0, kappa_val); } From 7024faeb2a6620488b9942c287624179b58dcc36 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sat, 20 Jan 2024 18:30:43 +0200 Subject: [PATCH 2/2] Update von mises test to use AD framework --- test/unit/math/mix/prob/von_mises_test.cpp | 122 ++++----------------- 1 file changed, 21 insertions(+), 101 deletions(-) diff --git a/test/unit/math/mix/prob/von_mises_test.cpp b/test/unit/math/mix/prob/von_mises_test.cpp index eed79db687f..6d959523201 100644 --- a/test/unit/math/mix/prob/von_mises_test.cpp +++ b/test/unit/math/mix/prob/von_mises_test.cpp @@ -1,103 +1,23 @@ #include -#include -#include -#include - -std::vector test_von_mises_lpdf(double y, double mu, double kappa) { - using stan::math::var; - using stan::math::von_mises_lpdf; - - var y_var = y; - var mu_var = mu; - var kappa_var = kappa; - - std::vector x; - x.push_back(y_var); - x.push_back(mu_var); - x.push_back(kappa_var); - - var logp = von_mises_lpdf(y_var, mu_var, kappa_var); - std::vector grad; - logp.grad(x, grad); - return grad; -} - -TEST(ProbAgradDistributionsVonMises, derivatives) { - using stan::math::fvar; - using stan::math::von_mises_lpdf; - - std::vector grad = test_von_mises_lpdf(0, 1, 0); - - fvar lp = von_mises_lpdf(0, 1, fvar(0, 1)); - EXPECT_FLOAT_EQ(grad[2], lp.d()); - - fvar kappa1(0, 1); - EXPECT_FLOAT_EQ(von_mises_lpdf(0, 1, kappa1).val_, -1.8378770664093453390819); - EXPECT_FLOAT_EQ(von_mises_lpdf(0, 1, kappa1).d_, 0.54030230586813976501); - - fvar> kappa2(0); - EXPECT_NO_THROW(von_mises_lpdf(0, 1, kappa2)); -} - -TEST(ProbAgradDistributionsVonMises, FvarVar_1stDeriv) { - using stan::math::fvar; - using stan::math::var; - using stan::math::von_mises_lpdf; - - fvar y_(0, 1); - double mu(1); - double kappa(0); - - fvar logp = von_mises_lpdf(y_, mu, kappa); - - std::vector y{y_.val_}; - std::vector g; - logp.val_.grad(y, g); - EXPECT_FLOAT_EQ(0, g[0]); -} - -TEST(ProbAgradDistributionsVonMises, FvarVar_2ndDeriv1) { - using stan::math::fvar; - using stan::math::var; - using stan::math::von_mises_lpdf; - - double y_(0); - fvar mu(1, 1); - double kappa(0); - fvar logp = von_mises_lpdf(y_, mu, kappa); - - std::vector y{mu.val_}; - std::vector g; - logp.d_.grad(y, g); - EXPECT_FLOAT_EQ(0, g[0]); -} - -TEST(ProbAgradDistributionsVonMises, FvarVar_2ndDeriv2) { - using stan::math::fvar; - using stan::math::var; - using stan::math::von_mises_lpdf; - - double y_(0); - double mu(1); - fvar kappa(0, 1); - fvar logp = von_mises_lpdf(y_, mu, kappa); - - std::vector y{kappa.val_}; - std::vector g; - logp.d_.grad(y, g); - // Fails: g[0] is -nan - // EXPECT_FLOAT_EQ(0, g[0]); -} - -// This test once failed sanitizer checks -- nothing explicitly tested in the -// code itself -TEST(ProbAgradDistributionsVonMises, sanitizer_error_fixed) { - using stan::math::var; - double y = boost::math::constants::third_pi(); - double mu = boost::math::constants::sixth_pi(); - std::vector kappa = {0.5}; - - auto lp = stan::math::von_mises_lpdf(y, mu, kappa); - - lp.grad(); +#include + +TEST(mathMixScalFun, von_mises_lpdf) { + auto f = [](const auto& y, const auto& mu, const auto& kappa) { + return stan::math::von_mises_lpdf(y, mu, kappa); + }; + + Eigen::VectorXd y(2); + y << 1.0, 2.0; + Eigen::VectorXd mu(2); + mu << 1.0, 0.5; + Eigen::VectorXd kappa(2); + kappa << 1.0, 0.5; + + stan::test::expect_ad(f, y[0], mu[0], kappa[0]); + stan::test::expect_ad(f, y[0], mu, kappa); + stan::test::expect_ad(f, y, mu[0], kappa); + stan::test::expect_ad(f, y, mu, kappa[0]); + stan::test::expect_ad(f, y[0], mu[0], kappa); + stan::test::expect_ad(f, y, mu[0], kappa[0]); + stan::test::expect_ad(f, y[0], mu, kappa[0]); }