Changeset 2445 in project


Ignore:
Timestamp:
11/15/06 18:43:59 (15 years ago)
Author:
Kon Lovett
Message:

Bug fixes.

Location:
srfi-27
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • srfi-27/srfi-27-distributions.scm

    r2444 r2445  
    3838
    3939; (could include mathh-constants here)
     40
    4041(define-constant PI     3.1415926535897932384626433832795028841972)
    4142(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) )
    4252
    4353;;; Vector Support
     
    180190                (error 'make-random-exponentials "mu must be a real in [0.0 1.0]" mu))
    181191        (values
    182                 (lambda ()
    183                         (* mu (- (log (rand)))))
     192                (if (= 1.0 mu)
     193                        (lambda () (- (log (rand))))
     194                        (lambda () (* mu (- (log (rand))))))
    184195                (lambda ()
    185196                        (values mu rand))) )
     
    222233                (error 'make-random-triangles "s must be a real" s))
    223234        (unless (and (real? l) (< s l))
    224                 (error 'make-random-triangles "l must be a real in (a +inf)" l))
     235                (error 'make-random-triangles "l must be a real in (s +inf)" l))
    225236        (unless (and (real? m) (<= s m l))
    226237                (error 'make-random-triangles "m must be a real in [s l]" m))
     
    248259                        (lambda ()
    249260                                ; FIXME O(mu) but O(log(mu)) desired for >> mu
    250                                 (do ([m 0 (+ 1 m)]
     261                                (do ([m 0 (fx+ 1 m)]
    251262                                                 [prod (rand) (* prod (rand))])
    252263                                                ([<= prod emu] m)))
     
    260271                (error 'make-random-bernoullis "p must be a real in [0.0 1.0]" p))
    261272        (values
    262                 (lambda ()
    263                         (if (zero? p)
    264                                 #f
    265                                 (<= (rand) p)))
     273                (cond
     274                        [(zero? p) (lambda () #f)]
     275                        [(= 1.0 p) (lambda () #t)]
     276                        [else (lambda () (<= (rand) p))])
    266277                (lambda ()
    267278                        (values p rand))) )
    268279
    269280;;; Binomials
     281
     282(define-constant LARGE-FIXNUM 1073741823)
    270283
    271284(define (make-random-binomials #!optional (t 1) (p 0.5) (rand (make-uniform-random-reals)))
     
    274287        (unless (and (real? p) (<= 0.0 p 1.0))
    275288                (error 'make-random-binomials "p must be a real in [0.0 1.0]" p))
    276         (let ([bernoullis (make-random-bernoullis rand p)])
     289        (let ([bernoullis (make-random-bernoullis p rand)])
    277290                (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))))))
     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))))
    285301                        (lambda ()
    286302                                (values t p rand)))) )
     
    343359                                (* theta (- (log (rand)))))
    344360                        ;else general case
    345                         (let ([norms (make-random-normals rand)]
     361                        (let ([norms (make-random-normals 0.0 1.0 rand)]
    346362                                                [unis
    347363                                                        (if (< alpha 1.0)
    348                                                                 (let ([1/alpha (/ 1.0 alpha)])
     364                                                                (let ([inv-alpha (inv alpha)])
    349365                                                                        (lambda ()
    350                                                                                 (expt (rand) 1/alpha)))
     366                                                                                (expt (rand) inv-alpha)))
    351367                                                                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)))])
     368                                (let* ([d (- (or (and (< alpha 1.0) (+ 1.0 alpha)) alpha) one-third)]
     369                                                         [c (inv (sqrt (* 9.0 d)))])
    354370                                        (lambda ()
    355371                                                (* theta
     
    384400                (error 'make-random-paretos "xmin must be a positive real" xmin))
    385401        (values
    386                 (make-random-exponentials (/ 1.0 (+ xmin ((make-random-gammas alpha (/ 1.0 xmin))))))
     402                (let ([gs (make-random-gammas alpha (inv xmin) rand)])
     403                        (make-random-exponentials 1.0 (lambda () (inv (+ xmin (gs))))))
    387404                (lambda () (values alpha xmin rand))) )
    388405
    389406;;; LŽvys
     407
     408;; See Stable Distributions - John P. Nolan, Formula 1.12
    390409
    391410(define (make-random-levys #!optional (gamma 1.0) (delta 0.0) (rand (make-uniform-random-reals)))
     
    398417                        (lambda ()
    399418                                (let ([r (rand)])
    400                                         (/ 1.0 (* r r))))
     419                                        (inv (* r r))))
    401420                        (lambda ()
    402421                                (let ([r (rand)])
    403                                         (+ delta (* gamma (/ 1.0 (* r r)))))))
     422                                        (+ delta (* gamma (inv (* r r)))))))
    404423                (lambda () (values gamma delta rand))) )
    405424
     
    412431(define (random-hollow-sphere! vec #!optional (mu 0.0) (sigma 1.0) (rand (make-uniform-random-reals)))
    413432        (random-normal-vector! vec rand mu sigma)
    414         (vector*-scale! vec (/ 1.0 (sqrt (vector*-sum-squares vec)))) )
     433        (vector*-scale! vec (inv (sqrt (vector*-sum-squares vec)))) )
    415434
    416435(define (random-solid-sphere! vec #!optional (mu 0.0) (sigma 1.0) (rand (make-uniform-random-reals)))
    417436        (random-hollow-sphere! vec rand mu sigma)
    418         (vector*-scale! vec (expt (rand) (/ 1.0 (vector*-length vec)))) )
     437        (vector*-scale! vec (expt (rand) (inv (vector*-length vec)))) )
  • srfi-27/srfi-27-eggdoc.scm

    r2444 r2445  
    3838                (author (url "mailto:klovett@pacbell.net" "Kon Lovett"))
    3939                (history
    40                         (version "1.3" "Added Levy distribution")
     40                        (version "1.3" "Added Levy distribution [per Benedikt Rosenau], dist. bug fixes & special cases")
    4141                        (version "1.2" "Added Pareto & Erlang distributions")
    4242                        (version "1.1" "Changed entropy source port close timeout handler")
Note: See TracChangeset for help on using the changeset viewer.