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

Last change on this file since 33848 was 33848, checked in by Kon Lovett, 4 years ago

re-flow, add "catalog" supp

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