dists.lsp

   1 ;;; -*- mode: lisp -*-
   2 ;;; Copyright (c) 2005--2007, by A.J. Rossini <blindglobe@gmail.com>
   3 ;;; See COPYRIGHT file for any additional restrictions (BSD license).
   4 ;;; Since 1991, ANSI was finally finished.  Edited for ANSI Common Lisp.
   5
   6 ;;;; dists -- Lisp-Stat interface to basic probability distribution routines
   7 ;;;;
   8 ;;;; Copyright (c) 1991, by Luke Tierney. Permission is granted for
   9 ;;;; unrestricted use.
  10
  11 ;;;;
  12 ;;;; Package Setup
  13 ;;;;
  14
  15 (defpackage :lisp-stat-probability
  16     (:use :common-lisp
  17           :lisp-stat-ffi-int
  18           :lisp-stat-macros)
  19     (:export
  20      log-gamma
  21      uniform-rand
  22      normal-cdf normal-quant normal-dens normal-rand bivnorm-cdf
  23      cauchy-cdf cauchy-quant cauchy-dens cauchy-rand
  24      gamma-cdf gamma-quant gamma-dens gamma-rand
  25      chisq-cdf chisq-quant chisq-dens chisq-rand
  26      beta-cdf beta-quant beta-dens beta-rand
  27      t-cdf t-quant t-dens t-rand
  28      f-cdf f-quant f-dens f-rand
  29      poisson-cdf poisson-quant poisson-pmf poisson-rand
  30      binomial-cdf binomial-quant binomial-pmf binomial-rand))
  31
  32 (in-package :lisp-stat-probability)
  33
  34 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  35 ;;;
  36 ;;;  CFFI support for Probability Distributions
  37 ;;;
  38 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  39
  40 ;;;
  41 ;;; C-callable uniform generator
  42 ;;;
  43
  44 (cffi:defcfun ("register_uni" register-uni)
  45     :void (f :pointer))
  46 (cffi:defcallback ccl-uni :int () (ccl-store-double (random 1.0)) 0)
  47 (register-uni (cffi:callback ccl-uni))
  48
  49 (defun one-uniform-rand () (random 1.0))
  50
  51 ;;;;
  52 ;;;; Log-gamma function
  53 ;;;;
  54
  55 (cffi:defcfun ("ccl_gamma" ccl-base-log-gamma)
  56     :double (x :double))
  57 (defun base-log-gamma (x)
  58   (ccl-base-log-gamma (float x 1d0)))
  59
  60 ;;;;
  61 ;;;; Normal distribution
  62 ;;;;
  63
  64 (cffi:defcfun ("ccl_normalcdf" ccl-base-normal-cdf)
  65     :double (x :double))
  66 (defun base-normal-cdf (x)
  67   (ccl-base-normal-cdf (float x 1d0)))
  68
  69 (cffi:defcfun ("ccl_normalquant" ccl-base-normal-quant)
  70     :double (x :double))
  71 (defun base-normal-quant (x)
  72   (ccl-base-normal-quant (float x 1d0)))
  73
  74 (cffi:defcfun ("ccl_normaldens" ccl-base-normal-dens)
  75         :double (x :double))
  76 (defun base-normal-dens (x)
  77   (ccl-base-normal-dens (float x 1d0)))
  78
  79 (cffi:defcfun ("ccl_normalrand" one-normal-rand)
  80     :float)
  81
  82 (cffi:defcfun ("ccl_bnormcdf" ccl-base-bivnorm-cdf)
  83     :double (x :double) (y :double) (z :double))
  84 (defun base-bivnorm-cdf (x y z)
  85   (ccl-base-bivnorm-cdf (float x 1d0) (float y 1d0) (float z 1d0)))
  86
  87 ;;;;
  88 ;;;; Cauchy distribution
  89 ;;;;
  90
  91 (cffi:defcfun ("ccl_cauchycdf" ccl-base-cauchy-cdf)
  92             :double (x :double))
  93 (defun base-cauchy-cdf (x)
  94   (ccl-base-cauchy-cdf (float x 1d0)))
  95
  96 (cffi:defcfun ("ccl_cauchyquant" ccl-base-cauchy-quant)
  97             :double (x :double))
  98 (defun base-cauchy-quant (x)
  99   (ccl-base-cauchy-quant (float x 1d0)))
 100
 101 (cffi:defcfun ("ccl_cauchydens" ccl-base-cauchy-dens)
 102             :double (x :double))
 103 (defun base-cauchy-dens (x)
 104   (ccl-base-cauchy-dens (float x 1d0)))
 105
 106 (cffi:defcfun ("ccl_cauchyrand" one-cauchy-rand)
 107     :double)
 108
 109 ;;;;
 110 ;;;; Gamma distribution
 111 ;;;;
 112
 113 (cffi:defcfun ("ccl_gammacdf" ccl-base-gamma-cdf)
 114             :double (x :double) (y :double))
 115 (defun base-gamma-cdf (x y)
 116   (ccl-base-gamma-cdf (float x 1d0) (float y 1d0)))
 117
 118 (cffi:defcfun ("ccl_gammaquant" ccl-base-gamma-quant)
 119             :double (x :double) (y :double))
 120 (defun base-gamma-quant (x y)
 121   (ccl-base-gamma-quant (float x 1d0) (float y 1d0)))
 122
 123 (cffi:defcfun ("ccl_gammadens" ccl-base-gamma-dens)
 124             :double (x :double) (y :double))
 125 (defun base-gamma-dens (x y)
 126   (ccl-base-gamma-dens (float x 1d0) (float y 1d0)))
 127
 128 (cffi:defcfun ("ccl_gammarand" ccl-gamma-rand)
 129             :double (x :double))
 130 (defun one-gamma-rand (x)
 131   (ccl-gamma-rand (float x 1d0)))
 132
 133 ;;;;
 134 ;;;; Chi-square distribution
 135 ;;;;
 136
 137 (cffi:defcfun ("ccl_chisqcdf" ccl-base-chisq-cdf)
 138             :double (x :double) (y :double))
 139 (defun base-chisq-cdf (x y)
 140   (ccl-base-chisq-cdf (float x 1d0) (float y 1d0)))
 141
 142 (cffi:defcfun ("ccl_chisqquant" ccl-base-chisq-quant)
 143             :double (x :double) (y :double))
 144 (defun base-chisq-quant (x y)
 145   (ccl-base-chisq-quant (float x 1d0) (float y 1d0)))
 146
 147 (cffi:defcfun ("ccl_chisqdens" ccl-base-chisq-dens)
 148             :double (x :double) (y :double))
 149 (defun base-chisq-dens (x y)
 150   (ccl-base-chisq-dens (float x 1d0) (float y 1d0)))
 151
 152 (cffi:defcfun ("ccl_chisqrand" ccl-chisq-rand)
 153             :double (x :double))
 154 (defun one-chisq-rand (x)
 155   (ccl-chisq-rand (float x 1d0)))
 156
 157 ;;;;
 158 ;;;; Beta distribution
 159 ;;;;
 160
 161 (cffi:defcfun ("ccl_betacdf" ccl-base-beta-cdf)
 162     :double (x :double) (y :double) (z :double))
 163 (defun base-beta-cdf (x y z)
 164   (ccl-base-beta-cdf (float x 1d0) (float y 1d0) (float z 1d0)))
 165
 166 (cffi:defcfun ("ccl_betaquant" ccl-base-beta-quant)
 167     :double (x :double) (y :double) (z :double))
 168 (defun base-beta-quant (x y z)
 169   (ccl-base-beta-quant (float x 1d0) (float y 1d0) (float z 1d0)))
 170
 171 (cffi:defcfun ("ccl_betadens" ccl-base-beta-dens)
 172     :double (x :double) (y :double) (z :double))
 173 (defun base-beta-dens (x y z)
 174   (ccl-base-beta-dens (float x 1d0) (float y 1d0) (float z 1d0)))
 175
 176 (cffi:defcfun ("ccl_betarand" ccl-beta-rand)
 177     :double (x :double) (y :double))
 178 (defun one-beta-rand (x y)
 179   (ccl-beta-rand (float x 1d0) (float y 1d0)))
 180
 181 ;;;;
 182 ;;;; t distribution
 183 ;;;;
 184
 185 (cffi:defcfun ("ccl_tcdf" ccl-base-t-cdf)
 186     :double (x :double) (y :double))
 187 (defun base-t-cdf (x y)
 188   (ccl-base-t-cdf (float x 1d0) (float y 1d0)))
 189
 190 (cffi:defcfun ("ccl_tquant" ccl-base-t-quant)
 191     :double (x :double) (y :double))
 192 (defun base-t-quant (x y)
 193   (ccl-base-t-quant (float x 1d0) (float y 1d0)))
 194
 195 (cffi:defcfun ("ccl_tdens" ccl-base-t-dens)
 196     :double (x :double) (y :double))
 197 (defun base-t-dens (x y)
 198   (ccl-base-t-dens (float x 1d0) (float y 1d0)))
 199
 200 (cffi:defcfun ("ccl_trand" ccl-t-rand)
 201     :double (x :double))
 202 (defun one-t-rand (x)
 203   (ccl-t-rand (float x 1d0)))
 204
 205 ;;;;
 206 ;;;; F distribution
 207 ;;;;
 208
 209 (cffi:defcfun ("ccl_fcdf" ccl-base-f-cdf)
 210     :double (x :double) (y :double) (z :double))
 211 (defun base-f-cdf (x y z)
 212   (ccl-base-f-cdf (float x 1d0) (float y 1d0) (float z 1d0)))
 213
 214 (cffi:defcfun ("ccl_fquant" ccl-base-f-quant)
 215     :double (x :double) (y :double) (z :double))
 216 (defun base-f-quant (x y z)
 217   (ccl-base-f-quant (float x 1d0) (float y 1d0) (float z 1d0)))
 218
 219 (cffi:defcfun ("ccl_fdens" ccl-base-f-dens)
 220     :double (x :double) (y :double) (z :double))
 221 (defun base-f-dens (x y z)
 222   (ccl-base-f-dens (float x 1d0) (float y 1d0) (float z 1d0)))
 223
 224 (cffi:defcfun ("ccl_frand" ccl-f-rand)
 225     :double (x :double) (y :double))
 226 (defun one-f-rand (x y) (ccl-f-rand (float x 1d0) (float y 1d0)))
 227
 228 ;;;;
 229 ;;;; Poisson distribution
 230 ;;;;
 231
 232 (cffi:defcfun ("ccl_poissoncdf" ccl-base-poisson-cdf)
 233     :double (x :double) (y :double))
 234 (defun base-poisson-cdf (x y)
 235   (ccl-base-poisson-cdf (float x 1d0) (float y 1d0)))
 236
 237 (cffi:defcfun ("ccl_poissonquant" ccl-base-poisson-quant)
 238     :int (x :double) (y :double))
 239 (defun base-poisson-quant (x y)
 240   (ccl-base-poisson-quant (float x 1d0) (float y 1d0)))
 241
 242 (cffi:defcfun ("ccl_poissonpmf" ccl-base-poisson-pmf)
 243     :double (x :int) (y :double))
 244 (defun base-poisson-pmf (x y)
 245   (ccl-base-poisson-pmf x (float y 1d0)))
 246
 247 (cffi:defcfun ("ccl_poissonrand" ccl-poisson-rand)
 248     :int (x :double))
 249 (defun one-poisson-rand (x)
 250   (ccl-poisson-rand (float x 1d0)))
 251
 252 ;;;;
 253 ;;;; Binomial distribution
 254 ;;;;
 255
 256 (cffi:defcfun ("ccl_binomialcdf" ccl-base-binomial-cdf)
 257     :double (x :double) (y :int) (z :double))
 258 (defun base-binomial-cdf (x y z)
 259   (ccl-base-binomial-cdf (float x 1d0) y (float z 1d0)))
 260
 261 (cffi:defcfun ("ccl_binomialquant" ccl-base-binomial-quant)
 262     :int (x :double) (y :int) (z :double))
 263 (defun base-binomial-quant (x y z)
 264   (ccl-base-binomial-quant (float x 1d0) y (float z 1d0)))
 265
 266 (cffi:defcfun ("ccl_binomialpmf" ccl-base-binomial-pmf)
 267     :double (x :int) (y :int) (z :double))
 268 (defun base-binomial-pmf (x y z)
 269   (ccl-base-binomial-pmf x y (float z 1d0)))
 270
 271 (cffi:defcfun ("ccl_binomialrand" ccl-binomial-rand)
 272     :int (x :int) (y :double))
 273 (defun one-binomial-rand (x y)
 274   (ccl-binomial-rand x (float y 1d0)))
 275
 276
 277
 278
 279
 280 ;;; definitions though macros
 281
 282 (defmacro defbaserand (name onefun &rest args)
 283   `(defun ,name (n ,@args)
 284      (let ((result nil))
 285        (dotimes (i n result)
 286                 (declare (fixnum i) (inline ,onefun))
 287                 (setf result (cons (,onefun ,@args) result))))))
 288
 289 (defbaserand base-uniform-rand one-uniform-rand)
 290 (defbaserand base-normal-rand one-normal-rand)
 291 (defbaserand base-cauchy-rand one-cauchy-rand)
 292 (defbaserand base-gamma-rand one-gamma-rand a)
 293 (defbaserand base-chisq-rand one-chisq-rand df)
 294 (defbaserand base-beta-rand one-beta-rand a b)
 295 (defbaserand base-t-rand one-t-rand df)
 296 (defbaserand base-f-rand one-f-rand ndf ddf)
 297 (defbaserand base-poisson-rand one-poisson-rand a)
 298 (defbaserand base-binomial-rand one-binomial-rand a b)
 299
 300 (make-rv-function log-gamma base-log-gamma x)
 301
 302 (make-rv-function uniform-rand base-uniform-rand n)
 303
 304 (make-rv-function normal-cdf base-normal-cdf x)
 305 (make-rv-function normal-quant base-normal-quant p)
 306 (make-rv-function normal-dens base-normal-dens x)
 307 (make-rv-function normal-rand base-normal-rand n)
 308 (make-rv-function bivnorm-cdf base-bivnorm-cdf x y r)
 309
 310 (make-rv-function cauchy-cdf base-cauchy-cdf x)
 311 (make-rv-function cauchy-quant base-cauchy-quant p)
 312 (make-rv-function cauchy-dens base-cauchy-dens x)
 313 (make-rv-function cauchy-rand base-cauchy-rand n)
 314
 315 (make-rv-function gamma-cdf base-gamma-cdf x a)
 316 (make-rv-function gamma-quant base-gamma-quant p a)
 317 (make-rv-function gamma-dens base-gamma-dens x a)
 318 (make-rv-function gamma-rand base-gamma-rand n a)
 319
 320 (make-rv-function chisq-cdf base-chisq-cdf x df)
 321 (make-rv-function chisq-quant base-chisq-quant p df)
 322 (make-rv-function chisq-dens base-chisq-dens x df)
 323 (make-rv-function chisq-rand base-chisq-rand n df)
 324
 325 (make-rv-function beta-cdf base-beta-cdf x a b)
 326 (make-rv-function beta-quant base-beta-quant p a b)
 327 (make-rv-function beta-dens base-beta-dens x a b)
 328 (make-rv-function beta-rand base-beta-rand n a b)
 329
 330 (make-rv-function t-cdf base-t-cdf x df)
 331 (make-rv-function t-quant base-t-quant p df)
 332 (make-rv-function t-dens base-t-dens x df)
 333 (make-rv-function t-rand base-t-rand n df)
 334
 335 (make-rv-function f-cdf base-f-cdf x ndf ddf)
 336 (make-rv-function f-quant base-f-quant p ndf ddf)
 337 (make-rv-function f-dens base-f-dens x ndf ddf)
 338 (make-rv-function f-rand base-f-rand n ndf ddf)
 339
 340 (make-rv-function poisson-cdf base-poisson-cdf x a)
 341 (make-rv-function poisson-quant base-poisson-quant p a)
 342 (make-rv-function poisson-pmf base-poisson-pmf x a)
 343 (make-rv-function poisson-rand base-poisson-rand n a)
 344
 345 (make-rv-function binomial-cdf base-binomial-cdf x a b)
 346 (make-rv-function binomial-quant base-binomial-quant p a b)
 347 (make-rv-function binomial-pmf base-binomial-pmf x a b)
 348 (make-rv-function binomial-rand base-binomial-rand n a b)
 349
 350 ;;;;
 351 ;;;; Documentation
 352 ;;;;
 353
 354 (setf (documentation 'bivnorm-cdf 'function)
 355 "Args: (x y r)
 356 Returns the value of the standard bivariate normal distribution function
 357 with correlation R at (X, Y). Vectorized.")
 358
 359 (setf (documentation 'normal-cdf 'function)
 360 "Args: (x)
 361 Returns the value of the standard normal distribution function at X.
 362 Vectorized.")
 363
 364 (setf (documentation 'beta-cdf 'function)
 365 "Args: (x alpha beta)
 366 Returns the value of the Beta(ALPHA, BETA) distribution function at X.
 367 Vectorized.")
 368
 369 (setf (documentation 'gamma-cdf 'function)
 370 "Args: (x alpha)
 371 Returns the value of the Gamma(alpha, 1) distribution function at X.
 372 Vectorized.")
 373
 374 (setf (documentation 'chisq-cdf 'function)
 375 "Args: (x df)
 376 Returns the value of the Chi-Square(DF) distribution function at X. Vectorized.")
 377
 378 (setf (documentation 't-cdf 'function)
 379 "Args: (x df)
 380 Returns the value of the T(DF) distribution function at X. Vectorized.")
 381
 382 (setf (documentation 'f-cdf 'function)
 383 "Args: (x ndf ddf)
 384 Returns the value of the F(NDF, DDF) distribution function at X. Vectorized.")
 385
 386 (setf (documentation 'cauchy-cdf 'function)
 387 "Args: (x)
 388 Returns the value of the standard Cauchy distribution function at X.
 389 Vectorized.")
 390
 391 (setf (documentation 'log-gamma 'function)
 392 "Args: (x)
 393 Returns the log gamma function of X. Vectorized.")
 394
 395 (setf (documentation 'normal-quant 'function)
 396 "Args (p)
 397 Returns the P-th quantile of the standard normal distribution. Vectorized.")
 398
 399 (setf (documentation 'cauchy-quant 'function)
 400 "Args (p)
 401 Returns the P-th quantile(s) of the standard Cauchy distribution. Vectorized.")
 402
 403 (setf (documentation 'beta-quant 'function)
 404 "Args: (p alpha beta)
 405 Returns the P-th quantile of the Beta(ALPHA, BETA) distribution. Vectorized.")
 406
 407 (setf (documentation 'gamma-quant 'function)
 408 "Args: (p alpha)
 409 Returns the P-th quantile of the Gamma(ALPHA, 1) distribution. Vectorized.")
 410
 411 (setf (documentation 'chisq-quant 'function)
 412 "Args: (p df)
 413 Returns the P-th quantile of the Chi-Square(DF) distribution. Vectorized.")
 414
 415 (setf (documentation 't-quant 'function)
 416 "Args: (p df)
 417 Returns the P-th quantile of the T(DF) distribution. Vectorized.")
 418
 419 (setf (documentation 'f-quant 'function)
 420 "Args: (p ndf ddf)
 421 Returns the P-th quantile of the F(NDF, DDF) distribution. Vectorized.")
 422
 423 (setf (documentation 'normal-dens 'function)
 424 "Args: (x)
 425 Returns the density at X of the standard normal distribution. Vectorized.")
 426
 427 (setf (documentation 'cauchy-dens 'function)
 428 "Args: (x)
 429 Returns the density at X of the standard Cauchy distribution. Vectorized.")
 430
 431 (setf (documentation 'beta-dens 'function)
 432 "Args: (x alpha beta)
 433 Returns the density at X of the Beta(ALPHA, BETA) distribution. Vectorized.")
 434
 435 (setf (documentation 'gamma-dens 'function)
 436 "Args: (x alpha)
 437 Returns the density at X of the Gamma(ALPHA, 1) distribution. Vectorized.")
 438
 439 (setf (documentation 'chisq-dens 'function)
 440 "Args: (x alpha)
 441 Returns the density at X of the Chi-Square(DF) distribution. Vectorized.")
 442
 443 (setf (documentation 't-dens 'function)
 444 "Args: (x alpha)
 445 Returns the density at X of the T(DF) distribution. Vectorized.")
 446
 447 (setf (documentation 'f-dens 'function)
 448 "Args: (x ndf ddf)
 449 Returns the density at X of the F(NDF, DDF) distribution. Vectorized.")
 450
 451 (setf (documentation 'uniform-rand 'function)
 452 "Args: (n)
 453 Returns a list of N uniform random variables from the range (0, 1).
 454 Vectorized.")
 455
 456 (setf (documentation 'normal-rand 'function)
 457 "Args: (n)
 458 Returns a list of N standard normal random numbers. Vectorized.")
 459
 460 (setf (documentation 'cauchy-rand 'function)
 461 "Args: (n)
 462 Returns a list of N standard Cauchy random numbers. Vectorized.")
 463
 464 (setf (documentation 't-rand 'function)
 465 "Args: (n df)
 466 Returns a list of N T(DF) random variables. Vectorized.")
 467
 468 (setf (documentation 'f-rand 'function)
 469 "Args: (n ndf ddf)
 470 Returns a list of N F(NDF, DDF) random variables. Vectorized.")
 471
 472 (setf (documentation 'gamma-rand 'function)
 473 "Args: (n a)
 474 Returns a list of N Gamma(A, 1) random variables. Vectorized.")
 475
 476 (setf (documentation 'chisq-rand 'function)
 477 "Args: (n df)
 478 Returns a list of N Chi-Square(DF) random variables. Vectorized.")
 479
 480 (setf (documentation 'beta-rand 'function)
 481 "Args: (n a b)
 482 Returns a list of N beta(A, B) random variables. Vectorized.")
 483
 484 (setf (documentation 'binomial-cdf 'function)
 485 "Args (x n p)
 486 Returns value of the Binomial(N, P) distribution function at X. Vectorized.")
 487
 488 (setf (documentation 'poisson-cdf 'function)
 489 "Args (x mu)
 490 Returns value of the Poisson(MU) distribution function at X. Vectorized.")
 491
 492 (setf (documentation 'binomial-pmf 'function)
 493 "Args (k n p)
 494 Returns value of the Binomial(N, P) pmf function at integer K. Vectorized.")
 495
 496 (setf (documentation 'poisson-pmf 'function)
 497 "Args (k mu)
 498 Returns value of the Poisson(MU) pmf function at integer K. Vectorized.")
 499
 500 (setf (documentation 'binomial-quant 'function)
 501 "Args: (x n p)
 502 Returns x-th quantile (left continuous inverse) of Binomial(N, P) cdf.
 503 Vectorized.")
 504
 505 (setf (documentation 'poisson-quant 'function)
 506 "Args: (x mu)
 507 Returns x-th quantile (left continuous inverse) of Poisson(MU) cdf.
 508 Vectorized.")
 509
 510 (setf (documentation 'binomial-rand 'function)
 511 "Args: (k n p)
 512 Returns list of K draws from the Binomial(N, P) distribution. Vectorized.")
 513
 514 (setf (documentation 'poisson-rand 'function)
 515 "Args: (k mu)
 516 Returns list of K draws from the Poisson(MU) distribution. Vectorized.")