Revision 1109 trunk/src/df1b2-separable/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.0-eps;
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())/
... This diff was truncated because it exceeds the maximum size that can be displayed.

Also available in: Unified diff