25 #include "bias_config.h"
26 #include <Base/Common/W32Compat.hh>
27 #include <MathAlgo/GSL.hh>
28 #include <Base/Debug/Debug.hh>
29 #include <Base/Debug/Error.hh>
30 #include <Base/Debug/Exception.hh>
31 #include <Base/Common/CompareFloatingPoint.hh>
37 # include <gsl/gsl_cdf.h>
38 # include <gsl/gsl_randist.h>
39 # include <gsl/gsl_sf_gamma.h>
52 BEXCEPTION(
"ChiSquareProbDensFun: invalid argument, x must be > 0: "<<x);
54 BEXCEPTION(
"ChiSquareProbDensFun: invalid argument, deg_fred must be > 0 "
56 double res = gsl_ran_chisq_pdf(x, deg_freedom);
57 BIASASSERT(!isnan(res) && !isinf(res));
60 BIASERR(
"recompile BIAS with gsl support enabled");
70 double res = gsl_cdf_chisq_P(x, deg_freedom);
71 BIASASSERT(!isnan(res) && !isinf(res));
74 BIASERR(
"recompile BIAS with gsl support enabled");
84 double res = gsl_cdf_chisq_Pinv(x, deg_freedom);
85 BIASASSERT(!isnan(res) && !isinf(res));
88 BIASERR(
"recompile BIAS with gsl support enabled");
98 double res = gsl_sf_gamma (x);
99 BIASASSERT(!isnan(res) && !isinf(res));
102 BIASERR(
"recompile BIAS with gsl support enabled");
112 double res = gsl_sf_gammainv (x);
113 BIASASSERT(!isnan(res) && !isinf(res));
116 BIASERR(
"recompile BIAS with gsl support enabled");
126 if ((a+b)<=GSL_SF_GAMMA_XMAX){
127 double res = gsl_ran_beta_pdf(x, a, b);
128 BIASASSERT(!isnan(res) && !isinf(res));
131 BEXCEPTION(
"a or b out of range: a="<<a<<
" b="<<b<<endl);
134 BIASERR(
"recompile BIAS with gsl support enabled");
145 double res = gsl_cdf_beta_P(x, a, b);
147 if (isnan(res) || isinf(res)){
148 cout <<
"a: "<<a<<
" b: "<<b<<
" x: "<<x<<endl;
150 cout <<
"a: "<<a<<
" b: "<<b<<
" x: "<<x<<endl;
152 BIASASSERT(!isnan(res) && !isinf(res));
155 BIASERR(
"recompile BIAS with gsl support enabled");
166 double res = gsl_cdf_gaussian_P(x, sigma);
168 if (isnan(res) || isinf(res)){
169 cout <<
"sigma: "<<sigma<<
" x: "<<x<<
" res: "<<res<<endl;
171 cout <<
"sigma: "<<sigma<<
" x: "<<x<<
" res: "<<res<<endl;
173 BIASASSERT(!isnan(res) && !isinf(res));
176 BIASERR(
"recompile BIAS with gsl support enabled");
187 if ((a+b)<=GSL_SF_GAMMA_XMAX){
188 double res = gsl_ran_beta_pdf(x, a, b);
189 BIASASSERT(!isnan(res) && !isinf(res));
191 }
else if (
Equal(a,b)){
198 double sigma_sq = (a*b)/(a+b)/(a+b)/(a+b+1);
199 double res = 1.0/sqrt(sigma_sq * 2.0 * M_PI) *
200 exp(-(dist*dist)/(2.0 * sigma_sq));
203 cout <<
"approximating beta distribution by gaussian: alpha = beta = "<<a
204 <<
" x = "<<x<<
" sigma: "<<sqrt(sigma_sq)<<
" result: "<<res<<endl;
208 BIASERR(
"a or b to big, cannot compute beta distribution\n");
213 BIASERR(
"recompile BIAS with gsl support enabled");
BIASMathAlgo_EXPORT double BetaDistribution(double a, double b, double x)
BIASMathAlgo_EXPORT double ChiSquareProbDensFun(double x, int deg_freedom)
BIASMathAlgo_EXPORT double CulmulativeNormalDistribution(double x, double sigma)
culmulative probability function of normal distribution
BIASMathAlgo_EXPORT double InvChiSquareCulmProbFun(double x, int deg_freedom)
BIASMathAlgo_EXPORT double Gamma(double x)
BIASMathAlgo_EXPORT double InvGamma(double x)
BIASMathAlgo_EXPORT double CulmBetaDistribution(double a, double b, double x)
BIASMathAlgo_EXPORT double ChiSquareCulmProbFun(double x, int deg_freedom)
bool Equal(const T left, const T right, const T eps)
comparison function for floating point values See http://www.boost.org/libs/test/doc/components/test_...
BIASMathAlgo_EXPORT double BetaDistributionApprox(double a, double b, double x)