PamUnfold/src/RanGen.cpp

#include <RanGen.h>

using namespace std;

ClassImp(RanGen);

RanGen::RanGen(Int_t seed) : TRandom3(seed){
}

RanGen::~RanGen(){};

Double_t RanGen::Gamma(Int_t n, Double_t theta){

  Double_t value = 0;

  for(Int_t i=0; i<n; i++) value -= TMath::Log( this->Uniform() );

  return value;

}


vector<Double_t> RanGen::Dirichlet(Int_t n, vector<Int_t> x){

  Double_t sum = 0;
  vector<Double_t> gammas(n);
  vector<Double_t> value(n);

  for(Int_t i=0; i<n; i++){
    gammas[i] = Gamma(x[i], 1);
    sum += gammas[i];
  }

  for(Int_t i=0; i<n; i++){
    if(sum)
      value[i] = gammas[i]/sum;
    else
      value[i] = 0;
  }

  return value;

}


Int_t RanGen::Binomial(Int_t n, Double_t p){
  unsigned int i, a, b, k = 0;
  while (n > 10)        /* This parameter is tunable */
    {
      double X;
      a = 1 + (n / 2);
      b = 1 + n - a;
      
      X = this->Beta((Int_t) a, (Int_t) b);
      
      if (X >= p)
        {
          n = a - 1;
          p /= X;
        }
      else
        {
          k += a;
          n = b - 1;
          p = (p - X) / (1 - X);
        }
    }
  
  for (i = 0; i < n; i++)
    {
      double u = this->Uniform();
      if (u < p)
        k++;
    }
  
  return k;
  
}

Double_t RanGen::Beta(const Int_t a, const Int_t b){
  double x1 = this->Gamma(a, 1.0);
  double x2 = this->Gamma(b, 1.0);
  
  return x1 / (x1 + x2);
}

vector<Int_t> RanGen::Multinomial(Int_t N, vector<Double_t> p){

/*  From the GSL Library:
    The multinomial distribution has the form
   
                                       N!           n_1  n_2      n_K
    prob(n_1, n_2, ... n_K) = -------------------- p_1  p_2  ... p_K
                              (n_1! n_2! ... n_K!) 
 
    where n_1, n_2, ... n_K are nonnegative integers, sum_{k=1,K} n_k = N,
    and p = (p_1, p_2, ..., p_K) is a probability distribution. 
 
    Random variates are generated using the conditional binomial method.
    This scales well with N and does not require a setup step.
 
    Ref: 
    C.S. David, The computer generation of multinomial random variates,
    Comp. Stat. Data Anal. 16 (1993) 205-217
 */

  size_t K = p.size();
  vector<Int_t> n(K);

  size_t k;
  double norm = 0.0;
  double sum_p = 0.0;
  
  unsigned int sum_n = 0;
  
  /* p[k] may contain non-negative weights that do not sum to 1.0.
   * Even a probability distribution will not exactly sum to 1.0
   * due to rounding errors. 
   */
  
  for (k = 0; k < K; k++)
    {
      norm += p[k];
    }
  
  for (k = 0; k < K; k++)
    {
      if (p[k] > 0.0)
        {
          n[k] = this->Binomial(N - sum_n, p[k] / (norm - sum_p));
        }
      else
        {
          n[k] = 0;
        }
      
      sum_p += p[k];
      sum_n += n[k];
    }
  return n;
  
}
1	#include <RanGen.h>
2
3	using namespace std;
4
5	ClassImp(RanGen);
6
7	RanGen::RanGen(Int_t seed) : TRandom3(seed){
8	}
9
10	RanGen::~RanGen(){};
11
12	Double_t RanGen::Gamma(Int_t n, Double_t theta){
13
14	Double_t value = 0;
15
16	for(Int_t i=0; i<n; i++) value -= TMath::Log( this->Uniform() );
17
18	return value;
19
20	}
21
22
23	vector<Double_t> RanGen::Dirichlet(Int_t n, vector<Int_t> x){
24
25	Double_t sum = 0;
26	vector<Double_t> gammas(n);
27	vector<Double_t> value(n);
28
29	for(Int_t i=0; i<n; i++){
30	gammas[i] = Gamma(x[i], 1);
31	sum += gammas[i];
32	}
33
34	for(Int_t i=0; i<n; i++){
35	if(sum)
36	value[i] = gammas[i]/sum;
37	else
38	value[i] = 0;
39	}
40
41	return value;
42
43	}
44
45
46	Int_t RanGen::Binomial(Int_t n, Double_t p){
47	unsigned int i, a, b, k = 0;
48	while (n > 10) /* This parameter is tunable */
49	{
50	double X;
51	a = 1 + (n / 2);
52	b = 1 + n - a;
53
54	X = this->Beta((Int_t) a, (Int_t) b);
55
56	if (X >= p)
57	{
58	n = a - 1;
59	p /= X;
60	}
61	else
62	{
63	k += a;
64	n = b - 1;
65	p = (p - X) / (1 - X);
66	}
67	}
68
69	for (i = 0; i < n; i++)
70	{
71	double u = this->Uniform();
72	if (u < p)
73	k++;
74	}
75
76	return k;
77
78	}
79
80	Double_t RanGen::Beta(const Int_t a, const Int_t b){
81	double x1 = this->Gamma(a, 1.0);
82	double x2 = this->Gamma(b, 1.0);
83
84	return x1 / (x1 + x2);
85	}
86
87	vector<Int_t> RanGen::Multinomial(Int_t N, vector<Double_t> p){
88
89	/* From the GSL Library:
90	The multinomial distribution has the form
91
92	N! n_1 n_2 n_K
93	prob(n_1, n_2, ... n_K) = -------------------- p_1 p_2 ... p_K
94	(n_1! n_2! ... n_K!)
95
96	where n_1, n_2, ... n_K are nonnegative integers, sum_{k=1,K} n_k = N,
97	and p = (p_1, p_2, ..., p_K) is a probability distribution.
98
99	Random variates are generated using the conditional binomial method.
100	This scales well with N and does not require a setup step.
101
102	Ref:
103	C.S. David, The computer generation of multinomial random variates,
104	Comp. Stat. Data Anal. 16 (1993) 205-217
105	*/
106
107	size_t K = p.size();
108	vector<Int_t> n(K);
109
110	size_t k;
111	double norm = 0.0;
112	double sum_p = 0.0;
113
114	unsigned int sum_n = 0;
115
116	/* p[k] may contain non-negative weights that do not sum to 1.0.
117	* Even a probability distribution will not exactly sum to 1.0
118	* due to rounding errors.
119	*/
120
121	for (k = 0; k < K; k++)
122	{
123	norm += p[k];
124	}
125
126	for (k = 0; k < K; k++)
127	{
128	if (p[k] > 0.0)
129	{
130	n[k] = this->Binomial(N - sum_n, p[k] / (norm - sum_p));
131	}
132	else
133	{
134	n[k] = 0;
135	}
136
137	sum_p += p[k];
138	sum_n += n[k];
139	}
140	return n;
141
142	}