1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
|
// $Id: betaDistribution.cpp 3985 2008-05-11 11:00:44Z adido $
#include "betaDistribution.h"
#include "gammaUtilities.h"
#include "betaUtilities.h"
#include "errorMsg.h"
#include "logFile.h"
#include <cmath>
betaDistribution::betaDistribution()
{
_alpha = 0.0;
_beta = 0.0;
_boundary.resize(0,0);
_rates.resize(0,0);
_ratesProb.resize(0,0);
_globalRate = 1;//??? 0.5 or 1
_discretizationType = MEDIAN;
}
// note that the order of initalization makes a diffrence.
betaDistribution::betaDistribution(const betaDistribution& other) :
_boundary(other._boundary),
_alpha(other._alpha),
_beta(other._beta),
_rates(other._rates),
_ratesProb(other._ratesProb),
_globalRate(other._globalRate),
_discretizationType(other._discretizationType){
}
betaDistribution::betaDistribution(MDOUBLE alpha,MDOUBLE beta,int in_number_of_categories,discretizationType in_discretizationType) :distribution(){
_globalRate=1.0;
_discretizationType = in_discretizationType;
setBetaParameters(in_number_of_categories,alpha,beta);
}
betaDistribution::~betaDistribution() {
_boundary.clear();
_rates.clear();
_ratesProb.clear();
}
void betaDistribution::setAlpha(MDOUBLE in_alpha) {
if (in_alpha == _alpha)
return;
setBetaParameters(categories(), in_alpha, _beta);
}
void betaDistribution::setBeta(MDOUBLE in_beta) {
if (in_beta == _beta)
return;
setBetaParameters( categories(), _alpha, in_beta);
}
void betaDistribution::setDiscretizationType(discretizationType in_discretizationType) {
if (in_discretizationType == _discretizationType)
return;
_discretizationType = in_discretizationType;
if (categories() > 1)
fill_rates();
}
void betaDistribution::change_number_of_categories(int in_number_of_categories) {
if (in_number_of_categories == categories())
return;
setBetaParameters( in_number_of_categories, _alpha, _beta);
}
void betaDistribution::setBetaParameters(int in_number_of_categories, MDOUBLE in_alpha, MDOUBLE in_beta) {
if ((in_alpha == _alpha) && (in_beta == _beta) && (in_number_of_categories == categories()))
return;
if (in_alpha < MINIMUM_ALPHA_PARAM)
in_alpha = MINIMUM_ALPHA_PARAM;// when alpha is very small there are underflaw problems
if (in_beta < MINIMUM_ALPHA_PARAM)
in_beta = MINIMUM_ALPHA_PARAM;// when beta is very small there are underflaw problems
_alpha = in_alpha;
_beta = in_beta;
_rates.clear();
_rates.resize(in_number_of_categories);
_ratesProb.clear();
_ratesProb.resize(in_number_of_categories, 1.0/in_number_of_categories);
_boundary.clear();
_boundary.resize(in_number_of_categories+1);
if (in_number_of_categories==1) {
_rates[0] = 1.0;
return;
}
if (categories() > 1) {
fill_rates();
return ;
}
}
int betaDistribution::fill_rates() {
fill_boundaries();
int i;
//LOG(5,<<endl<<" alpha = "<<_alpha<<" beta = "<< _beta<<endl);
//for (i=0; i<=categories(); ++i) cout<<endl<<_boundary[i];
//LOG(5,<<"\n====== the r categories are =====\n");
for (i=0; i<categories(); ++i) {
if (_discretizationType == MEAN)
_rates[i]=computeAverage_r(_boundary[i], _boundary[i+1], _alpha, _beta, categories());
else //_discretizationType == MEDIAN
_rates[i] =inverseCDFBeta(_alpha, _beta,static_cast<MDOUBLE>(i*2 +1)/(2*categories()));
//LOG(5,<<_rates[i]<<endl);
}
//LOG(5,<<endl<<_alpha<<endl);
return 0;
}
int betaDistribution::fill_boundaries() {
int i;
//LOG(5,<<endl<<"========BOUNDARY============="<<endl);
for (i=1; i<categories(); ++i)
{
_boundary[i]=inverseCDFBeta(_alpha, _beta,static_cast<MDOUBLE>(i)/categories());
//LOG(5,<<"_boundary[ "<<i<<"] ="<<_boundary[i]<<endl);
}
_boundary[0]=0;
_boundary[i]=1;
return 0;
}
const MDOUBLE betaDistribution::getCumulativeProb(const MDOUBLE x) const
{//
//since r~gamma(alpha, beta) then beta*r~ gamma(alpha,1)=gammp
//here we assume alpha=beta
return incompleteBeta(_alpha,_beta,x);
}
|
