EigenRand  0.3.0
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 ur_base = UniformRealGen<float>;
62  using PacketType = decltype(reinterpret_to_float(std::declval<Packet>()));
63 
64  auto mean = gamma.template packetOp<PacketType>(rng);
65  auto res = pset1<Packet>(0);
66  PacketType val = pset1<PacketType>(1), pne_mean = pexp(pnegate(mean));
67  if (pmovemask(pcmplt(pset1<PacketType>(12), mean)) == 0)
68  {
69  while (1)
70  {
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(padd(mean, pset1<PacketType>(1))));
85  while (1)
86  {
87  PacketType fres, yx, psin, pcos;
88  psincos(pmul(ppi, ur.template packetOp<PacketType>(rng)), psin, pcos);
89  yx = pdiv(psin, pcos);
90  fres = ptruncate(padd(pmul(psqrt_tmean, yx), mean));
91 
92  auto p1 = pmul(padd(pmul(yx, yx), pset1<PacketType>(1)), pset1<PacketType>(0.9));
93  auto p2 = pexp(psub(psub(pmul(fres, plog_mean), plgamma(padd(fres, pset1<PacketType>(1)))), pg1));
94 
95  auto c1 = pcmple(pset1<PacketType>(0), fres);
96  auto c2 = pcmple(ur.template packetOp<PacketType>(rng), pmul(p1, p2));
97 
98  auto cands = fres;
99  bool full = false;
100  gamma.cache_rest_cnt = cm.compress_append(cands, pand(c1, c2),
101  gamma.template get<PacketType>(), gamma.cache_rest_cnt, full);
102  if (full) return pcast<PacketType, Packet>(cands);
103  }
104  }
105  }
106  };
107 
108  template<typename Derived, typename Urng>
109  using NegativeBinomialType = CwiseNullaryOp<internal::scalar_rng_adaptor<NegativeBinomialGen<typename Derived::Scalar>, typename Derived::Scalar, Urng, true>, const Derived>;
110 
125  template<typename Derived, typename Urng>
126  inline const NegativeBinomialType<Derived, Urng>
127  negativeBinomial(Index rows, Index cols, Urng&& urng, typename Derived::Scalar trials = 1, double p = 0.5)
128  {
129  return {
130  rows, cols, { std::forward<Urng>(urng), NegativeBinomialGen<typename Derived::Scalar>{trials, p} }
131  };
132  }
133 
147  template<typename Derived, typename Urng>
148  inline const NegativeBinomialType<Derived, Urng>
149  negativeBinomialLike(Derived& o, Urng&& urng, typename Derived::Scalar trials = 1, double p = 0.5)
150  {
151  return {
152  o.rows(), o.cols(), { std::forward<Urng>(urng), NegativeBinomialGen<typename Derived::Scalar>{trials, p} }
153  };
154  }
155  }
156 }
157 #endif
Eigen::Rand::GenBase
Base class of all univariate random generators.
Definition: Basic.h:33
Eigen::Rand::GammaGen< float >
Eigen::Rand::negativeBinomial
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:127
Eigen::Rand::NegativeBinomialGen::NegativeBinomialGen
NegativeBinomialGen(_Scalar _trials=1, double _p=0.5)
Construct a new negative binomial generator.
Definition: GammaPoisson.h:43
Eigen::Rand::UniformRealGen< float >
Eigen::Rand::gamma
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:1118
Eigen::Rand::NegativeBinomialGen
Generator of integers on a negative binomial distribution.
Definition: GammaPoisson.h:30
Eigen::Rand::negativeBinomialLike
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:149