Skip to content

Commit d130b18

Browse files
yashiknoyashikno
committed
Merge commit 'b95dd5a46e8bdfdf984d11413dc86bc98306a3d8' into HEAD
2 parents 576711f + b95dd5a commit d130b18

8 files changed

Lines changed: 355 additions & 0 deletions

File tree

‎stan/math/prim/prob.hpp‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -202,6 +202,9 @@
202202
#include <stan/math/prim/prob/multi_student_t_log.hpp>
203203
#include <stan/math/prim/prob/multi_student_t_lpdf.hpp>
204204
#include <stan/math/prim/prob/multi_student_t_rng.hpp>
205+
#include <stan/math/prim/prob/multinomial_logit_log.hpp>
206+
#include <stan/math/prim/prob/multinomial_logit_lpmf.hpp>
207+
#include <stan/math/prim/prob/multinomial_logit_rng.hpp>
205208
#include <stan/math/prim/prob/multinomial_log.hpp>
206209
#include <stan/math/prim/prob/multinomial_lpmf.hpp>
207210
#include <stan/math/prim/prob/multinomial_rng.hpp>

‎stan/math/prim/prob/multinomial_logit_log.hpp‎

Lines changed: 34 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,34 @@
1+
#ifndef STAN_MATH_PRIM_PROB_MULTINOMIAL_LOGIT_LOG_HPP
2+
#define STAN_MATH_PRIM_PROB_MULTINOMIAL_LOGIT_LOG_HPP
3+
4+
#include <stan/math/prim/meta.hpp>
5+
#include <stan/math/prim/fun/Eigen.hpp>
6+
#include <stan/math/prim/prob/multinomial_logit_lpmf.hpp>
7+
#include <vector>
8+
9+
namespace stan {
10+
namespace math {
11+
12+
/** \ingroup multivar_dists
13+
* @deprecated use <code>multinomial_logit_lpmf</code>
14+
*/
15+
template <bool propto, typename T_beta, typename T_prob = scalar_type_t<T_beta>,
16+
require_eigen_col_vector_t<T_beta>* = nullptr>
17+
return_type_t<T_prob> multinomial_logit_log(const std::vector<int>& ns,
18+
const T_beta& beta) {
19+
return multinomial_logit_lpmf<propto, T_beta>(ns, beta);
20+
}
21+
22+
/** \ingroup multivar_dists
23+
* @deprecated use <code>multinomial_logit_lpmf</code>
24+
*/
25+
template <typename T_beta, typename T_prob = scalar_type_t<T_beta>,
26+
require_eigen_col_vector_t<T_beta>* = nullptr>
27+
return_type_t<T_prob> multinomial_logit_log(const std::vector<int>& ns,
28+
const T_beta& beta) {
29+
return multinomial_logit_lpmf<false>(ns, beta);
30+
}
31+
32+
} // namespace math
33+
} // namespace stan
34+
#endif

‎stan/math/prim/prob/multinomial_logit_lpmf.hpp‎

Lines changed: 60 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,60 @@
1+
#ifndef STAN_MATH_PRIM_PROB_MULTINOMIAL_LOGIT_LPMF_HPP
2+
#define STAN_MATH_PRIM_PROB_MULTINOMIAL_LOGIT_LPMF_HPP
3+
4+
#include <stan/math/prim/meta.hpp>
5+
#include <stan/math/prim/err.hpp>
6+
#include <stan/math/prim/fun/lgamma.hpp>
7+
#include <stan/math/prim/fun/log_sum_exp.hpp>
8+
#include <vector>
9+
10+
namespace stan {
11+
namespace math {
12+
13+
/** \ingroup multivar_dists
14+
* Multinomial log PMF in log parametrization.
15+
* Multinomial(ns| softmax(beta))
16+
*
17+
* @param ns Array of outcome counts
18+
* @param beta Vector of unnormalized log probabilities
19+
* @return log probability
20+
*/
21+
template <bool propto, typename T_beta, typename T_prob = scalar_type_t<T_beta>,
22+
require_eigen_col_vector_t<T_beta>* = nullptr>
23+
return_type_t<T_prob> multinomial_logit_lpmf(const std::vector<int>& ns,
24+
const T_beta& beta) {
25+
static const char* function = "multinomial_logit_lpmf";
26+
check_nonnegative(function, "Number of trials variable", ns);
27+
check_finite(function, "log-probabilities parameter", beta);
28+
check_size_match(function, "Size of number of trials variable", ns.size(),
29+
"rows of log-probabilities parameter", beta.rows());
30+
31+
return_type_t<T_prob> lp(0.0);
32+
33+
decltype(auto) ns_map = as_array_or_scalar(ns);
34+
35+
if (include_summand<propto>::value) {
36+
lp += lgamma(1 + ns_map.sum()) - lgamma(1 + ns_map).sum();
37+
}
38+
39+
if (include_summand<propto, T_prob>::value) {
40+
decltype(auto) beta_ref = to_ref(beta);
41+
T_prob alpha = log_sum_exp(beta_ref);
42+
for (unsigned int i = 0; i < ns.size(); ++i) {
43+
if (ns[i] != 0)
44+
lp += ns[i] * (beta_ref[i] - alpha);
45+
}
46+
}
47+
48+
return lp;
49+
}
50+
51+
template <typename T_beta, typename T_prob = scalar_type_t<T_beta>,
52+
require_eigen_col_vector_t<T_beta>* = nullptr>
53+
return_type_t<T_prob> multinomial_logit_lpmf(const std::vector<int>& ns,
54+
const T_beta& beta) {
55+
return multinomial_logit_lpmf<false>(ns, beta);
56+
}
57+
58+
} // namespace math
59+
} // namespace stan
60+
#endif

