EigenRand  0.4.0-alpha
Multinomial.h
Go to the documentation of this file.
1 
12 #ifndef EIGENRAND_MVDISTS_MULTINOMIAL_H
13 #define EIGENRAND_MVDISTS_MULTINOMIAL_H
14 
15 namespace Eigen
16 {
17  namespace Rand
18  {
25  template<typename _Scalar = int32_t, Index Dim = -1>
26  class MultinomialGen : public MvVecGenBase<MultinomialGen<_Scalar, Dim>, _Scalar, Dim>
27  {
28  static_assert(std::is_same<_Scalar, int32_t>::value, "`MultinomialGen` needs integral types.");
29  _Scalar trials;
30  Matrix<double, Dim, 1> probs;
31  DiscreteGen<_Scalar> discrete;
32  public:
40  template<typename WeightTy>
41  MultinomialGen(_Scalar _trials, const MatrixBase<WeightTy>& _weights)
42  : trials{ _trials }, probs{ _weights.template cast<double>() }, discrete(probs.data(), probs.data() + probs.size())
43  {
44  eigen_assert(_weights.cols() == 1);
45  for (Index i = 0; i < probs.size(); ++i)
46  {
47  eigen_assert(probs[i] >= 0);
48  }
49  probs /= probs.sum();
50  }
51 
52  MultinomialGen(const MultinomialGen&) = default;
53  MultinomialGen(MultinomialGen&&) = default;
54 
55  MultinomialGen& operator=(const MultinomialGen&) = default;
56  MultinomialGen& operator=(MultinomialGen&&) = default;
57 
58  Index dims() const { return probs.rows(); }
59 
60  template<typename Urng>
61  inline Matrix<_Scalar, Dim, -1> generate(Urng&& urng, Index samples)
62  {
63  const Index dim = probs.size();
64  Matrix<_Scalar, Dim, -1> ret(dim, samples);
65  //if (trials < 2500)
66  {
67  for (Index j = 0; j < samples; ++j)
68  {
69  ret.col(j) = generate(urng);
70  }
71  }
72  /*else
73  {
74  ret.row(0) = binomial<Matrix<_Scalar, -1, 1>>(samples, 1, urng, trials, probs[0]).eval().transpose();
75  for (Index j = 0; j < samples; ++j)
76  {
77  double rest_p = 1 - probs[0];
78  _Scalar t = trials - ret(0, j);
79  for (Index i = 1; i < dim - 1; ++i)
80  {
81  ret(i, j) = binomial<Matrix<_Scalar, 1, 1>>(1, 1, urng, t, probs[i] / rest_p)(0);
82  t -= ret(i, j);
83  rest_p -= probs[i];
84  }
85  ret(dim - 1, j) = 0;
86  }
87  ret.row(dim - 1).setZero();
88  ret.row(dim - 1).array() = trials - ret.colwise().sum().array();
89  }*/
90  return ret;
91  }
92 
93  template<typename Urng>
94  inline Matrix<_Scalar, Dim, 1> generate(Urng&& urng)
95  {
96  const Index dim = probs.size();
97  Matrix<_Scalar, Dim, 1> ret(dim);
98  //if (trials < 2500)
99  {
100  ret.setZero();
101  auto d = discrete.template generate<Matrix<_Scalar, -1, 1>>(trials, 1, urng).eval();
102  for (Index i = 0; i < trials; ++i)
103  {
104  ret[d[i]] += 1;
105  }
106  }
107  /*else
108  {
109  double rest_p = 1;
110  _Scalar t = trials;
111  for (Index i = 0; i < dim - 1; ++i)
112  {
113  ret[i] = binomial<Matrix<_Scalar, 1, 1>>(1, 1, urng, t, probs[i] / rest_p)(0);
114  t -= ret[i];
115  rest_p -= probs[i];
116  }
117  ret[dim - 1] = 0;
118  ret[dim - 1] = trials - ret.sum();
119  }*/
120  return ret;
121  }
122  };
123 
134  template<typename IntTy, typename WeightTy>
135  inline auto makeMultinomialGen(IntTy trials, const MatrixBase<WeightTy>& probs)
137  {
139  }
140 
147  template<typename _Scalar, Index Dim = -1>
148  class DirichletGen : public MvVecGenBase<DirichletGen<_Scalar, Dim>, _Scalar, Dim>
149  {
150  static_assert(std::is_floating_point<_Scalar>::value, "`DirichletGen` needs floating point types.");
151 
152  Matrix<_Scalar, Dim, 1> alpha;
153  std::vector<GammaGen<_Scalar>> gammas;
154  public:
161  template<typename AlphaTy>
162  DirichletGen(const MatrixBase<AlphaTy>& _alpha)
163  : alpha{ _alpha }
164  {
165  eigen_assert(_alpha.cols() == 1);
166  for (Index i = 0; i < alpha.size(); ++i)
167  {
168  eigen_assert(alpha[i] > 0);
169  gammas.emplace_back(alpha[i]);
170  }
171  }
172 
173  DirichletGen(const DirichletGen&) = default;
174  DirichletGen(DirichletGen&&) = default;
175 
176  Index dims() const { return alpha.rows(); }
177 
178  template<typename Urng>
179  inline Matrix<_Scalar, Dim, -1> generate(Urng&& urng, Index samples)
180  {
181  const Index dim = alpha.size();
182  Matrix<_Scalar, Dim, -1> ret(dim, samples);
183  Matrix<_Scalar, -1, 1> tmp(samples);
184  for (Index i = 0; i < dim; ++i)
185  {
186  tmp = gammas[i].generateLike(tmp, urng);
187  ret.row(i) = tmp.transpose();
188  }
189  ret.array().rowwise() /= ret.array().colwise().sum();
190  return ret;
191  }
192 
193  template<typename Urng>
194  inline Matrix<_Scalar, Dim, 1> generate(Urng&& urng)
195  {
196  const Index dim = alpha.size();
197  Matrix<_Scalar, Dim, 1> ret(dim);
198  for (Index i = 0; i < dim; ++i)
199  {
200  ret[i] = gammas[i].template generate<Matrix<_Scalar, 1, 1>>(1, 1, urng)(0);
201  }
202  ret /= ret.sum();
203  return ret;
204  }
205  };
206 
215  template<typename AlphaTy>
216  inline auto makeDirichletGen(const MatrixBase<AlphaTy>& alpha)
217  -> DirichletGen<typename MatrixBase<AlphaTy>::Scalar, MatrixBase<AlphaTy>::RowsAtCompileTime>
218  {
219  return DirichletGen<typename MatrixBase<AlphaTy>::Scalar, MatrixBase<AlphaTy>::RowsAtCompileTime>{ alpha };
220  }
221  }
222 }
223 
224 #endif
Generator of reals on a Dirichlet distribution.
Definition: Multinomial.h:149
DirichletGen(const MatrixBase< AlphaTy > &_alpha)
Construct a new Dirichlet generator.
Definition: Multinomial.h:162
Generator of integers on the interval [0, n), where the probability of each individual integer i is p...
Definition: Discrete.h:329
Generator of real vectors on a multinomial distribution.
Definition: Multinomial.h:27
MultinomialGen(_Scalar _trials, const MatrixBase< WeightTy > &_weights)
Construct a new multinomial generator.
Definition: Multinomial.h:41
Base class of all multivariate random vector generators.
Definition: Basic.h:84
auto makeMultinomialGen(IntTy trials, const MatrixBase< WeightTy > &probs) -> MultinomialGen< IntTy, MatrixBase< WeightTy >::RowsAtCompileTime >
helper function constructing Eigen::Rand::MultinomialGen
Definition: Multinomial.h:135
auto makeDirichletGen(const MatrixBase< AlphaTy > &alpha) -> DirichletGen< typename MatrixBase< AlphaTy >::Scalar, MatrixBase< AlphaTy >::RowsAtCompileTime >
helper function constructing Eigen::Rand::DirichletGen
Definition: Multinomial.h:216
const DiscreteType< Derived, Urng > discrete(Index rows, Index cols, Urng &&urng, RealIter first, RealIter last)
generates random integers on the interval [0, n), where the probability of each individual integer i ...
Definition: Discrete.h:1302