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