src/sf-dpq.c

   1 #include <gnumeric-config.h>
   2 #include <sf-dpq.h>
   3 #include <mathfunc.h>
   4
   5 #define give_log log_p
   6 #define R_D__0  (log_p ? gnm_ninf : 0.0)
   7 #define R_D__1  (log_p ? 0.0 : 1.0)
   8 #define R_DT_0  (lower_tail ? R_D__0 : R_D__1)
   9 #define R_DT_1  (lower_tail ? R_D__1 : R_D__0)
  10 #define M_1_SQRT_2PI    GNM_const(0.398942280401432677939946059934)  /* 1/sqrt(2pi) */
  11 #define M_SQRT_2PI GNM_const(2.506628274631000502415765284811045253006986740609938316629923)
  12
  13 /* ------------------------------------------------------------------------- */
  14
  15 #undef DEBUG_pfuncinverter
  16
  17 /**
  18  * pfuncinverter:
  19  * @p:
  20  * @shape:
  21  * @lower_tail:
  22  * @log_p:
  23  * @xlow:
  24  * @xhigh:
  25  * @x0:
  26  * @pfunc: (scope call):
  27  * @dpfunc_dx: (scope call):
  28  *
  29  * Returns:
  30  **/
  31 gnm_float
  32 pfuncinverter (gnm_float p, const gnm_float shape[],
  33                gboolean lower_tail, gboolean log_p,
  34                gnm_float xlow, gnm_float xhigh, gnm_float x0,
  35                GnmPFunc pfunc, GnmDPFunc dpfunc_dx)
  36 {
  37         gboolean have_xlow = gnm_finite (xlow);
  38         gboolean have_xhigh = gnm_finite (xhigh);
  39         gnm_float exlow, exhigh;
  40         gnm_float x = 0, e = 0, px;
  41         int i;
  42
  43         g_return_val_if_fail (pfunc != NULL, gnm_nan);
  44
  45         if (log_p ? (p > 0) : (p < 0 || p > 1))
  46                 return gnm_nan;
  47
  48         if (p == R_DT_0) return xlow;
  49         if (p == R_DT_1) return xhigh;
  50
  51         exlow = R_DT_0 - p;
  52         exhigh = R_DT_1 - p;
  53         if (!lower_tail) {
  54                 exlow = -exlow;
  55                 exhigh = -exhigh;
  56         }
  57
  58 #ifdef DEBUG_pfuncinverter
  59         g_printerr ("p=%.15g\n", p);
  60 #endif
  61
  62         for (i = 0; i < 100; i++) {
  63                 if (i == 0) {
  64                         if (x0 > xlow && x0 < xhigh)
  65                                 /* Use supplied guess.  */
  66                                 x = x0;
  67                         else if (have_xlow && x0 <= xlow)
  68                                 x = xlow + have_xhigh ? (xhigh - xlow) / 100 : 1;
  69                         else if (have_xhigh && x0 >= xhigh)
  70                                 x = xhigh - have_xlow ? (xhigh - xlow) / 100 : 1;
  71                         else
  72                                 x = 0;  /* Whatever */
  73                 } else if (i == 1) {
  74                         /*
  75                          * Under the assumption that the initial guess was
  76                          * good, pick a nearby point that is hopefully on
  77                          * the other side.  If we already have both sides,
  78                          * just bisect.
  79                          */
  80                         if (have_xlow && have_xhigh)
  81                                 x = (xlow + xhigh) / 2;
  82                         else if (have_xlow)
  83                                 x = xlow * 1.1;
  84                         else
  85                                 x = xhigh / 1.1;
  86                 } else if (have_xlow && have_xhigh) {
  87                         switch (i % 8) {
  88                         case 0:
  89                                 x = xhigh - (xhigh - xlow) *
  90                                         (exhigh / (exhigh - exlow));
  91                                 break;
  92                         case 4:
  93                                 /* Half-way in log-space.  */
  94                                 if (xlow >= 0 && xhigh >= 0)
  95                                         x = gnm_sqrt (MAX (GNM_MIN, xlow)) * gnm_sqrt (xhigh);
  96                                 else if (xlow <= 0 && xhigh <= 0)
  97                                         x = -gnm_sqrt (-xlow) * gnm_sqrt (MAX (GNM_MIN, -xhigh));
  98                                 else
  99                                         x = 0;
 100                                 break;
 101                         case 2:
 102                                 x = (xhigh + 1000 * xlow) / 1001;
 103                                 break;
 104                         case 6:
 105                                 x = (1000 * xhigh + xlow) / 1001;
 106                                 break;
 107                         default:
 108                                 x = (xhigh + xlow) / 2;
 109                         }
 110                 } else if (have_xlow) {
 111                         /* Agressively seek right in search of xhigh.  */
 112                         x = (xlow < 1) ? 1 : (2 * i) * xlow;
 113                 } else {
 114                         /* Agressively seek left in search of xlow.  */
 115                         x = (xhigh > -1) ? -1 : (2 * i) * xhigh;
 116                 }
 117
 118         newton_retry:
 119                 if ((have_xlow && x <= xlow) || (have_xhigh && x >= xhigh))
 120                         continue;
 121
 122                 px = pfunc (x, shape, lower_tail, log_p);
 123                 e = px - p;
 124                 if (!lower_tail) e = -e;
 125
 126 #ifdef DEBUG_pfuncinverter
 127                 g_printerr ("%3d:  x=%.15g  e=%.15g  l=%.15g  h=%.15g\n",
 128                             i, x, e, xlow, xhigh);
 129 #endif
 130
 131                 if (e == 0)
 132                         goto done;
 133                 else if (e > 0) {
 134                         xhigh = x;
 135                         exhigh = e;
 136                         have_xhigh = TRUE;
 137                 } else if (e < 0) {
 138                         xlow = x;
 139                         exlow = e;
 140                         have_xlow = TRUE;
 141                 } else {
 142                         /* We got a NaN.  */
 143                 }
 144
 145                 if (have_xlow && have_xhigh) {
 146                         gnm_float prec = (xhigh - xlow) /
 147                                 (gnm_abs (xlow) + gnm_abs (xhigh));
 148                         if (prec < GNM_EPSILON * 4) {
 149                                 x = (xhigh + xlow) / 2;
 150                                 e = pfunc (x, shape, lower_tail, log_p) - p;
 151                                 if (!lower_tail) e = -e;
 152                                 goto done;
 153                         }
 154
 155                         if (dpfunc_dx && i % 3 < 2 && (i == 0 || prec < 0.05)) {
 156                                 gnm_float d = dpfunc_dx (x, shape, log_p);
 157                                 if (log_p) d = gnm_exp (d - px);
 158 #ifdef DEBUG_pfuncinverter
 159                                 g_printerr ("Newton: d=%-.14g\n", d);
 160 #endif
 161                                 if (d) {
 162                                         /*
 163                                          * Deliberately overshoot a bit to help
 164                                          * with getting good points on both
 165                                          * sides of the root.
 166                                          */
 167                                         x = x - e / d * 1.000001;
 168                                         if (x > xlow && x < xhigh) {
 169 #ifdef DEBUG_pfuncinverter
 170                                                 g_printerr ("Newton ok\n");
 171 #endif
 172                                                 i++;
 173                                                 goto newton_retry;
 174                                         }
 175                                 } else {
 176 #ifdef DEBUG_pfuncinverter
 177                                                 g_printerr ("Newton d=0\n");
 178 #endif
 179                                 }
 180                         }
 181                 }
 182         }
 183
 184 #ifdef DEBUG_pfuncinverter
 185         g_printerr ("Failed to converge\n");
 186 #endif
 187
 188  done:
 189         /* Make sure to keep a lucky near-hit.  */
 190
 191         if (have_xhigh && gnm_abs (e) > exhigh)
 192                 e = exhigh, x = xhigh;
 193         if (have_xlow && gnm_abs (e) > -exlow)
 194                 e = exlow, x = xlow;
 195
 196 #ifdef DEBUG_pfuncinverter
 197         g_printerr ("--> %.15g\n\n", x);
 198 #endif
 199         return x;
 200 }
 201
 202 /**
 203  * discpfuncinverter:
 204  * @p:
 205  * @shape:
 206  * @lower_tail:
 207  * @log_p:
 208  * @xlow:
 209  * @xhigh:
 210  * @x0:
 211  * @pfunc: (scope call):
 212  *
 213  * Discrete pfuncs only.  (Specifically: only integer x are allowed).
 214  * Returns:
 215  */
 216 gnm_float
 217 discpfuncinverter (gnm_float p, const gnm_float shape[],
 218                    gboolean lower_tail, gboolean log_p,
 219                    gnm_float xlow, gnm_float xhigh, gnm_float x0,
 220                    GnmPFunc pfunc)
 221 {
 222         gboolean have_xlow = gnm_finite (xlow);
 223         gboolean have_xhigh = gnm_finite (xhigh);
 224         gboolean check_left = TRUE;
 225         gnm_float step;
 226         int i;
 227
 228         if (log_p ? (p > 0) : (p < 0 || p > 1))
 229                 return gnm_nan;
 230
 231         if (p == R_DT_0) return xlow;
 232         if (p == R_DT_1) return xhigh;
 233
 234         if (gnm_finite (x0) && x0 >= xlow && x0 <= xhigh)
 235                 ; /* Nothing -- guess is good.  */
 236         else if (have_xlow && have_xhigh)
 237                 x0 = (xlow + xhigh) / 2;
 238         else if (have_xhigh)
 239                 x0 = xhigh;
 240         else if (have_xlow)
 241                 x0 = xlow;
 242         else
 243                 x0 = 0;
 244         x0 = gnm_floor (x0 + 0.5);
 245         step = 1 + gnm_floor (gnm_abs (x0) * GNM_EPSILON);
 246
 247 #if 0
 248         g_printerr ("step=%.20g\n", step);
 249 #endif
 250
 251         for (i = 1; 1; i++) {
 252                 gnm_float ex0 = pfunc (x0, shape, lower_tail, log_p) - p;
 253 #if 0
 254                 g_printerr ("x=%.20g  e=%.20g\n", x0, ex0);
 255 #endif
 256                 if (!lower_tail) ex0 = -ex0;
 257
 258                 if (ex0 == 0)
 259                         return x0;
 260                 else if (ex0 < 0)
 261                         xlow = x0, have_xlow = TRUE, check_left = FALSE;
 262                 else if (ex0 > 0)
 263                         xhigh = x0, have_xhigh = TRUE, step = -gnm_abs (step);
 264
 265                 if (i > 1 && have_xlow && have_xhigh) {
 266                         gnm_float xmid = gnm_floor ((xlow + xhigh) / 2);
 267                         if (xmid - xlow < 0.5 ||
 268                             xmid - xlow < gnm_abs (xlow) * GNM_EPSILON) {
 269                                 if (check_left) {
 270                                         /*
 271                                          * The lower edge of the support might
 272                                          * have a probability higher than what
 273                                          * we are looking for.
 274                                          */
 275                                         gnm_float e = pfunc (xlow, shape, lower_tail, log_p) - p;
 276                                         if (!lower_tail) e = -e;
 277                                         if (e >= 0)
 278                                                 return xhigh = xlow;
 279                                 }
 280                                 return xhigh;
 281                         }
 282                         x0 = xmid;
 283                 } else {
 284                         gnm_float x1 = x0 + step;
 285
 286                         if (x1 == x0) {
 287                                 /* Probably infinite.  */
 288                                 return gnm_nan;
 289                         } else if (x1 >= xlow && x1 <= xhigh) {
 290                                 x0 = x1;
 291                                 step *= 2 * i;
 292                         } else {
 293                                 /* We went off the edge by walking too fast.  */
 294                                 gnm_float newstep = 1 + gnm_floor (gnm_abs (x0) * GNM_EPSILON);
 295                                 step = (step > 0) ? newstep : -newstep;
 296                                 x1 = x0 + step;
 297                                 if (x1 >= xlow && x1 <= xhigh)
 298                                         continue;
 299                                 /*
 300                                  * We don't seem to find a finite x on the
 301                                  * other side of the root.
 302                                  */
 303                                 return (step > 0) ? xhigh : xlow;
 304                         }
 305                 }
 306         }
 307
 308         g_assert_not_reached ();
 309 }
 310
 311 /* ------------------------------------------------------------------------ */
 312
 313 /**
 314  * dnorm:
 315  * @x: observation
 316  * @mu: mean of the distribution
 317  * @sigma: standard deviation of the distribution
 318  * @give_log: if %TRUE, log of the result will be returned instead
 319  *
 320  * Returns: density of the normal distribution.
 321  */
 322 gnm_float
 323 dnorm (gnm_float x, gnm_float mu, gnm_float sigma, gboolean give_log)
 324 {
 325         gnm_float x0;
 326
 327         if (gnm_isnan (x) || gnm_isnan (mu) || gnm_isnan (sigma))
 328                 return x + mu + sigma;
 329         if (sigma < 0)
 330                 return gnm_nan;
 331
 332         /* Center.  */
 333         x = x0 = gnm_abs (x - mu);
 334         x /= sigma;
 335
 336         if (give_log)
 337                 return -(M_LN_SQRT_2PI + 0.5 * x * x + gnm_log (sigma));
 338         else if (x > 3 + 2 * gnm_sqrt (gnm_log (GNM_MAX)))
 339                 /* Far into the tail; x > ~100 for long double  */
 340                 return 0;
 341         else if (x > 4) {
 342                 /*
 343                  * Split x into xh+xl such that:
 344                  * 1) xh*xh is exact
 345                  * 2) 0 <= xl <= 1/65536
 346                  * 3) 0 <= xl*xh < ~100/65536
 347                  */
 348                 gnm_float xh = gnm_floor (x * 65536) / 65536;  /* At most 24 bits */
 349                 gnm_float xl = (x0 - xh * sigma) / sigma;
 350                 return M_1_SQRT_2PI *
 351                         gnm_exp (-0.5 * (xh * xh)) *
 352                         gnm_exp (-xl * (0.5 * xl + xh)) /
 353                         sigma;
 354         } else
 355                 /* Near-center case.  */
 356                 return M_1_SQRT_2PI * expmx2h (x) / sigma;
 357 }
 358
 359 gnm_float
 360 pnorm2 (gnm_float x1, gnm_float x2)
 361 {
 362         if (gnm_isnan (x1) || gnm_isnan (x2))
 363                 return gnm_nan;
 364
 365         if (x1 > x2)
 366                 return 0 - pnorm2 (x2, x1);
 367
 368         /* A bunch of special cases:  */
 369         if (x1 == x2)
 370                 return 0.0;
 371         if (x1 == gnm_ninf)
 372                 return pnorm (x2, 0.0, 1.0, TRUE, FALSE);
 373         if (x2 == gnm_pinf)
 374                 return pnorm (x1, 0.0, 1.0, FALSE, FALSE);
 375         if (x1 == 0)
 376                 return gnm_erf (x2 / M_SQRT2gnum) / 2;
 377         if (x2 == 0)
 378                 return gnm_erf (x1 / -M_SQRT2gnum) / 2;
 379
 380         if (x1 <= 0 && x2 >= 0) {
 381                 /* The interval spans 0.  */
 382                 gnm_float p1 = pnorm2 (0, MIN (-x1, x2));
 383                 gnm_float p2 = pnorm2 (MIN (-x1, x2), MAX (-x1, x2));
 384                 return 2 * p1 + p2;
 385         } else if (x1 < 0) {
 386                 /* Both < 0 -- use symmetry */
 387                 return pnorm2 (-x2, -x1);
 388         } else {
 389                 /* Both >= 0 */
 390                 gnm_float p1C = pnorm (x1, 0.0, 1.0, FALSE, FALSE);
 391                 gnm_float p2C = pnorm (x2, 0.0, 1.0, FALSE, FALSE);
 392                 gnm_float raw = p1C - p2C;
 393                 gnm_float dx, d1, d2, ub, lb;
 394
 395                 if (gnm_abs (p1C - p2C) * 32 > gnm_abs (p1C + p2C))
 396                         return raw;
 397
 398                 /* dnorm is strictly decreasing in this area.  */
 399                 dx = x2 - x1;
 400                 d1 = dnorm (x1, 0.0, 1.0, FALSE);
 401                 d2 = dnorm (x2, 0.0, 1.0, FALSE);
 402                 ub = dx * d1;  /* upper bound */
 403                 lb = dx * d2;  /* lower bound */
 404
 405                 raw = MAX (raw, lb);
 406                 raw = MIN (raw, ub);
 407                 return raw;
 408         }
 409 }
 410
 411 /* ------------------------------------------------------------------------ */
 412
 413 /**
 414  * dlnorm:
 415  * @x: observation
 416  * @logmean: mean of the underlying normal distribution
 417  * @logsd: standard deviation of the underlying normal distribution
 418  * @give_log: if %TRUE, log of the result will be returned instead
 419  *
 420  * Returns: density of the log-normal distribution.
 421  */
 422 gnm_float
 423 dlnorm (gnm_float x, gnm_float logmean, gnm_float logsd, gboolean give_log)
 424 {
 425         void *state;
 426         GnmQuad qx, qlx, qs, qy, qt;
 427         static GnmQuad qsqrt2pi;
 428         gnm_float r;
 429
 430         if (gnm_isnan (x) || gnm_isnan (logmean) || gnm_isnan (logsd))
 431                 return x + logmean + logsd;
 432
 433         if (logsd <= 0)
 434                 return gnm_nan;
 435
 436         if (x <= 0)
 437                 return R_D__0;
 438
 439         state = gnm_quad_start ();
 440         if (qsqrt2pi.h == 0)
 441                 gnm_quad_sqrt (&qsqrt2pi, &gnm_quad_2pi);
 442         gnm_quad_init (&qx, x);
 443         gnm_quad_log (&qlx, &qx);
 444         gnm_quad_init (&qt, logmean);
 445         gnm_quad_sub (&qy, &qlx, &qt);
 446         gnm_quad_init (&qs, logsd);
 447         gnm_quad_div (&qy, &qy, &qs);
 448         gnm_quad_mul (&qy, &qy, &qy);
 449         qy.h *= -0.5; qy.l *= -0.5;
 450         gnm_quad_mul (&qt, &qs, &qx);
 451         gnm_quad_mul (&qt, &qt, &qsqrt2pi);
 452         if (give_log) {
 453                 gnm_quad_log (&qt, &qt);
 454                 gnm_quad_sub (&qy, &qy, &qt);
 455         } else {
 456                 gnm_quad_exp (&qy, NULL, &qy);
 457                 gnm_quad_div (&qy, &qy, &qt);
 458         }
 459         r = gnm_quad_value (&qy);
 460         gnm_quad_end (state);
 461
 462         return r;
 463 }
 464
 465 /* ------------------------------------------------------------------------ */
 466
 467 /**
 468  * plnorm:
 469  * @x: observation
 470  * @logmean: mean of the underlying normal distribution
 471  * @logsd: standard deviation of the underlying normal distribution
 472  * @lower_tail: if %TRUE, the lower tail of the distribution is considered.
 473  * @log_p: if %TRUE, log of the result will be returned instead
 474  *
 475  * Returns: cumulative density of the log-normal distribution.
 476  */
 477 gnm_float
 478 plnorm (gnm_float x, gnm_float logmean, gnm_float logsd, gboolean lower_tail, gboolean log_p)
 479 {
 480         if (gnm_isnan (x) || gnm_isnan (logmean) || gnm_isnan (logsd))
 481                 return x + logmean + logsd;
 482
 483         if (logsd <= 0)
 484                 return gnm_nan;
 485
 486         return (x > 0)
 487                 ? pnorm (gnm_log (x), logmean, logsd, lower_tail, log_p)
 488                 : R_D__0;
 489 }
 490
 491 /* ------------------------------------------------------------------------ */
 492
 493 /**
 494  * qlnorm:
 495  * @p: probability
 496  * @logmean: mean of the underlying normal distribution
 497  * @logsd: standard deviation of the underlying normal distribution
 498  * @lower_tail: if %TRUE, the lower tail of the distribution is considered.
 499  * @log_p: if %TRUE, @p is given as log probability
 500  *
 501  * Returns: the observation with cumulative probability @p for the
 502  * log-normal distribution.
 503  */
 504 gnm_float
 505 qlnorm (gnm_float p, gnm_float logmean, gnm_float logsd, gboolean lower_tail, gboolean log_p)
 506 {
 507         if (gnm_isnan (p) || gnm_isnan (logmean) || gnm_isnan (logsd))
 508                 return p + logmean + logsd;
 509
 510         if (log_p ? (p > 0) : (p < 0 || p > 1))
 511                 return gnm_nan;
 512
 513         return gnm_exp (qnorm (p, logmean, logsd, lower_tail, log_p));
 514 }
 515
 516 /* ------------------------------------------------------------------------ */
 517
 518 static gnm_float
 519 dcauchy1 (gnm_float x, const gnm_float shape[], gboolean give_log)
 520 {
 521         return dcauchy (x, shape[0], shape[1], give_log);
 522 }
 523
 524 static gnm_float
 525 pcauchy1 (gnm_float x, const gnm_float shape[], gboolean lower_tail, gboolean log_p)
 526 {
 527         return pcauchy (x, shape[0], shape[1], lower_tail, log_p);
 528 }
 529
 530 /**
 531  * qcauchy:
 532  * @p: probability
 533  * @location: center of distribution
 534  * @scale: scale parameter of the distribution
 535  * @lower_tail: if %TRUE, the lower tail of the distribution is considered.
 536  * @log_p: if %TRUE, @p is given as log probability
 537  *
 538  * Returns: the observation with cumulative probability @p for the
 539  * Cauchy distribution.
 540  */
 541
 542 gnm_float
 543 qcauchy (gnm_float p, gnm_float location, gnm_float scale,
 544          gboolean lower_tail, gboolean log_p)
 545 {
 546         gnm_float x;
 547
 548         if (gnm_isnan(p) || gnm_isnan(location) || gnm_isnan(scale))
 549                 return p + location + scale;
 550
 551         if (log_p ? (p > 0) : (p < 0 || p > 1))
 552                 return gnm_nan;
 553
 554         if (scale < 0 || !gnm_finite(scale)) return gnm_nan;
 555
 556         if (log_p) {
 557                 if (p > -1)
 558                         /* The "0" here is important for the p=0 case:  */
 559                         lower_tail = !lower_tail, p = 0 - gnm_expm1 (p);
 560                 else
 561                         p = gnm_exp (p);
 562                 log_p = FALSE;
 563         } else {
 564                 if (p > 0.5) {
 565                         p = 1 - p;
 566                         lower_tail = !lower_tail;
 567                 }
 568         }
 569         x = location + (lower_tail ? -scale : scale) * gnm_cotpi (p);
 570
 571         if (location != 0 && gnm_abs (x / location) < 0.25) {
 572                 /* Cancellation has occurred.  */
 573                 gnm_float shape[2];
 574                 shape[0] = location;
 575                 shape[1] = scale;
 576                 x = pfuncinverter (p, shape, lower_tail, log_p,
 577                                    gnm_ninf, gnm_pinf, x,
 578                                    pcauchy1, dcauchy1);
 579
 580         }
 581
 582         return x;
 583 }
 584
 585 /* ------------------------------------------------------------------------ */
 586
 587 static gnm_float
 588 phyper1 (gnm_float x, const gnm_float shape[],
 589          gboolean lower_tail, gboolean log_p)
 590 {
 591         return phyper (x, shape[0], shape[1], shape[2], lower_tail, log_p);
 592 }
 593
 594 gnm_float
 595 qhyper (gnm_float p, gnm_float NR, gnm_float NB, gnm_float n,
 596         gboolean lower_tail, gboolean log_p)
 597 {
 598         gnm_float y, shape[3];
 599         gnm_float N = NR + NB;
 600
 601         if (gnm_isnan (p) || gnm_isnan (N) || gnm_isnan (n))
 602                 return p + N + n;
 603         if(!gnm_finite (p) || !gnm_finite (N) ||
 604            NR < 0 || NB < 0 || n < 0 || n > N)
 605                 return gnm_nan;
 606
 607         shape[0] = NR;
 608         shape[1] = NB;
 609         shape[2] = n;
 610
 611         if (N > 2) {
 612                 gnm_float mu = n * NR / N;
 613                 gnm_float sigma =
 614                         gnm_sqrt (NR * NB * n * (N - n) / (N * N * (N - 1)));
 615                 gnm_float sigma_gamma =
 616                         (N - 2 * NR) * (N - 2 * n) / ((N - 2) * N);
 617
 618                 /* Cornish-Fisher expansion:  */
 619                 gnm_float z = qnorm (p, 0., 1., lower_tail, log_p);
 620                 y = mu + sigma * z + sigma_gamma * (z * z - 1) / 6;
 621         } else
 622                 y = 0;
 623
 624         return discpfuncinverter (p, shape, lower_tail, log_p,
 625                                   MAX (0, n - NB), MIN (n, NR), y,
 626                                   phyper1);
 627 }
 628
 629 /* ------------------------------------------------------------------------ */
 630 /**
 631  * drayleigh:
 632  * @x: observation
 633  * @scale: scale parameter
 634  * @give_log: if %TRUE, log of the result will be returned instead
 635  *
 636  * Returns: density of the Rayleigh distribution.
 637  */
 638
 639 gnm_float
 640 drayleigh (gnm_float x, gnm_float scale, gboolean give_log)
 641 {
 642         // This is tempting, but has lower precision since sqrt(2)
 643         // is inexact.
 644         //
 645         // return dweibull (x, 2, M_SQRT2gnum * scale, give_log);
 646
 647         if (scale <= 0)
 648                 return gnm_nan;
 649         if (x <= 0)
 650                 return R_D__0;
 651         else {
 652                 gnm_float p = dnorm (x, 0, scale, give_log);
 653                 return give_log
 654                         ? p + gnm_log (x / scale) + M_LN_SQRT_2PI
 655                         : p * x / scale * M_SQRT_2PI;
 656         }
 657 }
 658
 659 /* ------------------------------------------------------------------------ */
 660 /**
 661  * prayleigh:
 662  * @x: observation
 663  * @scale: scale parameter
 664  * @lower_tail: if %TRUE, the lower tail of the distribution is considered.
 665  * @log_p: if %TRUE, log of the result will be returned instead
 666  *
 667  * Returns: cumulative density of the Rayleigh distribution.
 668  */
 669
 670 gnm_float
 671 prayleigh (gnm_float x, gnm_float scale, gboolean lower_tail, gboolean log_p)
 672 {
 673         return pweibull (x, 2, M_SQRT2gnum * scale, lower_tail, log_p);
 674 }
 675
 676
 677 /**
 678  * qrayleigh:
 679  * @p: probability
 680  * @scale: scale parameter
 681  * @lower_tail: if %TRUE, the lower tail of the distribution is considered.
 682  * @log_p: if %TRUE, log of the result will be returned instead
 683  *
 684  * Returns: the observation with cumulative probability @p for the
 685  * Rayleigh distribution.
 686  */
 687 gnm_float
 688 qrayleigh (gnm_float p, gnm_float scale, gboolean lower_tail, gboolean log_p)
 689 {
 690         return qweibull (p, 2, M_SQRT2gnum * scale, lower_tail, log_p);
 691 }
 692
 693 /* ------------------------------------------------------------------------ */