1 #include <gnumeric-config.h>
2 #include <gnumeric.h>
3 #include <mathfunc.h>
4 #include <sf-dpq.h>
5 #include <sf-trig.h>
6 #include <sf-gamma.h>
7 #include "extra.h"
8
9 /* ------------------------------------------------------------------------- */
10
11 /* The skew-normal distribution. */
12
13 gnm_float
dsnorm(gnm_float x,gnm_float shape,gnm_float location,gnm_float scale,gboolean give_log)14 dsnorm (gnm_float x, gnm_float shape, gnm_float location, gnm_float scale, gboolean give_log)
15 {
16 if (gnm_isnan (x) || gnm_isnan (shape) || gnm_isnan (location) || gnm_isnan (scale))
17 return gnm_nan;
18
19 if (shape == 0.)
20 return dnorm (x, location, scale, give_log);
21 else if (give_log)
22 return M_LN2gnum + dnorm (x, location, scale, TRUE) + pnorm (shape * x, shape * location, scale, TRUE, TRUE);
23 else
24 return 2 * dnorm (x, location, scale, FALSE) * pnorm (shape * x, location/shape, scale, TRUE, FALSE);
25 }
26
27 gnm_float
psnorm(gnm_float x,gnm_float shape,gnm_float location,gnm_float scale,gboolean lower_tail,gboolean log_p)28 psnorm (gnm_float x, gnm_float shape, gnm_float location, gnm_float scale, gboolean lower_tail, gboolean log_p)
29 {
30 gnm_float result, h;
31
32 if (gnm_isnan (x) || gnm_isnan (shape) ||
33 gnm_isnan (location) || gnm_isnan (scale))
34 return gnm_nan;
35
36 if (shape == 0.)
37 return pnorm (x, location, scale, lower_tail, log_p);
38
39 /* Normalize */
40 h = (x - location) / scale;
41
42 /* Flip to a lower-tail problem. */
43 if (!lower_tail) {
44 h = -h;
45 shape = -shape;
46 lower_tail = !lower_tail;
47 }
48
49 if (gnm_abs (shape) < 10) {
50 gnm_float s = pnorm (h, 0, 1, lower_tail, FALSE);
51 gnm_float t = 2 * gnm_owent (h, shape);
52 result = s - t;
53 } else {
54 /*
55 * Make use of this result for Owen's T:
56 *
57 * T(h,a) = .5N(h) + .5N(ha) - N(h)N(ha) - T(ha,1/a)
58 */
59 gnm_float s = pnorm (h * shape, 0, 1, TRUE, FALSE);
60 gnm_float u = gnm_erf (h / M_SQRT2gnum);
61 gnm_float t = 2 * gnm_owent (h * shape, 1 / shape);
62 result = s * u + t;
63 }
64
65 /*
66 * Negatives can occur due to rounding errors and hopefully for no
67 * other reason.
68 */
69 result= CLAMP (result, 0.0, 1.0);
70
71 if (log_p)
72 return gnm_log (result);
73 else
74 return result;
75 }
76
77 static gnm_float
dsnorm1(gnm_float x,const gnm_float params[],gboolean log_p)78 dsnorm1 (gnm_float x, const gnm_float params[], gboolean log_p)
79 {
80 return dsnorm (x, params[0], params[1], params[2], log_p);
81 }
82
83 static gnm_float
psnorm1(gnm_float x,const gnm_float params[],gboolean lower_tail,gboolean log_p)84 psnorm1 (gnm_float x, const gnm_float params[],
85 gboolean lower_tail, gboolean log_p)
86 {
87 return psnorm (x, params[0], params[1], params[2], lower_tail, log_p);
88 }
89
90 gnm_float
qsnorm(gnm_float p,gnm_float shape,gnm_float location,gnm_float scale,gboolean lower_tail,gboolean log_p)91 qsnorm (gnm_float p, gnm_float shape, gnm_float location, gnm_float scale,
92 gboolean lower_tail, gboolean log_p)
93 {
94 gnm_float x0;
95 gnm_float params[3];
96
97 if (gnm_isnan (p) || gnm_isnan (shape) || gnm_isnan (location) || gnm_isnan (scale))
98 return gnm_nan;
99
100 if (shape == 0.)
101 return qnorm (p, location, scale, lower_tail, log_p);
102
103 if (!log_p && p > 0.9) {
104 /* We're far into the tail. Flip. */
105 p = 1 - p;
106 lower_tail = !lower_tail;
107 }
108
109 x0 = 0.0;
110 params[0] = shape;
111 params[1] = location;
112 params[2] = scale;
113 return pfuncinverter (p, params, lower_tail, log_p,
114 gnm_ninf, gnm_pinf, x0,
115 psnorm1, dsnorm1);
116 }
117
118 /* ------------------------------------------------------------------------- */
119
120 /* The skew-t distribution. */
121
122 gnm_float
dst(gnm_float x,gnm_float n,gnm_float shape,gboolean give_log)123 dst (gnm_float x, gnm_float n, gnm_float shape, gboolean give_log)
124 {
125 if (n <= 0 || gnm_isnan (x) || gnm_isnan (n) || gnm_isnan (shape))
126 return gnm_nan;
127
128 if (shape == 0.)
129 return dt (x, n, give_log);
130 else {
131 gnm_float pdf = dt (x, n, give_log);
132 gnm_float cdf = pt (shape * x * gnm_sqrt ((n + 1)/(x * x + n)),
133 n + 1, TRUE, give_log);
134 return give_log ? (M_LN2gnum + pdf + cdf) : (2. * pdf * cdf);
135 }
136 }
137
138 static gnm_float
gnm_atan_mpihalf(gnm_float x)139 gnm_atan_mpihalf (gnm_float x)
140 {
141 if (x > 0)
142 return gnm_acot (-x);
143 else
144 return gnm_atan (x) - (M_PIgnum / 2);
145 }
146
147 gnm_float
pst(gnm_float x,gnm_float n,gnm_float shape,gboolean lower_tail,gboolean log_p)148 pst (gnm_float x, gnm_float n, gnm_float shape, gboolean lower_tail, gboolean log_p)
149 {
150 gnm_float p;
151
152 if (n <= 0 || gnm_isnan (x) || gnm_isnan (n) || gnm_isnan (shape))
153 return gnm_nan;
154
155 if (shape == 0.)
156 return pt (x, n, lower_tail, log_p);
157
158 if (n > 100) {
159 /* Approximation */
160 return psnorm (x, shape, 0.0, 1.0, lower_tail, log_p);
161 }
162
163 /* Flip to a lower-tail problem. */
164 if (!lower_tail) {
165 x = -x;
166 shape = -shape;
167 lower_tail = !lower_tail;
168 }
169
170 /* Generic fallback. */
171 if (log_p)
172 gnm_log (pst (x, n, shape, TRUE, FALSE));
173
174 if (n != gnm_floor (n)) {
175 /* We would need numerical integration for this. */
176 return gnm_nan;
177 }
178
179 /*
180 * Use recurrence formula from "Recurrent relations for
181 * distributions of a skew-t and a linear combination of order
182 * statistics form a bivariate-t", Computational Statistics
183 * and Data Analysis volume 52, 2009 by Jamallizadeh,
184 * Khosravi, Balakrishnan.
185 *
186 * This brings us down to n==1 or n==2 for which explicit formulas
187 * are available.
188 */
189
190 p = 0;
191 while (n > 2) {
192 double a, lb, c, d, pv, v = n - 1;
193
194 d = v == 2
195 ? M_LN2gnum - gnm_log (M_PIgnum) + gnm_log (3) / 2
196 : (0.5 + M_LN2gnum / 2 - gnm_log (M_PIgnum) / 2 +
197 v / 2 * (gnm_log1p (-1 / (v - 1)) + gnm_log (v + 1)) -
198 0.5 * (gnm_log (v - 2) + gnm_log (v + 1)) +
199 stirlerr (v / 2 - 1) -
200 stirlerr ((v - 1) / 2));
201
202 a = v + 1 + x * x;
203 lb = (d - gnm_log (a) * v / 2);
204 c = pt (gnm_sqrt (v) * shape * x / gnm_sqrt (a), v, TRUE, FALSE);
205 pv = x * gnm_exp (lb) * c;
206 p += pv;
207
208 n -= 2;
209 x *= gnm_sqrt ((v - 1) / (v + 1));
210 }
211
212 g_return_val_if_fail (n == 1 || n == 2, gnm_nan);
213 if (n == 1) {
214 gnm_float p1;
215
216 p1 = (gnm_atan (x) + gnm_acos (shape / gnm_sqrt ((1 + shape * shape) * (1 + x * x)))) / M_PIgnum;
217 p += p1;
218 } else if (n == 2) {
219 gnm_float p2, f;
220
221 f = x / gnm_sqrt (2 + x * x);
222
223 p2 = (gnm_atan_mpihalf (shape) + f * gnm_atan_mpihalf (-shape * f)) / -M_PIgnum;
224
225 p += p2;
226 } else {
227 return gnm_nan;
228 }
229
230 /*
231 * Negatives can occur due to rounding errors and hopefully for no
232 * other reason.
233 */
234 p = CLAMP (p, 0.0, 1.0);
235
236 return p;
237 }
238
239
240 static gnm_float
dst1(gnm_float x,const gnm_float params[],gboolean log_p)241 dst1 (gnm_float x, const gnm_float params[], gboolean log_p)
242 {
243 return dst (x, params[0], params[1], log_p);
244 }
245
246 static gnm_float
pst1(gnm_float x,const gnm_float params[],gboolean lower_tail,gboolean log_p)247 pst1 (gnm_float x, const gnm_float params[],
248 gboolean lower_tail, gboolean log_p)
249 {
250 return pst (x, params[0], params[1], lower_tail, log_p);
251 }
252
253 gnm_float
qst(gnm_float p,gnm_float n,gnm_float shape,gboolean lower_tail,gboolean log_p)254 qst (gnm_float p, gnm_float n, gnm_float shape,
255 gboolean lower_tail, gboolean log_p)
256 {
257 gnm_float x0;
258 gnm_float params[2];
259
260 if (n <= 0 || gnm_isnan (p) || gnm_isnan (n) || gnm_isnan (shape))
261 return gnm_nan;
262
263 if (shape == 0.)
264 return qt (p, n, lower_tail, log_p);
265
266 if (!log_p && p > 0.9) {
267 /* We're far into the tail. Flip. */
268 p = 1 - p;
269 lower_tail = !lower_tail;
270 }
271
272 x0 = 0.0;
273 params[0] = n;
274 params[1] = shape;
275 return pfuncinverter (p, params, lower_tail, log_p,
276 gnm_ninf, gnm_pinf, x0,
277 pst1, dst1);
278 }
279
280 /* ------------------------------------------------------------------------- */
281
282 gnm_float
dgumbel(gnm_float x,gnm_float mu,gnm_float beta,gboolean give_log)283 dgumbel (gnm_float x, gnm_float mu, gnm_float beta, gboolean give_log)
284 {
285 gnm_float z, lp;
286
287 if (!(beta > 0) || gnm_isnan (mu) || gnm_isnan (beta) || gnm_isnan (x))
288 return gnm_nan;
289
290 z = (x - mu) / beta;
291 lp = -(z + gnm_exp (-z));
292 return give_log ? lp - gnm_log (beta) : gnm_exp (lp) / beta;
293 }
294
295 gnm_float
pgumbel(gnm_float x,gnm_float mu,gnm_float beta,gboolean lower_tail,gboolean log_p)296 pgumbel (gnm_float x, gnm_float mu, gnm_float beta, gboolean lower_tail, gboolean log_p)
297 {
298 gnm_float z, lp;
299
300 if (!(beta > 0) || gnm_isnan (mu) || gnm_isnan (beta) || gnm_isnan (x))
301 return gnm_nan;
302
303 z = (x - mu) / beta;
304 lp = -gnm_exp (-z);
305 if (lower_tail)
306 return log_p ? lp : gnm_exp (lp);
307 else
308 return log_p ? swap_log_tail (lp) : 0 - gnm_expm1 (lp);
309 }
310
311 gnm_float
qgumbel(gnm_float p,gnm_float mu,gnm_float beta,gboolean lower_tail,gboolean log_p)312 qgumbel (gnm_float p, gnm_float mu, gnm_float beta, gboolean lower_tail, gboolean log_p)
313 {
314 if (!(beta > 0) ||
315 gnm_isnan (mu) || gnm_isnan (beta) || gnm_isnan (p) ||
316 (log_p ? p > 0 : (p < 0 || p > 1)))
317 return gnm_nan;
318
319 if (log_p) {
320 if (!lower_tail)
321 p = swap_log_tail (p);
322 } else {
323 if (lower_tail)
324 p = gnm_log (p);
325 else
326 p = gnm_log1p (-p);
327 }
328
329 /* We're now in the log_p, lower_tail case. */
330 return mu - beta * gnm_log (-p);
331 }
332