Revision 1109 trunk/src/df1b2separable/df1b2invcumdbeta.cpp
df1b2invcumdbeta.cpp (revision 1109)  

15  15 
const df1b2variable& b,const df1b2variable& x); 
16  16 
df1b2variable inv_cumd_beta_stable(const df1b2variable& a, 
17  17 
const df1b2variable& b,const df1b2variable& x,double eps); 
18 


18  
19  19 
df3_three_variable gammln(const df3_three_variable& xx); 
20  20 
df3_three_variable betacf(const df3_three_variable& a, 
21  21 
const df3_three_variable& b,const df3_three_variable& x); 
...  ...  
25  25 
const df3_three_variable& b,const df3_three_variable& x,int maxit); 
26  26 
df3_three_variable betai(const df3_three_variable& a,const df3_three_variable& b, 
27  27 
double x,int maxit); 
28 


28  
29  29 
static double lnbeta(double a,double b) 
30  30 
{ 
31  31 
return gammln(a)+gammln(b)gammln(a+b); 
32  32 
} 
33 


33  
34  34 
static int sgn(double z) 
35  35 
{ 
36  36 
if (z>=0) 
...  ...  
38  38 
else 
39  39 
return 1; 
40  40 
} 
41 


41  
42  42 
double inv_cumd_beta(double _a,double _b,double _y); 
43  43 
double inv_cumd_beta_stable(double _a,double _b,double _y,double eps); 
44 


44  
45  45 
df1b2variable inv_cumd_beta_stable(const df1b2variable& _a, 
46  46 
const df1b2variable& _b,const df1b2variable& _y,double eps) 
47  47 
{ 
48  48 
ADUNCONST(df1b2variable,a); 
49  49 
ADUNCONST(df1b2variable,b); 
50  50 
ADUNCONST(df1b2variable,y); 
51 


51  
52  52 
double eps1=1.0eps; 
53  53 
// this gets the inverse without derivatives 
54  54 
double ca=value(a); 
55  55 
double cb=value(b); 
56  56 
double cx=inv_cumd_beta_stable(ca,cb,value(y),eps); 
57 


57  
58  58 
init_df3_three_variable vx(cx); 
59  59 
init_df3_three_variable va(_a); 
60  60 
init_df3_three_variable vb(_b); 
61 


61  
62  62 
// this gets the derivatives for the function itself 
63 


63  
64  64 
df3_three_variable z=(betai(va,vb,vx,25)betai(va,vb,eps,25))/ 
65  65 
(betai(va,vb,eps1,25)betai(va,vb,eps,25)); 
66 


66  
67  67 
// now solve for the derivatves of the inverse function 
68  68 
double F_x=1.0/(*z.get_u_x()); 
69  69 
double F_y=F_x* *z.get_u_y(); 
70  70 
double F_z=F_x* *z.get_u_z(); 
71 


71  
72  72 
double F_xx=F_x* *z.get_u_xx()/square(*z.get_u_x()); 
73  73  
74  74 
double F_xy=(F_xx * *z.get_u_x() * *z.get_u_y() + F_x * *z.get_u_xy())/ 
...  ...  
76  76  
77  77 
double F_xz=(F_xx * *z.get_u_x() * *z.get_u_z() + F_x * *z.get_u_xz())/ 
78  78 
*z.get_u_x(); 
79 
double F_yy=(F_xx * square(*z.get_u_y())


79 
double F_yy=(F_xx * square(*z.get_u_y()) 

80  80 
+2.0* F_xy* *z.get_u_y() 
81  81 
+F_x * *z.get_u_yy()); 
82  82  
83 
double F_yz=( F_xx * *z.get_u_y() * *z.get_u_z()


83 
double F_yz=( F_xx * *z.get_u_y() * *z.get_u_z() 

84  84 
+ F_x * *z.get_u_yz() 
85  85 
+ F_xy * *z.get_u_z() 
86  86 
+ F_xz * *z.get_u_y()); 
87  87  
88 
double F_zz=(F_xx * square(*z.get_u_z())


88 
double F_zz=(F_xx * square(*z.get_u_z()) 

89  89 
+2.0* F_xz* *z.get_u_z() 
90  90 
+F_x * *z.get_u_zz()); 
91  91  
...  ...  
95  95 
double F_xxy=(F_xxx * square(*z.get_u_x())* *z.get_u_y() 
96  96 
+ 2.0*F_xx* *z.get_u_x()* *z.get_u_xy() 
97  97 
+ F_xx* *z.get_u_xx() * *z.get_u_y() 
98 
+ F_xy * *z.get_u_xx()


98 
+ F_xy * *z.get_u_xx() 

99  99 
+ F_x * *z.get_u_xxy())/ 
Also available in: Unified diff