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.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
    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=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