htswanalysis/MACS/lib/gsl/gsl-1.11/randist/exppow.c

   1 /* randist/exppow.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2006, 2007 James Theiler, Brian Gough
   4  * Copyright (C) 2006 Giulio Bottazzi
   5  *
   6  * This program is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU General Public License as published by
   8  * the Free Software Foundation; either version 3 of the License, or (at
   9  * your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful, but
  12  * WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14  * General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU General Public License
  17  * along with this program; if not, write to the Free Software
  18  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  19  */
  20
  21 #include <config.h>
  22 #include <math.h>
  23 #include <gsl/gsl_math.h>
  24 #include <gsl/gsl_sf_gamma.h>
  25 #include <gsl/gsl_rng.h>
  26 #include <gsl/gsl_randist.h>
  27
  28 /* The exponential power probability distribution is
  29
  30    p(x) dx = (1/(2 a Gamma(1+1/b))) * exp(-|x/a|^b) dx
  31
  32    for -infty < x < infty. For b = 1 it reduces to the Laplace
  33    distribution.
  34
  35    The exponential power distribution is related to the gamma
  36    distribution by E = a * pow(G(1/b),1/b), where E is an exponential
  37    power variate and G is a gamma variate.
  38
  39    We use this relation for b < 1. For b >=1 we use rejection methods
  40    based on the laplace and gaussian distributions which should be
  41    faster.  For b>4 we revert to the gamma method.
  42
  43    See P. R. Tadikamalla, "Random Sampling from the Exponential Power
  44    Distribution", Journal of the American Statistical Association,
  45    September 1980, Volume 75, Number 371, pages 683-686.
  46
  47 */
  48
  49 double
  50 gsl_ran_exppow (const gsl_rng * r, const double a, const double b)
  51 {
  52   if (b < 1 || b > 4)
  53     {
  54       double u = gsl_rng_uniform (r);
  55       double v = gsl_ran_gamma (r, 1 / b, 1.0);
  56       double z = a * pow (v, 1 / b);
  57
  58       if (u > 0.5)
  59         {
  60           return z;
  61         }
  62       else
  63         {
  64           return -z;
  65         }
  66     }
  67   else if (b == 1)
  68     {
  69       /* Laplace distribution */
  70       return gsl_ran_laplace (r, a);
  71     }
  72   else if (b < 2)
  73     {
  74       /* Use laplace distribution for rejection method, from Tadikamalla */
  75
  76       double x, h, u;
  77
  78       double B = pow (1 / b, 1 / b);
  79
  80       do
  81         {
  82           x = gsl_ran_laplace (r, B);
  83           u = gsl_rng_uniform_pos (r);
  84           h = -pow (fabs (x), b) + fabs (x) / B - 1 + (1 / b);
  85         }
  86       while (log (u) > h);
  87
  88       return a * x;
  89     }
  90   else if (b == 2)
  91     {
  92       /* Gaussian distribution */
  93       return gsl_ran_gaussian (r, a / sqrt (2.0));
  94     }
  95   else
  96     {
  97       /* Use gaussian for rejection method, from Tadikamalla */
  98
  99       double x, h, u;
 100
 101       double B = pow (1 / b, 1 / b);
 102
 103       do
 104         {
 105           x = gsl_ran_gaussian (r, B);
 106           u = gsl_rng_uniform_pos (r);
 107           h = -pow (fabs (x), b) + (x * x) / (2 * B * B) + (1 / b) - 0.5;
 108         }
 109       while (log (u) > h);
 110
 111       return a * x;
 112     }
 113 }
 114
 115 double
 116 gsl_ran_exppow_pdf (const double x, const double a, const double b)
 117 {
 118   double p;
 119   double lngamma = gsl_sf_lngamma (1 + 1 / b);
 120   p = (1 / (2 * a)) * exp (-pow (fabs (x / a), b) - lngamma);
 121   return p;
 122 }