Hey all,
I need to calculate \log(\Phi(y)) in a numerically stable way (where \Phi is the CDF of a standard normal distribution), where y may be a very large or small number. What is the best way to do this in stan math? Currently, I am using
log(Phi_approx(y))
which seems to work well, but I am curious if there is a better way.
Thanks,
Eric.
Hello Eric,
The Stan Functions reference indeed recommends log(Phi_approx(y))
realstd_normal_lcdf(reals y)
The log of the cumulative standard normal distribution of y; std_normal_lcdf will underflow to −∞ for y below -37.5 and overflow to 0 for y above 8.25;log(Phi_approx(...))is more robust in the tails.
However, based on this Github issue, it seems that std_normal_lcdf has been improved at some point, actually making it more numerically stable than both log(Phi(y)) and log(Phi_approx(y)). So it seems that the functions reference is outdated.
So personally, I tend to use std_normal_lcdf(y).
Thanks Frank!
It appears that there’s a subtle difference between using log(Phi_approx(y)) and std_normal_lcdf(y) in C++.
My original code using Phi_approx (which compiles and works) looks like this:
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> first_term_lcdf = stan::math::log(
stan::math::Phi_approx(
-first_term_cdf_arg
)
);
where first_term_cdf_arg is of type Eigen::MatrixXd
However, when I change to using std_normal_lcdf(y) like so
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> first_term_lcdf =
stan::math::std_normal_lcdf(
-first_term_cdf_arg
);
I get the following compilation error:
softmax_approx.cpp:40:57: error: no viable conversion from 'return_type_t<CwiseUnaryOp<scalar_opposite_op<double>, const Matrix<double, -1, -1, 0, -1, -1>>>' (aka 'double') to 'Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>'
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> first_term_lcdf = stan::math::std_normal_lcdf(
^ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/Matrix.h:273:5: note: candidate constructor not viable: no known conversion from 'return_type_t<CwiseUnaryOp<scalar_opposite_op<double>, const Matrix<double, -1, -1, 0, -1, -1>>>' (aka 'double') to 'Matrix<double, -1, -1> &&' for 1st argument
Matrix(Matrix&& other) EIGEN_NOEXCEPT_IF(std::is_nothrow_move_constructible<Scalar>::value)
^
/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/Matrix.h:414:25: note: candidate constructor not viable: no known conversion from 'return_type_t<CwiseUnaryOp<scalar_opposite_op<double>, const Matrix<double, -1, -1, 0, -1, -1>>>' (aka 'double') to 'const Matrix<double, -1, -1> &' for 1st argument
EIGEN_STRONG_INLINE Matrix(const Matrix& other) : Base(other)
^
/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/Matrix.h:422:25: note: candidate template ignored: could not match 'EigenBase<OtherDerived>' against 'return_type_t<CwiseUnaryOp<scalar_opposite_op<double>, const Matrix<double, -1, -1, 0, -1, -1>>>' (aka 'double')
EIGEN_STRONG_INLINE Matrix(const EigenBase<OtherDerived> &other)
^
/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/Matrix.h:267:14: note: explicit constructor is not a candidate
explicit Matrix(internal::constructor_without_unaligned_array_assert)
^
/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/Matrix.h:321:34: note: explicit constructor is not a candidate
explicit EIGEN_STRONG_INLINE Matrix(const std::initializer_list<std::initializer_list<Scalar>>& list) : Base(list) {}
^
/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/Matrix.h:329:14: note: explicit constructor is not a candidate
explicit Matrix(const T& x)
^
/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/Matrix.h:435:14: note: explicit constructor is not a candidate
explicit Matrix(const RotationBase<OtherDerived,ColsAtCompileTime>& r);
^
In file included from softmax_approx.cpp:5:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/mix.hpp:24:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim.hpp:16:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/prob.hpp:256:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/prob/ordered_probit_log.hpp:5:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/prob/ordered_probit_lpmf.hpp:13:
/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/prob/std_normal_lcdf.hpp:45:28: error: implicit instantiation of undefined template 'stan::scalar_seq_view<Eigen::Matrix<double, -1, -1>>'
scalar_seq_view<T_y_ref> y_vec(y_ref);
^
softmax_approx.cpp:40:87: note: in instantiation of function template specialization 'stan::math::std_normal_lcdf<Eigen::CwiseUnaryOp<Eigen::internal::scalar_opposite_op<double>, const Eigen::Matrix<double, -1, -1>>, nullptr>' requested here
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> first_term_lcdf = stan::math::std_normal_lcdf(
^
/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/scalar_seq_view.hpp:18:7: note: template is declared here
class scalar_seq_view;
^
Any ideas on how to use this function correctly? Perhaps it needs a revision to have the same functionality as Phi_approx?
Maybe a good question for @Bob_Carpenter ?
Thanks!
Eric.
I know very little about Stan’s C++ code so unfortunately I can’t help you with the error messages you’re seeing. Hopefully one of the Stan developers can pitch in!
log(Phi_approx) will work elementwise over its input, but std_normal_lcdf will return the sum of all the terms, rather than a container
I see. @WardBrian is there a way to get the stability of std_normal_lcdf as an elementwise operation?
A good old fashioned double-for loop would do it, though you might lose out on some speed comparatively. If you’re not trying to get derivatives, it probably won’t be too bad, though
@WardBrian unfortunately I will need to use this to take derivatives (though above I just gave an example with a matrix of doubles).
Would there be any way to turn this into an elementwise operation?
If so, I would be happy to submit a PR. However, I’m not the most familiar with the internals of stan math so I might need a reference (e.g. another stan function source) for how to do this.
Thanks!
That would not be a small task, unfortunately. It’s related to the longstanding request (2015!) to also have elementwise versions of the PDFs (which also currently return a sum), as they internally use most of the same machinery. See [WIP] Allow distributions to return vectors by SteveBronder · Pull Request #2751 · stan-dev/math · GitHub and linked issues for some more detail