share/distrib/distrib.mac

   1 /*               COPYRIGHT NOTICE
   2
   3 Copyright (C) 2005-2017 Mario Rodriguez Riotorto
   4
   5 This program is free software; you can redistribute
   6 it and/or modify it under the terms of the
   7 GNU General Public License as published by
   8 the Free Software Foundation; either version 2
   9 of the License, or (at your option) any later version.
  10
  11 This program is distributed in the hope that it
  12 will be useful, but WITHOUT ANY WARRANTY;
  13 without even the implied warranty of MERCHANTABILITY
  14 or FITNESS FOR A PARTICULAR PURPOSE. See the
  15 GNU General Public License for more details at
  16 http://www.gnu.org/copyleft/gpl.html
  17
  18 For comments, suggestions and the like, feel free to contact the author at
  19
  20 To test:
  21 batch("rtest_distrib.mac", test) ;
  22
  23 riotorto AT yahoo DOT com
  24
  25 */
  26
  27
  28 /*             INTRODUCTION
  29
  30 This is a set of Maxima functions for univariate probability distributions,
  31 both continuous and discrete.
  32
  33 Continuous distributions:                   Discrete distributions:
  34    Normal              (*normal)              Binomial             (*binomial)
  35    Student             (*student_t)           Poisson              (*poisson)
  36    Chi^2               (*chi2)                Bernoulli            (*bernoulli)
  37    F                   (*f)                   Geometric            (*geometric)
  38    Exponential         (*exp)                 Discrete uniform     (*discrete_uniform)
  39    Lognormal           (*lognormal)           Hypergeometric       (*hypergeometric)
  40    Gamma               (*gamma)               Negative binomial    (*negative_binomial)
  41    Beta                (*beta)                Finite discrete      (*general_finite_discrete)
  42    Continuous uniform  (*continuous_uniform)
  43    Logistic            (*logistic)
  44    Pareto              (*pareto)
  45    Weibull             (*weibull)
  46    Rayleigh            (*rayleigh)
  47    Laplace             (*laplace)
  48    Cauchy              (*cauchy)
  49    Gumbel              (*gumbel)
  50    Noncentral Chi^2    (*noncentral_chi2)
  51    Noncentral Student  (*noncentral_student_t)
  52
  53 Functions:
  54    Density function            (pdf_*)
  55    Distribution function       (cdf_*)
  56    Quantile                    (quantile_*)
  57    Mean                        (mean_*)
  58    Variance                    (var_*)
  59    Standard deviation          (std_*)
  60    Skewness coefficient        (skewness_*)
  61    Kurtosis coefficient        (kurtosis_*)
  62    Random variate              (random_*)
  63
  64 For example,
  65    pdf_student_t(x,n) is the density function of the Student distribution
  66                    with n degrees of freedom
  67    std_pareto(a,b) is the standard deviation of the Pareto distribution
  68                    with parameters a and b
  69    kurtosis_poisson(m) is the kurtosis coefficient of the Poisson distribution
  70                    with mean m
  71
  72 Note: the Cauchy model has no moments, in this case only the density and
  73       the distribution functions, 'pdf_cauchy' and 'cdf_cauchy', are defined.
  74
  75 For questions, suggestions, bugs and the like, feel free
  76 to contact me at
  77
  78 riotorto @@@ yahoo DOT com
  79 */
  80
  81 put('distrib, 2, 'version) $
  82
  83 /* Sets the random state according to the computer clock time */
  84 set_random_state(make_random_state(true))$
  85
  86
  87 /* Loads numerical routines */
  88 load("numdistrib.lisp")$
  89
  90
  91
  92 /*         NORMAL (OR GAUSSIAN) DISTRIBUTION          */
  93
  94 pdf_normal(x,m,s) :=
  95    if maybe(s > 0) = false
  96      then  error("pdf_normal: standard deviation must be greater than zero")
  97      else  exp(-(x-m)^2/(2*s^2))/(sqrt(2*%pi)*s)$
  98
  99 cdf_normal(x,m,s) :=
 100    if maybe(s > 0) = false
 101      then error("cdf_normal: standard deviation must be greater than zero")
 102      else 1/2+erf((x-m)/(s*sqrt(2)))/2 $
 103
 104 /* R: qnorm(q,m,s) */
 105 quantile_normal(q,m,s) :=
 106    if maybe(s > 0 and q >= 0 and q <= 1) = false
 107      then error("quantile_normal: illegal parameters")
 108      else if sign(q) = 'zero
 109             then 'minf
 110           elseif sign(q-1) = 'zero
 111             then 'inf
 112             else m + sqrt(2)*s*inverse_erf(2*q-1) $
 113
 114 mean_normal(m,s) :=
 115    if maybe(s > 0) = false
 116      then error("mean_normal: standard deviation must be greater than zero")
 117      else m $
 118
 119 var_normal(m,s) :=
 120    if maybe(s > 0) = false
 121      then error("var_normal: standard deviation must be greater than zero")
 122      else s^2 $
 123
 124 std_normal(m,s) :=
 125    if maybe(s > 0) = false
 126      then error("std_normal: standard deviation must be greater than zero")
 127      else s $
 128
 129 skewness_normal(m,s) :=
 130    if maybe(s > 0) = false
 131      then error("skewness_normal: standard deviation must be greater than zero")
 132      else 0 $
 133
 134 kurtosis_normal(m,s) :=
 135    if maybe(s > 0) = false
 136      then error("kurtosis_normal: standard deviation must be greater than zero")
 137      else 0 $
 138
 139 random_normal(m,s,[num]) :=
 140    if maybe(s > 0) = false
 141      then error("random_normal: standard deviation must be greater than zero")
 142      else block([no],
 143             if length(num) = 0 then no: 0 else no: num[1],
 144             if integerp(no) and no >= 0
 145               then m + s * ?rndnormal(no)
 146               else error("random_normal: check sample size")) $
 147
 148
 149
 150 /*         STUDENT DISTRIBUTION          */
 151
 152 pdf_student_t(x,n) :=
 153    if maybe(n > 0) = false
 154      then error("pdf_student: number of degrees must be greater than zero")
 155      else gamma((n+1)/2) * (1+x*x/n)^(-(n+1)/2) / (sqrt(n*%pi) * gamma(n/2)) $
 156
 157 /* R: pt(x,n) */
 158 cdf_student_t(x,n) :=
 159    if maybe(n > 0) = false
 160      then error("cdf_student_t: number of degrees must be greater than zero")
 161      else (1+signum(x))/2 - signum(x) * beta_incomplete_regularized(n/2,1/2,n/(n+x^2)) / 2 $
 162
 163 /* R: qt(q,n) */
 164 quantile_student_t(q,n) :=
 165   if maybe(n > 0 and q >= 0 and q <= 1) = false
 166     then error("quantile_student_t: illegal parameters")
 167     else /* need numerical approximation */
 168          block([fq: float(q), fn: float(n), aux, sgn],
 169                if numberp(fq) and numberp(fn)
 170                  then (if fq = 0.0 then return('minf),
 171                        if fq = 1.0 then return('inf),
 172                        if fq = 0.5 then return(0),
 173                        if fq < 0.5
 174                          then (aux: 2*fq,
 175                                sgn: -1)
 176                          else (aux: 2*(1-fq),
 177                                sgn: 1),
 178                        sgn*sqrt(n*(1 / ?iibeta(aux,float(n/2),0.5)-1)))
 179                  else error("quantile_student: need numeric arguments for approximate procedure")) $
 180
 181 mean_student_t(n) :=
 182    if maybe(n > 0) = false
 183      then error("mean_student_t:: degrees of freedom must be greater than zero")
 184      else 0 $
 185
 186 var_student_t(n) :=
 187    if maybe(n > 2) = false
 188      then error("var_student_t: degrees of freedom must be greater than 2")
 189      else n / (n-2) $
 190
 191 std_student_t(n) :=
 192    if maybe(n > 2) = false
 193      then error("std_student_t: degrees of freedom must be greater than 2")
 194      else sqrt(n / (n-2)) $
 195
 196 skewness_student_t(n) :=
 197    if maybe(n > 3) = false
 198      then error("skewness_student_t: degrees of freedom must be greater than 3")
 199      else 0 $
 200
 201 kurtosis_student_t(n) :=
 202    if maybe(n > 4) = false
 203      then error("kurtosis_student_t: degrees of freedom must be greater than 4")
 204      else 6/(n-4) $
 205
 206 random_student_t(n,[num]) :=
 207    if maybe(n > 0) = false
 208      then error("random_student_t: degrees of freedom must be greater than zero")
 209      else block([no, fn: float(n)],
 210             if not numberp(fn)
 211               then error("random_student_t: need numeric arguments for approximate procedure"),
 212             if length(num) = 0 then no: 0 else no: num[1],
 213             if integerp(no) and no >= 0
 214               then ?rndstudent(fn,no)
 215               else error("random_student_t: check sample size") ) $
 216
 217
 218
 219 /*     NONCENTRAL STUDENT DISTRIBUTION      */
 220
 221 /* According to documentation on hgfred, sometimes it might be useful to load */
 222 /* package orthopoly. In R, dt(x,n,ncp)                                       */
 223 pdf_noncentral_student_t(x,n,ncp) :=
 224    if maybe(n > 0) = false
 225      then error("pdf_noncentral_student_t: degrees of freedom must be greater than 0")
 226    elseif sign(ncp) = 'zero
 227      then pdf_student_t(x,n)
 228      else n^(n/2) * factorial(n) * exp(-ncp^2/2) /
 229           (2^n * (n+x^2)^(n/2) * gamma(n/2)) *
 230           (sqrt(2)*ncp*x*hgfred([n/2+1],[3/2],ncp^2*x^2/(2*(n+x^2))) /  ((n+x^2) * gamma((n+1)/2))  +
 231            hgfred([(n+1)/2],[1/2],ncp^2*x^2/(2*(n+x^2))) / (sqrt(n+x^2) * gamma(n/2+1))) $
 232
 233 /* R: pt(x,n,ncp) */
 234 cdf_noncentral_student_t(x,n,ncp) :=
 235    if maybe(n > 0) = false
 236      then error("cdf_noncentral_student_t: degrees of freedom must be greater than 0")
 237    elseif sign(ncp) = 'zero
 238      then cdf_student_t(x,n)
 239      else /* numerical approximation */
 240           block([fx: float(x), fn: float(n), fncp: float(ncp)],
 241             if numberp(fx) and numberp(fn) and numberp(fncp)
 242               then ?cdfnt(fx,fn,fncp)
 243               else error("cdf_noncentral_student_t: need numeric arguments for approximate procedure")) $
 244
 245 /* R: qt(q,n,ncp) */
 246 quantile_noncentral_student_t(q,n,ncp) :=
 247   if maybe(n > 0 and q >= 0 and q <= 1) = false
 248     then error("quantile_noncentral_student_t: illegal parameters")
 249     else /* numerical approximation */
 250          block([fq: float(q), fn: float(n), fncp: float(ncp)],
 251             if numberp(fq) and numberp(fn)
 252               then (if fq = 0.0 then return('minf),
 253                     if fq = 1.0 then return('inf),
 254                     if numberp(fncp)
 255                       then if fncp=0.0
 256                              then return(quantile_student_t(fq,fn))
 257                              else ?qnct(fq, fn, fncp)
 258                       else error("quantile_noncentral_student_t: need numeric arguments for approximate procedure")) ) $
 259
 260 mean_noncentral_student_t(n,ncp) :=
 261    if maybe(n > 1) = false
 262      then error("mean_noncentral_student_t: degrees of freedom must be greater than 1")
 263      else ncp * sqrt(n/2) * gamma((n-1)/2) / gamma(n/2) $
 264
 265 var_noncentral_student_t(n,ncp) :=
 266    if maybe(n > 2) = false
 267      then error("var_noncentral_student_t: degrees of freedom must be greater than 2")
 268      else n*(1+ncp^2)/(n-2) - ncp^2*n/2 * (gamma((n-1)/2) / gamma(n/2))^2 $
 269
 270 std_noncentral_student_t(n,ncp) := sqrt(var_noncentral_student_t(n,ncp)) $
 271
 272 skewness_noncentral_student_t(n,ncp) :=
 273    if maybe(n > 3) = false
 274      then error("skewness_noncentral_student_t: degrees of freedom must be greater than 3")
 275      else block([m: mean_noncentral_student_t(n,ncp), v: var_noncentral_student_t(n,ncp)],
 276                  m / v^(3/2) * (n*(2*n-3+ncp^2)/((n-2)*(n-3)) - 2*v) ) $
 277
 278 kurtosis_noncentral_student_t(n,ncp) :=
 279    if maybe(n > 4) = false
 280      then error("kurtosis_noncentral_student_t: degrees of freedom must be greater than 4")
 281      else block([m: mean_noncentral_student_t(n,ncp), v: var_noncentral_student_t(n,ncp)],
 282                  (n^2*(3+6*ncp^2+ncp^4)/((n-2)*(n-4)) - m^2 * (n*((n+1)*ncp^2+3*(3*n-5))/((n-2)*(n-3)) - 3*v))/v^2 - 3) $
 283
 284 random_noncentral_student_t(n,ncp,[num]) :=
 285    if maybe(n > 0) = false
 286      then error("random_noncentral_student_t: degrees of freedom must be greater than zero")
 287      else block([no, fn: float(n), fncp: float(ncp)],
 288                  if numberp(fn) and numberp(fncp)
 289                    then (if length(num) = 0 then no: 0 else no: num[1],
 290                          if integerp(no) and no >= 0
 291                            then ?rndncstudent(fn,fncp,no)
 292                            else error("random_noncentral_student_t: check sample size"))
 293                    else error("random_noncentral_student_t: need numeric arguments for approximate procedure") ) $
 294
 295
 296
 297 /*         CHI-SQUARE (OR PEARSON'S)  DISTRIBUTION          */
 298 /*          chi2(n) is equivalent to gamma(n/2,2)           */
 299
 300 pdf_chi2(x,n) := pdf_gamma(x,n/2,2)$
 301
 302 cdf_chi2(x,n) := cdf_gamma(x,n/2,2)$
 303
 304 quantile_chi2(q,n) := quantile_gamma(q,n/2,2)$
 305
 306 mean_chi2(n) := mean_gamma(n/2,2)$
 307
 308 var_chi2(n) := var_gamma(n/2,2)$
 309
 310 std_chi2(n) := std_gamma(n/2,2)$
 311
 312 skewness_chi2(n) := skewness_gamma(n/2,2)$
 313
 314 kurtosis_chi2(n) := kurtosis_gamma(n/2,2)$
 315
 316 random_chi2(n,[num]) :=
 317    if maybe(n > 0) = false
 318      then error("random_chi2: degrees of freedom must be greater than zero")
 319      else block([no, fn: float(n)],
 320             if not numberp(fn)
 321               then error("random_chi2: need numeric arguments for approximate procedure"),
 322             if length(num) = 0 then no: 0 else no: num[1],
 323             if integerp(no) and no >= 0
 324               then ?rndchi2(fn,no)
 325               else error("random_chi2: check sample size")) $
 326
 327
 328
 329 /*     NONCENTRAL CHI-SQUARE DISTRIBUTION      */
 330
 331 /* R: dchisq(x,n,ncp) */
 332 pdf_noncentral_chi2(x,n,ncp) :=
 333    if maybe(n > 0) = false
 334      then error("pdf_noncentral_chi2: degrees of freedom must be greater than 0")
 335    elseif sign(ncp) = 'zero
 336      then pdf_chi2(x,n)
 337      else 1/2 * exp(-(x+ncp)/2) * (x/ncp)^(n/4-1/2) * bessel_i(n/2-1, sqrt(x*ncp)) * unit_step(x) $
 338
 339 /* R: pchisq(x,n,ncp) */
 340 cdf_noncentral_chi2(x,n,ncp) :=
 341    if maybe(n > 0) = false
 342      then error("cdf_noncentral_chi2: degrees of freedom must be greater than 0")
 343      elseif float(ncp)=0.0
 344        then cdf_chi2(x,n)
 345        elseif maybe(x >= 0) = false
 346          then 0
 347          else /* need numerical approximation */
 348               block([fx: float(x), fn: float(n), fncp: float(ncp)],
 349                     if numberp(fx) and numberp(fn) and numberp(fncp)
 350                       then ?cdfnchi2(fx, fn, fncp,
 351                                     1e-12,                /* maximum error */
 352                                     1.4210854715202e-14,  /* 8*DBL_EPSILON */
 353                                     1000000.0)            /* number of iterations */
 354                       else error("cdf_noncentral_chi2: need numeric arguments for approximate procedure") ) $
 355
 356 /* R: qchisq(q,n,ncp) */
 357 quantile_noncentral_chi2(q,n,ncp) :=
 358    if maybe(n > 0 and q >= 0 and q <= 1) = false
 359      then error("quantile_noncentral_chi2: degrees of freedom must be greater than 0")
 360      else /* numerical approximation */
 361           block([fq: float(q), fn: float(n), fncp: float(ncp)],
 362             if numberp(fq) and numberp(fn) and numberp(fncp)
 363               then (if fq = 0.0 then return(0),
 364                     if fq = 1.0 then return('inf),
 365                     if fncp=0.0 then return(quantile_chi2(q,n)),
 366                     ?qnchi2(fq, fn, fncp))
 367               else error("quantile_noncentral_chi2: need numeric arguments for approximate procedure")) $
 368
 369 mean_noncentral_chi2(n,ncp) :=
 370    if maybe(n > 0) = false
 371      then error("mean_noncentral_chi2: degrees of freedom must be greater than 0")
 372      else n + ncp $
 373
 374 var_noncentral_chi2(n,ncp) :=
 375    if maybe(n > 0) = false
 376      then error("var_noncentral_chi2: degrees of freedom must be greater than 0")
 377      else 2*(n + 2*ncp) $
 378
 379 std_noncentral_chi2(n,ncp) := sqrt(var_noncentral_chi2(n,ncp)) $
 380
 381 skewness_noncentral_chi2(n,ncp) :=
 382    if maybe(n > 0) = false
 383      then error("skewness_noncentral_chi2: degrees of freedom must be greater than 0")
 384      else 2^(3/2) * (n+3*ncp) /(n+2*ncp)^(3/2) $
 385
 386 kurtosis_noncentral_chi2(n,ncp) :=
 387    if maybe(n > 0) = false
 388      then error("kurtosis_noncentral_chi2: degrees of freedom must be greater than 0")
 389      else 12 * (n+4*ncp) /(n+2*ncp)^2 $
 390
 391 random_noncentral_chi2(n,ncp,[num]) :=
 392    if maybe(n > 0) = false
 393      then error("random_noncentral_chi2: degrees of freedom must be greater than 0")
 394      else block([no, fn: float(n), fncp: float(ncp)],
 395             if not numberp(fn) or not numberp(fncp)
 396               then error("random_noncentral_chi2: need numeric arguments for approximate procedure"),
 397             if length(num) = 0 then no: 0 else no: num[1],
 398             if integerp(no) and no >= 0
 399               then ?rndnchi2(fn,fncp,no)
 400               else error("random_noncentral_chi2: check sample size")) $
 401
 402
 403
 404 /*         F DISTRIBUTION          */
 405
 406 pdf_f(x,m,n) :=
 407    if maybe(n > 0 and m > 0) = false
 408      then error("pdf_f: degrees of freedom must be greater than 0")
 409      else gamma((m+n)/2)*(m/n)^(m/2)*x^(m/2-1)*(1+m*x/n)^(-(m+n)/2) / (gamma(m/2)*gamma(n/2)) * unit_step(x) $
 410
 411 /* R: pf(x,m,n) */
 412 cdf_f(x,m,n) :=
 413    if maybe(n > 0 and m > 0) = false
 414      then error("cdf_f: degrees of freedom must be greater than 0")
 415      else (1 - beta_incomplete_regularized(n/2, m/2, n/(n+m*x))) * unit_step(x) $
 416
 417 /* R: qf(q,m,n) */
 418 quantile_f(q,m,n) :=
 419    if maybe(n > 0 and m > 0 and q >= 0 and q <= 1) = false
 420      then error("quantile_f: check input parameters")
 421    else block([fq: float(q), fn: float(n), fm: float(m)],
 422           if numberp(fq) and numberp(fn) and numberp(fm)
 423             then(if fq = 0.0 then return(0),
 424                  if fq = 1.0 then return('inf),
 425                  fn * (1 / ?iibeta(1-fq, fn/2, fm/2) - 1) / fm)
 426             else error("quantile_f: need numeric arguments for approximate procedure")) $
 427
 428 mean_f(m,n) :=
 429    if maybe(n > 2 and m > 0) = false
 430      then error("mean_f: degrees of freedom must be m>0 and n>2")
 431      else n/(n-2) $
 432
 433 var_f(m,n) :=
 434    if maybe(n > 4 and m > 0) = false
 435      then error("var_f: degrees of freedom must be m>0 and n>4")
 436      else 2*n*n*(m+n-2)/(m*(n-2)*(n-2)*(n-4)) $
 437
 438 std_f(m,n) := sqrt(var_f(m,n)) $
 439
 440 skewness_f(m,n) :=
 441    if maybe(n > 6 and m > 0) = false
 442      then error("skewness_f: degrees of freedom must be m>0 and n>6")
 443      else (2*m+n-2)*sqrt(8*(n-4))/(sqrt(m*(m+n-2))*(n-6)) $
 444
 445 kurtosis_f(m,n) :=
 446    if maybe(n > 8 and m > 0) = false
 447      then error("kurtosis_f: degrees of freedom must be m>0 and n>8")
 448      else 12*((n-2)*(n-2)*(n-4)+m*(m+n-2)*(5*n-22)) / (m*(n-6)*(n-8)*(m+n-2)) $
 449
 450 random_f(m,n,[num]) :=
 451    if maybe(n > 0 and m > 0) = false
 452      then error("pdf_f: degrees of freedom must be greater than 0")
 453      else block([no, fm: float(m), fn: float(n)],
 454             if not numberp(fn) or not numberp(fm)
 455               then error("random_f: need numeric arguments for approximate procedure"),
 456             if length(num) = 0 then no: 0 else no: num[1],
 457             if integerp(no) and no >= 0
 458               then ?rndf(fm,fn,no)
 459               else error("random_f: check sample size")) $
 460
 461
 462
 463 /*         EXPONENTIAL DISTRIBUTION          */
 464 /*   exp(m) equivalent to Weibull(1,1/m)     */
 465
 466 pdf_exp(x,m) := pdf_weibull(x,1,1/m)$
 467
 468 cdf_exp(x,m) := cdf_weibull(x,1,1/m)$
 469
 470 quantile_exp(q,m) := quantile_weibull(q,1,1/m)$
 471
 472 mean_exp(m) := mean_weibull(1,1/m)$
 473
 474 var_exp(m) := var_weibull(1,1/m)$
 475
 476 std_exp(m) := std_weibull(1,1/m)$
 477
 478 skewness_exp(m) := skewness_weibull(1,1/m)$
 479
 480 kurtosis_exp(m) := kurtosis_weibull(1,1/m)$
 481
 482 random_exp(m,[num]) :=
 483    if maybe(m > 0) = false
 484      then error("random_exp: rate must be greater than 0")
 485      else block([no, fm: float(m)],
 486             if not numberp(fm)
 487               then error("random_exp: need numeric arguments for approximate procedure"),
 488             if length(num) = 0 then no: 0 else no: num[1],
 489             if integerp(no) and no >= 0
 490               then ?rndexp(fm,no)
 491               else error("random_exp: check sample size") ) $
 492
 493
 494
 495 /*         LOGNORMAL DISTRIBUTION          */
 496
 497 pdf_lognormal(x,m,s) :=
 498    if maybe(s > 0) = false
 499      then error("pdf_lognormal: parameter s must be greater than 0")
 500      else exp(-(log(x)-m)^2/(2*s^2))/(sqrt(2*%pi)*s*x) * unit_step(x) $
 501
 502 cdf_lognormal(x,m,s) :=
 503    if maybe(s > 0) = false
 504      then error("cdf_lognormal: parameter s must be greater than 0")
 505      else 1/2+erf((log(x)-m)/(s*sqrt(2)))/2 * unit_step(x) $
 506
 507 /* R: qlnorm(p,m,s)*/
 508 quantile_lognormal(q,m,s) :=
 509    if maybe(s > 0 and q >= 0 and q <= 1) = false
 510      then error("quantile_lognormal: check input parameters")
 511    elseif sign(q)='zero
 512      then 0
 513    elseif sign(q-1)='zero
 514      then 'inf
 515      else exp(m + sqrt(2)*s*inverse_erf(2*q-1)) $
 516
 517 mean_lognormal(m,s) :=
 518    if maybe(s > 0) = false
 519      then error("mean_lognormal: parameter s must be greater than 0")
 520      else exp(m+s^2/2) $
 521
 522 var_lognormal(m,s) :=
 523    if maybe(s > 0) = false
 524      then error("var_lognormal: parameter s must be greater than 0")
 525      else exp(2*m+s^2)*(exp(s^2)-1) $
 526
 527 std_lognormal(m,s) := sqrt(var_lognormal(m,s)) $
 528
 529 skewness_lognormal(m,s) :=
 530    if maybe(s > 0) = false
 531      then error("skewness_lognormal: parameter s must be greater than 0")
 532      else (exp(s^2)+2)*sqrt(exp(s^2)-1) $
 533
 534 kurtosis_lognormal(m,s) :=
 535    if maybe(s > 0) = false
 536      then error("kurtosis_lognormal: parameter s must be greater than 0")
 537      else exp(4*s^2)+2*exp(3*s^2)+3*exp(2*s^2)-3 $
 538
 539 random_lognormal(m,s,[num]) :=
 540    if maybe(s > 0) = false
 541      then error("random_lognormal: parameter s must be greater than 0")
 542      else block([no],
 543             if length(num) = 0 then no: 0 else no: num[1],
 544             if integerp(no) and no >= 0
 545               then exp(m + s * ?rndnormal(no))
 546               else error("random_lognormal: check sample size") ) $
 547
 548
 549
 550 /*         GAMMA DISTRIBUTION          */
 551
 552 /* R: dgamma(x,a,1/b) */
 553 pdf_gamma(x,a,b) :=
 554    if maybe(a > 0 and b > 0) = false
 555      then error("pdf_gamma: parameters a and b must be greater than 0")
 556      else x^(a-1)*exp(-x/b)/(b^a*gamma(a)) * unit_step(x) $
 557
 558 /* R: pgamma(x,a,1/b) */
 559 cdf_gamma(x,a,b) :=
 560    if maybe(a > 0 and b > 0) = false
 561      then error("cdf_gamma: parameters a and b must be greater than 0")
 562      else (1 - gamma_incomplete_regularized(a,x/b)) * unit_step(x) $
 563
 564 /* R: qgamma(q,a,1/b) */
 565 quantile_gamma(q,a,b) :=
 566    if maybe(a > 0 and b >0 and q >= 0 and q <= 1) = false
 567      then error("quantile_gamma: check input parameters")
 568    elseif sign(q)='zero
 569      then 0
 570    elseif sign(q-1)='zero
 571      then 'inf
 572      else /* approximate procedure */
 573           block([fq: float(q), fa: float(a)],
 574             if numberp(fq) and numberp(fa)
 575               then(if fq = 0.0 then return(0),
 576                    if fq = 1.0 then return('inf),
 577                    b * ?iigamma(fq, fa))
 578               else error("quantile_gamma: need numeric arguments for approximate procedure")) $
 579
 580 mean_gamma(a,b) :=
 581    if maybe(a > 0 and b > 0) = false
 582      then error("mean_gamma: parameters a and b must be greater than 0")
 583      else a*b $
 584
 585 var_gamma(a,b) :=
 586    if maybe(a > 0 and b > 0) = false
 587      then error("var_gamma: parameters a and b must be greater than 0")
 588      else a*b^2 $
 589
 590 std_gamma(a,b) := sqrt(var_gamma(a,b)) $
 591
 592 skewness_gamma(a,b) :=
 593    if maybe(a > 0 and b > 0) = false
 594      then error("skewness_gamma: parameters a and b must be greater than 0")
 595      else 2/sqrt(a) $
 596
 597 kurtosis_gamma(a,b) :=
 598    if maybe(a > 0 and b > 0) = false
 599      then error("kurtosis_gamma: parameters a and b must be greater than 0")
 600      else 6/a $
 601
 602 random_gamma(a,b,[num]) :=
 603    if maybe(a > 0 and b > 0) = false
 604      then error("random_gamma: parameters a and b must be greater than 0")
 605      else block([no, fa:float(a), fb:float(b)],
 606             if not numberp(fa) or not numberp(fb)
 607               then error("random_gamma: need numeric arguments for approximate procedure"),
 608             if length(num) = 0 then no: 0 else no: num[1],
 609             if integerp(no) and no >= 0
 610               then ?rndgamma(fa,fb,no)
 611               else error("random_gamma: check sample size") ) $
 612
 613
 614
 615 /*         BETA DISTRIBUTION          */
 616
 617 pdf_beta(x,a,b) :=
 618    if maybe(a > 0 and b > 0) = false
 619      then error("pdf_beta: parameters a and b must be greater than 0")
 620      else x^(a-1)*(1-x)^(b-1)/beta(a,b)*(unit_step(x)-unit_step(x-1)) $
 621
 622 /* R: pbeta(x,a,b) */
 623 cdf_beta(x,a,b) :=
 624    if maybe(a > 0 and b > 0) = false
 625      then error("cdf_beta: parameters a and b must be greater than 0")
 626      else beta_incomplete_regularized(a,b,x) * (unit_step(x)-unit_step(x-1)) + unit_step(x-1)$
 627
 628 /* R: qbeta(q,a,b) */
 629 quantile_beta(q,a,b) :=
 630    if maybe(a > 0 and b >0 and q >= 0 and q <= 1) = false
 631      then error("quantile_beta: check input parameters")
 632    elseif sign(q)='zero
 633      then 0
 634    elseif sign(q-1)='zero
 635      then 1
 636      else block([fq: float(q), fa: float(a), fb: float(b)],
 637           if numberp(fq) and numberp(fa) and numberp(fb)
 638             then(if fq = 0.0 then return(0),
 639                  if fq = 1.0 then return(1),
 640                  ?iibeta(fq, fa, fb))
 641             else error("quantile_beta: need numeric arguments for approximate procedure")) $
 642
 643 mean_beta(a,b) :=
 644    if maybe(a > 0 and b > 0) = false
 645      then error("mean_beta: parameters a and b must be greater than 0")
 646      else a/(a+b) $
 647
 648 var_beta(a,b) :=
 649    if maybe(a > 0 and b > 0) = false
 650      then error("var_beta: parameters a and b must be greater than 0")
 651      else a*b/((a+b)^2*(a+b+1)) $
 652
 653 std_beta(a,b) := sqrt(var_beta(a,b)) $
 654
 655 skewness_beta(a,b) :=
 656    if maybe(a > 0 and b > 0) = false
 657      then error("skewness_beta: parameters a and b must be greater than 0")
 658      else 2*(b-a)*sqrt(a+b+1)/(sqrt(a*b)*(a+b+2)) $
 659
 660 kurtosis_beta(a,b) :=
 661    if maybe(a > 0 and b > 0) = false
 662      then error("kurtosis_beta: parameters a and b must be greater than 0")
 663      else 3*(a+b+1)*(2*(a+b)^2+a*b*(a+b-6)) / (a*b*(a+b+2)*(a+b+3)) - 3 $
 664
 665 random_beta(a,b,[num]) :=
 666    if maybe(a > 0 and b > 0) = false
 667      then error("random_beta: parameters a and b must be greater than 0")
 668      else block([no, fa:float(a), fb:float(b)],
 669             if not numberp(fa) or not numberp(fb)
 670               then error("random_beta: need numeric arguments for approximate procedure"),
 671             if length(num) = 0 then no: 0 else no: num[1],
 672             if integerp(no) and no >= 0
 673               then ?rndbeta(fa,fb,no)
 674               else error("random_beta: check sample size") ) $
 675
 676
 677
 678 /*         CONTINUOUS UNIFORM DISTRIBUTION          */
 679
 680 pdf_continuous_uniform(x,a,b) :=
 681    if maybe(b - a > 0) = false
 682      then error("pdf_continuous_uniform: parameter b must be greater than a")
 683      else (unit_step(x-a)-unit_step(x-b))/(b-a) $
 684
 685 cdf_continuous_uniform(x,a,b) :=
 686    if maybe(b - a > 0) = false
 687      then error("cdf_continuous_uniform: parameter b must be greater than a")
 688      else (x-a)/(b-a)*(unit_step(x-a)-unit_step(x-b)) + unit_step(x-b) $
 689
 690 quantile_continuous_uniform(q,a,b) :=
 691    if maybe(b - a > 0 and q >= 0 and q <= 1) = false
 692      then error("quantile_continuous_uniform: check input parameters")
 693    elseif sign(q)='zero
 694      then a
 695    elseif sign(q-1)='zero
 696      then b
 697      else a + q * (b-a) $
 698
 699 mean_continuous_uniform(a,b) :=
 700    if maybe(b - a > 0) = false
 701      then error("mean_continuous_uniform: parameter b must be greater than a")
 702      else (a+b)/2 $
 703
 704 var_continuous_uniform(a,b) :=
 705    if maybe(b - a > 0) = false
 706      then error("var_continuous_uniform: parameter b must be greater than a")
 707      else (b-a)^2/12 $
 708
 709 std_continuous_uniform(a,b) := sqrt(var_continuous_uniform(a,b))$
 710
 711 skewness_continuous_uniform(a,b) :=
 712    if maybe(b - a > 0) = false
 713      then error("skewness_continuous_uniform: parameter b must be greater than a")
 714      else 0 $
 715
 716 kurtosis_continuous_uniform(a,b) :=
 717    if maybe(b - a > 0) = false
 718      then error("kurtosis_continuous_uniform: parameter b must be greater than a")
 719      else -6/5 $
 720
 721 /* This is a direct application of the maxima         */
 722 /* random function. Make describe(random) for details */
 723 random_continuous_uniform(a,b,[num]) :=
 724    if maybe(b - a > 0) = false
 725      then error("random_continuous_uniform: parameter b must be greater than a")
 726      else block([no, f:float(b-a), listarith:true],
 727             if not numberp(f)
 728               then error("random_continuous_uniform: need numeric arguments for approximate procedure"),
 729             if length(num) = 0 then no: 0 else no: num[1],
 730             if integerp(no)
 731               then if no = 0
 732                      then a + random(f)
 733                      else a + makelist(random(f),k,no)
 734               else error("random_continuous_uniform: check sample size") ) $
 735
 736
 737
 738 /*         LOGISTIC DISTRIBUTION          */
 739
 740 pdf_logistic(x,a,b) :=
 741    if maybe(b > 0) = false
 742      then error("pdf_logistic: parameter b must be greater than 0")
 743      else block([y: exp((a-x)/b)], y/(b*(1+y)^2)) $
 744
 745 cdf_logistic(x,a,b) :=
 746    if maybe(b > 0) = false
 747      then error("cdf_logistic: parameter b must be greater than 0")
 748      else 1/(1+exp((a-x)/b)) $
 749
 750 quantile_logistic(q,a,b) :=
 751    if maybe(b > 0 and q >= 0 and q <= 1) = false
 752      then error("quantile_logistic: check input parameters")
 753    elseif sign(q)='zero
 754      then 'minf
 755    elseif sign(q-1)='zero
 756      then 'inf
 757      else a - b * log(1/q-1) $
 758
 759 mean_logistic(a,b) :=
 760    if maybe(b > 0) = false
 761      then error("mean_logistic: parameter b must be greater than 0")
 762      else a $
 763
 764 var_logistic(a,b) :=
 765    if maybe(b > 0) = false
 766      then error("var_logistic: parameter b must be greater than 0")
 767      else b^2*%pi^2/3 $
 768
 769 std_logistic(a,b) := sqrt(var_logistic(a,b)) $
 770
 771 skewness_logistic(a,b) :=
 772    if maybe(b > 0) = false
 773      then error("skewness_logistic: parameter b must be greater than 0")
 774      else 0 $
 775
 776 kurtosis_logistic(a,b) :=
 777    if maybe(b > 0) = false
 778      then error("kurtosis_logistic: parameter b must be greater than 0")
 779      else 6/5 $
 780
 781 /* inverse method */
 782 random_logistic(a,b,[num]) :=
 783    if maybe(b > 0) = false
 784      then error("random_logistic: parameter b must be greater than 0")
 785      else block([no, fa:float(a), fb:float(b)],
 786             if not numberp(fa) or not numberp(fb)
 787               then error("random_logistic: need numeric arguments for approximate procedure"),
 788             if length(num) = 0 then no: 0 else no: num[1],
 789             if integerp(no)
 790               then if no = 0
 791                      then a - b * log(1/random(1.0) - 1.0)
 792                      else a - b * map('log,1/makelist(random(1.0),k,no) - 1.0)
 793               else error("random_logistic: check sample size")) $
 794
 795
 796
 797 /*         PARETO DISTRIBUTION          */
 798
 799 pdf_pareto(x,a,b) :=
 800    if maybe(a > 0 and b > 0) = false
 801      then error("pdf_pareto: parameters a and b must be greater than 0")
 802      else a*b^a*x^(-a-1) * unit_step(x-b) $
 803
 804 cdf_pareto(x,a,b) :=
 805    if maybe(a > 0 and b > 0) = false
 806      then error("pdf_pareto: parameters a and b must be greater than 0")
 807      else (1-(b/x)^a) * unit_step(x-b) $
 808
 809 quantile_pareto(q,a,b) :=
 810    if maybe(a > 0 and b > 0 and q >= 0 and q <= 1) = false
 811      then error("quantile_pareto: check input parameters")
 812    elseif sign(q)='zero
 813      then b
 814    elseif sign(q-1)='zero
 815      then 'inf
 816      else b / (1-q)^(1/a) $
 817
 818 mean_pareto(a,b) :=
 819    if maybe(a > 1) = false
 820      then error("mean_pareto: parameter a must be greater than 1")
 821    elseif maybe(b > 0) = false
 822      then error("mean_pareto: parameter b must be greater than 0")
 823      else a*b/(a-1) $
 824
 825 var_pareto(a,b) :=
 826    if maybe(a > 2) = false
 827      then error("var_pareto: parameter a must be greater than 2")
 828    elseif maybe(b > 0) = false
 829      then error("var_pareto: parameter b must be greater than 0")
 830      else a*b*b/((a-2)*(a-1)^2) $
 831
 832 std_pareto(a,b) :=
 833    if maybe(a > 2) = false
 834      then error("std_pareto: parameter a must be greater than 2")
 835    elseif maybe(b > 0) = false
 836      then error("std_pareto: parameter b must be greater than 0")
 837      else b*sqrt(a/(a-2))/(a-1) $
 838
 839 skewness_pareto(a,b) :=
 840    if maybe(a > 3) = false
 841      then error("skewness_pareto: parameter a must be greater than 3")
 842    elseif maybe(b > 0) = false
 843      then error("skewness_pareto: parameter b must be greater than 0")
 844      else 2*(a+1)*sqrt(a-2)/((a-3)*sqrt(a)) $
 845
 846 kurtosis_pareto(a,b) :=
 847    if maybe(a > 4) = false
 848      then error("kurtosis_pareto: parameter a must be greater than 4")
 849    elseif maybe(b > 0) = false
 850      then error("kurtosis_pareto: parameter b must be greater than 0")
 851      else (6*(a^3+a^2-6*a-2)) / (a*(a^2-7*a+12)) - 3 $
 852
 853 /* inverse method */
 854 random_pareto(a,b,[num]) :=
 855    if maybe(a > 0 and b > 0) = false
 856      then error("random_pareto: parameters a and b must be greater than 0")
 857      else block([no, fa:float(a), fb:float(b)],
 858             if not numberp(fa) or not numberp(fb)
 859               then error("random_pareto: need numeric arguments for approximate procedure"),
 860             if length(num) = 0 then no: 0 else no: num[1],
 861             if integerp(no)
 862               then if no = 0
 863                      then b / (1-random(1.0))^(1.0/a)
 864                      else b / (1-makelist(random(1.0),k,1,no))^(1.0/a)
 865               else error("random_pareto: check sample size")) $
 866
 867
 868
 869 /*         WEIBULL DISTRIBUTION          */
 870
 871 pdf_weibull(x,a,b) :=
 872    if maybe(a > 0 and b > 0) = false
 873      then error("pdf_weibull: parameters a and b must be greater than 0")
 874      else a/b*(x/b)^(a-1)*exp(-(x/b)^a) * unit_step(x) $
 875
 876 cdf_weibull(x,a,b) :=
 877    if maybe(a > 0 and b > 0) = false
 878      then error("cdf_weibull: parameters a and b must be greater than 0")
 879      else (1-exp(-(x/b)^a)) * unit_step(x) $
 880
 881 /* R: qweibull(q,a,b) */
 882 quantile_weibull(q,a,b) :=
 883    if maybe(a > 0 and b > 0 and q >= 0 and q <= 1) = false
 884      then error("quantile_weibull: check input parameters")
 885    elseif sign(q)='zero
 886      then 0
 887    elseif sign(q-1)='zero
 888      then 'inf
 889      else b * (-log(1-q))^(1/a) $
 890
 891 mean_weibull(a,b) :=
 892    if maybe(a > 0 and b > 0) = false
 893      then error("mean_weibull: parameters a and b must be greater than 0")
 894      else gamma(1/a+1)*b $
 895
 896 var_weibull(a,b) :=
 897    if maybe(a > 0 and b > 0) = false
 898      then error("var_weibull: parameters a and b must be greater than 0")
 899      else (gamma(2/a+1)-gamma(1/a+1)^2)*b^2 $
 900
 901 std_weibull(a,b) :=
 902    if maybe(a > 0 and b > 0) = false
 903      then error("std_weibull: parameters a and b must be greater than 0")
 904      else sqrt((gamma(2/a+1)-gamma(1/a+1)^2))*b $
 905
 906 skewness_weibull(a,b) :=
 907    if maybe(a > 0 and b > 0) = false
 908      then error("skewness_weibull: parameters a and b must be greater than 0")
 909      else (gamma(3/a+1)-3*gamma(1/a+1)*gamma(2/a+1)+2*gamma(1/a+1)^3)  /
 910           (gamma(2/a+1)-gamma(1/a+1)^2)^(3/2) $
 911
 912 kurtosis_weibull(a,b) :=
 913    if maybe(a > 0 and b > 0) = false
 914      then error("kurtosis_weibull: parameters a and b must be greater than 0")
 915      else (gamma(4/a+1)-4*gamma(1/a+1)*gamma(3/a+1)+
 916                        6*gamma(1/a+1)^2*gamma(2/a+1)-3*gamma(1/a+1)^4) /
 917           (gamma(2/a+1)-gamma(1/a+1)^2)^2 - 3 $
 918
 919 /* inverse method */
 920 random_weibull(a,b,[num]) :=
 921    if maybe(a > 0 and b > 0) = false
 922      then error("random_weibull: parameters a and b must be greater than 0")
 923      else block([no, fa:float(a), fb:float(b)],
 924             if not numberp(fa) or not numberp(fb)
 925               then error("random_weibull: need numeric arguments for approximate procedure"),
 926             if length(num) = 0 then no: 0 else no: num[1],
 927             if integerp(no)
 928               then if no = 0
 929                      then b * (-log(random(1.0)))^(1.0/a)
 930                      else b * (-map('log,makelist(random(1.0),k,1,no)))^(1.0/a)
 931               else error("random_weibull: check sample size")) $
 932
 933
 934
 935 /*          RAYLEIGH DISTRIBUTION              */
 936 /* Rayleigh(b) is equivalent to Weibull(2,1/b) */
 937
 938 pdf_rayleigh(x,b):=pdf_weibull(x,2,1/b)$
 939
 940 cdf_rayleigh(x,b):=cdf_weibull(x,2,1/b)$
 941
 942 quantile_rayleigh(q,b):=quantile_weibull(q,2,1/b)$
 943
 944 mean_rayleigh(b):=mean_weibull(2,1/b)$
 945
 946 var_rayleigh(b):=var_weibull(2,1/b)$
 947
 948 std_rayleigh(b):=std_weibull(2,1/b)$
 949
 950 skewness_rayleigh(b):=skewness_weibull(2,1/b)$
 951
 952 kurtosis_rayleigh(b):=kurtosis_weibull(2,1/b)$
 953
 954 /* inverse method */
 955 random_rayleigh(b,[num]) :=
 956    if maybe(b > 0) = false
 957      then error("random_rayleigh: parameter b must be greater than 0")
 958      else block([no, fb:float(b)],
 959             if not numberp(fb)
 960               then error("random_rayleigh: need numeric argument for approximate procedure"),
 961             if length(num) = 0 then no: 0 else no: num[1],
 962             if integerp(no)
 963               then if no = 0
 964                      then sqrt(-log(random(1.0))) / b
 965                      else sqrt(-map('log,makelist(random(1.0),k,1,no))) / b
 966               else error("random_rayleigh: check sample size")) $
 967
 968
 969
 970 /*         LAPLACE DISTRIBUTION          */
 971
 972 pdf_laplace(x,a,b) :=
 973    if maybe(b > 0) = false
 974      then error("pdf_laplace: parameter b must be greater than 0")
 975      else exp(-abs(x-a)/b)/(2*b) $
 976
 977 cdf_laplace(x,a,b) :=
 978    if maybe(b > 0) = false
 979      then error("pdf_laplace: parameter b must be greater than 0")
 980      else (1+signum(x-a)-signum(x-a)*exp(-abs(x-a)/b)) / 2  $
 981
 982 quantile_laplace(q,a,b) :=
 983    if maybe(b > 0 and q >= 0 and q <= 1) = false
 984      then error("quantile_laplace: check input parameters")
 985    elseif sign(q)='zero
 986      then 'minf
 987    elseif sign(q-1)='zero
 988      then 'inf
 989      else a - b * signum(2*q-1) * log(1 - signum(2*q-1) * (2*q-1))$
 990
 991 mean_laplace(a,b) :=
 992    if maybe(b > 0) = false
 993      then error("mean_laplace: parameter b must be greater than 0")
 994      else a $
 995
 996 var_laplace(a,b) :=
 997    if maybe(b > 0) = false
 998      then error("var_laplace: parameter b must be greater than 0")
 999      else 2*b*b $
1000
1001 std_laplace(a,b) :=
1002    if maybe(b > 0) = false
1003      then error("std_laplace: parameter b must be greater than 0")
1004      else sqrt(2)*b $
1005
1006 skewness_laplace(a,b) :=
1007    if maybe(b > 0) = false
1008      then error("skewness_laplace: parameter b must be greater than 0")
1009      else 0 $
1010
1011 kurtosis_laplace(a,b) :=
1012    if maybe(b > 0) = false
1013      then error("kurtosis_laplace: parameter b must be greater than 0")
1014      else 3 $
1015
1016 /* inverse method */
1017 random_laplace(a,b,[num]) :=
1018    if maybe(b > 0) = false
1019      then error("random_laplace: parameter b must be greater than 0")
1020      else block([no, fa:float(a), fb:float(b)],
1021             if not numberp(fa) or not numberp(fb)
1022               then error("random_laplace: need numeric arguments for approximate procedure"),
1023             if length(num) = 0 then no: 0 else no: num[1],
1024             if integerp(no)
1025               then if no = 0
1026                      then quantile_laplace(random(1.0),a,b)
1027                      else makelist(quantile_laplace(random(1.0),a,b),k,1,no)
1028               else error("random_laplace: check sample size")) $
1029
1030
1031
1032 /*         CAUCHY (OR LORENTZ) DISTRIBUTION          */
1033
1034 pdf_cauchy(x,a,b) :=
1035    if maybe(b > 0) = false
1036      then error("pdf_cauchy: parameter b must be greater than 0")
1037      else b/(%pi*((x-a)^2+b^2)) $
1038
1039 cdf_cauchy(x,a,b) :=
1040    if maybe(b > 0) = false
1041      then error("pdf_cauchy: parameter b must be greater than 0")
1042      else 1/2+atan((x-a)/b)/%pi $
1043
1044 quantile_cauchy(q,a,b) :=
1045    if maybe(b > 0 and q >= 0 and q <= 1) = false
1046      then error("quantile_cauchy: check input parameters")
1047    elseif sign(q)='zero
1048      then 'minf
1049    elseif sign(q-1)='zero
1050      then 'inf
1051      else a + b * tan(%pi * (q - 1/2)) $
1052
1053 /* Note: integrals for Cauchy moments are divergent */
1054
1055 /* inverse method */
1056 random_cauchy(a,b,[num]) :=
1057    if maybe(b > 0) = false
1058      then error("random_cauchy: parameter b must be greater than 0")
1059      else block([no, fa:float(a), fb:float(b)],
1060             if not numberp(fa) or not numberp(fb)
1061               then error("random_cauchy: need numeric arguments for approximate procedure"),
1062             if length(num) = 0 then no: 0 else no: num[1],
1063             if integerp(no)
1064               then if no = 0
1065                      then tan (float(%pi) * (random(1.0)-0.5))*fb+fa
1066                      else fa + fb * map('tan, 3.141592653589793 * (makelist(random(1.0),k,1,no)-0.5))
1067               else error("random_cauchy: check sample size")) $
1068
1069
1070
1071 /*         GUMBEL (OR EXTREME VALUE) DISTRIBUTION          */
1072
1073 pdf_gumbel(x,a,b) :=
1074    if maybe(b > 0) = false
1075      then error("pdf_gumbel: parameter b must be greater than 0")
1076      else exp((a-x)/b-exp((a-x)/b))/b $
1077
1078 cdf_gumbel(x,a,b) :=
1079    if maybe(b > 0) = false
1080      then error("cdf_gumbel: parameter b must be greater than 0")
1081      else exp(-exp((a-x)/b)) $
1082
1083 quantile_gumbel(q,a,b) :=
1084    if maybe(b > 0 and q >= 0 and q <= 1) = false
1085      then error("quantile_gumbel: check input parameters")
1086    elseif sign(q)='zero
1087      then 'minf
1088    elseif sign(q-1)='zero
1089      then 'inf
1090      else a - b * log(-log(q)) $
1091
1092 mean_gumbel(a,b) :=
1093    if maybe(b > 0) = false
1094      then error("mean_gumbel: parameter b must be greater than 0")
1095      else /* %gamma=Euler-Mascheroni constant */
1096           a + b*%gamma $
1097
1098 var_gumbel(a,b) :=
1099    if maybe(b > 0) = false
1100      then error("var_gumbel: parameter b must be greater than 0")
1101      else b*b*%pi*%pi/6 $
1102
1103 std_gumbel(a,b) :=
1104    if maybe(b > 0) = false
1105      then error("std_gumbel: parameter b must be greater than 0")
1106      else b*%pi/sqrt(6) $
1107
1108 skewness_gumbel(a,b) :=
1109    if maybe(b > 0) = false
1110      then error("std_gumbel: parameter b must be greater than 0")
1111      else 12*sqrt(6)*zeta(3)/%pi^3 $
1112
1113 kurtosis_gumbel(a,b) :=
1114    if maybe(b > 0) = false
1115      then error("kurtosis_gumbel: parameter b must be greater than 0")
1116      else 12/5 $
1117
1118 /* inverse method */
1119 random_gumbel(a,b,[num]) :=
1120    if maybe(b > 0) = false
1121      then error("random_gumbel: parameter b must be greater than 0")
1122      else block([no, fa:float(a), fb:float(b)],
1123             if not numberp(fa) or not numberp(fb)
1124               then error("random_gumbel: need numeric arguments for approximate procedure"),
1125             if length(num) = 0 then no: 0 else no: num[1],
1126             if integerp(no)
1127               then if no = 0
1128                      then a - b * log(-log(random(1.0)))
1129                      else a - b * map('log,-map('log,makelist(random(1.0),k,1,no)))
1130               else error("random_gumbel: check sample size")) $
1131
1132
1133
1134 /*         BINOMIAL DISTRIBUTION          */
1135
1136 /* R: dbinom(x,n,p) */
1137 pdf_binomial(x,n,p) :=
1138   if maybe(p >= 0 and p <= 1) = false
1139      then error("pdf_binomial: p must be a probability")
1140   /* when n is not an integer, R returns NaN, we throw an error */
1141   elseif maybe(n>=0)=false or numberp(n) and n-floor(n) > 0
1142      then error("pdf_binomial: n must be a positive integer")
1143   elseif sign(x) = 'neg or sign(n-x) = 'neg or numberp(x) and x-floor(x) > 0
1144      then 0
1145   elseif sign(p) = 'zero
1146      then kron_delta(x,0)
1147   elseif sign(1-p) = 'zero
1148      then kron_delta(x,n)
1149      else binomial(n,x)*p^x*(1-p)^(n-x) $
1150
1151 /* R: pbinom(x,n,p) */
1152 cdf_binomial(x,n,p):=
1153   if maybe(p >= 0 and p <= 1) = false
1154      then error("cdf_binomial: p must be a probability")
1155   elseif maybe(n>=0)=false or numberp(n) and n-floor(n) > 0
1156      then error("cdf_binomial: n must be a positive integer")
1157   elseif sign(x) = 'neg
1158      then 0
1159   elseif member(sign(x-n),['pos,'pz,'zero])
1160      then 1
1161      else beta_incomplete_regularized(n-floor(x),floor(x)+1,1-p) $
1162
1163 /* R: qbinom(q,n,p) */
1164 quantile_binomial(q,n,p):=
1165   if maybe(p >= 0 and p <= 1) = false
1166      then error("quantile_binomial: p must be a probability")
1167   elseif maybe(q >= 0 and q <= 1) = false
1168      then error("quantile_binomial: q must be a probability")
1169   elseif maybe(n>=0)=false or numberp(n) and n-floor(n) > 0
1170      then error("quantile_binomial: n must be a positive integer")
1171   elseif sign(q) = 'zero
1172      then 0
1173   elseif sign(1-q) = 'zero
1174      then n
1175   elseif float(q) <= float((1-p)^n)
1176      then 0
1177   elseif numberp(float(q)) and numberp(float(n)) and numberp(float(p))
1178     then /* partition method */
1179          block([a:0, b:n, m],
1180            while (b-a>1) do (
1181              m: 0.5*(a+b),
1182              if cdf_binomial(m,n,p) < q
1183                then a: m
1184                else b: m ),
1185            floor(b))
1186     else error("quantile_binomial: need numeric arguments for approximate procedure") $
1187
1188 mean_binomial(n,p) :=
1189   if maybe(p >= 0 and p <= 1) = false
1190      then error("mean_binomial: p must be a probability")
1191   elseif maybe(n>=0)=false or numberp(n) and n-floor(n) > 0
1192      then error("mean_binomial: n must be a positive integer")
1193      else n*p $
1194
1195 var_binomial(n,p) :=
1196   if maybe(p >= 0 and p <= 1) = false
1197      then error("var_binomial: p must be a probability")
1198   elseif maybe(n>=0)=false or numberp(n) and n-floor(n) > 0
1199      then error("var_binomial: n must be a positive integer")
1200      else n*p*(1-p) $
1201
1202 std_binomial(n,p) :=
1203   if maybe(p >= 0 and p <= 1) = false
1204      then error("std_binomial: p must be a probability")
1205   elseif maybe(n>=0)=false or numberp(n) and n-floor(n) > 0
1206      then error("std_binomial: n must be a positive integer")
1207      else sqrt(n*p*(1-p)) $
1208
1209 skewness_binomial(n,p) :=
1210   if maybe(p >= 0 and p <= 1) = false
1211      then error("std_binomial: p must be a probability")
1212   elseif maybe(n>=0)=false or numberp(n) and n-floor(n) > 0
1213      then error("std_binomial: n must be a positive integer")
1214      else (1-2*p)/sqrt(n*p*(1-p)) $
1215
1216 kurtosis_binomial(n,p) :=
1217   if maybe(p >= 0 and p <= 1) = false
1218      then error("kurtosis_binomial: p must be a probability")
1219   elseif maybe(n>=0)=false or numberp(n) and n-floor(n) > 0
1220      then error("kurtosis_binomial: n must be a positive integer")
1221      else (1-6*p*(1-p))/(n*p*(1-p)) $
1222
1223 random_binomial(n,p,[num]) :=
1224   if maybe(p >= 0 and p <= 1) = false
1225      then error("random_binomial: p must be a probability")
1226   elseif maybe(n>=0)=false or numberp(n) and n-floor(n) > 0
1227      then error("random_binomial: n must be a positive integer")
1228      else block([no, fn:float(n), fp:float(p)],
1229             if not numberp(fn) or not numberp(fp)
1230               then error("random_binomial: need numeric arguments for approximate procedure"),
1231             if length(num) = 0 then no: 0 else no: num[1],
1232             if integerp(no)
1233               then ?rndbinomial(n,fp,no)
1234               else error("random_binomial: check sample size")) $
1235
1236
1237
1238 /*         POISSON DISTRIBUTION          */
1239
1240 /* R: dpois(x,m) */
1241 pdf_poisson(x,m) :=
1242   if maybe(m >= 0) = false
1243      then error("pdf_poisson: m must be positive")
1244   elseif sign(m) = 'zero
1245      then kron_delta(0,x)
1246   elseif sign(x) = 'neg or numberp(x) and x-floor(x) > 0
1247      then 0
1248      else exp(-m)*m^x/x! $
1249
1250 /* R: ppois(x,m) */
1251 cdf_poisson(x,m):=
1252   if maybe(m >= 0) = false
1253      then error("cdf_poisson: m must be positive")
1254   elseif sign(x) = 'neg
1255      then 0
1256      else gamma_incomplete_regularized(floor(x)+1, m) $
1257
1258 /* R: qpois(q,m) */
1259 quantile_poisson(q,m):=
1260   if maybe(m >= 0) = false
1261      then error("quantile_poisson: m must be positive")
1262   elseif maybe(q >= 0 and q <= 1) = false
1263      then error("quantile_poisson: q must be a probability")
1264   elseif sign(m)='zero
1265      then 0
1266   elseif member(sign(q-exp(-m)),['nz,'neg,'zero])
1267      then 0
1268   elseif sign(1-q)='zero
1269      then 'inf
1270   elseif numberp(float(q)) and numberp(float(m))
1271     then /* partition method */
1272          block([a, b:1.0,mm],
1273            while (cdf_poisson(b,m) < q) do b : 2.0*b,
1274            a: b/2.0,
1275            while (b-a>1) do(
1276               mm: 0.5*(a+b),
1277               if cdf_poisson(mm,m) < q
1278                  then a: mm
1279                  else b: mm ),
1280            floor(b))
1281     else error("quantile_poisson: need numeric arguments for approximate procedure") $
1282
1283 mean_poisson(m):=
1284   if maybe(m >= 0) = false
1285      then error("mean_poisson: m must be positive")
1286      else m $
1287
1288 var_poisson(m):=
1289   if maybe(m >= 0) = false
1290      then error("var_poisson: m must be positive")
1291      else m $
1292
1293 std_poisson(m):=
1294   if maybe(m >= 0) = false
1295      then error("std_poisson: m must be positive")
1296      else sqrt(m) $
1297
1298 skewness_poisson(m):=
1299   if maybe(m >= 0) = false
1300      then error("skewness_poisson: m must be positive")
1301      else 1/sqrt(m) $
1302
1303 kurtosis_poisson(m):=
1304   if maybe(m >= 0) = false
1305      then error("kurtosis_poisson: m must be positive")
1306      else 1/m $
1307
1308 random_poisson(m,[num]) :=
1309   if maybe(m >= 0) = false
1310      then error("random_poisson: m must be positive")
1311      else block([no, fm:float(m)],
1312             if not numberp(fm)
1313               then error("random_poisson: need numeric argument for approximate procedure"),
1314             if length(num) = 0 then no: 0 else no: num[1],
1315             if integerp(no)
1316               then ?rndpoisson(fm,no)
1317               else error("random_poisson: check sample size")) $
1318
1319
1320
1321 /*           BERNOULLI DISTRIBUTION            */
1322 /* Bernoulli(p) is equivalent to binomial(1,p) */
1323
1324 pdf_bernoulli(x,p) :=
1325   if maybe(p >= 0 and p <= 1) = false
1326      then error("pdf_bernoulli: p must be a probability")
1327   elseif sign(x) = 'neg or sign(1-x) = 'neg or numberp(x) and x-floor(x) > 0
1328      then 0
1329   elseif sign(p) = 'zero
1330      then kron_delta(x,0)
1331   elseif sign(1-p) = 'zero
1332      then kron_delta(x,1)
1333      else p^x*(1-p)^(1-x) $
1334
1335 cdf_bernoulli(x,p):= cdf_binomial(x,1,p)$
1336
1337 quantile_bernoulli(q,p):= quantile_binomial(q,1,p)$
1338
1339 mean_bernoulli(p):=mean_binomial(1,p)$
1340
1341 var_bernoulli(p):=var_binomial(1,p)$
1342
1343 std_bernoulli(p):=std_binomial(1,p)$
1344
1345 skewness_bernoulli(p):=skewness_binomial(1,p)$
1346
1347 kurtosis_bernoulli(p):=kurtosis_binomial(1,p)$
1348
1349 /* This is a direct application of the maxima
1350    random function. Make describe(random) for details */
1351 random_bernoulli(p,[num]) :=
1352   if maybe(p >= 0 and p <= 1) = false
1353      then error("random_bernoulli: p must be a probability")
1354      else block([no, fp:float(p)],
1355             if not numberp(fp)
1356               then error("random_bernoulli: need numeric argument for approximate procedure"),
1357             if length(num) = 0 then no: 0 else no: num[1],
1358             if integerp(no)
1359               then if no = 0
1360                      then if random(1.0)<=fp then 1 else 0
1361                      else makelist(if random(1.0)<=fp then 1 else 0,k,1,no)
1362               else error("random_bernoulli: check sample size")) $
1363
1364
1365
1366 /*         GEOMETRIC (OR PASCAL) DISTRIBUTION          */
1367
1368 /* R: dgeom(x,p) */
1369 pdf_geometric(x,p) :=
1370   if maybe(p > 0 and p <= 1) = false
1371      then error("pdf_geometric: p must be a non zero probability")
1372   elseif sign(p-1) = 'zero
1373      then if sign(x) = 'zero then 1 else 0
1374   elseif sign(x) = 'neg or numberp(x) and x-floor(x) > 0
1375      then 0
1376      else p*(1-p)^x $
1377
1378 /* R: pgeom(q,p) */
1379 cdf_geometric(x,p) :=
1380   if maybe(p > 0 and p <= 1) = false
1381      then error("cdf_geometric: p must be a non zero probability")
1382   elseif sign(p-1) = 'zero
1383      then if member(sign(x), ['zero, 'pos, 'pz]) then 1 else 0
1384   elseif sign(x) = 'neg
1385      then 0
1386      else 1-(1-p)^(floor(x)+1) $
1387
1388 /* R: qgeom(q,p) */
1389 quantile_geometric(q,p) :=
1390   if maybe(p > 0 and p <= 1) = false
1391      then error("quantile_geometric: p must be a non zero probability")
1392   elseif maybe(q >= 0 and q <= 1) = false
1393      then error("quantile_geometric: q must be a probability")
1394   elseif sign(q-1) = 'zero
1395      then 'inf
1396   elseif sign(q) = 'zero
1397      then 0
1398   elseif sign(p-1) = 'zero
1399      then 1
1400      else ceiling(log(1-q)/log(1-p)-1) $
1401
1402 mean_geometric(p) :=
1403   if maybe(p > 0 and p <= 1) = false
1404      then error("mean_geometric: p must be a non zero probability")
1405      else 1/p-1 $
1406
1407 var_geometric(p) :=
1408   if maybe(p > 0 and p <= 1) = false
1409      then error("var_geometric: p must be a non zero probability")
1410      else (1-p)/p^2 $
1411
1412 std_geometric(p) :=
1413   if maybe(p > 0 and p <= 1) = false
1414      then error("std_geometric: p must be a non zero probability")
1415      else sqrt(1-p)/p $
1416
1417 skewness_geometric(p) :=
1418   if maybe(p > 0 and p <= 1) = false
1419      then error("skewness_geometric: p must be a non zero probability")
1420      else (2-p)/sqrt(1-p) $
1421
1422 kurtosis_geometric(p) :=
1423   if maybe(p > 0 and p <= 1) = false
1424      then error("kurtosis_geometric: p must be a non zero probability")
1425      else (p^2+6-6*p)/(1-p) $
1426
1427 random_geometric(p,[num]) :=
1428   if maybe(p > 0 and p <= 1) = false
1429      then error("random_geometric: p must be a non zero probability")
1430      else block([no, fp:float(p)],
1431             if not numberp(fp)
1432               then error("random_geometric: need numeric argument for approximate procedure"),
1433             if length(num) = 0 then no: 0 else no: num[1],
1434             if integerp(no)
1435               then ?rndgeo(fp,no)
1436               else error("random_geometric: check sample size")) $
1437
1438
1439
1440 /*         DISCRETE UNIFORM DISTRIBUTION          */
1441
1442 pdf_discrete_uniform(x,n) :=
1443   if maybe(n > 0) = false or numberp(n) and n-floor(n) > 0
1444     then error("pdf_discrete_uniform: n must be a positive integer")
1445   elseif sign(x-1) = 'neg or sign(n-x) = 'neg or numberp(x) and x-floor(x) > 0
1446     then 0
1447     else 1/n $
1448
1449 cdf_discrete_uniform(x,n) :=
1450   if maybe(n > 0) = false or numberp(n) and n-floor(n) > 0
1451     then error("cdf_discrete_uniform: n must be a positive integer")
1452   elseif sign(x-1) = 'neg
1453     then 0
1454   elseif member(sign(x-n), ['pos,'zero, 'pz])
1455     then 1
1456     else floor(x)/n $
1457
1458 quantile_discrete_uniform(q,n) :=
1459   if maybe(n > 0) = false or numberp(n) and n-floor(n) > 0
1460     then error("cdf_discrete_uniform: n must be a positive integer")
1461   elseif maybe(q >= 0 and q <= 1) = false
1462     then error("quantile_discrete_uniform: q must be a probability")
1463   elseif sign(q-1) = 'zero
1464     then n
1465   elseif sign(q) = 'zero
1466     then 1
1467     else ceiling(q*n) $
1468
1469 mean_discrete_uniform(n) :=
1470   if maybe(n > 0) = false or numberp(n) and n-floor(n) > 0
1471     then error("mean_discrete_uniform: n must be a positive integer")
1472     else (1+n)/2 $
1473
1474 var_discrete_uniform(n) :=
1475   if maybe(n > 0) = false or numberp(n) and n-floor(n) > 0
1476     then error("var_discrete_uniform: n must be a positive integer")
1477     else (n^2-1)/12 $
1478
1479 std_discrete_uniform(n) :=
1480   if maybe(n > 0) = false or numberp(n) and n-floor(n) > 0
1481     then error("std_discrete_uniform: n must be a positive integer")
1482     else sqrt((n^2-1)/12) $
1483
1484 skewness_discrete_uniform(n) :=
1485   if maybe(n > 0) = false or numberp(n) and n-floor(n) > 0
1486     then error("skewness_discrete_uniform: n must be a positive integer")
1487     else 0 $
1488
1489 kurtosis_discrete_uniform(n) :=
1490   if maybe(n > 0) = false or numberp(n) and n-floor(n) > 0
1491     then error("kurtosis_discrete_uniform: n must be a positive integer")
1492     else -6/5-12/(5*(n^2-1)) $
1493
1494 /* This is a direct application of the maxima
1495    random function. Make describe(random) for details */
1496 random_discrete_uniform(n,[num]) :=
1497   if maybe(n > 0) = false or numberp(n) and n-floor(n) > 0
1498      then error("random_discrete_uniform: n must be a positive integer")
1499      else block([no, fn:float(n)],
1500             if not numberp(fn)
1501               then error("random_discrete_uniform: need numeric argument for approximate procedure"),
1502             if length(num) = 0 then no: 0 else no: num[1],
1503             if integerp(no)
1504               then if no = 0
1505                      then 1+random(n)
1506                      else 1+makelist(random(n),k,1,no)
1507               else error("random_discrete_uniform: check sample size")) $
1508
1509
1510
1511 /*         HYPERGEOMETRIC DISTRIBUTION          */
1512
1513 /* R: dhyper(x, n1, n2, n) */
1514 pdf_hypergeometric(x,n1,n2,n) :=
1515   if maybe(n1 >= 0) = false or numberp(n1) and n1-floor(n1) > 0 or
1516      maybe(n2 >= 0) = false or numberp(n2) and n2-floor(n2) > 0 or
1517      maybe( n >= 0) = false or numberp(n)  and  n-floor(n)  > 0
1518     then error("pdf_hypergeometric: n1, n2, n must be a positive integers")
1519   elseif maybe (n <= n1 +n2) = false
1520     then error("pdf_hypergeometric: n must be less or equal than n1 + n2")
1521   elseif sign(x-max(0,n-n2))='neg or sign(min(n1,n)-x)='neg or numberp(x) and x-floor(x) > 0
1522     then 0
1523     else binomial(n1,x)*binomial(n2,n-x)/binomial(n1+n2,n) $
1524
1525 /* R: phyper(x, n1, n2, n) */
1526 cdf_hypergeometric(x,n1,n2,n) :=
1527   if maybe(n1 >= 0) = false or numberp(n1) and n1-floor(n1) > 0 or
1528      maybe(n2 >= 0) = false or numberp(n2) and n2-floor(n2) > 0 or
1529      maybe( n >= 0) = false or numberp(n)  and  n-floor(n)  > 0
1530     then error("pdf_hypergeometric: n1, n2, n must be a positive integers")
1531   elseif maybe (n <= n1 +n2) = false
1532     then error("cdf_hypergeometric: n must be less or equal than n1 + n2")
1533   elseif sign(x-max(0,n-n2))='neg
1534     then 0
1535   elseif member(sign(min(n1,n)-x),['neg,'zero,'nz])
1536     then 1
1537     else sum(binomial(n1,k)*binomial(n2,n-k) / binomial(n1+n2,n),k,0,floor(x)) $
1538
1539 quantile_hypergeometric(q, n1, n2, n) :=
1540   if maybe(n1 >= 0) = false or numberp(n1) and n1-floor(n1) > 0 or
1541      maybe(n2 >= 0) = false or numberp(n2) and n2-floor(n2) > 0 or
1542      maybe( n >= 0) = false or numberp(n)  and  n-floor(n)  > 0
1543     then error("quantile_hypergeometric: n1, n2, n must be a positive integers")
1544   elseif maybe (n <= n1 +n2) = false
1545     then error("quantile_hypergeometric: n must be less or equal than n1 + n2")
1546   elseif maybe(q >= 0 and q <= 1) = false
1547      then error("quantile_hypergeometric: q must be a probability")
1548   elseif sign(q-1) = 'zero
1549      then min(n1,n)
1550   elseif member(sign(q - pdf_hypergeometric(max(0, n-n2),n1,n2,n)), ['neg,'zero, 'nz])
1551      then max(0, n-n2)
1552   elseif numberp(float(q)) and numberp(float(n1)) and numberp(float(n2)) and numberp(float(n))
1553     then /* partition method */
1554          block([a: max(0, n-n2), b: min(n1,n),m],
1555            while (b-a>1) do (
1556              m: floor(0.5*(a+b)),
1557              if cdf_hypergeometric(m,n1,n2,n) < q
1558                then a: m
1559                else b: m),
1560            floor(b))
1561     else error("quantile_hypergeometric: need numeric arguments for approximate procedure") $
1562
1563 mean_hypergeometric(n1,n2,n) :=
1564   if maybe(n1 >= 0) = false or numberp(n1) and n1-floor(n1) > 0 or
1565      maybe(n2 >= 0) = false or numberp(n2) and n2-floor(n2) > 0 or
1566      maybe( n >= 0) = false or numberp(n)  and  n-floor(n)  > 0
1567     then error("mean_hypergeometric: n1, n2, n must be a positive integers")
1568   elseif maybe (n <= n1 +n2) = false
1569     then error("mean_hypergeometric: n must be less or equal than n1 + n2")
1570     else n*n1/(n1+n2) $
1571
1572 var_hypergeometric(n1,n2,n) :=
1573   if maybe(n1 >= 0) = false or numberp(n1) and n1-floor(n1) > 0 or
1574      maybe(n2 >= 0) = false or numberp(n2) and n2-floor(n2) > 0 or
1575      maybe( n >= 0) = false or numberp(n)  and  n-floor(n)  > 0
1576     then error("var_hypergeometric: n1, n2, n must be a positive integers")
1577   elseif maybe (n <= n1 +n2) = false
1578     then error("var_hypergeometric: n must be less or equal than n1 + n2")
1579     else block([t:n1+n2], n*n1*n2*(t-n)/(t*t*(t-1))) $
1580
1581 std_hypergeometric(n1,n2,n) :=
1582   if maybe(n1 >= 0) = false or numberp(n1) and n1-floor(n1) > 0 or
1583      maybe(n2 >= 0) = false or numberp(n2) and n2-floor(n2) > 0 or
1584      maybe( n >= 0) = false or numberp(n)  and  n-floor(n)  > 0
1585     then error("std_hypergeometric: n1, n2, n must be a positive integers")
1586   elseif maybe (n <= n1 +n2) = false
1587     then error("std_hypergeometric: n must be less or equal than n1 + n2")
1588     else block([t:n1+n2], sqrt(n*n1*n2*(t-n)/(t-1))/t) $
1589
1590 skewness_hypergeometric(n1,n2,n) :=
1591   if maybe(n1 >= 0) = false or numberp(n1) and n1-floor(n1) > 0 or
1592      maybe(n2 >= 0) = false or numberp(n2) and n2-floor(n2) > 0 or
1593      maybe( n >= 0) = false or numberp(n)  and  n-floor(n)  > 0
1594     then error("skewness_hypergeometric: n1, n2, n must be a positive integers")
1595   elseif maybe (n <= n1 +n2) = false
1596     then error("skewness_hypergeometric: n must be less or equal than n1 + n2")
1597     else block([t:n1+n2],  (n2-n1)*(t-2*n)*sqrt((t-1)/(n*n1*n2*(t-n)))/(t-2)) $
1598
1599 kurtosis_hypergeometric(n1,n2,n) :=
1600   if maybe(n1 >= 0) = false or numberp(n1) and n1-floor(n1) > 0 or
1601      maybe(n2 >= 0) = false or numberp(n2) and n2-floor(n2) > 0 or
1602      maybe( n >= 0) = false or numberp(n)  and  n-floor(n)  > 0
1603     then error("kurtosis_hypergeometric: n1, n2, n must be a positive integers")
1604   elseif maybe (n <= n1 +n2) = false
1605     then error("kurtosis_hypergeometric: n must be less or equal than n1 + n2")
1606     else block([t:n1+n2],
1607            t*t*(t-1)/((t-2)*(t-3)*n*n1*n2*(t-n)) *
1608            (t*(t+1)-6*n*(t-n)+3*n1*n2*(t*t*(n-2)-t*n*n+6*n*(t-n))/t^2) - 3) $
1609
1610 random_hypergeometric(n1,n2,n,[num]) :=
1611   if maybe(n1 >= 0) = false or numberp(n1) and n1-floor(n1) > 0 or
1612      maybe(n2 >= 0) = false or numberp(n2) and n2-floor(n2) > 0 or
1613      maybe( n >= 0) = false or numberp(n)  and  n-floor(n)  > 0
1614     then error("random_hypergeometric: n1, n2, n must be positive integers")
1615   elseif maybe (n <= n1 +n2) = false
1616     then error("random_hypergeometric: n must be less or equal than n1 + n2")
1617     else block([no],
1618            if not integerp(n1) or not integerp(n2) or not integerp(n)
1619              then error("random_hypergeometric: need numeric argument for approximate procedure"),
1620            if length(num) = 0 then no: 0 else no: num[1],
1621            if integerp(no)
1622              then ?rndhypergeo(n1,n2,n,no)
1623              else error("random_hypergeometric: check sample size")) $
1624
1625
1626
1627 /*         NEGATIVE BINOMIAL DISTRIBUTION          */
1628
1629 /* R: dnbinom(x, n, p) */
1630 pdf_negative_binomial(x,n,p) :=
1631   if maybe(p > 0 and p <= 1) = false
1632     then error("pdf_negative_binomial: p must be a positive probability")
1633   elseif maybe(n>0)=false or numberp(n) and n-floor(n) > 0
1634     then error("pdf_negative_binomial: n must be a positive integer")
1635   elseif sign(p-1) = 'zero or sign(x) = 'neg or numberp(x) and x-floor(x) > 0
1636     then 0
1637     else gamma(n+x)*p^n*(1-p)^x/(x!*gamma(n)) $
1638
1639 /* R: pnbinom(x, n, p) */
1640 cdf_negative_binomial(x,n,p) :=
1641   if maybe(p > 0 and p <= 1) = false
1642     then error("cdf_negative_binomial: p must be a positive probability")
1643   elseif maybe(n>0)=false or numberp(n) and n-floor(n) > 0
1644     then error("cdf_negative_binomial: n must be a positive integer")
1645   elseif sign(x) = 'neg
1646     then 0
1647     else beta_incomplete_regularized(n,floor(x)+1,p) $
1648
1649 /* R: qnbinom(q,n,p) */
1650 quantile_negative_binomial(q,n,p) :=
1651   if maybe(p > 0 and p <= 1) = false
1652     then error("quantile_negative_binomial: p must be a positive probability")
1653   elseif maybe(n>0)=false or numberp(n) and n-floor(n) > 0
1654     then error("quantile_negative_binomial: n must be a positive integer")
1655   elseif maybe(q >= 0 and q <= 1) = false
1656      then error("quantile_negative_binomial: q must be a probability")
1657   elseif member(sign(q-p^n), ['zero, 'neg, 'nz])
1658      then 0
1659   elseif sign(1-q) = 'zero
1660      then 'inf
1661   elseif numberp(float(q)) and numberp(float(n)) and numberp(float(p))
1662     then /* partition method */
1663          block([a, b: 1.0, m, fq: float(q)],
1664            while (float(beta_incomplete_regularized(n,floor(b)+1,p)) < fq) do b: b*2,
1665            a: b/2,
1666            while (b-a > 1) do(
1667               m: floor(0.5*(a+b)),
1668               if float(beta_incomplete_regularized(n,m+1,p)) < fq
1669                 then a: m
1670                 else b: m ),
1671            b)
1672     else error("quantile_negative_binomial: need numeric arguments for approximate procedure") $
1673
1674 mean_negative_binomial(n,p) :=
1675   if maybe(p > 0 and p <= 1) = false
1676     then error("mean_negative_binomial: p must be a positive probability")
1677   elseif maybe(n>0)=false or numberp(n) and n-floor(n) > 0
1678     then error("mean_negative_binomial: n must be a positive integer")
1679     else n*(1-p)/p $
1680
1681 var_negative_binomial(n,p) :=
1682   if maybe(p > 0 and p <= 1) = false
1683     then error("var_negative_binomial: p must be a positive probability")
1684   elseif maybe(n>0)=false or numberp(n) and n-floor(n) > 0
1685     then error("var_negative_binomial: n must be a positive integer")
1686     else n*(1-p)/p^2 $
1687
1688 std_negative_binomial(n,p) :=
1689   if maybe(p > 0 and p <= 1) = false
1690     then error("std_negative_binomial: p must be a positive probability")
1691   elseif maybe(n>0)=false or numberp(n) and n-floor(n) > 0
1692     then error("std_negative_binomial: n must be a positive integer")
1693     else sqrt(n*(1-p))/p $
1694
1695 skewness_negative_binomial(n,p) :=
1696   if maybe(p > 0 and p <= 1) = false
1697     then error("skewness_negative_binomial: p must be a positive probability")
1698   elseif maybe(n>0)=false or numberp(n) and n-floor(n) > 0
1699     then error("skewness_negative_binomial: n must be a positive integer")
1700     else (2-p)/sqrt(n*(1-p)) $
1701
1702 kurtosis_negative_binomial(n,p) :=
1703   if maybe(p > 0 and p <= 1) = false
1704     then error("kurtosis_negative_binomial: p must be a positive probability")
1705   elseif maybe(n>0)=false or numberp(n) and n-floor(n) > 0
1706     then error("kurtosis_negative_binomial: n must be a positive integer")
1707     else (p*p+6-6*p)/(n*(1-p)) $
1708
1709 random_negative_binomial(n,p,[num]) :=
1710   if maybe(p > 0 and p <= 1) = false
1711     then error("random_negative_binomial: p must be a positive probability")
1712   elseif maybe(n>0)=false or numberp(n) and n-floor(n) > 0
1713     then error("random_negative_binomial: n must be a positive integer")
1714     else block([no, fp:float(p), fn:float(n)],
1715            if not numberp(fp) or not numberp(fn)
1716              then error("random_negative_binomial: need numeric argument for approximate procedure"),
1717            if length(num) = 0 then no: 0 else no: num[1],
1718            if integerp(no)
1719              then ?rndnegbinom(fn,fp,no)
1720              else error("random_negative_binomial: check sample size")) $
1721
1722
1723
1724 /*         GENERAL FINITE DISCRETE MODEL          */
1725
1726 pdf_general_finite_discrete(x,v) :=
1727   if not listp(v) or length(v)=0 or every(lambda([z], maybe(z >= 0)), v) = false
1728     then error("pdf_general_finite_discrete: v must be a list of non negative expressions")
1729   elseif not numberp(float(x))
1730     then error("pdf_general_finite_discrete: x must be a number")
1731   elseif member(sign(x), ['neg, 'zero, 'nz]) or sign(length(v)-x) = 'neg or x-floor(x) > 0
1732     then 0
1733     else v[floor(x)] / sum(v[k],k,1,length(v)) $
1734
1735 cdf_general_finite_discrete(x,v) :=
1736   if not listp(v) or length(v)=0 or every(lambda([z], maybe(z >= 0)), v) = false
1737     then error("cdf_general_finite_discrete: v must be a list of non negative expressions")
1738   elseif not numberp(float(x))
1739     then error("cdf_general_finite_discrete: x must be a number")
1740   elseif sign(x-1) = 'neg
1741     then 0
1742   elseif member(sign(x - length(v)), ['pos,'zero, 'pz])
1743     then 1
1744     else sum(v[k],k,1,floor(x)) / sum(v[k],k,1,length(v))$
1745
1746 quantile_general_finite_discrete(q,v) :=
1747   if not listp(v) or length(v)=0 or
1748      every(lambda([z], numberp(float(z)) and maybe(z >= 0)), v) = false
1749     then error("cdf_general_finite_discrete: v must be a list of non negative numbers")
1750   elseif maybe(q >= 0 and q <= 1) = false
1751     then error("quantile_general_finite_discrete: q must be a probability")
1752   elseif sign(q-1) = 'zero
1753     then length(v)
1754   elseif sign(q) = 'zero
1755     then 1
1756     else block([s:0, p, k:1],
1757            p: makelist(s:s+i, i, v/apply("+",v)),
1758            while (q>p[k]) do k: k+1,
1759            k )  $
1760
1761 mean_general_finite_discrete(v) :=
1762   if not listp(v) or length(v)=0 or every(lambda([z], maybe(z >= 0)), v) = false
1763     then error("mean_general_finite_discrete: v must be a list of non negative expressions")
1764     else block([p],
1765           p: v / apply("+", v),
1766           makelist(k,k,1,length(v)) . p ) $
1767
1768 var_general_finite_discrete(v) :=
1769   if not listp(v) or length(v)=0 or every(lambda([z], maybe(z >= 0)), v) = false
1770     then error("var_general_finite_discrete: v must be a list of non negative expressions")
1771     else block([p,m],
1772            p: v / apply("+", v),
1773            m: makelist(k,k,1,length(v)) . p,
1774            (makelist(k,k,1,length(v)) - m)^2 . p ) $
1775
1776 std_general_finite_discrete(v) :=
1777   if not listp(v) or length(v)=0 or every(lambda([z], maybe(z >= 0)), v) = false
1778     then error("std_general_finite_discrete: v must be a list of non negative expressions")
1779     else block([p,m],
1780            p: v / apply("+", v),
1781            m: makelist(k,k,1,length(v)) . p,
1782            sqrt((makelist(k,k,1,length(v)) - m)^2 . p)) $
1783
1784 skewness_general_finite_discrete(v) :=
1785   if not listp(v) or length(v)=0 or every(lambda([z], maybe(z >= 0)), v) = false
1786     then error("skewness_general_finite_discrete: v must be a list of non negative expressions")
1787     else block([p,m],
1788            p: v / apply("+", v),
1789            m: makelist(k,k,1,length(v)) . p,
1790            (makelist(k,k,1,length(v)) - m)^3 . p / var_discrete_model(v)^(3/2)) $
1791
1792 kurtosis_general_finite_discrete(v) :=
1793   if not listp(v) or length(v)=0 or every(lambda([z], maybe(z >= 0)), v) = false
1794     then error("kurtosis_general_finite_discrete: v must be a list of non negative expressions")
1795     else block([p,m],
1796            p: v / apply("+", v),
1797            m: makelist(k,k,1,length(v)) . p,
1798            (makelist(k,k,1,length(v)) - m)^4 . p / var_discrete_model(v)^2 - 3) $
1799
1800 random_general_finite_discrete(v,[num]) :=
1801   if not listp(v) or length(v)=0 or
1802      every(lambda([z], numberp(float(z)) and maybe(z >= 0)), v) = false
1803     then error("random_general_finite_discrete: v must be a list of non negative numbers")
1804     else (if length(num) = 0 then no: 0 else no: num[1],
1805           if integerp(no)
1806             then block([fv: float(v), s: 0, p, k, r],
1807                    fv: fv / apply("+", fv),
1808                    p: makelist(s:s+k, k, fv),
1809                    if no = 0
1810                      then
1811                        (r: random(1.0),
1812                         k: 1,
1813                         while (r > p[k]) do k: k+1,
1814                         k)
1815                      else
1816                        makelist((r: random(1.0),
1817                                  k: 1,
1818                                  while (r > p[k]) do k: k+1,
1819                                  k),
1820                                 i, no) )
1821             else error("random_general_finite_discrete: check sample size")) $
1822
1823
1824 /* INVERSE GAMMA DISTRIBUTION
1825  *
1826  * Formulas from: https://en.wikipedia.org/wiki/Inverse-gamma_distribution
1827  *
1828  * Two real parameters, a > 0 (shape), and b > 0 (scale).
1829  *
1830  * Support: x > 0
1831  */
1832
1833 pdf_inverse_gamma (x, a, b) :=
1834     if maybe (a > 0) = false or maybe (b > 0) = false
1835         then error ("pdf_inverse_gamma: parameters must be positive; found:", a, b)
1836         else b^a/gamma(a) * x^(-a-1) * exp(-b/x) * unit_step(x);
1837
1838 cdf_inverse_gamma (x, a, b) :=
1839     if maybe (a > 0) = false or maybe (b > 0) = false
1840         then error ("cdf_inverse_gamma: parameters must be positive; found:", a, b)
1841         else gamma_incomplete(a, b/x) / gamma(a) * unit_step(x);
1842
1843 quantile_inverse_gamma (x, a, b) :=
1844     if maybe (a > 0) = false or maybe (b > 0) = false
1845         then error ("quantile_inverse_gamma: parameters must be positive; found:", a, b)
1846         else if equal(x, 0) then 0
1847                  elseif equal(x, 1) then inf
1848                  else 1/quantile_gamma(1 - x, a, 1/b);
1849
1850 mean_inverse_gamma (a, b) :=
1851     if maybe (a > 0) = false or maybe (b > 0) = false
1852         then error ("mean_inverse_gamma: parameters must be positive; found:", a, b)
1853     elseif maybe (a > 1) = false
1854         then error ("mean_inverse_gamma: mean undefined for shape parameter =", a)
1855         else b/(a - 1);
1856
1857 mode_inverse_gamma (a, b) :=
1858     if maybe (a > 0) = false or maybe (b > 0) = false
1859         then error ("mode_inverse_gamma: parameters must be positive; found:", a, b)
1860         else b/(a + 1);
1861
1862 var_inverse_gamma (a, b) :=
1863     if maybe (a > 0) = false or maybe (b > 0) = false
1864         then error ("var_inverse_gamma: parameters must be positive; found:", a, b)
1865     elseif maybe (a > 2) = false
1866         then error ("var_inverse_gamma: variance undefined for shape parameter =", a)
1867         else b^2/(a - 1)^2/(a - 2);
1868
1869 std_inverse_gamma (a, b) :=
1870     if maybe (a > 0) = false or maybe (b > 0) = false
1871         then error ("std_inverse_gamma: parameters must be positive; found:", a, b)
1872     elseif maybe (a > 2) = false
1873         then error ("std_inverse_gamma: standard deviation undefined for shape parameter =", a)
1874         else b/(a - 1)/sqrt(a - 2);
1875
1876 skewness_inverse_gamma (a, b) :=
1877     if maybe (a > 0) = false or maybe (b > 0) = false
1878         then error ("skewness_inverse_gamma: parameters must be positive; found:", a, b)
1879     elseif maybe (a > 3) = false
1880         then error ("skewness_inverse_gamma: skewness undefined for shape parameter =", a)
1881         else 4*sqrt(a - 2)/(a - 3);
1882
1883 kurtosis_inverse_gamma (a, b) :=
1884     if maybe (a > 0) = false or maybe (b > 0) = false
1885         then error ("kurtosis_inverse_gamma: parameters must be positive; found:", a, b)
1886     elseif maybe (a > 4) = false
1887         then error ("kurtosis_inverse_gamma: kurtosis undefined for shape parameter =", a)
1888         else 6*(5*a - 11)/(a - 3)/(a - 4);
1889
1890 random_inverse_gamma (a, b, [n]) :=
1891     if maybe (a > 0) = false or maybe (b > 0) = false
1892         then error ("random_inverse_gamma: parameters must be positive; found:", a, b)
1893         else block ([listarith: true],
1894                     1/random_gamma (a, 1/b, if n = [] then 1 else n[1]));