Changeset 34968 in project


Ignore:
Timestamp:
12/29/17 20:58:25 (9 months ago)
Author:
kon
Message:

atomized distributions module

Location:
release/4/srfi-27/trunk
Files:
14 added
3 edited

Legend:

Unmodified
Added
Removed
  • release/4/srfi-27/trunk/srfi-27-distributions.scm

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

    r34867 r34968  
    1515        (numbers "2.8")
    1616        (mathh "3.2.0")
    17   (dsssl-utils "2.1.0")
     17  #;(dsssl-utils "2.1.0")
    1818        #;(random-bsd "0.2"))
    1919 (test-depends test)
    2020 (files
     21 "srfi-27.setup" "srfi-27.meta" #;"srfi-27.release-info"
     22 "srfi-27.scm"
     23 "srfi-27-implementation"
     24 "random-source.scm"
    2125 ;"bsdrnd.scm"
    22  "composite-random-source.scm" "composite-entropy-source.scm"
    23  "mrg32k3a.scm" "srfi-27.meta" "srfi-27-uniform-random.scm" "mwc.scm"
    24  "srfi-27.setup" "entropy-procedure.scm" "srfi-27.release-info"
    25  "srfi-27-implementation" "srfi-27-numbers.scm"
    26  "entropy-unix.scm" "entropy-linux.scm" "entropy-source.scm"
    27  "entropy-clock.scm" "random-source.scm" "srfi-27-vector-support.scm"
    28  "srfi-27-distributions.scm" "srfi-27-vector.scm" "srfi-27.scm"
    29  "entropy-windows.scm" "moa.scm" "entropy-support.scm"
     26 "composite-random-source.scm"
     27 "mrg32k3a.scm" "mwc.scm" "moa.scm"
     28 "srfi-27-numbers.scm"
     29 "entropy-source.scm"
     30 "composite-entropy-source.scm"
     31 "entropy-unix.scm" "entropy-linux.scm" "entropy-windows.scm"
     32 "entropy-clock.scm" "entropy-procedure.scm" "entropy-port.scm"
     33 "entropy-support.scm"
     34 "srfi-27-uniform-random.scm"
     35 "srfi-27-distributions.scm"
     36 "srfi-27-bernoullis.scm"
     37 "srfi-27-binomials.scm"
     38 "srfi-27-cauchys.scm"
     39 "srfi-27-erlangs.scm"
     40 "srfi-27-exponentials.scm"
     41 "srfi-27-gammas.scm"
     42 "srfi-27-geometrics.scm"
     43 "srfi-27-levys.scm"
     44 "srfi-27-lognormals.scm"
     45 "srfi-27-normals.scm"
     46 "srfi-27-paretos.scm"
     47 "srfi-27-poissons.scm"
     48 "srfi-27-triangles.scm"
     49 "srfi-27-weibulls.scm"
     50 "srfi-27-vector.scm" "srfi-27-vector-support.scm"
    3051 "source-registration.scm" "tests/test-diehard.scm" "tests/test-confidence.scm"
    31  "tests/test-mrg32k3a.scm" "tests/run.scm" "entropy-port.scm") )
     52 "tests/test-mrg32k3a.scm"
     53 "tests/srfi-27-test.scm" "tests/run.scm") )
  • release/4/srfi-27/trunk/srfi-27.setup

    r34967 r34968  
    124124  #:compile-options `(-scrutinize ,@publoptn) )
    125125
     126(setup-shared-extension-module 'srfi-27-bernoullis (extension-version "3.3.0")
     127  #:inline? #t
     128  #:types? #t
     129  #:compile-options `(-scrutinize ,@publoptn) )
     130
     131(setup-shared-extension-module 'srfi-27-binomials (extension-version "3.3.0")
     132  #:inline? #t
     133  #:types? #t
     134  #:compile-options `(-scrutinize ,@publoptn) )
     135
     136(setup-shared-extension-module 'srfi-27-cauchys (extension-version "3.3.0")
     137  #:inline? #t
     138  #:types? #t
     139  #:compile-options `(-scrutinize ,@publoptn) )
     140
     141(setup-shared-extension-module 'srfi-27-normals (extension-version "3.3.0")
     142  #:inline? #t
     143  #:types? #t
     144  #:compile-options `(-scrutinize ,@publoptn) )
     145
     146;needs normals
     147(setup-shared-extension-module 'srfi-27-gammas (extension-version "3.3.0")
     148  #:inline? #t
     149  #:types? #t
     150  #:compile-options `(-scrutinize ,@publoptn) )
     151
     152;needs gammas
     153(setup-shared-extension-module 'srfi-27-erlangs (extension-version "3.3.0")
     154  #:inline? #t
     155  #:types? #t
     156  #:compile-options `(-scrutinize ,@publoptn) )
     157
     158(setup-shared-extension-module 'srfi-27-exponentials (extension-version "3.3.0")
     159  #:inline? #t
     160  #:types? #t
     161  #:compile-options `(-scrutinize ,@publoptn) )
     162
     163(setup-shared-extension-module 'srfi-27-geometrics (extension-version "3.3.0")
     164  #:inline? #t
     165  #:types? #t
     166  #:compile-options `(-scrutinize ,@publoptn) )
     167
     168(setup-shared-extension-module 'srfi-27-levys (extension-version "3.3.0")
     169  #:inline? #t
     170  #:types? #t
     171  #:compile-options `(-scrutinize ,@publoptn) )
     172
     173(setup-shared-extension-module 'srfi-27-lognormals (extension-version "3.3.0")
     174  #:inline? #t
     175  #:types? #t
     176  #:compile-options `(-scrutinize ,@publoptn) )
     177
     178;needs gammas exponentials
     179(setup-shared-extension-module 'srfi-27-paretos (extension-version "3.3.0")
     180  #:inline? #t
     181  #:types? #t
     182  #:compile-options `(-scrutinize ,@publoptn) )
     183
     184(setup-shared-extension-module 'srfi-27-poissons (extension-version "3.3.0")
     185  #:inline? #t
     186  #:types? #t
     187  #:compile-options `(-scrutinize ,@publoptn) )
     188
     189(setup-shared-extension-module 'srfi-27-triangles (extension-version "3.3.0")
     190  #:inline? #t
     191  #:types? #t
     192  #:compile-options `(-scrutinize ,@publoptn) )
     193
     194(setup-shared-extension-module 'srfi-27-weibulls (extension-version "3.3.0")
     195  #:inline? #t
     196  #:types? #t
     197  #:compile-options `(-scrutinize ,@publoptn) )
     198
    126199(setup-shared-extension-module 'srfi-27-distributions (extension-version "3.3.0")
    127200  #:inline? #t
Note: See TracChangeset for help on using the changeset viewer.