‎stan/math/prim/prob/multinomial_logit_rng.hpp‎

Lines changed: 52 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,52 @@
1+
#ifndef STAN_MATH_PRIM_PROB_MULTINOMIAL_LOGIT_RNG_HPP
2+
#define STAN_MATH_PRIM_PROB_MULTINOMIAL_LOGIT_RNG_HPP
3+
4+
#include <stan/math/prim/meta.hpp>
5+
#include <stan/math/prim/err.hpp>
6+
#include <stan/math/prim/fun/softmax.hpp>
7+
#include <stan/math/prim/prob/binomial_rng.hpp>
8+
#include <vector>
9+
10+
namespace stan {
11+
namespace math {
12+
13+
/** \ingroup multivar_dists
14+
* Return a draw from a Multinomial distribution given a
15+
* a vector of unnormalized log probabilities and a pseudo-random
16+
* number generator.
17+
*
18+
* @tparam RNG Type of pseudo-random number generator.
19+
* @param beta Vector of unnormalized log probabilities.
20+
* @param N Total count
21+
* @param rng Pseudo-random number generator.
22+
* @return Multinomial random variate
23+
*/
24+
template <class RNG, typename T_beta,
25+
require_eigen_col_vector_t<T_beta>* = nullptr>
26+
inline std::vector<int> multinomial_logit_rng(const T_beta& beta, int N,
27+
RNG& rng) {
28+
static const char* function = "multinomial_logit_rng";
29+
check_finite(function, "Log-probabilities parameter", beta);
30+
check_positive(function, "number of trials variables", N);
31+
32+
plain_type_t<T_beta> theta = softmax(beta);
33+
std::vector<int> result(theta.size(), 0);
34+
double mass_left = 1.0;
35+
int n_left = N;
36+
37+
for (int k = 0; n_left > 0 && k < theta.size(); ++k) {
38+
double p = theta[k] / mass_left;
39+
if (p > 1.0) {
40+
p = 1.0;
41+
}
42+
result[k] = binomial_rng(n_left, p, rng);
43+
n_left -= result[k];
44+
mass_left -= theta[k];
45+
}
46+
47+
return result;
48+
} // namespace math
49+
50+
} // namespace math
51+
} // namespace stan
52+
#endif

‎test/unit/math/mix/prob/multinomial_logit_test.cpp‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,13 @@
1+
#include <test/unit/math/test_ad.hpp>
2+
3+
TEST(mathMixScalFun, multinomialLogit) {
4+
std::vector<int> ns{0, 1, 2, 3};
5+
Eigen::VectorXd beta(4);
6+
beta << 0.1, 0.1, 0.5, 0.3;
7+
8+
auto f = [&ns](const auto& b) {
9+
return stan::math::multinomial_logit_lpmf(ns, b);
10+
};
11+
12+
stan::test::expect_ad(f, beta);
13+
}

