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