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

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

Added levy dist.

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