‎test/unit/math/prim/prob/multinomial_logit_log_test.cpp‎

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
1+
#include <stan/math/prim.hpp>
2+
#include <gtest/gtest.h>
3+
#include <vector>
4+
5+
TEST(ProbMultinomialLogit, log_matches_lpmf) {
6+
using stan::math::multinomial_logit_log;
7+
using stan::math::multinomial_logit_lpmf;
8+
9+
std::vector<int> ns;
10+
ns.push_back(1);
11+
ns.push_back(2);
12+
ns.push_back(3);
13+
Eigen::Matrix<double, Eigen::Dynamic, 1> theta(3, 1);
14+
theta << log(0.2), log(0.3), log(0.5);
15+
EXPECT_FLOAT_EQ((multinomial_logit_lpmf(ns, theta)),
16+
(multinomial_logit_log(ns, theta)));
17+
EXPECT_FLOAT_EQ((multinomial_logit_lpmf<true>(ns, theta)),
18+
(multinomial_logit_log<true>(ns, theta)));
19+
EXPECT_FLOAT_EQ((multinomial_logit_lpmf<false>(ns, theta)),
20+
(multinomial_logit_log<false>(ns, theta)));
21+
}

‎test/unit/math/prim/prob/multinomial_logit_test.cpp‎

Lines changed: 126 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,126 @@
1+
#include <stan/math/prim.hpp>
2+
#include <gtest/gtest.h>
3+
#include <boost/random/mersenne_twister.hpp>
4+
#include <boost/math/distributions.hpp>
5+
#include <limits>
6+
#include <vector>
7+
8+
using Eigen::Dynamic;
9+
using Eigen::Matrix;
10+
11+
TEST(ProbDistributionsMultinomialLogit, RNGSize) {
12+
boost::random::mt19937 rng;
13+
Matrix<double, Dynamic, 1> beta(5);
14+
beta << log(0.3), log(0.1), log(0.2), log(0.2), log(0.2);
15+
std::vector<int> sample = stan::math::multinomial_logit_rng(beta, 10, rng);
16+
EXPECT_EQ(5U, sample.size());
17+
}
18+
19+
TEST(ProbDistributionsMultinomialLogit, MultinomialLogit) {
20+
std::vector<int> ns;
21+
ns.push_back(1);
22+
ns.push_back(2);
23+
ns.push_back(3);
24+
Matrix<double, Dynamic, 1> beta(3, 1);
25+
beta << log(0.2), log(0.3), log(0.5);
26+
EXPECT_FLOAT_EQ(-2.002481, stan::math::multinomial_logit_log(ns, beta));
27+
}
28+
TEST(ProbDistributionsMultinomialLogit, Propto) {
29+
std::vector<int> ns;
30+
ns.push_back(1);
31+
ns.push_back(2);
32+
ns.push_back(3);
33+
Matrix<double, Dynamic, 1> beta(3, 1);
34+
beta << log(0.2), log(0.3), log(0.5);
35+
EXPECT_FLOAT_EQ(0.0, stan::math::multinomial_logit_log<true>(ns, beta));
36+
}
37+
38+
using stan::math::multinomial_logit_log;
39+
40+
TEST(ProbDistributionsMultinomialLogit, error) {
41+
double nan = std::numeric_limits<double>::quiet_NaN();
42+
double inf = std::numeric_limits<double>::infinity();
43+
44+
std::vector<int> ns;
45+
ns.push_back(1);
46+
ns.push_back(2);
47+
ns.push_back(3);
48+
Matrix<double, Dynamic, 1> beta(3, 1);
49+
beta << log(0.2), log(0.3), log(0.5);
50+
51+
EXPECT_NO_THROW(multinomial_logit_log(ns, beta));
52+
53+
ns[1] = 0;
54+
EXPECT_NO_THROW(multinomial_logit_log(ns, beta));
55+
ns[1] = -1;
56+
EXPECT_THROW(multinomial_logit_log(ns, beta), std::domain_error);
57+
ns[1] = 1;
58+
59+
beta(0) = nan;
60+
EXPECT_THROW(multinomial_logit_log(ns, beta), std::domain_error);
61+
beta(0) = inf;
62+
EXPECT_THROW(multinomial_logit_log(ns, beta), std::domain_error);
63+
beta(0) = -inf;
64+
EXPECT_THROW(multinomial_logit_log(ns, beta), std::domain_error);
65+
66+
beta(0) = 0.2;
67+
beta(1) = 0.3;
68+
beta(2) = 0.5;
69+
70+
ns.resize(2);
71+
EXPECT_THROW(multinomial_logit_log(ns, beta), std::invalid_argument);
72+
}
73+
74+
TEST(ProbDistributionsMultinomialLogit, zeros) {
75+
double result;
76+
std::vector<int> ns;
77+
ns.push_back(0);
78+
ns.push_back(1);
79+
ns.push_back(2);
80+
Matrix<double, Dynamic, 1> beta(3, 1);
81+
beta << log(0.2), log(0.3), log(0.5);
82+
83+
result = multinomial_logit_log(ns, beta);
84+
EXPECT_FALSE(std::isnan(result));
85+
86+
std::vector<int> ns2;
87+
ns2.push_back(0);
88+
ns2.push_back(0);
89+
ns2.push_back(0);
90+
91+
double result2 = multinomial_logit_log(ns2, beta);
92+
EXPECT_FLOAT_EQ(0.0, result2);
93+
}
94+
95+
TEST(ProbDistributionsMultinomialLogit, chiSquareGoodnessFitTest) {
96+
boost::random::mt19937 rng;
97+
int M = 10;
98+
int trials = 1000;
99+
int N = M * trials;
100+
101+
int K = 3;
102+
Matrix<double, Dynamic, 1> beta(K);
103+
beta << -1, 1, -10;
104+
Eigen::VectorXd theta = stan::math::softmax(beta);
105+
boost::math::chi_squared mydist(K - 1);
106+
107+
double expect[K];
108+
for (int i = 0; i < K; ++i)
109+
expect[i] = N * theta(i);
110+
111+
int bin[K];
112+
for (int i = 0; i < K; ++i)
113+
bin[i] = 0;
114+
115+
for (int count = 0; count < M; ++count) {
116+
std::vector<int> a = stan::math::multinomial_logit_rng(beta, trials, rng);
117+
for (int i = 0; i < K; ++i)
118+
bin[i] += a[i];
119+
}
120+
121+
double chi = 0;
122+
for (int j = 0; j < K; j++)
123+
chi += ((bin[j] - expect[j]) * (bin[j] - expect[j])) / expect[j];
124+
125+
EXPECT_TRUE(chi < quantile(complement(mydist, 1e-6)));
126+
}

