3 ;;; Copyright (c) 2005--2008, by A.J. Rossini <blindglobe@gmail.com>
4 ;;; See COPYRIGHT file for any additional restrictions (BSD license).
5 ;;; Since 1991, ANSI was finally finished. Edited for ANSI Common Lisp.
7 ;;; dists -- Lisp-Stat interface to basic probability distribution routines
9 ;;; Copyright (c) 1991, by Luke Tierney. Permission is granted for
12 (in-package :lisp-stat-probability
)
14 ;;; This stuff needs to be improved. We really could use something
15 ;;; like the R libraries, but of course in a better packaged manner.
17 ;;; Currently, there is a function for everything. Probably better to
18 ;;; simplify by thinking about a more generic approach, distribution
19 ;;; being specified by keyword.
23 "stupid dummy function, need to implement rng seeding tool."
26 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
28 ;;; CFFI support for Probability Distributions
30 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
33 ;;; C-callable uniform generator
36 (defcfun ("register_uni" register-uni
)
38 (defcallback ccl-uni
:int
() (ccl-store-double (random 1.0)) 0)
39 (register-uni (callback ccl-uni
))
41 (defun one-uniform-rand () (random 1.0))
44 ;;; Log-gamma function
47 (defcfun ("ccl_gamma" ccl-base-log-gamma
)
49 (defun base-log-gamma (x)
50 (ccl-base-log-gamma (float x
1d0
)))
53 ;;; Normal distribution
56 (defcfun ("ccl_normalcdf" ccl-base-normal-cdf
)
58 (defun base-normal-cdf (x)
59 (ccl-base-normal-cdf (float x
1d0
)))
61 (defcfun ("ccl_normalquant" ccl-base-normal-quant
)
63 (defun base-normal-quant (x)
64 (ccl-base-normal-quant (float x
1d0
)))
66 (defcfun ("ccl_normaldens" ccl-base-normal-dens
)
68 (defun base-normal-dens (x)
69 (ccl-base-normal-dens (float x
1d0
)))
71 (defcfun ("ccl_normalrand" one-normal-rand
)
74 (defcfun ("ccl_bnormcdf" ccl-base-bivnorm-cdf
)
75 :double
(x :double
) (y :double
) (z :double
))
76 (defun base-bivnorm-cdf (x y z
)
77 (ccl-base-bivnorm-cdf (float x
1d0
) (float y
1d0
) (float z
1d0
)))
80 ;;; Cauchy distribution
83 (defcfun ("ccl_cauchycdf" ccl-base-cauchy-cdf
)
85 (defun base-cauchy-cdf (x)
86 (ccl-base-cauchy-cdf (float x
1d0
)))
88 (defcfun ("ccl_cauchyquant" ccl-base-cauchy-quant
)
90 (defun base-cauchy-quant (x)
91 (ccl-base-cauchy-quant (float x
1d0
)))
93 (defcfun ("ccl_cauchydens" ccl-base-cauchy-dens
)
95 (defun base-cauchy-dens (x)
96 (ccl-base-cauchy-dens (float x
1d0
)))
98 (defcfun ("ccl_cauchyrand" one-cauchy-rand
)
102 ;;;; Gamma distribution
105 (defcfun ("ccl_gammacdf" ccl-base-gamma-cdf
)
106 :double
(x :double
) (y :double
))
107 (defun base-gamma-cdf (x y
)
108 (ccl-base-gamma-cdf (float x
1d0
) (float y
1d0
)))
110 (defcfun ("ccl_gammaquant" ccl-base-gamma-quant
)
111 :double
(x :double
) (y :double
))
112 (defun base-gamma-quant (x y
)
113 (ccl-base-gamma-quant (float x
1d0
) (float y
1d0
)))
115 (defcfun ("ccl_gammadens" ccl-base-gamma-dens
)
116 :double
(x :double
) (y :double
))
117 (defun base-gamma-dens (x y
)
118 (ccl-base-gamma-dens (float x
1d0
) (float y
1d0
)))
120 (defcfun ("ccl_gammarand" ccl-gamma-rand
)
122 (defun one-gamma-rand (x)
123 (ccl-gamma-rand (float x
1d0
)))
126 ;;;; Chi-square distribution
129 (defcfun ("ccl_chisqcdf" ccl-base-chisq-cdf
)
130 :double
(x :double
) (y :double
))
131 (defun base-chisq-cdf (x y
)
132 (ccl-base-chisq-cdf (float x
1d0
) (float y
1d0
)))
134 (defcfun ("ccl_chisqquant" ccl-base-chisq-quant
)
135 :double
(x :double
) (y :double
))
136 (defun base-chisq-quant (x y
)
137 (ccl-base-chisq-quant (float x
1d0
) (float y
1d0
)))
139 (defcfun ("ccl_chisqdens" ccl-base-chisq-dens
)
140 :double
(x :double
) (y :double
))
141 (defun base-chisq-dens (x y
)
142 (ccl-base-chisq-dens (float x
1d0
) (float y
1d0
)))
144 (defcfun ("ccl_chisqrand" ccl-chisq-rand
)
146 (defun one-chisq-rand (x)
147 (ccl-chisq-rand (float x
1d0
)))
150 ;;;; Beta distribution
153 (defcfun ("ccl_betacdf" ccl-base-beta-cdf
)
154 :double
(x :double
) (y :double
) (z :double
))
155 (defun base-beta-cdf (x y z
)
156 (ccl-base-beta-cdf (float x
1d0
) (float y
1d0
) (float z
1d0
)))
158 (defcfun ("ccl_betaquant" ccl-base-beta-quant
)
159 :double
(x :double
) (y :double
) (z :double
))
160 (defun base-beta-quant (x y z
)
161 (ccl-base-beta-quant (float x
1d0
) (float y
1d0
) (float z
1d0
)))
163 (defcfun ("ccl_betadens" ccl-base-beta-dens
)
164 :double
(x :double
) (y :double
) (z :double
))
165 (defun base-beta-dens (x y z
)
166 (ccl-base-beta-dens (float x
1d0
) (float y
1d0
) (float z
1d0
)))
168 (defcfun ("ccl_betarand" ccl-beta-rand
)
169 :double
(x :double
) (y :double
))
170 (defun one-beta-rand (x y
)
171 (ccl-beta-rand (float x
1d0
) (float y
1d0
)))
177 (defcfun ("ccl_tcdf" ccl-base-t-cdf
)
178 :double
(x :double
) (y :double
))
179 (defun base-t-cdf (x y
)
180 (ccl-base-t-cdf (float x
1d0
) (float y
1d0
)))
182 (defcfun ("ccl_tquant" ccl-base-t-quant
)
183 :double
(x :double
) (y :double
))
184 (defun base-t-quant (x y
)
185 (ccl-base-t-quant (float x
1d0
) (float y
1d0
)))
187 (defcfun ("ccl_tdens" ccl-base-t-dens
)
188 :double
(x :double
) (y :double
))
189 (defun base-t-dens (x y
)
190 (ccl-base-t-dens (float x
1d0
) (float y
1d0
)))
192 (defcfun ("ccl_trand" ccl-t-rand
)
194 (defun one-t-rand (x)
195 (ccl-t-rand (float x
1d0
)))
201 (defcfun ("ccl_fcdf" ccl-base-f-cdf
)
202 :double
(x :double
) (y :double
) (z :double
))
203 (defun base-f-cdf (x y z
)
204 (ccl-base-f-cdf (float x
1d0
) (float y
1d0
) (float z
1d0
)))
206 (defcfun ("ccl_fquant" ccl-base-f-quant
)
207 :double
(x :double
) (y :double
) (z :double
))
208 (defun base-f-quant (x y z
)
209 (ccl-base-f-quant (float x
1d0
) (float y
1d0
) (float z
1d0
)))
211 (defcfun ("ccl_fdens" ccl-base-f-dens
)
212 :double
(x :double
) (y :double
) (z :double
))
213 (defun base-f-dens (x y z
)
214 (ccl-base-f-dens (float x
1d0
) (float y
1d0
) (float z
1d0
)))
216 (defcfun ("ccl_frand" ccl-f-rand
)
217 :double
(x :double
) (y :double
))
218 (defun one-f-rand (x y
) (ccl-f-rand (float x
1d0
) (float y
1d0
)))
221 ;;;; Poisson distribution
224 (defcfun ("ccl_poissoncdf" ccl-base-poisson-cdf
)
225 :double
(x :double
) (y :double
))
226 (defun base-poisson-cdf (x y
)
227 (ccl-base-poisson-cdf (float x
1d0
) (float y
1d0
)))
229 (defcfun ("ccl_poissonquant" ccl-base-poisson-quant
)
230 :int
(x :double
) (y :double
))
231 (defun base-poisson-quant (x y
)
232 (ccl-base-poisson-quant (float x
1d0
) (float y
1d0
)))
234 (defcfun ("ccl_poissonpmf" ccl-base-poisson-pmf
)
235 :double
(x :int
) (y :double
))
236 (defun base-poisson-pmf (x y
)
237 (ccl-base-poisson-pmf x
(float y
1d0
)))
239 (defcfun ("ccl_poissonrand" ccl-poisson-rand
)
241 (defun one-poisson-rand (x)
242 (ccl-poisson-rand (float x
1d0
)))
245 ;;;; Binomial distribution
248 (defcfun ("ccl_binomialcdf" ccl-base-binomial-cdf
)
249 :double
(x :double
) (y :int
) (z :double
))
250 (defun base-binomial-cdf (x y z
)
251 (ccl-base-binomial-cdf (float x
1d0
) y
(float z
1d0
)))
253 (defcfun ("ccl_binomialquant" ccl-base-binomial-quant
)
254 :int
(x :double
) (y :int
) (z :double
))
255 (defun base-binomial-quant (x y z
)
256 (ccl-base-binomial-quant (float x
1d0
) y
(float z
1d0
)))
258 (defcfun ("ccl_binomialpmf" ccl-base-binomial-pmf
)
259 :double
(x :int
) (y :int
) (z :double
))
260 (defun base-binomial-pmf (x y z
)
261 (ccl-base-binomial-pmf x y
(float z
1d0
)))
263 (defcfun ("ccl_binomialrand" ccl-binomial-rand
)
264 :int
(x :int
) (y :double
))
265 (defun one-binomial-rand (x y
)
266 (ccl-binomial-rand x
(float y
1d0
)))
272 ;;; definitions though macros
274 (defmacro defbaserand
(name onefun
&rest args
)
275 `(defun ,name
(n ,@args
)
277 (dotimes (i n result
)
278 (declare (fixnum i
) (inline ,onefun
))
279 (setf result
(cons (,onefun
,@args
) result
))))))
281 (defbaserand base-uniform-rand one-uniform-rand
)
282 (defbaserand base-normal-rand one-normal-rand
)
283 (defbaserand base-cauchy-rand one-cauchy-rand
)
284 (defbaserand base-gamma-rand one-gamma-rand a
)
285 (defbaserand base-chisq-rand one-chisq-rand df
)
286 (defbaserand base-beta-rand one-beta-rand a b
)
287 (defbaserand base-t-rand one-t-rand df
)
288 (defbaserand base-f-rand one-f-rand ndf ddf
)
289 (defbaserand base-poisson-rand one-poisson-rand a
)
290 (defbaserand base-binomial-rand one-binomial-rand a b
)
292 (make-rv-function log-gamma base-log-gamma x
)
294 (make-rv-function uniform-rand base-uniform-rand n
)
296 (make-rv-function normal-cdf base-normal-cdf x
)
297 (make-rv-function normal-quant base-normal-quant p
)
298 (make-rv-function normal-dens base-normal-dens x
)
299 (make-rv-function normal-rand base-normal-rand n
)
300 (make-rv-function bivnorm-cdf base-bivnorm-cdf x y r
)
302 (make-rv-function cauchy-cdf base-cauchy-cdf x
)
303 (make-rv-function cauchy-quant base-cauchy-quant p
)
304 (make-rv-function cauchy-dens base-cauchy-dens x
)
305 (make-rv-function cauchy-rand base-cauchy-rand n
)
307 (make-rv-function gamma-cdf base-gamma-cdf x a
)
308 (make-rv-function gamma-quant base-gamma-quant p a
)
309 (make-rv-function gamma-dens base-gamma-dens x a
)
310 (make-rv-function gamma-rand base-gamma-rand n a
)
312 (make-rv-function chisq-cdf base-chisq-cdf x df
)
313 (make-rv-function chisq-quant base-chisq-quant p df
)
314 (make-rv-function chisq-dens base-chisq-dens x df
)
315 (make-rv-function chisq-rand base-chisq-rand n df
)
317 (make-rv-function beta-cdf base-beta-cdf x a b
)
318 (make-rv-function beta-quant base-beta-quant p a b
)
319 (make-rv-function beta-dens base-beta-dens x a b
)
320 (make-rv-function beta-rand base-beta-rand n a b
)
322 (make-rv-function t-cdf base-t-cdf x df
)
323 (make-rv-function t-quant base-t-quant p df
)
324 (make-rv-function t-dens base-t-dens x df
)
325 (make-rv-function t-rand base-t-rand n df
)
327 (make-rv-function f-cdf base-f-cdf x ndf ddf
)
328 (make-rv-function f-quant base-f-quant p ndf ddf
)
329 (make-rv-function f-dens base-f-dens x ndf ddf
)
330 (make-rv-function f-rand base-f-rand n ndf ddf
)
332 (make-rv-function poisson-cdf base-poisson-cdf x a
)
333 (make-rv-function poisson-quant base-poisson-quant p a
)
334 (make-rv-function poisson-pmf base-poisson-pmf x a
)
335 (make-rv-function poisson-rand base-poisson-rand n a
)
337 (make-rv-function binomial-cdf base-binomial-cdf x a b
)
338 (make-rv-function binomial-quant base-binomial-quant p a b
)
339 (make-rv-function binomial-pmf base-binomial-pmf x a b
)
340 (make-rv-function binomial-rand base-binomial-rand n a b
)
346 (setf (documentation 'bivnorm-cdf
'function
)
348 Returns the value of the standard bivariate normal distribution function
349 with correlation R at (X, Y). Vectorized.")
351 (setf (documentation 'normal-cdf
'function
)
353 Returns the value of the standard normal distribution function at X.
356 (setf (documentation 'beta-cdf
'function
)
357 "Args: (x alpha beta)
358 Returns the value of the Beta(ALPHA, BETA) distribution function at X.
361 (setf (documentation 'gamma-cdf
'function
)
363 Returns the value of the Gamma(alpha, 1) distribution function at X.
366 (setf (documentation 'chisq-cdf
'function
)
368 Returns the value of the Chi-Square(DF) distribution function at X. Vectorized.")
370 (setf (documentation 't-cdf
'function
)
372 Returns the value of the T(DF) distribution function at X. Vectorized.")
374 (setf (documentation 'f-cdf
'function
)
376 Returns the value of the F(NDF, DDF) distribution function at X. Vectorized.")
378 (setf (documentation 'cauchy-cdf
'function
)
380 Returns the value of the standard Cauchy distribution function at X.
383 (setf (documentation 'log-gamma
'function
)
385 Returns the log gamma function of X. Vectorized.")
387 (setf (documentation 'normal-quant
'function
)
389 Returns the P-th quantile of the standard normal distribution. Vectorized.")
391 (setf (documentation 'cauchy-quant
'function
)
393 Returns the P-th quantile(s) of the standard Cauchy distribution. Vectorized.")
395 (setf (documentation 'beta-quant
'function
)
396 "Args: (p alpha beta)
397 Returns the P-th quantile of the Beta(ALPHA, BETA) distribution. Vectorized.")
399 (setf (documentation 'gamma-quant
'function
)
401 Returns the P-th quantile of the Gamma(ALPHA, 1) distribution. Vectorized.")
403 (setf (documentation 'chisq-quant
'function
)
405 Returns the P-th quantile of the Chi-Square(DF) distribution. Vectorized.")
407 (setf (documentation 't-quant
'function
)
409 Returns the P-th quantile of the T(DF) distribution. Vectorized.")
411 (setf (documentation 'f-quant
'function
)
413 Returns the P-th quantile of the F(NDF, DDF) distribution. Vectorized.")
415 (setf (documentation 'normal-dens
'function
)
417 Returns the density at X of the standard normal distribution. Vectorized.")
419 (setf (documentation 'cauchy-dens
'function
)
421 Returns the density at X of the standard Cauchy distribution. Vectorized.")
423 (setf (documentation 'beta-dens
'function
)
424 "Args: (x alpha beta)
425 Returns the density at X of the Beta(ALPHA, BETA) distribution. Vectorized.")
427 (setf (documentation 'gamma-dens
'function
)
429 Returns the density at X of the Gamma(ALPHA, 1) distribution. Vectorized.")
431 (setf (documentation 'chisq-dens
'function
)
433 Returns the density at X of the Chi-Square(DF) distribution. Vectorized.")
435 (setf (documentation 't-dens
'function
)
437 Returns the density at X of the T(DF) distribution. Vectorized.")
439 (setf (documentation 'f-dens
'function
)
441 Returns the density at X of the F(NDF, DDF) distribution. Vectorized.")
443 (setf (documentation 'uniform-rand
'function
)
445 Returns a list of N uniform random variables from the range (0, 1).
448 (setf (documentation 'normal-rand
'function
)
450 Returns a list of N standard normal random numbers. Vectorized.")
452 (setf (documentation 'cauchy-rand
'function
)
454 Returns a list of N standard Cauchy random numbers. Vectorized.")
456 (setf (documentation 't-rand
'function
)
458 Returns a list of N T(DF) random variables. Vectorized.")
460 (setf (documentation 'f-rand
'function
)
462 Returns a list of N F(NDF, DDF) random variables. Vectorized.")
464 (setf (documentation 'gamma-rand
'function
)
466 Returns a list of N Gamma(A, 1) random variables. Vectorized.")
468 (setf (documentation 'chisq-rand
'function
)
470 Returns a list of N Chi-Square(DF) random variables. Vectorized.")
472 (setf (documentation 'beta-rand
'function
)
474 Returns a list of N beta(A, B) random variables. Vectorized.")
476 (setf (documentation 'binomial-cdf
'function
)
478 Returns value of the Binomial(N, P) distribution function at X. Vectorized.")
480 (setf (documentation 'poisson-cdf
'function
)
482 Returns value of the Poisson(MU) distribution function at X. Vectorized.")
484 (setf (documentation 'binomial-pmf
'function
)
486 Returns value of the Binomial(N, P) pmf function at integer K. Vectorized.")
488 (setf (documentation 'poisson-pmf
'function
)
490 Returns value of the Poisson(MU) pmf function at integer K. Vectorized.")
492 (setf (documentation 'binomial-quant
'function
)
494 Returns x-th quantile (left continuous inverse) of Binomial(N, P) cdf.
497 (setf (documentation 'poisson-quant
'function
)
499 Returns x-th quantile (left continuous inverse) of Poisson(MU) cdf.
502 (setf (documentation 'binomial-rand
'function
)
504 Returns list of K draws from the Binomial(N, P) distribution. Vectorized.")
506 (setf (documentation 'poisson-rand
'function
)
508 Returns list of K draws from the Poisson(MU) distribution. Vectorized.")