add yule_simon_rng - #3285
Conversation
Sorry, something went wrong.
Jenkins Console Log Machine information No LSB modules are available. Distributor ID: Ubuntu Description: Ubuntu 20.04.3 LTS Release: 20.04 Codename: focal CPU: G++: Clang: |
Sorry, something went wrong.
|
Rtools & ARM64 CI / rev tests (windows-11-arm) (pull_request) [----------] 11 tests from StanOde/cos_arg_test/4, where TypeParam = std::tuple<ode_adjoint_functor, ode_adjoint_functor>
[ RUN ] StanOde/cos_arg_test/4.y0_error
[ OK ] StanOde/cos_arg_test/4.y0_error (0 ms)
[ RUN ] StanOde/cos_arg_test/4.t0_error
[ OK ] StanOde/cos_arg_test/4.t0_error (0 ms)
[ RUN ] StanOde/cos_arg_test/4.ts_error
[ OK ] StanOde/cos_arg_test/4.ts_error (1 ms)
[ RUN ] StanOde/cos_arg_test/4.one_arg_error
[ OK ] StanOde/cos_arg_test/4.one_arg_error (1 ms)
[ RUN ] StanOde/cos_arg_test/4.two_arg_error
[ OK ] StanOde/cos_arg_test/4.two_arg_error (2 ms)
[ RUN ] StanOde/cos_arg_test/4.rhs_wrong_size_error
[ OK ] StanOde/cos_arg_test/4.rhs_wrong_size_error (0 ms)
[ RUN ] StanOde/cos_arg_test/4.error_name
[ OK ] StanOde/cos_arg_test/4.error_name (0 ms)
[ RUN ] StanOde/cos_arg_test/4.tol_error
[ OK ] StanOde/cos_arg_test/4.tol_error (0 ms)
[ RUN ] StanOde/cos_arg_test/4.value
[ OK ] StanOde/cos_arg_test/4.value (0 ms)
[ RUN ] StanOde/cos_arg_test/4.grad
[ OK ] StanOde/cos_arg_test/4.grad (6 ms)
[ RUN ] StanOde/cos_arg_test/4.tol_grad
./test/unit/math/rev/functor/test_fixture_ode_cos_scalar.hpp:1317: Failure
Expected equality of these values:
a.adj()
Which is: -0.43164793
-0.5 * (tsd[0] * tsd[0] - t0d * t0d) * y0d(0) * exp(-0.5 * ad * (tsd[0] * tsd[0] - t0d * t0d))
Which is: -0.43164769
[ FAILED ] StanOde/cos_arg_test/4.tol_grad, where TypeParam = std::tuple<ode_adjoint_functor, ode_adjoint_functor> (5 ms)
[----------] 11 tests from StanOde/cos_arg_test/4 (15 ms total)
[----------] Global test environment tear-down
[==========] 55 tests from 5 test suites ran. (37 ms total)
[ PASSED ] 54 tests.
[ FAILED ] 1 test, listed below:
[ FAILED ] StanOde/cos_arg_test/4.tol_grad, where TypeParam = std::tuple<ode_adjoint_functor, ode_adjoint_functor>
1 FAILED TEST
test\unit\math\rev\functor\cos_ode_typed_test.exe --gtest_output="xml:test/unit/math/rev/functor/cos_ode_typed_test.xml" failed
Seems a test on ode fails, but this might be unrelated to this PR? My mac can pass the rev tests. |
Sorry, something went wrong.
|
@lingium it passed on re-run, so I think it was just random noise |
Sorry, something went wrong.
SteveBronder
left a comment
There was a problem hiding this comment.
Two minor comments but imo looks good!
Sorry, something went wrong.
| scalar_seq_view<T_w> w_vec(w); | ||
| size_t size_w = stan::math::size(w); | ||
|
|
||
| VectorBuilder<true, int, T_alpha> output(size_w); | ||
| for (size_t n = 0; n < size_w; ++n) { | ||
| double p = stan::math::exp(-w_vec.val(n)); | ||
| double odds_ratio_p | ||
| = stan::math::exp(stan::math::log(p) - stan::math::log1m(p)); | ||
| output[n] = neg_binomial_rng(1.0, odds_ratio_p, rng) + 1; | ||
| } | ||
|
|
||
| return output.data(); |
There was a problem hiding this comment.
Since you use neg_binomial_rng I think we can use the vectorized version of everything here
| scalar_seq_view<T_w> w_vec(w); | |
| size_t size_w = stan::math::size(w); | |
| VectorBuilder<true, int, T_alpha> output(size_w); | |
| for (size_t n = 0; n < size_w; ++n) { | |
| double p = stan::math::exp(-w_vec.val(n)); | |
| double odds_ratio_p | |
| = stan::math::exp(stan::math::log(p) - stan::math::log1m(p)); | |
| output[n] = neg_binomial_rng(1.0, odds_ratio_p, rng) + 1; | |
| } | |
| return output.data(); | |
| auto p = stan::math::exp(-w); | |
| auto odds_ratio_p = stan::math::exp(stan::math::log(p) - stan::math::log1m(p)); | |
| return neg_binomial_rng(1.0, odds_ratio_p, rng) + 1; |
Sorry, something went wrong.
There was a problem hiding this comment.
Seems neg_binomial_rng(...) + 1; cannot pass compile, let's keep the loop version?
Sorry, something went wrong.
There was a problem hiding this comment.
How about like this?
/** \ingroup prob_dists
* Return a yule-simon random variate with the given shape parameter,
* using the given random number generator.
*
* alpha can be a scalar or a one-dimensional container.
*
* @tparam T_alpha A scalar, `std::vector` or Eigen vector
* @tparam RNG type of random number generator
*
* @param alpha (Sequence of) shape parameter(s)
* @param rng random number generator
* @return (Sequence of) yule-simon random variate(s)
* @throw std::domain_error if alpha is nonpositive
*/
template <typename T_alpha, typename RNG>
inline auto yule_simon_rng(T_alpha&& alpha, RNG& rng) {
static constexpr const char *function = "yule_simon_rng";
decltype(auto) alpha_ref = to_ref(std::forward<T_alpha>(alpha));
check_positive_finite(function, "Shape parameter", alpha_ref);
auto w = exponential_rng(std::forward<decltype(alpha_ref)>(alpha_ref), rng);
auto w_arr = as_array_or_scalar(w);
const auto p = stan::math::exp(-w_arr);
const auto odds_ratio_p
= stan::math::exp(stan::math::log(p) - stan::math::log1m(p));
if constexpr (is_stan_scalar_v<T_alpha>) {
return neg_binomial_rng(1.0, odds_ratio_p, rng) + 1;
} else {
return to_array_1d(as_array_or_scalar(neg_binomial_rng(1.0, std::move(odds_ratio_p), rng)) + 1);
}
}
Sorry, something went wrong.
There was a problem hiding this comment.
@lingium just pinging to see if you saw this
Sorry, something went wrong.
Jenkins Console Log Machine information No LSB modules are available. Distributor ID: Ubuntu Description: Ubuntu 20.04.3 LTS Release: 20.04 Codename: focal CPU: G++: Clang: |
Sorry, something went wrong.
Jenkins Console Log Machine information No LSB modules are available. Distributor ID: Ubuntu Description: Ubuntu 20.04.3 LTS Release: 20.04 Codename: focal CPU: G++: Clang: |
Sorry, something went wrong.
SteveBronder
left a comment
There was a problem hiding this comment.
Looks good!
Sorry, something went wrong.
|
@lingium I'd be happy to do the stanc3 pr for all these new functions if you don't mind writing the signatures you'd expect them to have out somewhere |
Sorry, something went wrong.
|
Hi @WardBrian , it'd be the best if you could help with stanc3, I'll do with the doc. The signatures are yule_simon_lpmfreal yule_simon_lpmf(int n, real alpha);
real yule_simon_lpmf(int n, vector alpha);
real yule_simon_lpmf(array[] int n, real alpha);
real yule_simon_lpmf(array[] int n, vector alpha);yule_simon_lcdfreal yule_simon_lcdf(int n, real alpha);
real yule_simon_lcdf(int n, vector alpha);
real yule_simon_lcdf(array[] int n, real alpha);
real yule_simon_lcdf(array[] int n, vector alpha);yule_simon_lccdfreal yule_simon_lccdf(int n, real alpha);
real yule_simon_lccdf(int n, vector alpha);
real yule_simon_lccdf(array[] int n, real alpha);
real yule_simon_lccdf(array[] int n, vector alpha);yule_simon_cdfreal yule_simon_cdf(int n, real alpha);
real yule_simon_cdf(int n, vector alpha);
real yule_simon_cdf(array[] int n, real alpha);
real yule_simon_cdf(array[] int n, vector alpha);yule_simon_rngint yule_simon_rng(real alpha);
array[] int yule_simon_rng(vector alpha); |
Sorry, something went wrong.
Summary
With this PR the rng of Yule-Simon distribution are added.
See issue #3239
For $$Y\sim \text{Yule-Simon}(\alpha)$$:
$$ \Pr(Y=y\mid \alpha)= \alpha \mathrm{B}(y,\alpha+1) \quad y\in{1,2,\dots} $$
$$\mathrm{B}$$ is the Beta function.
It has a Beta–Geometric mixture. To sample form $$Y$$, we do:
$$ \begin{align} Y\mid p &\sim \text{Geometric}(p)\\ p=e^{-W} &\sim \text{Beta}(\alpha,1)\\ W &\sim\text{Exp}(\alpha) \end{align} $$
Stan does not have Geometric yet so we use neg_binomial_rng(1.0, odds_ratio_p).
Tests
Test is written follow the guide.
Side Effects
No.
Release notes
yule_simon_rng is available if merged.
Checklist
Copyright holder: Zhi Ling
The copyright holder is typically you or your assignee, such as a university or company. By submitting this pull request, the copyright holder is agreeing to the 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/)
the basic tests are passing
the code is written in idiomatic C++ and changes are documented in the doxygen
the new changes are tested