‎test/unit/math/rev/prob/multinomial_logit_test.cpp‎

Lines changed: 46 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,46 @@
1+
#include <stan/math/rev.hpp>
2+
#include <test/unit/math/rev/util.hpp>
3+
#include <test/unit/math/rev/prob/expect_eq_diffs.hpp>
4+
#include <gtest/gtest.h>
5+
#include <string>
6+
#include <vector>
7+
8+
using stan::math::multinomial_logit_lpmf;
9+
10+
template <typename T_prob>
11+
void expect_propto(std::vector<int>& ns1, T_prob beta1, std::vector<int>& ns2,
12+
T_prob beta2, std::string message) {
13+
expect_eq_diffs(multinomial_logit_lpmf<false>(ns1, beta1),
14+
multinomial_logit_lpmf<false>(ns2, beta2),
15+
multinomial_logit_lpmf<true>(ns1, beta1),
16+
multinomial_logit_lpmf<true>(ns2, beta2), message);
17+
}
18+
19+
using Eigen::Dynamic;
20+
using Eigen::Matrix;
21+
using stan::math::var;
22+
23+
TEST(AgradDistributionsMultinomialLogit, Propto) {
24+
std::vector<int> ns;
25+
ns.push_back(1);
26+
ns.push_back(2);
27+
ns.push_back(3);
28+
Matrix<var, Dynamic, 1> beta1(3, 1);
29+
beta1 << log(0.3), log(0.5), log(0.2);
30+
Matrix<var, Dynamic, 1> beta2(3, 1);
31+
beta2 << log(0.1), log(0.2), log(0.7);
32+
33+
expect_propto(ns, beta1, ns, beta2, "var: beta");
34+
}
35+
36+
TEST(AgradDistributionsMultinomialLogit, check_varis_on_stack) {
37+
std::vector<int> ns;
38+
ns.push_back(1);
39+
ns.push_back(2);
40+
ns.push_back(3);
41+
Matrix<var, Dynamic, 1> beta(3, 1);
42+
beta << log(0.3), log(0.5), log(0.2);
43+
44+
test::check_varis_on_stack(multinomial_logit_lpmf<false>(ns, beta));
45+
test::check_varis_on_stack(multinomial_logit_lpmf<true>(ns, beta));
46+
}

0 commit comments

Comments
 (0)