EigenRand  0.4.0-alpha
GammaPoisson.h
Go to the documentation of this file.
1 
12 #ifndef EIGENRAND_DISTS_GAMMAPOISSON_H
13 #define EIGENRAND_DISTS_GAMMAPOISSON_H
14 
15 #include <memory>
16 #include <iterator>
17 #include <limits>
18 
19 namespace Eigen
20 {
21  namespace Rand
22  {
28  template<typename _Scalar>
29  class NegativeBinomialGen : public GenBase<NegativeBinomialGen<_Scalar>, _Scalar>
30  {
31  static_assert(std::is_same<_Scalar, int32_t>::value, "negativeBinomial needs integral types.");
33  GammaGen<float> gamma;
34  public:
35  using Scalar = _Scalar;
36 
43  NegativeBinomialGen(_Scalar _trials = 1, double _p = 0.5)
44  : gamma{ (float)_trials, (float)((1 - _p) / _p) }
45 
46  {
47  }
48 
49  template<typename Rng>
50  EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng)
51  {
52  using namespace Eigen::internal;
53  float v = gamma(rng);
54  return PoissonGen<_Scalar>{v}(rng);
55  }
56 
57  template<typename Packet, typename Rng>
58  EIGEN_STRONG_INLINE const Packet packetOp(Rng&& rng)
59  {
60  using namespace Eigen::internal;
61  using PacketType = decltype(reinterpret_to_float(std::declval<Packet>()));
62 
63  auto mean = gamma.template packetOp<PacketType>(rng);
64  auto res = pset1<Packet>(0);
65  PacketType val = pset1<PacketType>(1), pne_mean = pexp(pnegate(mean));
66  if (pmovemask(pcmplt(pset1<PacketType>(12), mean)) == 0)
67  {
68  for (int _i = 0; ; ++_i)
69  {
70  EIGENRAND_CHECK_INFINITY_LOOP();
71  val = pmul(val, ur.template packetOp<PacketType>(rng));
72  auto c = reinterpret_to_int(pcmplt(pne_mean, val));
73  if (pmovemask(c) == 0) break;
74  res = padd(res, pnegate(c));
75  }
76  return res;
77  }
78  else
79  {
80  auto& cm = Rand::detail::CompressMask<sizeof(Packet)>::get_inst();
81  const PacketType ppi = pset1<PacketType>(constant::pi),
82  psqrt_tmean = psqrt(pmul(pset1<PacketType>(2), mean)),
83  plog_mean = plog(mean),
84  pg1 = psub(pmul(mean, plog_mean), plgamma_approx(padd(mean, pset1<PacketType>(1))));
85  for (int _i = 0; ; ++_i)
86  {
87  EIGENRAND_CHECK_INFINITY_LOOP();
88  PacketType fres, yx, psin, pcos;
89  psincos(pmul(ppi, ur.template packetOp<PacketType>(rng)), psin, pcos);
90  yx = pdiv(psin, pcos);
91  fres = ptruncate(padd(pmul(psqrt_tmean, yx), mean));
92 
93  auto p1 = pmul(padd(pmul(yx, yx), pset1<PacketType>(1)), pset1<PacketType>(0.9));
94  auto p2 = pexp(psub(psub(pmul(fres, plog_mean), plgamma_approx(padd(fres, pset1<PacketType>(1)))), pg1));
95 
96  auto c1 = pcmple(pset1<PacketType>(0), fres);
97  auto c2 = pcmple(ur.template packetOp<PacketType>(rng), pmul(p1, p2));
98 
99  auto cands = fres;
100  bool full = false;
101  gamma.cache_rest_cnt = cm.compress_append(cands, pand(c1, c2),
102  gamma.template get<PacketType>(), gamma.cache_rest_cnt, full);
103  if (full) return pcast<PacketType, Packet>(cands);
104  }
105  }
106  }
107  };
108 
109  template<typename Derived, typename Urng>
110  using NegativeBinomialType = CwiseNullaryOp<internal::scalar_rng_adaptor<NegativeBinomialGen<typename Derived::Scalar>, typename Derived::Scalar, Urng, true>, const Derived>;
111 
126  template<typename Derived, typename Urng>
127  inline const NegativeBinomialType<Derived, Urng>
128  negativeBinomial(Index rows, Index cols, Urng&& urng, typename Derived::Scalar trials = 1, double p = 0.5)
129  {
130  return {
131  rows, cols, { std::forward<Urng>(urng), NegativeBinomialGen<typename Derived::Scalar>{trials, p} }
132  };
133  }
134 
148  template<typename Derived, typename Urng>
149  inline const NegativeBinomialType<Derived, Urng>
150  negativeBinomialLike(Derived& o, Urng&& urng, typename Derived::Scalar trials = 1, double p = 0.5)
151  {
152  return {
153  o.rows(), o.cols(), { std::forward<Urng>(urng), NegativeBinomialGen<typename Derived::Scalar>{trials, p} }
154  };
155  }
156  }
157 }
158 #endif
Base class of all univariate random generators.
Definition: Basic.h:33
Generator of integers on a negative binomial distribution.
Definition: GammaPoisson.h:30
NegativeBinomialGen(_Scalar _trials=1, double _p=0.5)
Construct a new negative binomial generator.
Definition: GammaPoisson.h:43
Generator of integers on a Poisson distribution.
Definition: Discrete.h:735
const GammaType< Derived, Urng > gamma(Index rows, Index cols, Urng &&urng, typename Derived::Scalar alpha=1, typename Derived::Scalar beta=1)
generates reals on a gamma distribution with arbitrary shape and scale parameter.
Definition: NormalExp.h:1158
const NegativeBinomialType< Derived, Urng > negativeBinomialLike(Derived &o, Urng &&urng, typename Derived::Scalar trials=1, double p=0.5)
generates reals on the negative binomial distribution.
Definition: GammaPoisson.h:150
const NegativeBinomialType< Derived, Urng > negativeBinomial(Index rows, Index cols, Urng &&urng, typename Derived::Scalar trials=1, double p=0.5)
generates reals on the negative binomial distribution.
Definition: GammaPoisson.h:128