source: project/srfi-27/srfi-27-distributions.scm @ 2445

Last change on this file since 2445 was 2445, checked in by Kon Lovett, 15 years ago

Bug fixes.

File size: 13.1 KB
Line 
1;;;; srfi-27-distributions.scm
2;;;; Kon Lovett, May '06
3
4(use srfi-4 srfi-27-structures srfi-27-parameters vector-lib structures misc-extn)
5
6(eval-when (compile)
7        (declare
8                (not usual-integrations
9                        number? real? integer? positive? zero? negative?
10                        = >= <= < >
11                        inexact->exact
12                        exp log sqrt floor tan expt
13                        + - * /)
14                (inline)
15                (export
16                        make-uniform-random-integers
17                        make-uniform-random-reals
18                        make-random-exponentials
19                        make-random-normals
20                        make-random-triangles
21                        make-random-poissons
22                        make-random-bernoullis
23                        make-random-binomials
24                        make-random-geometrics
25                        make-random-lognormals
26                        make-random-cauchys
27                        make-random-gammas
28                        make-random-erlangs
29                        make-random-paretos
30                        make-random-levys
31                        random-normal-vector!
32                        random-hollow-sphere!
33                        random-solid-sphere!)
34        )
35)
36
37;;;
38
39; (could include mathh-constants here)
40
41(define-constant PI     3.1415926535897932384626433832795028841972)
42(define-constant E      2.7182818284590452353602874713526624977572)
43
44(define-constant one-third 0.3333333333333333333333333333333333333333)
45
46;;;
47
48; (in case special processing needed near limits)
49
50(define-inline (inv n)
51        (/ 1.0 n) )
52
53;;; Vector Support
54
55(define (f64vector-map! proc vec)
56        (let ([len (f64vector-length vec)])
57                (let loop ([i 0])
58                        (if (fx= i len)
59                                vec
60                                (begin
61                                        (f64vector-set! vec (proc i (f64vector-ref vec i)) i)
62                                        (loop (fx+ 1)))))) )
63
64(define (f64vector-fold proc state vec)
65        (let ([len (f64vector-length vec)])
66                (let loop ([i 0] [state state])
67                        (if (fx= i len)
68                                state
69                                (loop (fx+ 1) (proc i state (f64vector-ref vec i)))))) )
70
71(define (vector*-length vec)
72        (cond
73                [(vector? vec)
74                        (vector-length vec)]
75                [(f64vector? vec)
76                        (f64vector-length vec)]
77                #;[(and (array? vec) (fx= 1 (array-rank vec)))
78                        (car (array-dimensions vec))]
79                [else
80                        (error "not a proper vector" vec)]) )
81
82(define (vector*-map! proc vec)
83        (cond
84                [(vector? vec)
85                        (vector-map! proc vec)]
86                [(f64vector? vec)
87                        (f64vector-map! proc vec)]
88                #;[(and (array? vec) (fx= 1 (array-rank vec)))
89                        (array-map! vec (cut proc #f <>))]
90                [else
91                        (error "not a proper vector" vec)]) )
92
93(define (vector*-fold proc seed vec)
94        (cond
95                [(vector? vec)
96                        (vector-fold proc seed vec)]
97                [(f64vector? vec)
98                        (f64vector-fold proc seed vec)]
99                #;[(and (array? vec) (fx= 1 (array-rank vec)))
100                        (array-fold (cut proc #f <> <>) seed vec)]
101                [else
102                        (error "not a proper vector" vec)])     )
103
104(define-inline (vector*-scale! vec const)
105        (vector*-map! (lambda (i elm) (* elm const)) vec) )
106
107(define-inline (vector*-sum-squares vec)
108        (vector*-fold (lambda (i sum elm) (+ sum (* elm elm))) 0 vec) )
109
110;;; Uniform random integers in [low high]
111
112(define (make-uniform-random-integers #!optional (high #f) (low 0) (unit 1) (src (current-random-source-structure)))
113        (let (
114                        [make-rand
115                                (lambda (stc src)
116                                        (in-structure stc
117                                                ((import random-source-make-integers random-source-maximum-range))
118                                                (unless high
119                                                        (set! high (- (random-source-maximum-range src) 1)))
120                                                (let ([rand (random-source-make-integers src)]
121                                                                        [range (/ (+ (- high low) 1) unit)])
122                                                        (values
123                                                                (cond
124                                                                        [(and (= 0 low) (= 1 unit))
125                                                                                (lambda ()
126                                                                                        (rand range))]
127                                                                        [(= 0 low)
128                                                                                (lambda ()
129                                                                                        (* (rand range) unit))]
130                                                                        [else
131                                                                                (lambda ()
132                                                                                        (+ low (* (rand range) unit)))])
133                                                                (lambda ()
134                                                                        (values high low unit src))))))])
135                (unless (or (number? high) (boolean? high))
136                        (set! src high)
137                        (set! high #f))
138                (unless (number? low)
139                        (set! src low)
140                        (set! low 0))
141                (unless (number? unit)
142                        (set! src unit)
143                        (set! unit 1))
144                (cond
145                        [(%random-source? src)
146                                (make-rand (current-random-source-structure) src)]
147                        [(random-source-structure? src)
148                                (in-structure src
149                                        ((import default-random-source))
150                                        (make-rand src default-random-source))]
151                        [else
152                                (error 'make-uniform-random-integers "unknown random source" src)])) )
153
154;;; Uniform random reals in (0.0 1.0)
155
156(define (make-uniform-random-reals #!optional (unit #f) (src (current-random-source-structure)))
157        (cond
158                [(number? unit)
159                        (set! unit (list unit))]
160                [(boolean? unit)
161                        (set! unit '())]
162                [else
163                        (set! src unit)
164                        (set! unit '())])
165        (cond
166                [(%random-source? src)
167                        (in-structure (current-random-source-structure)
168                                ((import random-source-make-reals))
169                                (values
170                                        (apply random-source-make-reals src unit)
171                                        (lambda ()
172                                                (values unit src))))]
173                [(random-source-structure? src)
174                        (in-structure src
175                                ((import random-source-make-reals default-random-source))
176                                (values
177                                        (apply random-source-make-reals default-random-source unit)
178                                        (lambda ()
179                                                (values unit src))))]
180                [else
181                        (error 'make-uniform-random-reals "unknown random source" src)]) )
182
183;;; Exponentials
184
185;; Knuth's "The Art of Computer Programming", Vol. II, 2nd ed.,
186;; Section 3.4.1.D.
187;;
188(define (make-random-exponentials #!optional (mu 1.0) (rand (make-uniform-random-reals)))
189        (unless (and (real? mu) (<= 0.0 mu 1.0))
190                (error 'make-random-exponentials "mu must be a real in [0.0 1.0]" mu))
191        (values
192                (if (= 1.0 mu)
193                        (lambda () (- (log (rand))))
194                        (lambda () (* mu (- (log (rand))))))
195                (lambda ()
196                        (values mu rand))) )
197
198;;; Normals
199
200;; Knuth's "The Art of Computer Programming", Vol. II, 2nd ed.,
201;; Algorithm P of Section 3.4.1.C.
202
203(define (make-random-normals #!optional (mu 0.0) (sigma 1.0) (rand (make-uniform-random-reals)))
204        (unless (real? mu)
205                (error 'make-random-normals "mu must be a real" mu))
206        (unless (and (real? sigma) (not (zero? sigma)))
207                (error 'make-random-normals "sigma must be a non-zero real" sigma))
208  (let ([next #f])
209                (values
210                        (lambda ()
211                                (if next
212                                        (let ([result next])
213                                                (set! next #f)
214                                                (+ mu (* sigma result)))
215                                        (let loop ()
216                                                (let* ([v1 (- (* 2 (rand)) 1.0)]
217                                                                         [v2 (- (* 2 (rand)) 1.0)]
218                                                                         [s (+ (* v1 v1) (* v2 v2))])
219                                                        (if (>= s 1.0)
220                                                                (loop)
221                                                                (let ([scale (sqrt (/ (* -2 (log s)) s))])
222                                                                        (set! next (* scale v2))
223                                                                        (+ mu (* sigma scale v1))))))))
224                        (lambda ()
225                                (values mu sigma rand)))) )
226
227;;; Triangles
228
229;; s - smallest, m - most probable, l - largest
230
231(define (make-random-triangles #!optional (s 0.0) (m 0.5) (l 1.0) (rand (make-uniform-random-reals)))
232        (unless (real? s)
233                (error 'make-random-triangles "s must be a real" s))
234        (unless (and (real? l) (< s l))
235                (error 'make-random-triangles "l must be a real in (s +inf)" l))
236        (unless (and (real? m) (<= s m l))
237                (error 'make-random-triangles "m must be a real in [s l]" m))
238        (let ([d1 (- m s)]
239                                [d2 (- l s)]
240                                [d3 (sqrt (- l m))])
241                (let ([q1 (/ d1 d2)]
242                                        [p1 (sqrt (* d1 d2))])
243                        (values
244                                (lambda ()
245                                        (let ([u (rand)])
246                                                (if (<= u q1)
247                                                        (+ s (* p1 (sqrt u)))
248                                                        (- l (* d3 (sqrt (- (* d2 u) d1)))))))
249                                (lambda ()
250                                        (values s m l rand))))) )
251
252;;; Poisson
253
254(define (make-random-poissons #!optional (mu 1.0) (rand (make-uniform-random-reals)))
255        (unless (and (real? mu) (not (negative? mu)))
256                (error 'make-random-poissons "mu must be a non-negative real" mu))
257        (let ([emu (exp (- mu))])
258                (values
259                        (lambda ()
260                                ; FIXME O(mu) but O(log(mu)) desired for >> mu
261                                (do ([m 0 (fx+ 1 m)]
262                                                 [prod (rand) (* prod (rand))])
263                                                ([<= prod emu] m)))
264                        (lambda ()
265                                (values mu rand)))) )
266
267;;; Bernoulli
268
269(define (make-random-bernoullis #!optional (p 0.5) (rand (make-uniform-random-reals)))
270        (unless (and (real? p) (<= 0.0 p 1.0))
271                (error 'make-random-bernoullis "p must be a real in [0.0 1.0]" p))
272        (values
273                (cond
274                        [(zero? p) (lambda () #f)]
275                        [(= 1.0 p) (lambda () #t)]
276                        [else (lambda () (<= (rand) p))])
277                (lambda ()
278                        (values p rand))) )
279
280;;; Binomials
281
282(define-constant LARGE-FIXNUM 1073741823)
283
284(define (make-random-binomials #!optional (t 1) (p 0.5) (rand (make-uniform-random-reals)))
285        (unless (and (integer? t) (not (negative? t)))
286                (error 'make-random-binomials "t must be a non-negative integer" t))
287        (unless (and (real? p) (<= 0.0 p 1.0))
288                (error 'make-random-binomials "p must be a real in [0.0 1.0]" p))
289        (let ([bernoullis (make-random-bernoullis p rand)])
290                (values
291                        ;FIXME O(t) but O(log(t)) desired for >> t
292                        (if (< t LARGE-FIXNUM)
293                                (lambda ()
294                                        (do ([i 0 (fx+ 1 i)]
295                                                         [n 0 (if (bernoullis) (fx+ 1 n) n)])
296                                                        ([fx>= i t] n)))
297                                (lambda ()
298                                        (do ([i 0 (+ 1 i)]
299                                                         [n 0 (if (bernoullis) (+ 1 n) n)])
300                                                        ([>= i t] n))))
301                        (lambda ()
302                                (values t p rand)))) )
303
304;;; Geometrics
305
306(define (make-random-geometrics #!optional (p 0.5) (rand (make-uniform-random-reals)))
307        (unless (and (real? p) (<= 0.0 p 1.0))
308                (error 'make-random-geometrics "p must be a real in [0.0 1.0]" p))
309        (let ([log-p (log p)])
310                (values
311                        (lambda ()
312                                (+ 1 (inexact->exact (floor (/ (log (- 1.0 (rand))) log-p)))))
313                        (lambda ()
314                                (values p rand)))) )
315
316;;; Lognormals
317
318(define (make-random-lognormals #!optional (mu 1.0) (sigma 1.0) (rand (make-uniform-random-reals)))
319        (unless (and (real? mu) (not (zero? mu)))
320                (error 'make-random-lognormals "mu must be a non-zero real" mu))
321        (unless (and (real? sigma) (not (negative? sigma)))
322                (error 'make-random-lognormals "sigma must be a non-negative real" sigma))
323        (let ([normals (make-random-normals rand)]
324                                [nmu (log (* mu (/ mu (sqrt (+ (* sigma sigma) (* mu mu))))))]
325                                [nsigma (sqrt (log (+ 1 (* sigma (/ sigma mu mu)))))])
326                (values
327                        (lambda ()
328                                (exp (+ nmu (* (normals) nsigma))))
329                        (lambda ()
330                                (values mu sigma rand)))) )
331
332;;; Cauchys
333
334(define (make-random-cauchys #!optional (median 0.0) (sigma 1.0) (rand (make-uniform-random-reals)))
335        (unless (real? median)
336                (error 'make-random-cauchys "median must be a real" median))
337        (unless (and (real? sigma) (positive? sigma))
338                (error 'make-random-cauchys "sigma must be a positive real" sigma))
339        (values
340                (lambda ()
341                        (+ median (* sigma (tan (* PI (- (rand) 0.5))))))
342                (lambda ()
343                        (values median sigma rand))) )
344
345;;; Gammas
346
347;; "A Simple Method for Generating Gamma Variables", George Marsaglia & Wai Wan Tsang,
348;; ACM Transactions on Mathematical Software, Vol. 26, No. 3, September 2000, Pages 363Ð372.
349
350(define (make-random-gammas #!optional (alpha 1.0) (theta 1.0) (rand (make-uniform-random-reals)))
351        (unless (and (real? alpha) (positive? alpha))
352                (error 'make-random-gammas "alpha must be a positive real" alpha))
353        (unless (and (real? theta) (positive? theta))
354                (error 'make-random-gammas "theta must be a positive real" theta))
355        (values
356                (if (= 1.0 alpha)
357                        ;then special case
358                        (lambda ()
359                                (* theta (- (log (rand)))))
360                        ;else general case
361                        (let ([norms (make-random-normals 0.0 1.0 rand)]
362                                                [unis
363                                                        (if (< alpha 1.0)
364                                                                (let ([inv-alpha (inv alpha)])
365                                                                        (lambda ()
366                                                                                (expt (rand) inv-alpha)))
367                                                                rand)])
368                                (let* ([d (- (or (and (< alpha 1.0) (+ 1.0 alpha)) alpha) one-third)]
369                                                         [c (inv (sqrt (* 9.0 d)))])
370                                        (lambda ()
371                                                (* theta
372                                                        (let loop ([x (norms)])
373                                                                (let ([v (+ 1.0 (* c x))])
374                                                                        (if (and (< 0.0 v)
375                                                                                                                (let ([v (* v v v)]
376                                                                                                                                        [u (unis)]
377                                                                                                                                        [x^2 (* x x)])
378                                                                                                                        (or (< u (- 1.0 (* 0.0331 x^2 x^2)))
379                                                                                                                                        (< (log u) (+ (* 0.5 x^2) (* d (- 1.0 (+ v (log v)))))))))
380                                                                                (* d v)
381                                                                                (loop (norms))))))))))
382                (lambda ()
383                        (values alpha theta rand))) )
384
385;;; Erlangs
386
387(define (make-random-erlangs #!optional (alpha 1) (theta 1.0) (rand (make-uniform-random-reals)))
388        (unless (and (integer? alpha) (positive? alpha))
389                (error 'make-random-erlangs "alpha must be a positive integer" alpha))
390        (unless (and (real? theta) (positive? theta))
391                (error 'make-random-erlangs "theta must be a positive real" theta))
392        (make-random-gammas alpha theta rand) )
393
394;;; Pareto
395
396(define (make-random-paretos #!optional (alpha 1) (xmin 1.0) (rand (make-uniform-random-reals)))
397        (unless (and (real? alpha) (positive? alpha))
398                (error 'make-random-paretos "alpha must be a positive real" alpha))
399        (unless (and (real? xmin) (positive? xmin))
400                (error 'make-random-paretos "xmin must be a positive real" xmin))
401        (values
402                (let ([gs (make-random-gammas alpha (inv xmin) rand)])
403                        (make-random-exponentials 1.0 (lambda () (inv (+ xmin (gs))))))
404                (lambda () (values alpha xmin rand))) )
405
406;;; LŽvys
407
408;; See Stable Distributions - John P. Nolan, Formula 1.12
409
410(define (make-random-levys #!optional (gamma 1.0) (delta 0.0) (rand (make-uniform-random-reals)))
411        (unless (and (real? delta) (not (negative? delta)))
412                (error 'make-random-levys "delta must be a non-negative real" delta))
413        (unless (and (real? gamma) (positive? gamma))
414                (error 'make-random-levys "gamma must be a positive real" gamma))
415        (values
416                (if (and (= 1.0 gamma) (= 0.0 delta))
417                        (lambda ()
418                                (let ([r (rand)])
419                                        (inv (* r r))))
420                        (lambda ()
421                                (let ([r (rand)])
422                                        (+ delta (* gamma (inv (* r r)))))))
423                (lambda () (values gamma delta rand))) )
424
425;;; Normal vectors
426
427(define (random-normal-vector! vec #!optional (mu 0.0) (sigma 1.0) (rand (make-uniform-random-reals)))
428        (let ([norms (if (number? mu) (make-random-normals rand mu sigma) rand)])
429                (vector*-map! (lambda (i elm) (norms)) vec)) )
430
431(define (random-hollow-sphere! vec #!optional (mu 0.0) (sigma 1.0) (rand (make-uniform-random-reals)))
432        (random-normal-vector! vec rand mu sigma)
433        (vector*-scale! vec (inv (sqrt (vector*-sum-squares vec)))) )
434
435(define (random-solid-sphere! vec #!optional (mu 0.0) (sigma 1.0) (rand (make-uniform-random-reals)))
436        (random-hollow-sphere! vec rand mu sigma)
437        (vector*-scale! vec (expt (rand) (inv (vector*-length vec)))) )
Note: See TracBrowser for help on using the repository browser.