17 |
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA. |
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA. |
18 |
*/ |
*/ |
19 |
|
|
20 |
|
#include <config.h> |
21 |
|
#include <math.h> |
22 |
|
#include <gsl/gsl_cdf.h> |
23 |
|
#include <gsl/gsl_math.h> |
24 |
|
#include <gsl/gsl_randist.h> |
25 |
|
#include <gsl/gsl_sf_gamma.h> |
26 |
|
|
27 |
|
#include <stdio.h> |
28 |
|
|
29 |
|
double |
30 |
|
gsl_cdf_tdist_Pinv (const double p, const double nu) |
31 |
|
{ |
32 |
|
double x; |
33 |
|
|
34 |
|
if (p == 1.0) |
35 |
|
{ |
36 |
|
return GSL_POSINF; |
37 |
|
} |
38 |
|
else if (p == 0.0) |
39 |
|
{ |
40 |
|
return GSL_NEGINF; |
41 |
|
} |
42 |
|
|
43 |
|
if (p < 0.05) |
44 |
|
{ |
45 |
|
|
46 |
|
} |
47 |
|
else if (p > 0.95) |
48 |
|
{ |
49 |
|
|
50 |
|
} |
51 |
|
else |
52 |
|
{ |
53 |
|
double xg = gsl_cdf_ugaussian_Pinv (p); |
54 |
|
double x0 = xg + xg*(xg*xg+1.0)/(4.0*nu) ; |
55 |
|
x = x0; |
56 |
|
} |
57 |
|
|
58 |
|
{ |
59 |
|
double dp, phi; |
60 |
|
|
61 |
|
start: |
62 |
|
dp = p - gsl_cdf_tdist_P (x, nu); |
63 |
|
phi = gsl_ran_tdist_pdf (x, nu); |
64 |
|
|
65 |
|
if (dp == 0.0) |
66 |
|
goto end; |
67 |
|
|
68 |
|
{ |
69 |
|
double step = dp / phi; |
70 |
|
|
71 |
|
printf("x=%.18e step=%.18e\n", x, step); |
72 |
|
x += step; |
73 |
|
|
74 |
|
if (fabs (step) > 1e-10 * fabs(x)) |
75 |
|
goto start; |
76 |
|
} |
77 |
|
} |
78 |
|
|
79 |
|
end: |
80 |
|
|
81 |
|
return x; |
82 |
|
} |
83 |
|
|
84 |
#if 0 |
#if 0 |
85 |
/* |
/* |
86 |
* Invert the T distribution. Uses a method as shown in |
* Invert the T distribution. Uses a method as shown in |
143 |
return p; |
return p; |
144 |
} |
} |
145 |
|
|
|
double |
|
|
gsl_cdf_ut_Pinv (double prob, double nu) |
|
|
{ |
|
146 |
double initial_result; |
double initial_result; |
147 |
double result; |
double result; |
148 |
double d; |
double d; |