blob 32a466424731d37f76ea76b8b1bba0c5771da93a
1 ;;; -*- mode: lisp -*-
2 ;;; Copyright (c) 2005--2007, by A.J. Rossini <blindglobe@gmail.com>
4 ;;; Since 1991, ANSI was finally finished. Edited for ANSI Common Lisp.
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.
11 ;;; This stuff needs to be improved. We really could use something
12 ;;; like the R libraries, but of course in a better packaged manner.
14 ;;; Currently, there is a function for everything. Probably better to
15 ;;; simplify by thinking about a more generic approach, distribution
16 ;;; being specified by keyword.
18 ;;;
19 ;;; Package Setup
20 ;;;
22 (in-package :cl-user)
24 (defpackage :lisp-stat-probability
25 (:use :common-lisp
26 :cffi
27 :lisp-stat-ffi-int
28 :lisp-stat-macros)
29 (:export log-gamma set-seed
30 uniform-rand
31 normal-cdf normal-quant normal-dens normal-rand
32 bivnorm-cdf
33 cauchy-cdf cauchy-quant cauchy-dens cauchy-rand
34 gamma-cdf gamma-quant gamma-dens gamma-rand
35 chisq-cdf chisq-quant chisq-dens chisq-rand
36 beta-cdf beta-quant beta-dens beta-rand
37 t-cdf t-quant t-dens t-rand
38 f-cdf f-quant f-dens f-rand
39 poisson-cdf poisson-quant poisson-pmf poisson-rand
40 binomial-cdf binomial-quant binomial-pmf binomial-rand))
42 (in-package :lisp-stat-probability)
45 (defun set-seed (x)
46 "stupid dummy function, need to implement rng seeding tool."
47 (values x))
49 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
50 ;;;
51 ;;; CFFI support for Probability Distributions
52 ;;;
53 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
55 ;;;
56 ;;; C-callable uniform generator
57 ;;;
59 (defcfun ("register_uni" register-uni)
60 :void (f :pointer))
61 (defcallback ccl-uni :int () (ccl-store-double (random 1.0)) 0)
62 (register-uni (callback ccl-uni))
64 (defun one-uniform-rand () (random 1.0))
66 ;;;
67 ;;; Log-gamma function
68 ;;;
70 (defcfun ("ccl_gamma" ccl-base-log-gamma)
71 :double (x :double))
72 (defun base-log-gamma (x)
73 (ccl-base-log-gamma (float x 1d0)))
75 ;;;
76 ;;; Normal distribution
77 ;;;
79 (defcfun ("ccl_normalcdf" ccl-base-normal-cdf)
80 :double (x :double))
81 (defun base-normal-cdf (x)
82 (ccl-base-normal-cdf (float x 1d0)))
84 (defcfun ("ccl_normalquant" ccl-base-normal-quant)
85 :double (x :double))
86 (defun base-normal-quant (x)
87 (ccl-base-normal-quant (float x 1d0)))
89 (defcfun ("ccl_normaldens" ccl-base-normal-dens)
90 :double (x :double))
91 (defun base-normal-dens (x)
92 (ccl-base-normal-dens (float x 1d0)))
94 (defcfun ("ccl_normalrand" one-normal-rand)
95 :float)
97 (defcfun ("ccl_bnormcdf" ccl-base-bivnorm-cdf)
98 :double (x :double) (y :double) (z :double))
99 (defun base-bivnorm-cdf (x y z)
100 (ccl-base-bivnorm-cdf (float x 1d0) (float y 1d0) (float z 1d0)))
103 ;;; Cauchy distribution
106 (defcfun ("ccl_cauchycdf" ccl-base-cauchy-cdf)
107 :double (x :double))
108 (defun base-cauchy-cdf (x)
109 (ccl-base-cauchy-cdf (float x 1d0)))
111 (defcfun ("ccl_cauchyquant" ccl-base-cauchy-quant)
112 :double (x :double))
113 (defun base-cauchy-quant (x)
114 (ccl-base-cauchy-quant (float x 1d0)))
116 (defcfun ("ccl_cauchydens" ccl-base-cauchy-dens)
117 :double (x :double))
118 (defun base-cauchy-dens (x)
119 (ccl-base-cauchy-dens (float x 1d0)))
121 (defcfun ("ccl_cauchyrand" one-cauchy-rand)
122 :double)
124 ;;;;
125 ;;;; Gamma distribution
126 ;;;;
128 (defcfun ("ccl_gammacdf" ccl-base-gamma-cdf)
129 :double (x :double) (y :double))
130 (defun base-gamma-cdf (x y)
131 (ccl-base-gamma-cdf (float x 1d0) (float y 1d0)))
133 (defcfun ("ccl_gammaquant" ccl-base-gamma-quant)
134 :double (x :double) (y :double))
135 (defun base-gamma-quant (x y)
136 (ccl-base-gamma-quant (float x 1d0) (float y 1d0)))
139 :double (x :double) (y :double))
140 (defun base-gamma-dens (x y)
141 (ccl-base-gamma-dens (float x 1d0) (float y 1d0)))
143 (defcfun ("ccl_gammarand" ccl-gamma-rand)
144 :double (x :double))
145 (defun one-gamma-rand (x)
146 (ccl-gamma-rand (float x 1d0)))
148 ;;;;
149 ;;;; Chi-square distribution
150 ;;;;
152 (defcfun ("ccl_chisqcdf" ccl-base-chisq-cdf)
153 :double (x :double) (y :double))
154 (defun base-chisq-cdf (x y)
155 (ccl-base-chisq-cdf (float x 1d0) (float y 1d0)))
157 (defcfun ("ccl_chisqquant" ccl-base-chisq-quant)
158 :double (x :double) (y :double))
159 (defun base-chisq-quant (x y)
160 (ccl-base-chisq-quant (float x 1d0) (float y 1d0)))
162 (defcfun ("ccl_chisqdens" ccl-base-chisq-dens)
163 :double (x :double) (y :double))
164 (defun base-chisq-dens (x y)
165 (ccl-base-chisq-dens (float x 1d0) (float y 1d0)))
167 (defcfun ("ccl_chisqrand" ccl-chisq-rand)
168 :double (x :double))
169 (defun one-chisq-rand (x)
170 (ccl-chisq-rand (float x 1d0)))
172 ;;;;
173 ;;;; Beta distribution
174 ;;;;
176 (defcfun ("ccl_betacdf" ccl-base-beta-cdf)
177 :double (x :double) (y :double) (z :double))
178 (defun base-beta-cdf (x y z)
179 (ccl-base-beta-cdf (float x 1d0) (float y 1d0) (float z 1d0)))
181 (defcfun ("ccl_betaquant" ccl-base-beta-quant)
182 :double (x :double) (y :double) (z :double))
183 (defun base-beta-quant (x y z)
184 (ccl-base-beta-quant (float x 1d0) (float y 1d0) (float z 1d0)))
187 :double (x :double) (y :double) (z :double))
188 (defun base-beta-dens (x y z)
189 (ccl-base-beta-dens (float x 1d0) (float y 1d0) (float z 1d0)))
191 (defcfun ("ccl_betarand" ccl-beta-rand)
192 :double (x :double) (y :double))
193 (defun one-beta-rand (x y)
194 (ccl-beta-rand (float x 1d0) (float y 1d0)))
196 ;;;;
197 ;;;; t distribution
198 ;;;;
200 (defcfun ("ccl_tcdf" ccl-base-t-cdf)
201 :double (x :double) (y :double))
202 (defun base-t-cdf (x y)
203 (ccl-base-t-cdf (float x 1d0) (float y 1d0)))
205 (defcfun ("ccl_tquant" ccl-base-t-quant)
206 :double (x :double) (y :double))
207 (defun base-t-quant (x y)
208 (ccl-base-t-quant (float x 1d0) (float y 1d0)))
210 (defcfun ("ccl_tdens" ccl-base-t-dens)
211 :double (x :double) (y :double))
212 (defun base-t-dens (x y)
213 (ccl-base-t-dens (float x 1d0) (float y 1d0)))
215 (defcfun ("ccl_trand" ccl-t-rand)
216 :double (x :double))
217 (defun one-t-rand (x)
218 (ccl-t-rand (float x 1d0)))
220 ;;;;
221 ;;;; F distribution
222 ;;;;
224 (defcfun ("ccl_fcdf" ccl-base-f-cdf)
225 :double (x :double) (y :double) (z :double))
226 (defun base-f-cdf (x y z)
227 (ccl-base-f-cdf (float x 1d0) (float y 1d0) (float z 1d0)))
229 (defcfun ("ccl_fquant" ccl-base-f-quant)
230 :double (x :double) (y :double) (z :double))
231 (defun base-f-quant (x y z)
232 (ccl-base-f-quant (float x 1d0) (float y 1d0) (float z 1d0)))
234 (defcfun ("ccl_fdens" ccl-base-f-dens)
235 :double (x :double) (y :double) (z :double))
236 (defun base-f-dens (x y z)
237 (ccl-base-f-dens (float x 1d0) (float y 1d0) (float z 1d0)))
239 (defcfun ("ccl_frand" ccl-f-rand)
240 :double (x :double) (y :double))
241 (defun one-f-rand (x y) (ccl-f-rand (float x 1d0) (float y 1d0)))
243 ;;;;
244 ;;;; Poisson distribution
245 ;;;;
247 (defcfun ("ccl_poissoncdf" ccl-base-poisson-cdf)
248 :double (x :double) (y :double))
249 (defun base-poisson-cdf (x y)
250 (ccl-base-poisson-cdf (float x 1d0) (float y 1d0)))
252 (defcfun ("ccl_poissonquant" ccl-base-poisson-quant)
253 :int (x :double) (y :double))
254 (defun base-poisson-quant (x y)
255 (ccl-base-poisson-quant (float x 1d0) (float y 1d0)))
257 (defcfun ("ccl_poissonpmf" ccl-base-poisson-pmf)
258 :double (x :int) (y :double))
259 (defun base-poisson-pmf (x y)
260 (ccl-base-poisson-pmf x (float y 1d0)))
262 (defcfun ("ccl_poissonrand" ccl-poisson-rand)
263 :int (x :double))
264 (defun one-poisson-rand (x)
265 (ccl-poisson-rand (float x 1d0)))
267 ;;;;
268 ;;;; Binomial distribution
269 ;;;;
271 (defcfun ("ccl_binomialcdf" ccl-base-binomial-cdf)
272 :double (x :double) (y :int) (z :double))
273 (defun base-binomial-cdf (x y z)
274 (ccl-base-binomial-cdf (float x 1d0) y (float z 1d0)))
276 (defcfun ("ccl_binomialquant" ccl-base-binomial-quant)
277 :int (x :double) (y :int) (z :double))
278 (defun base-binomial-quant (x y z)
279 (ccl-base-binomial-quant (float x 1d0) y (float z 1d0)))
281 (defcfun ("ccl_binomialpmf" ccl-base-binomial-pmf)
282 :double (x :int) (y :int) (z :double))
283 (defun base-binomial-pmf (x y z)
284 (ccl-base-binomial-pmf x y (float z 1d0)))
286 (defcfun ("ccl_binomialrand" ccl-binomial-rand)
287 :int (x :int) (y :double))
288 (defun one-binomial-rand (x y)
289 (ccl-binomial-rand x (float y 1d0)))
295 ;;; definitions though macros
297 (defmacro defbaserand (name onefun &rest args)
298 `(defun ,name (n ,@args)
299 (let ((result nil))
300 (dotimes (i n result)
301 (declare (fixnum i) (inline ,onefun))
302 (setf result (cons (,onefun ,@args) result))))))
304 (defbaserand base-uniform-rand one-uniform-rand)
305 (defbaserand base-normal-rand one-normal-rand)
306 (defbaserand base-cauchy-rand one-cauchy-rand)
307 (defbaserand base-gamma-rand one-gamma-rand a)
308 (defbaserand base-chisq-rand one-chisq-rand df)
309 (defbaserand base-beta-rand one-beta-rand a b)
310 (defbaserand base-t-rand one-t-rand df)
311 (defbaserand base-f-rand one-f-rand ndf ddf)
312 (defbaserand base-poisson-rand one-poisson-rand a)
313 (defbaserand base-binomial-rand one-binomial-rand a b)
315 (make-rv-function log-gamma base-log-gamma x)
317 (make-rv-function uniform-rand base-uniform-rand n)
319 (make-rv-function normal-cdf base-normal-cdf x)
320 (make-rv-function normal-quant base-normal-quant p)
321 (make-rv-function normal-dens base-normal-dens x)
322 (make-rv-function normal-rand base-normal-rand n)
323 (make-rv-function bivnorm-cdf base-bivnorm-cdf x y r)
325 (make-rv-function cauchy-cdf base-cauchy-cdf x)
326 (make-rv-function cauchy-quant base-cauchy-quant p)
327 (make-rv-function cauchy-dens base-cauchy-dens x)
328 (make-rv-function cauchy-rand base-cauchy-rand n)
330 (make-rv-function gamma-cdf base-gamma-cdf x a)
331 (make-rv-function gamma-quant base-gamma-quant p a)
332 (make-rv-function gamma-dens base-gamma-dens x a)
333 (make-rv-function gamma-rand base-gamma-rand n a)
335 (make-rv-function chisq-cdf base-chisq-cdf x df)
336 (make-rv-function chisq-quant base-chisq-quant p df)
337 (make-rv-function chisq-dens base-chisq-dens x df)
338 (make-rv-function chisq-rand base-chisq-rand n df)
340 (make-rv-function beta-cdf base-beta-cdf x a b)
341 (make-rv-function beta-quant base-beta-quant p a b)
342 (make-rv-function beta-dens base-beta-dens x a b)
343 (make-rv-function beta-rand base-beta-rand n a b)
345 (make-rv-function t-cdf base-t-cdf x df)
346 (make-rv-function t-quant base-t-quant p df)
347 (make-rv-function t-dens base-t-dens x df)
348 (make-rv-function t-rand base-t-rand n df)
350 (make-rv-function f-cdf base-f-cdf x ndf ddf)
351 (make-rv-function f-quant base-f-quant p ndf ddf)
352 (make-rv-function f-dens base-f-dens x ndf ddf)
353 (make-rv-function f-rand base-f-rand n ndf ddf)
355 (make-rv-function poisson-cdf base-poisson-cdf x a)
356 (make-rv-function poisson-quant base-poisson-quant p a)
357 (make-rv-function poisson-pmf base-poisson-pmf x a)
358 (make-rv-function poisson-rand base-poisson-rand n a)
360 (make-rv-function binomial-cdf base-binomial-cdf x a b)
361 (make-rv-function binomial-quant base-binomial-quant p a b)
362 (make-rv-function binomial-pmf base-binomial-pmf x a b)
363 (make-rv-function binomial-rand base-binomial-rand n a b)
365 ;;;;
366 ;;;; Documentation
367 ;;;;
369 (setf (documentation 'bivnorm-cdf 'function)
370 "Args: (x y r)
371 Returns the value of the standard bivariate normal distribution function
372 with correlation R at (X, Y). Vectorized.")
374 (setf (documentation 'normal-cdf 'function)
375 "Args: (x)
376 Returns the value of the standard normal distribution function at X.
377 Vectorized.")
379 (setf (documentation 'beta-cdf 'function)
380 "Args: (x alpha beta)
381 Returns the value of the Beta(ALPHA, BETA) distribution function at X.
382 Vectorized.")
384 (setf (documentation 'gamma-cdf 'function)
385 "Args: (x alpha)
386 Returns the value of the Gamma(alpha, 1) distribution function at X.
387 Vectorized.")
389 (setf (documentation 'chisq-cdf 'function)
390 "Args: (x df)
391 Returns the value of the Chi-Square(DF) distribution function at X. Vectorized.")
393 (setf (documentation 't-cdf 'function)
394 "Args: (x df)
395 Returns the value of the T(DF) distribution function at X. Vectorized.")
397 (setf (documentation 'f-cdf 'function)
398 "Args: (x ndf ddf)
399 Returns the value of the F(NDF, DDF) distribution function at X. Vectorized.")
401 (setf (documentation 'cauchy-cdf 'function)
402 "Args: (x)
403 Returns the value of the standard Cauchy distribution function at X.
404 Vectorized.")
406 (setf (documentation 'log-gamma 'function)
407 "Args: (x)
408 Returns the log gamma function of X. Vectorized.")
410 (setf (documentation 'normal-quant 'function)
411 "Args (p)
412 Returns the P-th quantile of the standard normal distribution. Vectorized.")
414 (setf (documentation 'cauchy-quant 'function)
415 "Args (p)
416 Returns the P-th quantile(s) of the standard Cauchy distribution. Vectorized.")
418 (setf (documentation 'beta-quant 'function)
419 "Args: (p alpha beta)
420 Returns the P-th quantile of the Beta(ALPHA, BETA) distribution. Vectorized.")
422 (setf (documentation 'gamma-quant 'function)
423 "Args: (p alpha)
424 Returns the P-th quantile of the Gamma(ALPHA, 1) distribution. Vectorized.")
426 (setf (documentation 'chisq-quant 'function)
427 "Args: (p df)
428 Returns the P-th quantile of the Chi-Square(DF) distribution. Vectorized.")
430 (setf (documentation 't-quant 'function)
431 "Args: (p df)
432 Returns the P-th quantile of the T(DF) distribution. Vectorized.")
434 (setf (documentation 'f-quant 'function)
435 "Args: (p ndf ddf)
436 Returns the P-th quantile of the F(NDF, DDF) distribution. Vectorized.")
438 (setf (documentation 'normal-dens 'function)
439 "Args: (x)
440 Returns the density at X of the standard normal distribution. Vectorized.")
442 (setf (documentation 'cauchy-dens 'function)
443 "Args: (x)
444 Returns the density at X of the standard Cauchy distribution. Vectorized.")
446 (setf (documentation 'beta-dens 'function)
447 "Args: (x alpha beta)
448 Returns the density at X of the Beta(ALPHA, BETA) distribution. Vectorized.")
450 (setf (documentation 'gamma-dens 'function)
451 "Args: (x alpha)
452 Returns the density at X of the Gamma(ALPHA, 1) distribution. Vectorized.")
454 (setf (documentation 'chisq-dens 'function)
455 "Args: (x alpha)
456 Returns the density at X of the Chi-Square(DF) distribution. Vectorized.")
458 (setf (documentation 't-dens 'function)
459 "Args: (x alpha)
460 Returns the density at X of the T(DF) distribution. Vectorized.")
462 (setf (documentation 'f-dens 'function)
463 "Args: (x ndf ddf)
464 Returns the density at X of the F(NDF, DDF) distribution. Vectorized.")
466 (setf (documentation 'uniform-rand 'function)
467 "Args: (n)
468 Returns a list of N uniform random variables from the range (0, 1).
469 Vectorized.")
471 (setf (documentation 'normal-rand 'function)
472 "Args: (n)
473 Returns a list of N standard normal random numbers. Vectorized.")
475 (setf (documentation 'cauchy-rand 'function)
476 "Args: (n)
477 Returns a list of N standard Cauchy random numbers. Vectorized.")
479 (setf (documentation 't-rand 'function)
480 "Args: (n df)
481 Returns a list of N T(DF) random variables. Vectorized.")
483 (setf (documentation 'f-rand 'function)
484 "Args: (n ndf ddf)
485 Returns a list of N F(NDF, DDF) random variables. Vectorized.")
487 (setf (documentation 'gamma-rand 'function)
488 "Args: (n a)
489 Returns a list of N Gamma(A, 1) random variables. Vectorized.")
491 (setf (documentation 'chisq-rand 'function)
492 "Args: (n df)
493 Returns a list of N Chi-Square(DF) random variables. Vectorized.")
495 (setf (documentation 'beta-rand 'function)
496 "Args: (n a b)
497 Returns a list of N beta(A, B) random variables. Vectorized.")
499 (setf (documentation 'binomial-cdf 'function)
500 "Args (x n p)
501 Returns value of the Binomial(N, P) distribution function at X. Vectorized.")
503 (setf (documentation 'poisson-cdf 'function)
504 "Args (x mu)
505 Returns value of the Poisson(MU) distribution function at X. Vectorized.")
507 (setf (documentation 'binomial-pmf 'function)
508 "Args (k n p)
509 Returns value of the Binomial(N, P) pmf function at integer K. Vectorized.")
511 (setf (documentation 'poisson-pmf 'function)
512 "Args (k mu)
513 Returns value of the Poisson(MU) pmf function at integer K. Vectorized.")
515 (setf (documentation 'binomial-quant 'function)
516 "Args: (x n p)
517 Returns x-th quantile (left continuous inverse) of Binomial(N, P) cdf.
518 Vectorized.")
520 (setf (documentation 'poisson-quant 'function)
521 "Args: (x mu)
522 Returns x-th quantile (left continuous inverse) of Poisson(MU) cdf.
523 Vectorized.")
525 (setf (documentation 'binomial-rand 'function)
526 "Args: (k n p)
527 Returns list of K draws from the Binomial(N, P) distribution. Vectorized.")
529 (setf (documentation 'poisson-rand 'function)
530 "Args: (k mu)
531 Returns list of K draws from the Poisson(MU) distribution. Vectorized.")