plugins/fn-r/extra.c

   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
  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
  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
  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
  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
  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
 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
 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
 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
 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
 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
 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
 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
 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
 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 }