## Revision 601

1
/*
2
* \$Id\$
3
* Author: Unknown
4
*/
5
/**
6
* \file
7
* Description not yet available.
8
*/
9
10

11
/**
12
* Description not yet available.
13
* \param
14
*/
15
static double lnbeta(double a,double b)
16
{
17
return gammln(a)+gammln(b)-gammln(a+b);
18
}
19

20
/**
21
* Description not yet available.
22
* \param
23
*/
24
double inv_cumd_beta_stable(double a,double b,double y,double eps)
25
{
26
double eps1=1.0-eps;
27

28
int icount=0;
29
if (y<0.0 || y>1.0 || a<=0.0 || b<=0.0 )
30
{
31
cerr << "Illegal value in inv_cumd_beta" << endl;
32
33
}
34

35
double mu=a/(a+b);
36
double x=mu;
37
if (x<=eps)
38
x=1.1*eps;
39

40
if(x>=eps1)
41
x=eps1-0.1*eps;
42

43
double s=log(x/(1.0-x));
44
double lower=log(eps/eps1);
45
double upper=-lower; // bracket the minimum
46

47
double bet=exp(lnbeta(a,b)); // normalizing constant
48
double denom=1.0/(incbet(a,b,eps1)-incbet(a,b,eps));
49
double finit=incbet(a,b,x)-incbet(a,b,eps)*denom;
50

51
if (finit>y)
52
{
53
upper=s;
54
}
55
else if (finit<y)
56
{
57
lower=s;
58
}
59

60
double d=0.0;
61
do
62
{
63
double f,dx; // der of x wrt s
64
x=1.0/(1.0+exp(-s));  //transform from s to x
65

66
f=(incbet(a,b,x)-incbet(a,b,eps))*denom;
67

68
dx=exp(-s)/square(1+exp(-s)); // der of x wrt s
69

70
d=y-f;
71

72
if (d<0)   // x is too large
73
{
74
if (s<upper)
75
upper=s;
76
}
77
else if (d>0) // x is too small
78
{
79
if (s>lower)
80
{
81
lower=s;
82
}
83
}
84

85
double xa1=pow(x,a-1);
86
double xb1;
87
xb1=pow(1-x,b-1);
88

89
double fp=(xa1*xb1/bet)*dx*denom; // derivative of cumulative dist fun
90
// wrt x
91
double h=d/fp;
92

93
double stry=s+h;
94
if (h<0.0)
95
{
96
if (stry<lower)
... This diff was truncated because it exceeds the maximum size that can be displayed.

Also available in: Unified diff