source: project/release/4/srfi-27/trunk/srfi-27-distributions.scm @ 22477

Last change on this file since 22477 was 22477, checked in by Kon Lovett, 10 years ago

reciprocal not inverse.

File size: 12.2 KB
Line 
1;;;; srfi-27-distributions.scm
2;;;; Kon Lovett, May '06
3
4; Chicken Generic Arithmetic! (could use fp routines)
5
6(module srfi-27-distributions
7
8  (;export
9    make-random-exponentials
10    make-random-normals
11    make-random-triangles
12    make-random-poissons
13    make-random-bernoullis
14    make-random-binomials
15    make-random-geometrics
16    make-random-lognormals
17    make-random-cauchys
18    make-random-gammas
19    make-random-erlangs
20    make-random-paretos
21    make-random-levys
22    make-random-weibulls)
23
24  (import
25    scheme
26    chicken
27    (only type-errors
28      error-argument-type)
29    (only type-checks
30      check-procedure
31      check-cardinal-integer
32      check-real
33      check-open-interval
34      check-closed-interval)
35    (only srfi-27-uniform-random
36      make-uniform-random-reals))
37
38  (require-library
39    type-errors type-checks
40    srfi-27-uniform-random)
41
42;;; Chicken Generic Arithmetic Argument Checks
43
44(define (check-nonzero-real loc obj #!optional argnam)
45  (unless (and (real? obj) (not (zero? obj)))
46    (error-argument-type loc obj "nonzero-real" argnam)) )
47
48(define (check-nonnegative-real loc obj #!optional argnam)
49  (unless (and (real? obj) (not (negative? obj)))
50    (error-argument-type loc obj "nonnegative-real" argnam)) )
51
52(define (check-positive-real loc obj #!optional argnam)
53  (unless (and (real? obj) (positive? obj))
54    (error-argument-type loc obj "positive-real" argnam)) )
55
56(define (check-real-open-interval loc obj mn mx #!optional argnam)
57  (check-real loc obj argnam)
58  (check-open-interval loc obj mn mx argnam) )
59
60(define (check-real-closed-interval loc obj mn mx #!optional argnam)
61  (check-real loc obj argnam)
62  (check-closed-interval loc obj mn mx argnam) )
63
64#;
65(define (check-real-precision loc obj #!optional argnam)
66  (check-real-open-interval loc obj 0 1 argnam) )
67
68(define (check-real-unit loc obj #!optional argnam)
69  (check-real-closed-interval loc obj 0 1 argnam) )
70
71;;;
72
73(define-constant PI     3.1415926535897932384626433832795028841972)
74(define-constant FP1/3  0.3333333333333333333333333333333333333333)
75
76;;;
77
78; (in case special processing needed near limits TBD)
79(define-inline (*reciprocal n) (/ 1.0 n))
80(define-inline (*-reciprocal n) (/ -1.0 n))
81
82(define (fxadd1 n) (fx+ 1 n))
83
84;;; Normal distribution
85
86;; Knuth's "The Art of Computer Programming", Vol. II, 2nd ed.,
87;; Algorithm P of Section 3.4.1.C.
88
89(define (*make-random-normals mu sigma randoms)
90  (let ((next #f))
91    (lambda ()
92      (if next
93          (let ((result next))
94            (set! next #f)
95            (+ mu (* sigma result)))
96          (let loop ()
97            (let* ((v1 (- (* 2.0 (randoms)) 1.0))
98                   (v2 (- (* 2.0 (randoms)) 1.0))
99                   (s (+ (* v1 v1) (* v2 v2))))
100              (if (<= 1.0 s)
101                  (loop)
102                  (let ((scale (sqrt (/ (* -2.0 (log s)) s))))
103                    (set! next (* scale v2))
104                    (+ mu (* sigma scale v1))))))))) )
105
106(define (make-random-normals #!key (mu 0.0) (sigma 1.0) (randoms (make-uniform-random-reals)))
107  (check-real 'make-random-normals mu 'mu)
108  (check-nonzero-real 'make-random-normals sigma 'sigma)
109  (check-procedure 'make-random-normals randoms 'randoms)
110  (values
111    (*make-random-normals mu sigma randoms)
112    (lambda () (values mu sigma randoms))) )
113
114;;; Exponential distribution
115
116;; Knuth's "The Art of Computer Programming", Vol. II, 2nd ed.,
117;; Section 3.4.1.D.
118
119(define (*make-random-exponentials mu randoms)
120  (if (= 1.0 mu)
121      (lambda () (- (log (randoms))))
122      (lambda () (* mu (- (log (randoms)))))) )
123
124(define (make-random-exponentials #!key (mu 1.0) (randoms (make-uniform-random-reals)))
125  (check-real-unit 'make-random-exponentials mu 'mu)
126  (check-procedure 'make-random-exponentials randoms 'randoms)
127  (values
128    (*make-random-exponentials mu randoms)
129    (lambda () (values mu randoms))) )
130
131;;; Triangle distribution
132
133;; s - smallest, m - most probable, l - largest
134
135(define (*make-random-triangles s m l randoms)
136  (let ((d1 (- m s))
137        (d2 (- l s))
138        (d3 (sqrt (- l m))))
139    (let ((q1 (/ d1 d2))
140          (p1 (sqrt (* d1 d2))))
141      (lambda ()
142        (let ((u (randoms)))
143          (if (<= u q1)
144              (+ s (* p1 (sqrt u)))
145              (- l (* d3 (sqrt (- (* d2 u) d1))))))))) )
146
147(define (make-random-triangles #!key (s 0.0) (m 0.5) (l 1.0) (randoms (make-uniform-random-reals)))
148  (check-real 'make-random-triangles s 's)
149  (check-real 'make-random-triangles m 'm)
150  (check-real 'make-random-triangles l 'l)
151  (check-real-open-interval 'make-random-triangles l s +inf.0 'l)
152  (check-real-closed-interval 'make-random-triangles m s l 'm)
153  (check-procedure 'make-random-triangles randoms 'randoms)
154  (values
155    (*make-random-triangles s m l randoms)
156    (lambda () (values s m l randoms))) )
157
158;;; Poisson distribution
159
160(define (*make-random-poissons mu randoms)
161  (let ((emu (exp (- mu))))
162    (lambda ()
163      ; FIXME O(mu) but O(log(mu)) desired for >> mu
164      (do ((m 0 (fxadd1 m))
165           (prod (randoms) (* prod (randoms))))
166          ((<= prod emu) m)))) )
167
168(define (make-random-poissons #!key (mu 1.0) (randoms (make-uniform-random-reals)))
169  (check-nonnegative-real 'make-random-poissons mu 'mu)
170  (check-procedure 'make-random-poissons randoms 'randoms)
171  (values
172    (*make-random-poissons mu randoms)
173    (lambda () (values mu randoms))) )
174
175;;; Bernoulli distribution
176
177(define (*make-random-bernoullis p randoms)
178  (cond ((= 0.0 p) (lambda () #f))
179        ((= 1.0 p) (lambda () #t))
180        (else      (lambda () (<= (randoms) p)))) )
181
182(define (make-random-bernoullis #!key (p 0.5) (randoms (make-uniform-random-reals)))
183  (check-real-unit 'make-random-bernoullis p 'p)
184  (check-procedure 'make-random-bernoullis randoms 'randoms)
185  (values
186    (*make-random-bernoullis p randoms)
187    (lambda () (values p randoms))) )
188
189;;; Binomial distribution
190
191(define (*make-random-binomials t p randoms)
192  (let ((bernoullis (*make-random-bernoullis p randoms)))
193    ;FIXME O(t) but O(log(t)) desired for >> t
194    (if (fixnum? t)
195        (lambda ()
196          (do ((i 0 (fxadd1 i))
197               (n 0 (if (bernoullis) (fxadd1 n) n)))
198              ((fx<= t i) n)))
199        (lambda ()
200          (do ((i 0 (add1 i))
201               (n 0 (if (bernoullis) (add1 n) n)))
202              ((<= t i) n))))) )
203
204(define (make-random-binomials #!key (t 1) (p 0.5) (randoms (make-uniform-random-reals)))
205  (check-cardinal-integer 'make-random-binomials t 't)
206  (check-real-unit 'make-random-binomials p 'p)
207  (check-procedure 'make-random-binomials randoms 'randoms)
208  (values
209    (*make-random-binomials t p randoms)
210    (lambda () (values t p randoms))) )
211
212;;; Geometric distribution
213
214(define (*make-random-geometrics p randoms)
215  (let ((log-p (log p)))
216    (lambda () (+ 1 (inexact->exact (floor (/ (log (- 1.0 (randoms))) log-p)))))) )
217
218(define (make-random-geometrics #!key (p 0.5) (randoms (make-uniform-random-reals)))
219  (check-real-unit 'make-random-geometrics p 'p)
220  (check-procedure 'make-random-geometrics randoms 'randoms)
221  (values
222    (*make-random-geometrics p randoms)
223    (lambda () (values p randoms))) )
224
225;;; Lognormal distribution
226
227(define (*make-random-lognormals mu sigma randoms)
228  (let ((normals (*make-random-normals 0.0 1.0 randoms))
229        (nmu (log (* mu (/ mu (sqrt (+ (* sigma sigma) (* mu mu)))))))
230        (nsigma (sqrt (log (+ 1.0 (* sigma (/ sigma mu mu)))))) )
231    (lambda () (exp (+ nmu (* (normals) nsigma))))) )
232
233(define (make-random-lognormals #!key (mu 1.0) (sigma 1.0) (randoms (make-uniform-random-reals)))
234  (check-nonzero-real 'make-random-lognormals mu 'mu)
235  (check-nonnegative-real 'make-random-lognormals sigma 'sigma)
236  (check-procedure 'make-random-lognormals randoms 'randoms)
237  (values
238    (*make-random-lognormals mu sigma randoms)
239    (lambda () (values mu sigma randoms))) )
240
241;;; Cauchy distribution
242
243(define (*make-random-cauchys median sigma randoms)
244  (lambda () (+ median (* sigma (tan (* PI (- (randoms) 0.5)))))) )
245
246(define (make-random-cauchys #!key (median 0.0) (sigma 1.0) (randoms (make-uniform-random-reals)))
247  (check-real 'make-random-cauchys median 'median)
248  (check-positive-real 'make-random-cauchys sigma 'sigma)
249  (check-procedure 'make-random-cauchys randoms 'randoms)
250  (values
251    (*make-random-cauchys median sigma randoms)
252    (lambda () (values median sigma randoms))) )
253
254;;; Gamma distribution
255
256;; "A Simple Method for Generating Gamma Variables", George Marsaglia & Wai Wan Tsang,
257;; ACM Transactions on Mathematical Software, Vol. 26, No. 3, September 2000, Pages 363 372.
258
259(define (*make-random-gammas alpha theta randoms)
260  (if (= 1.0 alpha)
261      ; then special case
262      (lambda () (* theta (- (log (randoms)))) )
263      ; else general case
264      (let ((norms (*make-random-normals 0.0 1.0 randoms))
265            (unis
266              (if (< alpha 1.0)
267                  (let ((inv-alpha (*reciprocal alpha)))
268                    (lambda () (expt (randoms) inv-alpha) ) )
269                  randoms)))
270        (let* ((d (- (or (and (< alpha 1.0) (+ 1.0 alpha)) alpha) FP1/3))
271               (c (*reciprocal (sqrt (* 9.0 d)))))
272          (lambda ()
273            (* theta
274               (let loop ()
275                 (let* ((x (norms))
276                        (v (+ 1.0 (* c x))))
277                   (if (and (< 0.0 v)
278                            (let ((v (* v v v))
279                                  (u (unis))
280                                  (x^2 (* x x)))
281                              (or (< u (- 1.0 (* 0.0331 x^2 x^2)))
282                                  (< (log u) (+ (* 0.5 x^2) (* d (- 1.0 (+ v (log v)))))))))
283                       (* d v)
284                       (loop) ) ) ) ) ) ) ) ) )
285
286(define (make-random-gammas #!key (alpha 1.0) (theta 1.0) (randoms (make-uniform-random-reals)))
287  (check-positive-real 'make-random-gammas alpha 'alpha)
288  (check-positive-real 'make-random-gammas theta 'theta)
289  (check-procedure 'make-random-gammas randoms 'randoms)
290  (values
291    (*make-random-gammas alpha theta randoms)
292    (lambda () (values alpha theta randoms))) )
293
294;;; Erlang distribution
295
296(define (*make-random-erlangs alpha theta randoms)
297  (*make-random-gammas (exact->inexact alpha) (exact->inexact theta) randoms) )
298
299(define (make-random-erlangs #!key (alpha 1) (theta 1.0) (randoms (make-uniform-random-reals)))
300  (check-positive-real 'make-random-erlangs alpha 'alpha)
301  (check-positive-real 'make-random-erlangs theta 'theta)
302  (check-procedure 'make-random-erlangs randoms 'randoms)
303  (values
304    (*make-random-erlangs alpha theta randoms)
305    (lambda () (values alpha theta randoms))) )
306
307;;; Pareto distribution
308
309(define (*make-random-paretos alpha xmin randoms)
310  (let ((gammas (*make-random-gammas alpha (*reciprocal xmin) randoms)))
311    (*make-random-exponentials 1.0 (lambda () (*reciprocal (+ xmin (gammas)))))) )
312
313(define (make-random-paretos #!key (alpha 1.0) (xmin 1.0) (randoms (make-uniform-random-reals)))
314  (check-positive-real 'make-random-paretos alpha 'alpha)
315  (check-positive-real 'make-random-paretos xmin 'xmin)
316  (check-procedure 'make-random-paretos randoms 'randoms)
317  (values
318    (*make-random-paretos alpha xmin randoms)
319    (lambda () (values alpha xmin randoms))) )
320
321;;; Levy distribution
322
323;; See Stable Distributions - John P. Nolan, Formula 1.12
324
325(define (*make-random-levys gamma delta randoms)
326  (if (and (= 1.0 gamma) (= 0.0 delta))
327      (lambda () (let ((r (randoms))) (*reciprocal (* r r))))
328      (lambda () (let ((r (randoms))) (+ delta (* gamma (*reciprocal (* r r))))))) )
329
330(define (make-random-levys #!key (gamma 1.0) (delta 0.0) (randoms (make-uniform-random-reals)))
331  (check-nonnegative-real 'make-random-levys delta 'delta)
332  (check-positive-real 'make-random-levys gamma 'gamma)
333  (check-procedure 'make-random-levys randoms 'randoms)
334  (values
335    (*make-random-levys gamma delta randoms)
336    (lambda () (values gamma delta randoms))) )
337
338;;; Weibull distribution
339
340(define (*make-random-weibulls shape scale randoms)
341  (let ((invscale (*-reciprocal scale))
342        (invshape (*reciprocal shape)) )
343    (lambda () (expt (* invscale (log (- 1.0 (randoms)))) invshape)) ) )
344
345(define (make-random-weibulls #!key (shape 1.0) (scale 1.0) (randoms (make-uniform-random-reals)))
346  (check-positive-real 'make-random-weibulls shape 'shape)
347  (check-positive-real 'make-random-weibulls scale 'scale)
348  (check-procedure 'make-random-weibulls randoms 'randoms)
349  (values
350    (*make-random-weibulls shape scale randoms)
351    (lambda () (values shape scale randoms))) )
352
353) ;module srfi-27-distributions
Note: See TracBrowser for help on using the repository browser.