Revision 601
trunk/src/linad99/ccumdbetainv_incbet.cpp (revision 601)  

1 
/* 

2 
* $Id$ 

3 
* Author: Unknown 

4 
*/ 

5 
/** 

6 
* \file 

7 
* Description not yet available. 

8 
*/ 

9 
#include <admodel.h> 

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.0eps; 

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 
ad_exit(1); 

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=eps10.1*eps; 

42  
43 
double s=log(x/(1.0x)); 

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=yf; 

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,a1); 

86 
double xb1; 

87 
xb1=pow(1x,b1); 

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) 
Also available in: Unified diff