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	mayorov	1.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			}