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
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.
24 (defpackage :lisp-stat-probability
29 (:export log-gamma set-seed
31 normal-cdf normal-quant normal-dens normal-rand
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
)
46 "stupid dummy function, need to implement rng seeding tool."
50 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
52 ;;; CFFI support for Probability Distributions
54 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
57 ;;; C-callable uniform generator
60 (defcfun ("register_uni" register-uni
)
62 (defcallback ccl-uni
:int
() (ccl-store-double (random 1.0)) 0)
63 (register-uni (callback ccl-uni
))
65 (defun one-uniform-rand () (random 1.0))
68 ;;; Log-gamma function
71 (defcfun ("ccl_gamma" ccl-base-log-gamma
)
73 (defun base-log-gamma (x)
74 (ccl-base-log-gamma (float x
1d0
)))
77 ;;; Normal distribution
80 (defcfun ("ccl_normalcdf" ccl-base-normal-cdf
)
82 (defun base-normal-cdf (x)
83 (ccl-base-normal-cdf (float x
1d0
)))
85 (defcfun ("ccl_normalquant" ccl-base-normal-quant
)
87 (defun base-normal-quant (x)
88 (ccl-base-normal-quant (float x
1d0
)))
90 (defcfun ("ccl_normaldens" ccl-base-normal-dens
)
92 (defun base-normal-dens (x)
93 (ccl-base-normal-dens (float x
1d0
)))
95 (defcfun ("ccl_normalrand" one-normal-rand
)
98 (defcfun ("ccl_bnormcdf" ccl-base-bivnorm-cdf
)
99 :double
(x :double
) (y :double
) (z :double
))
100 (defun base-bivnorm-cdf (x y z
)
101 (ccl-base-bivnorm-cdf (float x
1d0
) (float y
1d0
) (float z
1d0
)))
104 ;;; Cauchy distribution
107 (defcfun ("ccl_cauchycdf" ccl-base-cauchy-cdf
)
109 (defun base-cauchy-cdf (x)
110 (ccl-base-cauchy-cdf (float x
1d0
)))
112 (defcfun ("ccl_cauchyquant" ccl-base-cauchy-quant
)
114 (defun base-cauchy-quant (x)
115 (ccl-base-cauchy-quant (float x
1d0
)))
117 (defcfun ("ccl_cauchydens" ccl-base-cauchy-dens
)
119 (defun base-cauchy-dens (x)
120 (ccl-base-cauchy-dens (float x
1d0
)))
122 (defcfun ("ccl_cauchyrand" one-cauchy-rand
)
126 ;;;; Gamma distribution
129 (defcfun ("ccl_gammacdf" ccl-base-gamma-cdf
)
130 :double
(x :double
) (y :double
))
131 (defun base-gamma-cdf (x y
)
132 (ccl-base-gamma-cdf (float x
1d0
) (float y
1d0
)))
134 (defcfun ("ccl_gammaquant" ccl-base-gamma-quant
)
135 :double
(x :double
) (y :double
))
136 (defun base-gamma-quant (x y
)
137 (ccl-base-gamma-quant (float x
1d0
) (float y
1d0
)))
139 (defcfun ("ccl_gammadens" ccl-base-gamma-dens
)
140 :double
(x :double
) (y :double
))
141 (defun base-gamma-dens (x y
)
142 (ccl-base-gamma-dens (float x
1d0
) (float y
1d0
)))
144 (defcfun ("ccl_gammarand" ccl-gamma-rand
)
146 (defun one-gamma-rand (x)
147 (ccl-gamma-rand (float x
1d0
)))
150 ;;;; Chi-square distribution
153 (defcfun ("ccl_chisqcdf" ccl-base-chisq-cdf
)
154 :double
(x :double
) (y :double
))
155 (defun base-chisq-cdf (x y
)
156 (ccl-base-chisq-cdf (float x
1d0
) (float y
1d0
)))
158 (defcfun ("ccl_chisqquant" ccl-base-chisq-quant
)
159 :double
(x :double
) (y :double
))
160 (defun base-chisq-quant (x y
)
161 (ccl-base-chisq-quant (float x
1d0
) (float y
1d0
)))
163 (defcfun ("ccl_chisqdens" ccl-base-chisq-dens
)
164 :double
(x :double
) (y :double
))
165 (defun base-chisq-dens (x y
)
166 (ccl-base-chisq-dens (float x
1d0
) (float y
1d0
)))
168 (defcfun ("ccl_chisqrand" ccl-chisq-rand
)
170 (defun one-chisq-rand (x)
171 (ccl-chisq-rand (float x
1d0
)))
174 ;;;; Beta distribution
177 (defcfun ("ccl_betacdf" ccl-base-beta-cdf
)
178 :double
(x :double
) (y :double
) (z :double
))
179 (defun base-beta-cdf (x y z
)
180 (ccl-base-beta-cdf (float x
1d0
) (float y
1d0
) (float z
1d0
)))
182 (defcfun ("ccl_betaquant" ccl-base-beta-quant
)
183 :double
(x :double
) (y :double
) (z :double
))
184 (defun base-beta-quant (x y z
)
185 (ccl-base-beta-quant (float x
1d0
) (float y
1d0
) (float z
1d0
)))
187 (defcfun ("ccl_betadens" ccl-base-beta-dens
)
188 :double
(x :double
) (y :double
) (z :double
))
189 (defun base-beta-dens (x y z
)
190 (ccl-base-beta-dens (float x
1d0
) (float y
1d0
) (float z
1d0
)))
192 (defcfun ("ccl_betarand" ccl-beta-rand
)
193 :double
(x :double
) (y :double
))
194 (defun one-beta-rand (x y
)
195 (ccl-beta-rand (float x
1d0
) (float y
1d0
)))
201 (defcfun ("ccl_tcdf" ccl-base-t-cdf
)
202 :double
(x :double
) (y :double
))
203 (defun base-t-cdf (x y
)
204 (ccl-base-t-cdf (float x
1d0
) (float y
1d0
)))
206 (defcfun ("ccl_tquant" ccl-base-t-quant
)
207 :double
(x :double
) (y :double
))
208 (defun base-t-quant (x y
)
209 (ccl-base-t-quant (float x
1d0
) (float y
1d0
)))
211 (defcfun ("ccl_tdens" ccl-base-t-dens
)
212 :double
(x :double
) (y :double
))
213 (defun base-t-dens (x y
)
214 (ccl-base-t-dens (float x
1d0
) (float y
1d0
)))
216 (defcfun ("ccl_trand" ccl-t-rand
)
218 (defun one-t-rand (x)
219 (ccl-t-rand (float x
1d0
)))
225 (defcfun ("ccl_fcdf" ccl-base-f-cdf
)
226 :double
(x :double
) (y :double
) (z :double
))
227 (defun base-f-cdf (x y z
)
228 (ccl-base-f-cdf (float x
1d0
) (float y
1d0
) (float z
1d0
)))
230 (defcfun ("ccl_fquant" ccl-base-f-quant
)
231 :double
(x :double
) (y :double
) (z :double
))
232 (defun base-f-quant (x y z
)
233 (ccl-base-f-quant (float x
1d0
) (float y
1d0
) (float z
1d0
)))
235 (defcfun ("ccl_fdens" ccl-base-f-dens
)
236 :double
(x :double
) (y :double
) (z :double
))
237 (defun base-f-dens (x y z
)
238 (ccl-base-f-dens (float x
1d0
) (float y
1d0
) (float z
1d0
)))
240 (defcfun ("ccl_frand" ccl-f-rand
)
241 :double
(x :double
) (y :double
))
242 (defun one-f-rand (x y
) (ccl-f-rand (float x
1d0
) (float y
1d0
)))
245 ;;;; Poisson distribution
248 (defcfun ("ccl_poissoncdf" ccl-base-poisson-cdf
)
249 :double
(x :double
) (y :double
))
250 (defun base-poisson-cdf (x y
)
251 (ccl-base-poisson-cdf (float x
1d0
) (float y
1d0
)))
253 (defcfun ("ccl_poissonquant" ccl-base-poisson-quant
)
254 :int
(x :double
) (y :double
))
255 (defun base-poisson-quant (x y
)
256 (ccl-base-poisson-quant (float x
1d0
) (float y
1d0
)))
258 (defcfun ("ccl_poissonpmf" ccl-base-poisson-pmf
)
259 :double
(x :int
) (y :double
))
260 (defun base-poisson-pmf (x y
)
261 (ccl-base-poisson-pmf x
(float y
1d0
)))
263 (defcfun ("ccl_poissonrand" ccl-poisson-rand
)
265 (defun one-poisson-rand (x)
266 (ccl-poisson-rand (float x
1d0
)))
269 ;;;; Binomial distribution
272 (defcfun ("ccl_binomialcdf" ccl-base-binomial-cdf
)
273 :double
(x :double
) (y :int
) (z :double
))
274 (defun base-binomial-cdf (x y z
)
275 (ccl-base-binomial-cdf (float x
1d0
) y
(float z
1d0
)))
277 (defcfun ("ccl_binomialquant" ccl-base-binomial-quant
)
278 :int
(x :double
) (y :int
) (z :double
))
279 (defun base-binomial-quant (x y z
)
280 (ccl-base-binomial-quant (float x
1d0
) y
(float z
1d0
)))
282 (defcfun ("ccl_binomialpmf" ccl-base-binomial-pmf
)
283 :double
(x :int
) (y :int
) (z :double
))
284 (defun base-binomial-pmf (x y z
)
285 (ccl-base-binomial-pmf x y
(float z
1d0
)))
287 (defcfun ("ccl_binomialrand" ccl-binomial-rand
)
288 :int
(x :int
) (y :double
))
289 (defun one-binomial-rand (x y
)
290 (ccl-binomial-rand x
(float y
1d0
)))
296 ;;; definitions though macros
298 (defmacro defbaserand
(name onefun
&rest args
)
299 `(defun ,name
(n ,@args
)
301 (dotimes (i n result
)
302 (declare (fixnum i
) (inline ,onefun
))
303 (setf result
(cons (,onefun
,@args
) result
))))))
305 (defbaserand base-uniform-rand one-uniform-rand
)
306 (defbaserand base-normal-rand one-normal-rand
)
307 (defbaserand base-cauchy-rand one-cauchy-rand
)
308 (defbaserand base-gamma-rand one-gamma-rand a
)
309 (defbaserand base-chisq-rand one-chisq-rand df
)
310 (defbaserand base-beta-rand one-beta-rand a b
)
311 (defbaserand base-t-rand one-t-rand df
)
312 (defbaserand base-f-rand one-f-rand ndf ddf
)
313 (defbaserand base-poisson-rand one-poisson-rand a
)
314 (defbaserand base-binomial-rand one-binomial-rand a b
)
316 (make-rv-function log-gamma base-log-gamma x
)
318 (make-rv-function uniform-rand base-uniform-rand n
)
320 (make-rv-function normal-cdf base-normal-cdf x
)
321 (make-rv-function normal-quant base-normal-quant p
)
322 (make-rv-function normal-dens base-normal-dens x
)
323 (make-rv-function normal-rand base-normal-rand n
)
324 (make-rv-function bivnorm-cdf base-bivnorm-cdf x y r
)
326 (make-rv-function cauchy-cdf base-cauchy-cdf x
)
327 (make-rv-function cauchy-quant base-cauchy-quant p
)
328 (make-rv-function cauchy-dens base-cauchy-dens x
)
329 (make-rv-function cauchy-rand base-cauchy-rand n
)
331 (make-rv-function gamma-cdf base-gamma-cdf x a
)
332 (make-rv-function gamma-quant base-gamma-quant p a
)
333 (make-rv-function gamma-dens base-gamma-dens x a
)
334 (make-rv-function gamma-rand base-gamma-rand n a
)
336 (make-rv-function chisq-cdf base-chisq-cdf x df
)
337 (make-rv-function chisq-quant base-chisq-quant p df
)
338 (make-rv-function chisq-dens base-chisq-dens x df
)
339 (make-rv-function chisq-rand base-chisq-rand n df
)
341 (make-rv-function beta-cdf base-beta-cdf x a b
)
342 (make-rv-function beta-quant base-beta-quant p a b
)
343 (make-rv-function beta-dens base-beta-dens x a b
)
344 (make-rv-function beta-rand base-beta-rand n a b
)
346 (make-rv-function t-cdf base-t-cdf x df
)
347 (make-rv-function t-quant base-t-quant p df
)
348 (make-rv-function t-dens base-t-dens x df
)
349 (make-rv-function t-rand base-t-rand n df
)
351 (make-rv-function f-cdf base-f-cdf x ndf ddf
)
352 (make-rv-function f-quant base-f-quant p ndf ddf
)
353 (make-rv-function f-dens base-f-dens x ndf ddf
)
354 (make-rv-function f-rand base-f-rand n ndf ddf
)
356 (make-rv-function poisson-cdf base-poisson-cdf x a
)
357 (make-rv-function poisson-quant base-poisson-quant p a
)
358 (make-rv-function poisson-pmf base-poisson-pmf x a
)
359 (make-rv-function poisson-rand base-poisson-rand n a
)
361 (make-rv-function binomial-cdf base-binomial-cdf x a b
)
362 (make-rv-function binomial-quant base-binomial-quant p a b
)
363 (make-rv-function binomial-pmf base-binomial-pmf x a b
)
364 (make-rv-function binomial-rand base-binomial-rand n a b
)
370 (setf (documentation 'bivnorm-cdf
'function
)
372 Returns the value of the standard bivariate normal distribution function
373 with correlation R at (X, Y). Vectorized.")
375 (setf (documentation 'normal-cdf
'function
)
377 Returns the value of the standard normal distribution function at X.
380 (setf (documentation 'beta-cdf
'function
)
381 "Args: (x alpha beta)
382 Returns the value of the Beta(ALPHA, BETA) distribution function at X.
385 (setf (documentation 'gamma-cdf
'function
)
387 Returns the value of the Gamma(alpha, 1) distribution function at X.
390 (setf (documentation 'chisq-cdf
'function
)
392 Returns the value of the Chi-Square(DF) distribution function at X. Vectorized.")
394 (setf (documentation 't-cdf
'function
)
396 Returns the value of the T(DF) distribution function at X. Vectorized.")
398 (setf (documentation 'f-cdf
'function
)
400 Returns the value of the F(NDF, DDF) distribution function at X. Vectorized.")
402 (setf (documentation 'cauchy-cdf
'function
)
404 Returns the value of the standard Cauchy distribution function at X.
407 (setf (documentation 'log-gamma
'function
)
409 Returns the log gamma function of X. Vectorized.")
411 (setf (documentation 'normal-quant
'function
)
413 Returns the P-th quantile of the standard normal distribution. Vectorized.")
415 (setf (documentation 'cauchy-quant
'function
)
417 Returns the P-th quantile(s) of the standard Cauchy distribution. Vectorized.")
419 (setf (documentation 'beta-quant
'function
)
420 "Args: (p alpha beta)
421 Returns the P-th quantile of the Beta(ALPHA, BETA) distribution. Vectorized.")
423 (setf (documentation 'gamma-quant
'function
)
425 Returns the P-th quantile of the Gamma(ALPHA, 1) distribution. Vectorized.")
427 (setf (documentation 'chisq-quant
'function
)
429 Returns the P-th quantile of the Chi-Square(DF) distribution. Vectorized.")
431 (setf (documentation 't-quant
'function
)
433 Returns the P-th quantile of the T(DF) distribution. Vectorized.")
435 (setf (documentation 'f-quant
'function
)
437 Returns the P-th quantile of the F(NDF, DDF) distribution. Vectorized.")
439 (setf (documentation 'normal-dens
'function
)
441 Returns the density at X of the standard normal distribution. Vectorized.")
443 (setf (documentation 'cauchy-dens
'function
)
445 Returns the density at X of the standard Cauchy distribution. Vectorized.")
447 (setf (documentation 'beta-dens
'function
)
448 "Args: (x alpha beta)
449 Returns the density at X of the Beta(ALPHA, BETA) distribution. Vectorized.")
451 (setf (documentation 'gamma-dens
'function
)
453 Returns the density at X of the Gamma(ALPHA, 1) distribution. Vectorized.")
455 (setf (documentation 'chisq-dens
'function
)
457 Returns the density at X of the Chi-Square(DF) distribution. Vectorized.")
459 (setf (documentation 't-dens
'function
)
461 Returns the density at X of the T(DF) distribution. Vectorized.")
463 (setf (documentation 'f-dens
'function
)
465 Returns the density at X of the F(NDF, DDF) distribution. Vectorized.")
467 (setf (documentation 'uniform-rand
'function
)
469 Returns a list of N uniform random variables from the range (0, 1).
472 (setf (documentation 'normal-rand
'function
)
474 Returns a list of N standard normal random numbers. Vectorized.")
476 (setf (documentation 'cauchy-rand
'function
)
478 Returns a list of N standard Cauchy random numbers. Vectorized.")
480 (setf (documentation 't-rand
'function
)
482 Returns a list of N T(DF) random variables. Vectorized.")
484 (setf (documentation 'f-rand
'function
)
486 Returns a list of N F(NDF, DDF) random variables. Vectorized.")
488 (setf (documentation 'gamma-rand
'function
)
490 Returns a list of N Gamma(A, 1) random variables. Vectorized.")
492 (setf (documentation 'chisq-rand
'function
)
494 Returns a list of N Chi-Square(DF) random variables. Vectorized.")
496 (setf (documentation 'beta-rand
'function
)
498 Returns a list of N beta(A, B) random variables. Vectorized.")
500 (setf (documentation 'binomial-cdf
'function
)
502 Returns value of the Binomial(N, P) distribution function at X. Vectorized.")
504 (setf (documentation 'poisson-cdf
'function
)
506 Returns value of the Poisson(MU) distribution function at X. Vectorized.")
508 (setf (documentation 'binomial-pmf
'function
)
510 Returns value of the Binomial(N, P) pmf function at integer K. Vectorized.")
512 (setf (documentation 'poisson-pmf
'function
)
514 Returns value of the Poisson(MU) pmf function at integer K. Vectorized.")
516 (setf (documentation 'binomial-quant
'function
)
518 Returns x-th quantile (left continuous inverse) of Binomial(N, P) cdf.
521 (setf (documentation 'poisson-quant
'function
)
523 Returns x-th quantile (left continuous inverse) of Poisson(MU) cdf.
526 (setf (documentation 'binomial-rand
'function
)
528 Returns list of K draws from the Binomial(N, P) distribution. Vectorized.")
530 (setf (documentation 'poisson-rand
'function
)
532 Returns list of K draws from the Poisson(MU) distribution. Vectorized.")