44 |
|
|
45 |
|
|
46 |
double |
double |
47 |
gsl_cdf_tdist_Pinv (const double p, const double nu) |
gsl_cdf_tdist_Pinv (const double P, const double nu) |
48 |
{ |
{ |
49 |
double x, ptail; |
double x, ptail; |
50 |
|
|
51 |
if (p == 1.0) |
if (P == 1.0) |
52 |
{ |
{ |
53 |
return GSL_POSINF; |
return GSL_POSINF; |
54 |
} |
} |
55 |
else if (p == 0.0) |
else if (P == 0.0) |
56 |
{ |
{ |
57 |
return GSL_NEGINF; |
return GSL_NEGINF; |
58 |
} |
} |
59 |
|
|
60 |
if (nu == 1.0) |
if (nu == 1.0) |
61 |
{ |
{ |
62 |
x = tan (M_PI * (p - 0.5)); |
x = tan (M_PI * (P - 0.5)); |
63 |
} |
} |
64 |
else if (nu == 2.0) |
else if (nu == 2.0) |
65 |
{ |
{ |
66 |
double a = 2 * p - 1; |
double a = 2 * P - 1; |
67 |
x = a / sqrt (2 * (1 - a * a)); |
x = a / sqrt (2 * (1 - a * a)); |
68 |
} |
} |
69 |
|
|
70 |
ptail = (p < 0.5) ? p : 1 - p; |
ptail = (P < 0.5) ? P : 1 - P; |
71 |
|
|
72 |
if (sqrt (M_PI * nu / 2) * ptail > pow (0.05, nu / 2)) |
if (sqrt (M_PI * nu / 2) * ptail > pow (0.05, nu / 2)) |
73 |
{ |
{ |
74 |
double xg = gsl_cdf_ugaussian_Pinv (p); |
double xg = gsl_cdf_ugaussian_Pinv (P); |
75 |
x = inv_cornish_fisher (xg, nu); |
x = inv_cornish_fisher (xg, nu); |
76 |
} |
} |
77 |
else |
else |
80 |
|
|
81 |
double beta = gsl_sf_beta (0.5, nu / 2); |
double beta = gsl_sf_beta (0.5, nu / 2); |
82 |
|
|
83 |
if (p < 0.5) |
if (P < 0.5) |
84 |
{ |
{ |
85 |
x = -sqrt (nu) * pow (beta * nu * p, -1.0 / nu); |
x = -sqrt (nu) * pow (beta * nu * P, -1.0 / nu); |
86 |
} |
} |
87 |
else |
else |
88 |
{ |
{ |
89 |
x = sqrt (nu) * pow (beta * nu * (1 - p), -1.0 / nu); |
x = sqrt (nu) * pow (beta * nu * (1 - P), -1.0 / nu); |
90 |
} |
} |
91 |
|
|
92 |
/* Correct nu -> nu/(1+nu/x^2) in the leading term to account |
/* Correct nu -> nu/(1+nu/x^2) in the leading term to account |
98 |
} |
} |
99 |
|
|
100 |
{ |
{ |
101 |
double dp, phi; |
double dP, phi; |
102 |
|
|
103 |
start: |
start: |
104 |
dp = p - gsl_cdf_tdist_P (x, nu); |
dP = P - gsl_cdf_tdist_P (x, nu); |
105 |
phi = gsl_ran_tdist_pdf (x, nu); |
phi = gsl_ran_tdist_pdf (x, nu); |
106 |
|
|
107 |
if (dp == 0.0) |
if (dP == 0.0) |
108 |
goto end; |
goto end; |
109 |
|
|
110 |
{ |
{ |
111 |
double lambda = dp / phi; |
double lambda = dP / phi; |
112 |
double step0 = lambda; |
double step0 = lambda; |
113 |
double step1 = ((nu + 1) * x / (x * x + nu)) * (lambda * lambda / 4.0); |
double step1 = ((nu + 1) * x / (x * x + nu)) * (lambda * lambda / 4.0); |
114 |
|
|
119 |
step += step1; |
step += step1; |
120 |
} |
} |
121 |
|
|
122 |
if (p > 0.5 && x + step < 0) |
if (P > 0.5 && x + step < 0) |
123 |
x /= 2; |
x /= 2; |
124 |
else if (p < 0.5 && x + step > 0) |
else if (P < 0.5 && x + step > 0) |
125 |
x /= 2; |
x /= 2; |
126 |
else |
else |
127 |
x += step; |
x += step; |
137 |
} |
} |
138 |
|
|
139 |
double |
double |
140 |
gsl_cdf_tdist_Qinv (const double p, const double nu) |
gsl_cdf_tdist_Qinv (const double Q, const double nu) |
141 |
{ |
{ |
142 |
double x, ptail; |
double x, qtail; |
143 |
|
|
144 |
if (p == 0.0) |
if (Q == 0.0) |
145 |
{ |
{ |
146 |
return GSL_POSINF; |
return GSL_POSINF; |
147 |
} |
} |
148 |
else if (p == 1.0) |
else if (Q == 1.0) |
149 |
{ |
{ |
150 |
return GSL_NEGINF; |
return GSL_NEGINF; |
151 |
} |
} |
152 |
|
|
153 |
if (nu == 1.0) |
if (nu == 1.0) |
154 |
{ |
{ |
155 |
x = tan (M_PI * (0.5 - p)); |
x = tan (M_PI * (0.5 - Q)); |
156 |
} |
} |
157 |
else if (nu == 2.0) |
else if (nu == 2.0) |
158 |
{ |
{ |
159 |
double a = 2 * (1 - p) - 1; |
double a = 2 * (1 - Q) - 1; |
160 |
x = a / sqrt (2 * (1 - a * a)); |
x = a / sqrt (2 * (1 - a * a)); |
161 |
} |
} |
162 |
|
|
163 |
ptail = (p < 0.5) ? p : 1 - p; |
qtail = (Q < 0.5) ? Q : 1 - Q; |
164 |
|
|
165 |
if (sqrt (M_PI * nu / 2) * ptail > pow (0.05, nu / 2)) |
if (sqrt (M_PI * nu / 2) * qtail > pow (0.05, nu / 2)) |
166 |
{ |
{ |
167 |
double xg = gsl_cdf_ugaussian_Qinv (p); |
double xg = gsl_cdf_ugaussian_Qinv (Q); |
168 |
x = inv_cornish_fisher (xg, nu); |
x = inv_cornish_fisher (xg, nu); |
169 |
} |
} |
170 |
else |
else |
173 |
|
|
174 |
double beta = gsl_sf_beta (0.5, nu / 2); |
double beta = gsl_sf_beta (0.5, nu / 2); |
175 |
|
|
176 |
if (p < 0.5) |
if (Q < 0.5) |
177 |
{ |
{ |
178 |
x = sqrt (nu) * pow (beta * nu * p, -1.0 / nu); |
x = sqrt (nu) * pow (beta * nu * Q, -1.0 / nu); |
179 |
} |
} |
180 |
else |
else |
181 |
{ |
{ |
182 |
x = -sqrt (nu) * pow (beta * nu * (1 - p), -1.0 / nu); |
x = -sqrt (nu) * pow (beta * nu * (1 - Q), -1.0 / nu); |
183 |
} |
} |
184 |
|
|
185 |
/* Correct nu -> nu/(1+nu/x^2) in the leading term to account |
/* Correct nu -> nu/(1+nu/x^2) in the leading term to account |
191 |
} |
} |
192 |
|
|
193 |
{ |
{ |
194 |
double dp, phi; |
double dQ, phi; |
195 |
|
|
196 |
start: |
start: |
197 |
dp = -(p - gsl_cdf_tdist_Q (x, nu)); |
dQ = Q - gsl_cdf_tdist_Q (x, nu); |
198 |
phi = gsl_ran_tdist_pdf (x, nu); |
phi = gsl_ran_tdist_pdf (x, nu); |
199 |
|
|
200 |
if (dp == 0.0) |
if (dQ == 0.0) |
201 |
goto end; |
goto end; |
202 |
|
|
203 |
{ |
{ |
204 |
double lambda = dp / phi; |
double lambda = - dQ / phi; |
205 |
double step0 = lambda; |
double step0 = lambda; |
206 |
double step1 = ((nu + 1) * x / (x * x + nu)) * (lambda * lambda / 4.0); |
double step1 = ((nu + 1) * x / (x * x + nu)) * (lambda * lambda / 4.0); |
207 |
|
|
212 |
step += step1; |
step += step1; |
213 |
} |
} |
214 |
|
|
215 |
if (p < 0.5 && x + step < 0) |
if (Q < 0.5 && x + step < 0) |
216 |
x /= 2; |
x /= 2; |
217 |
else if (p > 0.5 && x + step > 0) |
else if (Q > 0.5 && x + step > 0) |
218 |
x /= 2; |
x /= 2; |
219 |
else |
else |
220 |
x += step; |
x += step; |