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

Vectorized unbounded continuous random number generators (Issue 45) #622

Merged
merged 17 commits into from
Nov 16, 2017
Merged

Vectorized unbounded continuous random number generators (Issue 45) #622

Show file tree
Hide file tree
Changes from 7 commits
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
72eb950
Changed vector_rng_test_helper tools to work with anywhere between 1-…
bbbales2 Sep 10, 2017
fa7187b
Changed int to size_t to avoid signed comparison warning (Issue 45)
bbbales2 Sep 18, 2017
270c96a
Merge remote-tracking branch 'stan/develop' into feature/issue-45-vec…
bbbales2 Sep 26, 2017
9a85e8c
Merge remote-tracking branch 'stan/develop' into feature/issue-45-vec…
bbbales2 Sep 29, 2017
ba6c624
Moved around & symbols and if statements. Avoided cast w/ templating
bbbales2 Sep 30, 2017
275889c
This commit addresses Bob's comments (from like September 30 or somet…
bbbales2 Oct 1, 2017
d2281d1
Added int and std::vector<int> tests (Issue 45)
bbbales2 Oct 5, 2017
427d316
Merge branch 'develop' into feature/issue-45-vectorized_rngs_unbounde…
mitzimorris Oct 24, 2017
8e83e4f
Cleaned up spacing, changed docs to make interface clearer (Issue 45)
bbbales2 Oct 26, 2017
ee338ff
Fixed cpplint issues in test (Issue 45)
bbbales2 Oct 26, 2017
25daa7c
Merge remote-tracking branch 'stan/develop' into feature/issue-45-vec…
bbbales2 Oct 26, 2017
9bf8a0e
Swapped std::log for log in double_exponential test to get things wor…
bbbales2 Oct 28, 2017
aa207eb
Made a couple fixes to get Jenkins tests to pass (errors not reproduc…
bbbales2 Oct 28, 2017
4517668
Fix cpplint (Issue 45)
bbbales2 Oct 28, 2017
eff2de9
Fixed references to std::log and std::abs in double_exponential_rng. …
bbbales2 Oct 29, 2017
ddb60d9
Merge remote-tracking branch 'stan/develop' into feature/issue-45-vec…
bbbales2 Nov 6, 2017
47cdfc9
Moved around variable declarations below argument checks (Issue 45)
bbbales2 Nov 14, 2017
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
53 changes: 34 additions & 19 deletions stan/math/prim/scal/prob/cauchy_rng.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,51 +1,66 @@
#ifndef STAN_MATH_PRIM_SCAL_PROB_CAUCHY_RNG_HPP
#define STAN_MATH_PRIM_SCAL_PROB_CAUCHY_RNG_HPP

#include <boost/random/cauchy_distribution.hpp>
#include <boost/random/variate_generator.hpp>
#include <stan/math/prim/scal/err/check_consistent_sizes.hpp>
#include <stan/math/prim/scal/err/check_finite.hpp>
#include <stan/math/prim/scal/err/check_not_nan.hpp>
#include <stan/math/prim/scal/err/check_positive_finite.hpp>
#include <stan/math/prim/scal/fun/constants.hpp>
#include <stan/math/prim/scal/fun/log1p.hpp>
#include <stan/math/prim/scal/fun/square.hpp>
#include <stan/math/prim/scal/fun/value_of.hpp>
#include <stan/math/prim/scal/meta/include_summand.hpp>
#include <stan/math/prim/scal/meta/max_size.hpp>
#include <stan/math/prim/scal/meta/scalar_seq_view.hpp>
#include <stan/math/prim/scal/meta/VectorBuilder.hpp>
#include <boost/random/cauchy_distribution.hpp>
#include <boost/random/variate_generator.hpp>
#include <string>

namespace stan {
namespace math {

/**
* Return a pseudorandom Cauchy variate for the given location and scale
* using the specified random number generator.
*
* mu and sigma can each be either scalars or vectors. All vector inputs
* must be the same length.
*
* @tparam T_loc Type of location parameter
* @tparam T_scale Type of scale parameter
* @tparam RNG type of random number generator
* @param mu location parameter
* @param sigma positive scale parameter
* @param mu (Sequence of) location parameter(s)
* @param sigma (Sequence of) scale parameter(s)
* @param rng random number generator
* @return Cauchy random variate
* @throw std::domain_error if mu is infinite or sigma is nonpositive
* @throw std::invalid_argument if vector arguments are not the same length
*/
template <class RNG>
inline double
cauchy_rng(double mu,
double sigma,
template <typename T_loc, typename T_scale, class RNG>
inline typename VectorBuilder<true, double, T_loc, T_scale>::type
cauchy_rng(const T_loc& mu,
const T_scale& sigma,
RNG& rng) {
using boost::variate_generator;
using boost::random::cauchy_distribution;

static const std::string function = "cauchy_rng";

scalar_seq_view<T_loc> mu_vec(mu);
scalar_seq_view<T_scale> sigma_vec(sigma);

check_finite(function, "Location parameter", mu);
Copy link
Contributor

Choose a reason for hiding this comment

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

Checks should go first before any work. In general, try to put exceptional/odd/edge-case branches first so that when you're by them, you can scan the main branch of the program together. Would you mind fixing this on all of these? Sorry for not catchign that earlier.

check_positive_finite(function, "Scale parameter", sigma);
check_consistent_sizes(function,
"Location parameter", mu,
"Scale Parameter", sigma);

variate_generator<RNG&, cauchy_distribution<> >
cauchy_rng(rng, cauchy_distribution<>(mu, sigma));
return cauchy_rng();
}
size_t N = max_size(mu, sigma);

VectorBuilder<true, double, T_loc, T_scale> output(N);

for (size_t n = 0; n < N; n++) {
variate_generator<RNG&, cauchy_distribution<> >
cauchy_rng(rng, cauchy_distribution<>(mu_vec[n], sigma_vec[n]));
output[n] = cauchy_rng();
}

return output.data();
}
}
}
#endif
61 changes: 36 additions & 25 deletions stan/math/prim/scal/prob/double_exponential_rng.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,58 +1,69 @@
#ifndef STAN_MATH_PRIM_SCAL_PROB_DOUBLE_EXPONENTIAL_RNG_HPP
#define STAN_MATH_PRIM_SCAL_PROB_DOUBLE_EXPONENTIAL_RNG_HPP

#include <boost/random/uniform_01.hpp>
#include <boost/random/variate_generator.hpp>
#include <stan/math/prim/scal/err/check_consistent_sizes.hpp>
#include <stan/math/prim/scal/err/check_finite.hpp>
#include <stan/math/prim/scal/err/check_not_nan.hpp>
#include <stan/math/prim/scal/err/check_positive_finite.hpp>
#include <stan/math/prim/scal/fun/value_of.hpp>
#include <stan/math/prim/scal/fun/log1m.hpp>

#include <stan/math/prim/scal/fun/constants.hpp>
#include <stan/math/prim/scal/meta/include_summand.hpp>
#include <stan/math/prim/scal/fun/sign.hpp>
#include <stan/math/prim/scal/meta/max_size.hpp>
#include <stan/math/prim/scal/meta/scalar_seq_view.hpp>
#include <stan/math/prim/scal/meta/VectorBuilder.hpp>
#include <boost/random/uniform_01.hpp>
#include <boost/random/variate_generator.hpp>
#include <string>

namespace stan {
namespace math {

/**
* Return a pseudorandom double exponential variate with the given location
* and scale using the specified random number generator.
*
* mu and sigma can each be either scalars or vectors. All vector inputs
* must be the same length.
*
* @tparam T_loc Type of location parameter
* @tparam T_scale Type of scale parameter
* @tparam RNG class of random number generator
* @param mu location parameter
* @param sigma positive scale parameter
* @param mu (Sequence of) location parameter(s)
* @param sigma (Sequence of) scale parameter(s)
* @param rng random number generator
* @return double exponential random variate
* @throw std::domain_error if mu is infinite or sigma is nonpositive
* @throw std::invalid_argument if vector arguments are not the same length
*/
template <class RNG>
inline double
double_exponential_rng(double mu,
double sigma,
template <typename T_loc, typename T_scale, class RNG>
inline typename VectorBuilder<true, double, T_loc, T_scale>::type
double_exponential_rng(const T_loc& mu,
const T_scale& sigma,
RNG& rng) {
static const std::string function = "double_exponential_rng";

using boost::variate_generator;
using boost::random::uniform_01;
using std::log;
using std::abs;

scalar_seq_view<T_loc> mu_vec(mu);
scalar_seq_view<T_scale> sigma_vec(sigma);

check_finite(function, "Location parameter", mu);
check_positive_finite(function, "Scale parameter", sigma);
check_consistent_sizes(function,
"Location parameter", mu,
"Scale Parameter", sigma);

size_t N = max_size(mu, sigma);

VectorBuilder<true, double, T_loc, T_scale> output(N);

variate_generator<RNG&, uniform_01<> > rng_unit_01(rng, uniform_01<>());
for (size_t n = 0; n < N; n++) {
double laplaceRN = rng_unit_01();
Copy link
Contributor

Choose a reason for hiding this comment

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

[optional fix]

What I was suggesting in the original comment was something like this

double z = rng_unit_01();
double sign = bernoulli(0.5) ? 1 : -1;
output[n] = mu_vec[n] - sign * sigma_vec[n] * log(z)

A slightly more efficient version would be

double z = uniform(-1, 1);
output[n] = mu_vec[n] - sign(z) * sigma_vec[n] * log(abs(z))

For sign(), see: http://www.boost.org/doc/libs/1_62_0/libs/math/doc/html/math_toolkit/sign_functions.html

For uniform(-1, 1), see: http://www.boost.org/doc/libs/1_60_0/boost/random/uniform_real_distribution.hpp

For bernoulli(0.5), see: http://www.boost.org/doc/libs/1_60_0/boost/random/bernoulli_distribution.hpp

I got there through the version here with algebra, but it turns out to be pretty close to the Wikipedia version, only using uniform(-1, 1) and the recognition that 1 - uniform(-1, 1) has the same distribution as uniform(-1, 1) and the same sign randomization.

https://en.wikipedia.org/wiki/Laplace_distribution#Generating_random_variables_according_to_the_Laplace_distribution

double a = (0.5 - laplaceRN > 0) ? 1.0 : -1.0;
output[n] = mu_vec[n] -
sigma_vec[n] * a * log1m(2 * abs(0.5 - laplaceRN));
}

variate_generator<RNG&, uniform_01<> >
rng_unit_01(rng, uniform_01<>());
double a = 0;
double laplaceRN = rng_unit_01();
if (0.5 - laplaceRN > 0)
a = 1.0;
else if (0.5 - laplaceRN < 0)
a = -1.0;
return mu - sigma * a * log1m(2 * abs(0.5 - laplaceRN));
return output.data();
}
}
}
Expand Down
64 changes: 48 additions & 16 deletions stan/math/prim/scal/prob/exp_mod_normal_rng.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,38 +3,70 @@

#include <stan/math/prim/scal/err/check_consistent_sizes.hpp>
#include <stan/math/prim/scal/err/check_finite.hpp>
#include <stan/math/prim/scal/err/check_not_nan.hpp>
#include <stan/math/prim/scal/err/check_positive_finite.hpp>
#include <stan/math/prim/scal/fun/constants.hpp>
#include <stan/math/prim/scal/meta/length.hpp>
#include <stan/math/prim/scal/meta/max_size.hpp>
#include <stan/math/prim/scal/meta/VectorBuilder.hpp>
#include <stan/math/prim/scal/prob/normal_rng.hpp>
#include <stan/math/prim/scal/prob/exponential_rng.hpp>
#include <stan/math/prim/scal/meta/include_summand.hpp>
#include <stan/math/prim/scal/fun/value_of.hpp>
#include <stan/math/prim/scal/prob/normal_rng.hpp>
#include <boost/random/normal_distribution.hpp>
#include <boost/random/variate_generator.hpp>
#include <string>

namespace stan {
namespace math {

template <class RNG>
inline double
exp_mod_normal_rng(double mu,
double sigma,
double lambda,
/**
* Return a pseudorandom Exponentially modified normal variate for the
Copy link
Contributor

Choose a reason for hiding this comment

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

[optional fix]

I'd rather not capitalize anything other than proper names and names of distributions (not named after people) don't count.

* given location, scale, and inverse scale using the specified random
* number generator.
*
* mu, sigma, and lambda can each be either scalars or vectors. All vector
Copy link
Contributor

Choose a reason for hiding this comment

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

Can't they be scalars, standard vectors, or Eigen vectors?

* inputs must be the same length.
*
* @tparam T_loc Type of location parameter
* @tparam T_scale Type of scale parameter
* @tparam T_inv_scale Type of inverse scale parameter
* @tparam RNG type of random number generator
* @param mu (Sequence of) location parameter(s)
* @param sigma (Sequence of) scale parameter(s)
* @param lambda (Sequence of) inverse scale parameter(s)
* @param rng random number generator
* @return Exponentially modified normal random variate
* @throw std::domain_error if mu is infinite, sigma is nonpositive,
* or lambda is nonpositive
* @throw std::invalid_argument if vector arguments are not the same length
Copy link
Contributor

Choose a reason for hiding this comment

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

So I think this should say "container arguments" where you have "vector arguments".

*/
template <typename T_loc, typename T_scale, typename T_inv_scale, class RNG>
inline
typename VectorBuilder<true, double, T_loc, T_scale, T_inv_scale>::type
exp_mod_normal_rng(const T_loc& mu,
const T_scale& sigma,
const T_inv_scale& lambda,
RNG& rng) {
static const std::string function = "exp_mod_normal_rng";

scalar_seq_view<T_loc> mu_vec(mu);
scalar_seq_view<T_scale> sigma_vec(sigma);
scalar_seq_view<T_inv_scale> lambda_vec(lambda);

check_finite(function, "Location parameter", mu);
check_positive_finite(function, "Inv_scale parameter", lambda);
check_positive_finite(function, "Scale parameter", sigma);
check_positive_finite(function, "Inv_scale parameter", lambda);
check_consistent_sizes(function,
"Location parameter", mu,
"Scale Parameter", sigma,
"Inv_scale Parameter", lambda);

return normal_rng(mu, sigma, rng)
+ exponential_rng(lambda, rng);
}
size_t N = max_size(mu, sigma, lambda);

Copy link
Contributor

Choose a reason for hiding this comment

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

[optional fix]

I'd remove all the vertical space that's not absolutely necessary for readability. And that includes braces around for loops iwth single statements. Vertical space is precious!

VectorBuilder<true, double, T_loc, T_scale, T_inv_scale> output(N);

for (size_t n = 0; n < N; n++) {
Copy link
Contributor

Choose a reason for hiding this comment

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

[optional fix]

we prefer ++n unless you need the result of n++ (again, this is in Google's style guide, but not enforced by cpplint)

output[n] = normal_rng(mu_vec[n], sigma_vec[n], rng) +
Copy link
Contributor

Choose a reason for hiding this comment

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

This isn't going to be too much of an efficiency issue, but generally you don't want to call your public-facing functions internally because they're going to be doing redundant argument checks. It'll wind up creating more objects, too. Overally, though, this is probably OK like this and not worth optimizing now.

exponential_rng(lambda_vec[n], rng);
Copy link
Contributor

Choose a reason for hiding this comment

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

[optional fix]

I promised Sean not to require changes not legislated by cpplint, but I really think opertors should go on the following line for readability. You read code by scanning the structure down the left and it's lost with operators at theend of lines.

}

return output.data();
}
}
}
#endif
52 changes: 33 additions & 19 deletions stan/math/prim/scal/prob/gumbel_rng.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,50 +3,64 @@

#include <stan/math/prim/scal/err/check_consistent_sizes.hpp>
#include <stan/math/prim/scal/err/check_finite.hpp>
#include <stan/math/prim/scal/err/check_not_nan.hpp>
#include <stan/math/prim/scal/err/check_positive.hpp>
#include <stan/math/prim/scal/meta/length.hpp>
#include <stan/math/prim/scal/err/check_positive_finite.hpp>
#include <stan/math/prim/scal/meta/max_size.hpp>
#include <stan/math/prim/scal/meta/scalar_seq_view.hpp>
#include <stan/math/prim/scal/meta/VectorBuilder.hpp>
#include <stan/math/prim/scal/meta/return_type.hpp>
#include <stan/math/prim/scal/fun/constants.hpp>
#include <stan/math/prim/scal/meta/include_summand.hpp>
#include <stan/math/prim/scal/fun/value_of.hpp>
#include <boost/random/uniform_01.hpp>
#include <boost/random/variate_generator.hpp>
#include <string>

namespace stan {
namespace math {

/**
* Return a pseudorandom Gumbel variate with the given location and scale
* using the specified random number generator.
*
* mu and sigma can each be either scalars or vectors. All vector inputs
* must be the same length.
*
* @tparam T_loc Type of location parameter
* @tparam T_scale Type of scale parameter
* @tparam RNG type of random number generator
* @param mu location parameter
* @param beta positive scale parameter
* @param mu (Sequence of) location parameter(s)
* @param beta (Sequence of) scale parameter(s)
* @param rng random number generator
* @return Gumbel random variate
* @throw std::domain_error if mu is infinite or beta is nonpositive.
* @throw std::invalid_argument if vector arguments are not the same length
*/
template <class RNG>
inline double
gumbel_rng(double mu,
double beta,
template <typename T_loc, typename T_scale, class RNG>
inline typename VectorBuilder<true, double, T_loc, T_scale>::type
gumbel_rng(const T_loc& mu,
const T_scale& beta,
Copy link
Contributor

Choose a reason for hiding this comment

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

[optional fix]

pack as many arguments onto a line as possible. Remember, vertical space is precious. (again, in Google style guide, not in cpplint)

RNG& rng) {
using boost::variate_generator;
using boost::uniform_01;

static const std::string function = "gumbel_rng";

check_finite(function, "Location parameter", mu);
check_positive(function, "Scale parameter", beta);
check_positive_finite(function, "Scale parameter", beta);
check_consistent_sizes(function,
"Location parameter", mu,
"Scale Parameter", beta);

variate_generator<RNG&, uniform_01<> >
uniform01_rng(rng, uniform_01<>());
return mu - beta * std::log(-std::log(uniform01_rng()));
}
scalar_seq_view<T_loc> mu_vec(mu);
scalar_seq_view<T_scale> beta_vec(beta);

size_t N = max_size(mu, beta);

VectorBuilder<true, double, T_loc, T_scale> output(N);

variate_generator<RNG&, uniform_01<> > uniform01_rng(rng, uniform_01<>());
for (size_t n = 0; n < N; n++) {
output[n] = mu_vec[n] -
Copy link
Contributor

Choose a reason for hiding this comment

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

[optional]

Same comments as last time on vertical space and operators.

beta_vec[n] * std::log(-std::log(uniform01_rng()));
}

return output.data();
}
}
}
#endif
Loading