Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Feature/matrix input to mvn [WIP] #2540

Closed
wants to merge 2 commits into from
Closed

Feature/matrix input to mvn [WIP] #2540

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
aa929c9
initial implementation
adamhaber Jul 18, 2021
03de55a
Merge branch 'stan-dev:develop' into feature/matrix-input-to-mvn
adamhaber Jul 18, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Next Next commit
initial implementation
  • Loading branch information
@adamhaber
adamhaber committed Jul 18, 2021
commit aa929c989303e702b6bee7fbbb536d43d729cdb4
21 changes: 20 additions & 1 deletion stan/math/prim/prob/multi_normal_lpdf.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,8 @@
namespace stan {
namespace math {

template <bool propto, typename T_y, typename T_loc, typename T_covar>
template <bool propto, typename T_y, typename T_loc, typename T_covar,
require_all_vector_t<T_y, T_loc>* = nullptr>
return_type_t<T_y, T_loc, T_covar> multi_normal_lpdf(const T_y& y,
const T_loc& mu,
const T_covar& Sigma) {
Expand Down Expand Up @@ -102,6 +103,24 @@ return_type_t<T_y, T_loc, T_covar> multi_normal_lpdf(const T_y& y,
return lp;
}

template <bool propto, typename T_y, typename T_loc, typename T_covar,
require_all_not_vector_t<T_y, T_loc>* = nullptr>
return_type_t<T_y, T_loc, T_covar> multi_normal_lpdf(const T_y& y,
const T_loc& mu,
const T_covar& Sigma) {
using lp_type = return_type_t<T_y, T_loc, T_covar>;

const size_t N = y.rows();
lp_type lp(0.0);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you take away the using lp_type... statement on 111 and declare lp as that type? It's easier to read when that's together.


for (size_t n = 0; n < N; ++n) {
auto current_y = y.row(n);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this code safe? What happens if y and mu have different lengths? Won't this just read off the end?

I'd like to see this fixed before it gets into the code base. Perhaps use a check_matching_sizes or something similar?

auto current_mu = mu.row(n);
lp += multi_normal_lpdf<false>(current_y, current_mu, Sigma);
}
return lp;
}

template <typename T_y, typename T_loc, typename T_covar>
inline return_type_t<T_y, T_loc, T_covar> multi_normal_lpdf(
const T_y& y, const T_loc& mu, const T_covar& Sigma) {
Expand Down
19 changes: 18 additions & 1 deletion stan/math/prim/prob/multi_normal_rng.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ namespace math {
* std::invalid_argument if the length of (each) mu is not equal to the number
* of rows and columns in S
*/
template <typename T_loc, class RNG>
template <typename T_loc, class RNG, require_vector_t<T_loc>* = nullptr>
inline typename StdVectorBuilder<true, Eigen::VectorXd, T_loc>::type
multi_normal_rng(const T_loc& mu,
const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>& S,
Expand Down Expand Up @@ -80,6 +80,23 @@ multi_normal_rng(const T_loc& mu,
return output.data();
}

template <typename T_loc, class RNG, require_not_vector_t<T_loc>* = nullptr>
inline typename Eigen::MatrixXd multi_normal_rng(
const T_loc& mu,
const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>& S, RNG& rng) {
using Eigen::MatrixXd;

const size_t N = mu.rows();
const size_t M = mu.cols();

MatrixXd output(N, M);
for (int n = 0; n < N; ++n) {
auto current_mu = mu.row(n);
output.row(n) = multi_normal_rng(current_mu, S, rng);
}
return output;
}

} // namespace math
} // namespace stan
#endif
16 changes: 16 additions & 0 deletions test/unit/math/prim/prob/multi_normal_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,22 @@ TEST(ProbDistributionsMultiNormal, NotVectorized) {
EXPECT_NO_THROW(stan::math::multi_normal_rng(mu, Sigma, rng));
}

TEST(ProbDistributionsMultiNormal, MatrixNotVectorized) {
using Eigen::Dynamic;
using Eigen::Matrix;
using std::vector;

boost::random::mt19937 rng;
Matrix<double, Dynamic, Dynamic> y(2, 3);
y << 2.0, -2.0, 11.0, 2.0, -2.0, 11.0;
Matrix<double, Dynamic, Dynamic> mu(2, 3);
mu << 1.0, -1.0, 3.0, 1.0, -1.0, 3.0;
Matrix<double, Dynamic, Dynamic> Sigma(3, 3);
Sigma << 9.0, -3.0, 0.0, -3.0, 4.0, 0.0, 0.0, 0.0, 5.0;
EXPECT_FLOAT_EQ(-23.47816, stan::math::multi_normal_lpdf(y, mu, Sigma));
EXPECT_NO_THROW(stan::math::multi_normal_rng(mu, Sigma, rng));
}

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please add a test where y and mu have different number of rows. This will trigger a failure and we should know what that is through a unit test.

TEST(ProbDistributionsMultiNormal, Vectorized) {
using Eigen::Dynamic;
using Eigen::Matrix;
Expand Down