1 |
/* specfunc/beta.c |
/* cdf/cdf_beta.c |
2 |
* |
* |
3 |
* Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman |
* Copyright (C) 2003 Brian Gough. |
4 |
* |
* |
5 |
* This program is free software; you can redistribute it and/or modify |
* This program is free software; you can redistribute it and/or modify |
6 |
* it under the terms of the GNU General Public License as published by |
* it under the terms of the GNU General Public License as published by |
14 |
* |
* |
15 |
* You should have received a copy of the GNU General Public License |
* You should have received a copy of the GNU General Public License |
16 |
* along with this program; if not, write to the Free Software |
* along with this program; if not, write to the Free Software |
17 |
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. |
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA. |
|
*/ |
|
|
|
|
|
/* Author: G. Jungman |
|
|
* Tweaked by Jason Stover for the cdf module. |
|
18 |
*/ |
*/ |
19 |
|
|
20 |
#include <config.h> |
#include <config.h> |
21 |
|
#include <math.h> |
22 |
|
#include <gsl/gsl_cdf.h> |
23 |
|
#include <gsl/gsl_sf_gamma.h> |
24 |
#include <gsl/gsl_math.h> |
#include <gsl/gsl_math.h> |
|
#include <gsl/gsl_errno.h> |
|
25 |
|
|
26 |
#include "exp.h" |
#include "beta_inc.c" |
|
#include "log.h" |
|
|
#include "psi.h" |
|
|
#include "gamma.h" |
|
|
#include "error.h" |
|
27 |
|
|
28 |
int |
double |
29 |
gsl_cdf_lnbeta_e(const double x, const double y, gsl_cdf_result * result) |
gsl_cdf_beta_P (double x, double a, double b) |
30 |
{ |
{ |
31 |
/* CHECK_POINTER(result) */ |
double P; |
32 |
|
|
33 |
if(x <= 0.0 || y <= 0.0) { |
if (x <= 0.0 || x > 1.0) |
34 |
DOMAIN_ERROR(result); |
{ |
35 |
} |
return 0.0; |
|
else { |
|
|
const double max = GSL_MAX(x,y); |
|
|
const double min = GSL_MIN(x,y); |
|
|
const double rat = min/max; |
|
|
|
|
|
if(rat < 0.2) { |
|
|
/* min << max, so be careful |
|
|
* with the subtraction |
|
|
*/ |
|
|
double lnpre_val; |
|
|
double lnpre_err; |
|
|
double lnpow_val; |
|
|
double lnpow_err; |
|
|
double t1, t2, t3; |
|
|
gsl_cdf_result lnopr; |
|
|
gsl_cdf_result gsx, gsy, gsxy; |
|
|
gsl_cdf_gammastar_e(x, &gsx); |
|
|
gsl_cdf_gammastar_e(y, &gsy); |
|
|
gsl_cdf_gammastar_e(x+y, &gsxy); |
|
|
gsl_cdf_log_1plusx_e(rat, &lnopr); |
|
|
lnpre_val = log(gsx.val*gsy.val/gsxy.val * M_SQRT2*M_SQRTPI); |
|
|
lnpre_err = gsx.err/gsx.val + gsy.err/gsy.val + gsxy.err/gsxy.val; |
|
|
t1 = min*log(rat); |
|
|
t2 = 0.5*log(min); |
|
|
t3 = (x+y-0.5)*lnopr.val; |
|
|
lnpow_val = t1 - t2 - t3; |
|
|
lnpow_err = GSL_DBL_EPSILON * (fabs(t1) + fabs(t2) + fabs(t3)); |
|
|
lnpow_err += fabs(x+y-0.5) * lnopr.err; |
|
|
result->val = lnpre_val + lnpow_val; |
|
|
result->err = lnpre_err + lnpow_err; |
|
|
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
|
|
return GSL_SUCCESS; |
|
36 |
} |
} |
|
else { |
|
|
gsl_cdf_result lgx, lgy, lgxy; |
|
|
int stat_gx = gsl_cdf_lngamma_e(x, &lgx); |
|
|
int stat_gy = gsl_cdf_lngamma_e(y, &lgy); |
|
|
int stat_gxy = gsl_cdf_lngamma_e(x+y, &lgxy); |
|
|
result->val = lgx.val + lgy.val - lgxy.val; |
|
|
result->err = lgx.err + lgy.err + lgxy.err; |
|
|
result->err += GSL_DBL_EPSILON * (fabs(lgx.val) + fabs(lgy.val) + fabs(lgxy.val)); |
|
|
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
|
|
return GSL_ERROR_SELECT_3(stat_gx, stat_gy, stat_gxy); |
|
|
} |
|
|
} |
|
|
} |
|
37 |
|
|
38 |
|
P = beta_inc_AXPY (1.0, 0.0, a, b, x); |
39 |
|
|
40 |
int |
return P; |
|
gsl_cdf_beta_e(const double x, const double y, gsl_cdf_result * result) |
|
|
{ |
|
|
if(x < 50.0 && y < 50.0) { |
|
|
gsl_cdf_result gx, gy, gxy; |
|
|
gsl_cdf_gamma_e(x, &gx); |
|
|
gsl_cdf_gamma_e(y, &gy); |
|
|
gsl_cdf_gamma_e(x+y, &gxy); |
|
|
result->val = (gx.val*gy.val)/gxy.val; |
|
|
result->err = gx.err * gy.val/gxy.val; |
|
|
result->err += gy.err * gx.val/gxy.val; |
|
|
result->err += (gx.val*gy.val)/(gxy.val*gxy.val) * gxy.err; |
|
|
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
|
|
return GSL_SUCCESS; |
|
|
} |
|
|
else { |
|
|
gsl_cdf_result lb; |
|
|
int stat_lb = gsl_cdf_lnbeta_e(x, y, &lb); |
|
|
if(stat_lb == GSL_SUCCESS) { |
|
|
return gsl_cdf_exp_err_e(lb.val, lb.err, result); |
|
|
} |
|
|
else { |
|
|
result->val = 0.0; |
|
|
result->err = 0.0; |
|
|
return stat_lb; |
|
|
} |
|
|
} |
|
41 |
} |
} |
42 |
|
|
43 |
|
double |
44 |
|
gsl_cdf_beta_Q (double x, double a, double b) |
45 |
|
{ |
46 |
|
double P; |
47 |
|
|
48 |
/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/ |
if (x < 0.0 || x >= 1.0) |
49 |
|
{ |
50 |
#include "eval.h" |
return 0.0; |
51 |
|
} |
52 |
|
|
53 |
double gsl_cdf_lnbeta(const double x, const double y) |
P = beta_inc_AXPY (-1.0, 1.0, a, b, x); |
|
{ |
|
|
EVAL_RESULT(gsl_cdf_lnbeta_e(x, y, &result)); |
|
|
} |
|
54 |
|
|
55 |
double gsl_cdf_beta(const double x, const double y) |
return P; |
|
{ |
|
|
EVAL_RESULT(gsl_cdf_beta_e(x, y, &result)); |
|
56 |
} |
} |