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.
6 ;;; dists -- Lisp-Stat interface to basic probability distribution routines
8 ;;; Copyright (c) 1991, by Luke Tierney. Permission is granted for
17 (defpackage :lisp-stat-probability
24 normal-cdf normal-quant normal-dens normal-rand
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
))
35 (in-package :lisp-stat-probability
)
37 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
39 ;;; CFFI support for Probability Distributions
41 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
44 ;;; C-callable uniform generator
47 (cffi:defcfun
("register_uni" register-uni
)
49 (cffi:defcallback ccl-uni
:int
() (ccl-store-double (random 1.0)) 0)
50 (register-uni (cffi:callback ccl-uni
))
52 (defun one-uniform-rand () (random 1.0))
55 ;;; Log-gamma function
58 (cffi:defcfun
("ccl_gamma" ccl-base-log-gamma
)
60 (defun base-log-gamma (x)
61 (ccl-base-log-gamma (float x
1d0
)))
64 ;;; Normal distribution
67 (cffi:defcfun
("ccl_normalcdf" ccl-base-normal-cdf
)
69 (defun base-normal-cdf (x)
70 (ccl-base-normal-cdf (float x
1d0
)))
72 (cffi:defcfun
("ccl_normalquant" ccl-base-normal-quant
)
74 (defun base-normal-quant (x)
75 (ccl-base-normal-quant (float x
1d0
)))
77 (cffi:defcfun
("ccl_normaldens" ccl-base-normal-dens
)
79 (defun base-normal-dens (x)
80 (ccl-base-normal-dens (float x
1d0
)))
82 (cffi:defcfun
("ccl_normalrand" one-normal-rand
)
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
)))
91 ;;; Cauchy distribution
94 (cffi:defcfun
("ccl_cauchycdf" ccl-base-cauchy-cdf
)
96 (defun base-cauchy-cdf (x)
97 (ccl-base-cauchy-cdf (float x
1d0
)))
99 (cffi:defcfun
("ccl_cauchyquant" ccl-base-cauchy-quant
)
101 (defun base-cauchy-quant (x)
102 (ccl-base-cauchy-quant (float x
1d0
)))
104 (cffi:defcfun
("ccl_cauchydens" ccl-base-cauchy-dens
)
106 (defun base-cauchy-dens (x)
107 (ccl-base-cauchy-dens (float x
1d0
)))
109 (cffi:defcfun
("ccl_cauchyrand" one-cauchy-rand
)
113 ;;;; Gamma distribution
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
)))
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
)))
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
)))
131 (cffi:defcfun
("ccl_gammarand" ccl-gamma-rand
)
133 (defun one-gamma-rand (x)
134 (ccl-gamma-rand (float x
1d0
)))
137 ;;;; Chi-square distribution
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
)))
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
)))
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
)))
155 (cffi:defcfun
("ccl_chisqrand" ccl-chisq-rand
)
157 (defun one-chisq-rand (x)
158 (ccl-chisq-rand (float x
1d0
)))
161 ;;;; Beta distribution
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
)))
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
)))
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
)))
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
)))
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
)))
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
)))
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
)))
203 (cffi:defcfun
("ccl_trand" ccl-t-rand
)
205 (defun one-t-rand (x)
206 (ccl-t-rand (float x
1d0
)))
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
)))
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
)))
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
)))
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
)))
232 ;;;; Poisson distribution
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
)))
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
)))
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
)))
250 (cffi:defcfun
("ccl_poissonrand" ccl-poisson-rand
)
252 (defun one-poisson-rand (x)
253 (ccl-poisson-rand (float x
1d0
)))
256 ;;;; Binomial distribution
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
)))
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
)))
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
)))
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
)))
283 ;;; definitions though macros
285 (defmacro defbaserand
(name onefun
&rest args
)
286 `(defun ,name
(n ,@args
)
288 (dotimes (i n result
)
289 (declare (fixnum i
) (inline ,onefun
))
290 (setf result
(cons (,onefun
,@args
) result
))))))
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
)
303 (make-rv-function log-gamma base-log-gamma x
)
305 (make-rv-function uniform-rand base-uniform-rand n
)
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
)
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
)
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
)
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
)
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
)
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
)
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
)
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
)
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
)
357 (setf (documentation 'bivnorm-cdf
'function
)
359 Returns the value of the standard bivariate normal distribution function
360 with correlation R at (X, Y). Vectorized.")
362 (setf (documentation 'normal-cdf
'function
)
364 Returns the value of the standard normal distribution function at X.
367 (setf (documentation 'beta-cdf
'function
)
368 "Args: (x alpha beta)
369 Returns the value of the Beta(ALPHA, BETA) distribution function at X.
372 (setf (documentation 'gamma-cdf
'function
)
374 Returns the value of the Gamma(alpha, 1) distribution function at X.
377 (setf (documentation 'chisq-cdf
'function
)
379 Returns the value of the Chi-Square(DF) distribution function at X. Vectorized.")
381 (setf (documentation 't-cdf
'function
)
383 Returns the value of the T(DF) distribution function at X. Vectorized.")
385 (setf (documentation 'f-cdf
'function
)
387 Returns the value of the F(NDF, DDF) distribution function at X. Vectorized.")
389 (setf (documentation 'cauchy-cdf
'function
)
391 Returns the value of the standard Cauchy distribution function at X.
394 (setf (documentation 'log-gamma
'function
)
396 Returns the log gamma function of X. Vectorized.")
398 (setf (documentation 'normal-quant
'function
)
400 Returns the P-th quantile of the standard normal distribution. Vectorized.")
402 (setf (documentation 'cauchy-quant
'function
)
404 Returns the P-th quantile(s) of the standard Cauchy distribution. Vectorized.")
406 (setf (documentation 'beta-quant
'function
)
407 "Args: (p alpha beta)
408 Returns the P-th quantile of the Beta(ALPHA, BETA) distribution. Vectorized.")
410 (setf (documentation 'gamma-quant
'function
)
412 Returns the P-th quantile of the Gamma(ALPHA, 1) distribution. Vectorized.")
414 (setf (documentation 'chisq-quant
'function
)
416 Returns the P-th quantile of the Chi-Square(DF) distribution. Vectorized.")
418 (setf (documentation 't-quant
'function
)
420 Returns the P-th quantile of the T(DF) distribution. Vectorized.")
422 (setf (documentation 'f-quant
'function
)
424 Returns the P-th quantile of the F(NDF, DDF) distribution. Vectorized.")
426 (setf (documentation 'normal-dens
'function
)
428 Returns the density at X of the standard normal distribution. Vectorized.")
430 (setf (documentation 'cauchy-dens
'function
)
432 Returns the density at X of the standard Cauchy distribution. Vectorized.")
434 (setf (documentation 'beta-dens
'function
)
435 "Args: (x alpha beta)
436 Returns the density at X of the Beta(ALPHA, BETA) distribution. Vectorized.")
438 (setf (documentation 'gamma-dens
'function
)
440 Returns the density at X of the Gamma(ALPHA, 1) distribution. Vectorized.")
442 (setf (documentation 'chisq-dens
'function
)
444 Returns the density at X of the Chi-Square(DF) distribution. Vectorized.")
446 (setf (documentation 't-dens
'function
)
448 Returns the density at X of the T(DF) distribution. Vectorized.")
450 (setf (documentation 'f-dens
'function
)
452 Returns the density at X of the F(NDF, DDF) distribution. Vectorized.")
454 (setf (documentation 'uniform-rand
'function
)
456 Returns a list of N uniform random variables from the range (0, 1).
459 (setf (documentation 'normal-rand
'function
)
461 Returns a list of N standard normal random numbers. Vectorized.")
463 (setf (documentation 'cauchy-rand
'function
)
465 Returns a list of N standard Cauchy random numbers. Vectorized.")
467 (setf (documentation 't-rand
'function
)
469 Returns a list of N T(DF) random variables. Vectorized.")
471 (setf (documentation 'f-rand
'function
)
473 Returns a list of N F(NDF, DDF) random variables. Vectorized.")
475 (setf (documentation 'gamma-rand
'function
)
477 Returns a list of N Gamma(A, 1) random variables. Vectorized.")
479 (setf (documentation 'chisq-rand
'function
)
481 Returns a list of N Chi-Square(DF) random variables. Vectorized.")
483 (setf (documentation 'beta-rand
'function
)
485 Returns a list of N beta(A, B) random variables. Vectorized.")
487 (setf (documentation 'binomial-cdf
'function
)
489 Returns value of the Binomial(N, P) distribution function at X. Vectorized.")
491 (setf (documentation 'poisson-cdf
'function
)
493 Returns value of the Poisson(MU) distribution function at X. Vectorized.")
495 (setf (documentation 'binomial-pmf
'function
)
497 Returns value of the Binomial(N, P) pmf function at integer K. Vectorized.")
499 (setf (documentation 'poisson-pmf
'function
)
501 Returns value of the Poisson(MU) pmf function at integer K. Vectorized.")
503 (setf (documentation 'binomial-quant
'function
)
505 Returns x-th quantile (left continuous inverse) of Binomial(N, P) cdf.
508 (setf (documentation 'poisson-quant
'function
)
510 Returns x-th quantile (left continuous inverse) of Poisson(MU) cdf.
513 (setf (documentation 'binomial-rand
'function
)
515 Returns list of K draws from the Binomial(N, P) distribution. Vectorized.")
517 (setf (documentation 'poisson-rand
'function
)
519 Returns list of K draws from the Poisson(MU) distribution. Vectorized.")