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
15 (defpackage :lisp-stat-probability
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
))
32 (in-package :lisp-stat-probability
)
34 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
36 ;;; CFFI support for Probability Distributions
38 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
41 ;;; C-callable uniform generator
44 (cffi:defcfun
("register_uni" register-uni
)
46 (cffi:defcallback ccl-uni
:int
() (ccl-store-double (random 1.0)) 0)
47 (register-uni (cffi:callback ccl-uni
))
49 (defun one-uniform-rand () (random 1.0))
52 ;;;; Log-gamma function
55 (cffi:defcfun
("ccl_gamma" ccl-base-log-gamma
)
57 (defun base-log-gamma (x)
58 (ccl-base-log-gamma (float x
1d0
)))
61 ;;;; Normal distribution
64 (cffi:defcfun
("ccl_normalcdf" ccl-base-normal-cdf
)
66 (defun base-normal-cdf (x)
67 (ccl-base-normal-cdf (float x
1d0
)))
69 (cffi:defcfun
("ccl_normalquant" ccl-base-normal-quant
)
71 (defun base-normal-quant (x)
72 (ccl-base-normal-quant (float x
1d0
)))
74 (cffi:defcfun
("ccl_normaldens" ccl-base-normal-dens
)
76 (defun base-normal-dens (x)
77 (ccl-base-normal-dens (float x
1d0
)))
79 (cffi:defcfun
("ccl_normalrand" one-normal-rand
)
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
)))
88 ;;;; Cauchy distribution
91 (cffi:defcfun
("ccl_cauchycdf" ccl-base-cauchy-cdf
)
93 (defun base-cauchy-cdf (x)
94 (ccl-base-cauchy-cdf (float x
1d0
)))
96 (cffi:defcfun
("ccl_cauchyquant" ccl-base-cauchy-quant
)
98 (defun base-cauchy-quant (x)
99 (ccl-base-cauchy-quant (float x
1d0
)))
101 (cffi:defcfun
("ccl_cauchydens" ccl-base-cauchy-dens
)
103 (defun base-cauchy-dens (x)
104 (ccl-base-cauchy-dens (float x
1d0
)))
106 (cffi:defcfun
("ccl_cauchyrand" one-cauchy-rand
)
110 ;;;; Gamma distribution
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
)))
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
)))
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
)))
128 (cffi:defcfun
("ccl_gammarand" ccl-gamma-rand
)
130 (defun one-gamma-rand (x)
131 (ccl-gamma-rand (float x
1d0
)))
134 ;;;; Chi-square distribution
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
)))
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
)))
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
)))
152 (cffi:defcfun
("ccl_chisqrand" ccl-chisq-rand
)
154 (defun one-chisq-rand (x)
155 (ccl-chisq-rand (float x
1d0
)))
158 ;;;; Beta distribution
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
)))
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
)))
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
)))
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
)))
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
)))
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
)))
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
)))
200 (cffi:defcfun
("ccl_trand" ccl-t-rand
)
202 (defun one-t-rand (x)
203 (ccl-t-rand (float x
1d0
)))
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
)))
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
)))
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
)))
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
)))
229 ;;;; Poisson distribution
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
)))
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
)))
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
)))
247 (cffi:defcfun
("ccl_poissonrand" ccl-poisson-rand
)
249 (defun one-poisson-rand (x)
250 (ccl-poisson-rand (float x
1d0
)))
253 ;;;; Binomial distribution
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
)))
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
)))
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
)))
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
)))
280 ;;; definitions though macros
282 (defmacro defbaserand
(name onefun
&rest args
)
283 `(defun ,name
(n ,@args
)
285 (dotimes (i n result
)
286 (declare (fixnum i
) (inline ,onefun
))
287 (setf result
(cons (,onefun
,@args
) result
))))))
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
)
300 (make-rv-function log-gamma base-log-gamma x
)
302 (make-rv-function uniform-rand base-uniform-rand n
)
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
)
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
)
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
)
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
)
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
)
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
)
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
)
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
)
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
)
354 (setf (documentation 'bivnorm-cdf
'function
)
356 Returns the value of the standard bivariate normal distribution function
357 with correlation R at (X, Y). Vectorized.")
359 (setf (documentation 'normal-cdf
'function
)
361 Returns the value of the standard normal distribution function at X.
364 (setf (documentation 'beta-cdf
'function
)
365 "Args: (x alpha beta)
366 Returns the value of the Beta(ALPHA, BETA) distribution function at X.
369 (setf (documentation 'gamma-cdf
'function
)
371 Returns the value of the Gamma(alpha, 1) distribution function at X.
374 (setf (documentation 'chisq-cdf
'function
)
376 Returns the value of the Chi-Square(DF) distribution function at X. Vectorized.")
378 (setf (documentation 't-cdf
'function
)
380 Returns the value of the T(DF) distribution function at X. Vectorized.")
382 (setf (documentation 'f-cdf
'function
)
384 Returns the value of the F(NDF, DDF) distribution function at X. Vectorized.")
386 (setf (documentation 'cauchy-cdf
'function
)
388 Returns the value of the standard Cauchy distribution function at X.
391 (setf (documentation 'log-gamma
'function
)
393 Returns the log gamma function of X. Vectorized.")
395 (setf (documentation 'normal-quant
'function
)
397 Returns the P-th quantile of the standard normal distribution. Vectorized.")
399 (setf (documentation 'cauchy-quant
'function
)
401 Returns the P-th quantile(s) of the standard Cauchy distribution. Vectorized.")
403 (setf (documentation 'beta-quant
'function
)
404 "Args: (p alpha beta)
405 Returns the P-th quantile of the Beta(ALPHA, BETA) distribution. Vectorized.")
407 (setf (documentation 'gamma-quant
'function
)
409 Returns the P-th quantile of the Gamma(ALPHA, 1) distribution. Vectorized.")
411 (setf (documentation 'chisq-quant
'function
)
413 Returns the P-th quantile of the Chi-Square(DF) distribution. Vectorized.")
415 (setf (documentation 't-quant
'function
)
417 Returns the P-th quantile of the T(DF) distribution. Vectorized.")
419 (setf (documentation 'f-quant
'function
)
421 Returns the P-th quantile of the F(NDF, DDF) distribution. Vectorized.")
423 (setf (documentation 'normal-dens
'function
)
425 Returns the density at X of the standard normal distribution. Vectorized.")
427 (setf (documentation 'cauchy-dens
'function
)
429 Returns the density at X of the standard Cauchy distribution. Vectorized.")
431 (setf (documentation 'beta-dens
'function
)
432 "Args: (x alpha beta)
433 Returns the density at X of the Beta(ALPHA, BETA) distribution. Vectorized.")
435 (setf (documentation 'gamma-dens
'function
)
437 Returns the density at X of the Gamma(ALPHA, 1) distribution. Vectorized.")
439 (setf (documentation 'chisq-dens
'function
)
441 Returns the density at X of the Chi-Square(DF) distribution. Vectorized.")
443 (setf (documentation 't-dens
'function
)
445 Returns the density at X of the T(DF) distribution. Vectorized.")
447 (setf (documentation 'f-dens
'function
)
449 Returns the density at X of the F(NDF, DDF) distribution. Vectorized.")
451 (setf (documentation 'uniform-rand
'function
)
453 Returns a list of N uniform random variables from the range (0, 1).
456 (setf (documentation 'normal-rand
'function
)
458 Returns a list of N standard normal random numbers. Vectorized.")
460 (setf (documentation 'cauchy-rand
'function
)
462 Returns a list of N standard Cauchy random numbers. Vectorized.")
464 (setf (documentation 't-rand
'function
)
466 Returns a list of N T(DF) random variables. Vectorized.")
468 (setf (documentation 'f-rand
'function
)
470 Returns a list of N F(NDF, DDF) random variables. Vectorized.")
472 (setf (documentation 'gamma-rand
'function
)
474 Returns a list of N Gamma(A, 1) random variables. Vectorized.")
476 (setf (documentation 'chisq-rand
'function
)
478 Returns a list of N Chi-Square(DF) random variables. Vectorized.")
480 (setf (documentation 'beta-rand
'function
)
482 Returns a list of N beta(A, B) random variables. Vectorized.")
484 (setf (documentation 'binomial-cdf
'function
)
486 Returns value of the Binomial(N, P) distribution function at X. Vectorized.")
488 (setf (documentation 'poisson-cdf
'function
)
490 Returns value of the Poisson(MU) distribution function at X. Vectorized.")
492 (setf (documentation 'binomial-pmf
'function
)
494 Returns value of the Binomial(N, P) pmf function at integer K. Vectorized.")
496 (setf (documentation 'poisson-pmf
'function
)
498 Returns value of the Poisson(MU) pmf function at integer K. Vectorized.")
500 (setf (documentation 'binomial-quant
'function
)
502 Returns x-th quantile (left continuous inverse) of Binomial(N, P) cdf.
505 (setf (documentation 'poisson-quant
'function
)
507 Returns x-th quantile (left continuous inverse) of Poisson(MU) cdf.
510 (setf (documentation 'binomial-rand
'function
)
512 Returns list of K draws from the Binomial(N, P) distribution. Vectorized.")
514 (setf (documentation 'poisson-rand
'function
)
516 Returns list of K draws from the Poisson(MU) distribution. Vectorized.")