From 0936d8c20b60f9621dee38b95f2f0f2e23bc8e5d Mon Sep 17 00:00:00 2001 From: Neil Toronto Date: Wed, 21 Nov 2012 21:16:35 -0700 Subject: [PATCH] Reworked distribution API, finally happy with it (as happy as I can be without being able to partially instantiate polymorphic parent struct types) Added docs for math/distributions (about 75% finished) Started docs for math/array (very incomplete) --- .../private/distributions/bernoulli-dist.rkt | 31 +- .../math/private/distributions/beta-dist.rkt | 64 ++- .../private/distributions/binomial-dist.rkt | 23 +- .../private/distributions/cauchy-dist.rkt | 30 +- .../math/private/distributions/delta-dist.rkt | 14 +- .../private/distributions/discrete-dist.rkt | 37 +- .../private/distributions/dist-struct.rkt | 163 +++--- .../distributions/exponential-dist.rkt | 25 +- .../math/private/distributions/gamma-dist.rkt | 26 +- .../private/distributions/geometric-dist.rkt | 46 +- .../distributions/impl/binomial-random.rkt | 110 ++-- .../distributions/impl/gamma-random.rkt | 115 ++-- .../distributions/impl/normal-random.rkt | 20 +- .../distributions/impl/poisson-random.rkt | 65 ++- .../private/distributions/logistic-dist.rkt | 28 +- .../private/distributions/normal-dist.rkt | 23 +- .../private/distributions/poisson-dist.rkt | 20 +- .../private/distributions/triangle-dist.rkt | 49 +- .../private/distributions/truncated-dist.rkt | 47 +- .../private/distributions/uniform-dist.rkt | 37 +- collects/math/private/distributions/utils.rkt | 49 +- collects/math/private/flonum/flonum-log.rkt | 3 +- collects/math/private/utils.rkt | 60 +- collects/math/scribblings/math-array.scrbl | 92 ++++ .../math/scribblings/math-distributions.scrbl | 517 ++++++++++++++++++ .../math/scribblings/math-number-theory.scrbl | 2 +- collects/math/scribblings/math.scrbl | 2 + collects/math/scribblings/utils.rkt | 7 +- 28 files changed, 1284 insertions(+), 421 deletions(-) create mode 100644 collects/math/scribblings/math-array.scrbl create mode 100644 collects/math/scribblings/math-distributions.scrbl diff --git a/collects/math/private/distributions/bernoulli-dist.rkt b/collects/math/private/distributions/bernoulli-dist.rkt index e523ca9728..d50c094ec5 100644 --- a/collects/math/private/distributions/bernoulli-dist.rkt +++ b/collects/math/private/distributions/bernoulli-dist.rkt @@ -3,14 +3,16 @@ (require racket/performance-hint racket/promise "../../flonum.rkt" + "../../vector.rkt" + "../unsafe.rkt" "dist-struct.rkt" "utils.rkt") (provide flbernoulli-pdf flbernoulli-cdf flbernoulli-inv-cdf - flbernoulli-random - Bernoulli-Dist bernoulli-dist bernoulli-dist? bernoulli-dist-prob) + flbernoulli-sample + Bernoulli-Dist bernoulli-dist bernoulli-dist-prob) (: flbernoulli-pdf (Flonum Flonum Any -> Flonum)) (define (flbernoulli-pdf q k log?) @@ -46,17 +48,17 @@ [else (cond [log? (if (p . fl<= . (fllog1p (- q))) 0.0 1.0)] [else (if (p . fl<= . (fl- 1.0 q)) 0.0 1.0)])])) -(: flbernoulli-random (Flonum -> Flonum)) -(define (flbernoulli-random q) - (cond [(not (flprobability? q)) +nan.0] - [else (if ((random) . > . q) 0.0 1.0)])) +(: flbernoulli-sample (Flonum Integer -> FlVector)) +(define (flbernoulli-sample q n) + (cond [(n . < . 0) (raise-argument-error 'flbernoulli-sample "Natural" 1 q n)] + [(not (flprobability? q)) (build-flvector n (λ (_) +nan.0))] + [else (build-flvector n (λ (_) (if ((random) . > . q) 0.0 1.0)))])) +(define-real-dist: bernoulli-dist Bernoulli-Dist + bernoulli-dist-struct ([prob : Flonum])) (begin-encourage-inline - (define-distribution-type: Bernoulli-Dist (Ordered-Dist Real Flonum) - bernoulli-dist ([prob : Flonum])) - (: bernoulli-dist (case-> (-> Bernoulli-Dist) (Real -> Bernoulli-Dist))) (define (bernoulli-dist [q 0.5]) @@ -67,9 +69,12 @@ (flbernoulli-cdf q (fl k) log? 1-p?))) (define inv-cdf (opt-lambda: ([p : Real] [log? : Any #f] [1-p? : Any #f]) (flbernoulli-inv-cdf q (fl p) log? 1-p?))) - (define (random) (flbernoulli-random q)) - (make-bernoulli-dist pdf random cdf inv-cdf - 0.0 1.0 (delay (if (q . fl<= . 0.5) 0.0 1.0)) - q))) + (define sample (case-lambda: + [() (unsafe-flvector-ref (flbernoulli-sample q 1) 0)] + [([n : Integer]) (flvector->list (flbernoulli-sample q n))])) + (bernoulli-dist-struct + pdf sample cdf inv-cdf + 0.0 1.0 (delay (if (q . fl<= . 0.5) 0.0 1.0)) + q))) ) diff --git a/collects/math/private/distributions/beta-dist.rkt b/collects/math/private/distributions/beta-dist.rkt index 35178da2f3..50e4bec761 100644 --- a/collects/math/private/distributions/beta-dist.rkt +++ b/collects/math/private/distributions/beta-dist.rkt @@ -1,8 +1,11 @@ #lang typed/racket/base -(require racket/performance-hint +(require racket/fixnum + racket/performance-hint racket/promise "../../flonum.rkt" + "../../vector.rkt" + "../unsafe.rkt" "../functions/beta.rkt" "../functions/incomplete-beta.rkt" "impl/beta-pdf.rkt" @@ -14,8 +17,8 @@ (provide flbeta-pdf flbeta-cdf flbeta-inv-cdf - flbeta-random - Beta-Dist beta-dist beta-dist? beta-dist-alpha beta-dist-beta) + flbeta-sample + Beta-Dist beta-dist beta-dist-alpha beta-dist-beta) (: flbeta-pdf (Flonum Flonum Flonum Any -> Flonum)) (define (flbeta-pdf a b x log?) @@ -32,33 +35,38 @@ [log? (fllog-beta-inc a b x 1-p? #t)] [else (flbeta-inc a b x 1-p? #t)])) -(: flbeta-random (Flonum Flonum -> Flonum)) -(define (flbeta-random a b) - (define x (standard-flgamma-random a)) - (define y (standard-flgamma-random b)) - (fl/ x (fl+ x y))) +(: flbeta-sample (Flonum Flonum Integer -> FlVector)) +(define (flbeta-sample a b n) + (cond [(n . < . 0) (raise-argument-error 'flbeta-sample "Natural" 2 a b n)] + [else + (define xs (flgamma-sample a 1.0 n)) + (define ys (flgamma-sample b 1.0 n)) + (for ([i (in-range n)]) + (define x (unsafe-flvector-ref xs i)) + (define y (unsafe-flvector-ref ys i)) + (unsafe-flvector-set! xs i (fl/ x (fl+ x y)))) + xs])) + +(define-real-dist: beta-dist Beta-Dist + beta-dist-struct ([alpha : Flonum] [beta : Flonum])) (begin-encourage-inline - (define-distribution-type: Beta-Dist (Ordered-Dist Real Flonum) - beta-dist ([alpha : Flonum] [beta : Flonum])) - - (: beta-dist (case-> (-> Beta-Dist) - (Real Real -> Beta-Dist))) - (define beta-dist - (case-lambda - [() (beta-dist 1.0 1.0)] - [(a b) - (let ([a (fl a)] [b (fl b)]) - (define pdf (opt-lambda: ([x : Real] [log? : Any #f]) - (flbeta-pdf a b (fl x) log?))) - (define cdf (opt-lambda: ([x : Real] [log? : Any #f] [1-p? : Any #f]) - (flbeta-cdf a b (fl x) log? 1-p?))) - (define inv-cdf (opt-lambda: ([p : Real] [log? : Any #f] [1-p? : Any #f]) - (flbeta-inv-cdf a b (fl p) log? 1-p?))) - (define (random) (flbeta-random a b)) - (make-beta-dist pdf random cdf inv-cdf - 0.0 1.0 (delay (flbeta-inv-cdf a b 0.5 #f #f)) - a b))])) + (: beta-dist (Real Real -> Beta-Dist)) + (define (beta-dist a b) + (let ([a (fl a)] [b (fl b)]) + (define pdf (opt-lambda: ([x : Real] [log? : Any #f]) + (flbeta-pdf a b (fl x) log?))) + (define cdf (opt-lambda: ([x : Real] [log? : Any #f] [1-p? : Any #f]) + (flbeta-cdf a b (fl x) log? 1-p?))) + (define inv-cdf (opt-lambda: ([p : Real] [log? : Any #f] [1-p? : Any #f]) + (flbeta-inv-cdf a b (fl p) log? 1-p?))) + (define sample (case-lambda: + [() (unsafe-flvector-ref (flbeta-sample a b 1) 0)] + [([n : Integer]) (flvector->list (flbeta-sample a b n))])) + (beta-dist-struct + pdf sample cdf inv-cdf + 0.0 1.0 (delay (flbeta-inv-cdf a b 0.5 #f #f)) + a b))) ) diff --git a/collects/math/private/distributions/binomial-dist.rkt b/collects/math/private/distributions/binomial-dist.rkt index 4f4ccdccdb..01997ded23 100644 --- a/collects/math/private/distributions/binomial-dist.rkt +++ b/collects/math/private/distributions/binomial-dist.rkt @@ -3,6 +3,8 @@ (require racket/performance-hint racket/promise "../../flonum.rkt" + "../../vector.rkt" + "../unsafe.rkt" "../functions/incomplete-beta.rkt" (prefix-in impl: "impl/binomial-pdf.rkt") "impl/binomial-random.rkt" @@ -13,8 +15,8 @@ (provide flbinomial-pdf flbinomial-cdf flbinomial-inv-cdf - flbinomial-random - Binomial-Dist binomial-dist binomial-dist? binomial-dist-count binomial-dist-prob) + flbinomial-sample + Binomial-Dist binomial-dist binomial-dist-count binomial-dist-prob) (: flbinomial-pdf (Flonum Flonum Flonum Any -> Flonum)) (define (flbinomial-pdf n q k log?) @@ -51,11 +53,11 @@ 0.0 n (flmax 0.0 (flmin n z)))])) +(define-real-dist: binomial-dist Binomial-Dist + binomial-dist-struct ([count : Flonum] [prob : Flonum])) + (begin-encourage-inline - (define-distribution-type: Binomial-Dist (Ordered-Dist Real Flonum) - binomial-dist ([count : Flonum] [prob : Flonum])) - (: binomial-dist (case-> (-> Binomial-Dist) (Real -> Binomial-Dist) (Real Real -> Binomial-Dist))) @@ -67,9 +69,12 @@ (flbinomial-cdf n q (fl k) log? 1-p?))) (define inv-cdf (opt-lambda: ([p : Real] [log? : Any #f] [1-p? : Any #f]) (flbinomial-inv-cdf n q (fl p) log? 1-p?))) - (define (random) (flbinomial-random n q)) - (make-binomial-dist pdf random cdf inv-cdf - 0.0 +inf.0 (delay (flfloor (fl* n q))) - n q))) + (define sample (case-lambda: + [() (unsafe-flvector-ref (flbinomial-sample n q 1) 0)] + [([m : Integer]) (flvector->list (flbinomial-sample n q m))])) + (binomial-dist-struct + pdf sample cdf inv-cdf + 0.0 +inf.0 (delay (flfloor (fl* n q))) + n q))) ) diff --git a/collects/math/private/distributions/cauchy-dist.rkt b/collects/math/private/distributions/cauchy-dist.rkt index aa09a052b4..82294e32a9 100644 --- a/collects/math/private/distributions/cauchy-dist.rkt +++ b/collects/math/private/distributions/cauchy-dist.rkt @@ -4,14 +4,16 @@ racket/promise "../../flonum.rkt" "../../base.rkt" + "../../vector.rkt" + "../unsafe.rkt" "dist-struct.rkt" "utils.rkt") (provide flcauchy-pdf flcauchy-cdf flcauchy-inv-cdf - flcauchy-random - Cauchy-Dist cauchy-dist cauchy-dist? cauchy-dist-center cauchy-dist-scale) + flcauchy-sample + Cauchy-Dist cauchy-dist cauchy-dist-mode cauchy-dist-scale) (: flcauchy-pdf (Float Float Float Any -> Float)) (define flcauchy-pdf @@ -54,29 +56,37 @@ (: flcauchy-inv-cdf (Float Float Float Any Any -> Float)) (define flcauchy-inv-cdf (make-symmetric-location-scale-flinv-cdf standard-flcauchy-inv-cdf)) -(: flcauchy-random (Float Float -> Float)) -(define flcauchy-random (make-symmetric-location-scale-flrandom standard-flcauchy-inv-cdf)) +(: flcauchy-sample-single (Flonum Flonum -> Flonum)) +(define flcauchy-sample-single + (make-symmetric-location-scale-flrandom standard-flcauchy-inv-cdf)) + +(: flcauchy-sample (Float Float Integer -> FlVector)) +(define (flcauchy-sample c s n) + (cond [(n . < . 0) (raise-argument-error 'flcauchy-sample "Natural" 2 c s n)] + [else (build-flvector n (λ (_) (flcauchy-sample-single c s)))])) ;; =================================================================================================== ;; Distribution object +(define-real-dist: cauchy-dist Cauchy-Dist + cauchy-dist-struct ([mode : Float] [scale : Float])) + (begin-encourage-inline - (define-distribution-type: Cauchy-Dist (Ordered-Dist Real Flonum) - cauchy-dist ([center : Float] [scale : Float])) - (: cauchy-dist (case-> (-> Cauchy-Dist) (Real -> Cauchy-Dist) (Real Real -> Cauchy-Dist))) (define (cauchy-dist [c 0.0] [s 1.0]) - (let ([c (fl c)] [s (fl s)]) + (let ([c (fl c)] [s (fl s)]) (define pdf (opt-lambda: ([x : Real] [log? : Any #f]) (flcauchy-pdf c s (fl x) log?))) (define cdf (opt-lambda: ([x : Real] [log? : Any #f] [complement? : Any #f]) (flcauchy-cdf c s (fl x) log? complement?))) (define inv-cdf (opt-lambda: ([p : Real] [log? : Any #f] [complement? : Any #f]) (flcauchy-inv-cdf c s (fl p) log? complement?))) - (define (random) (flcauchy-random c s)) - (make-cauchy-dist pdf random cdf inv-cdf -inf.0 +inf.0 (delay c) c s))) + (define sample (case-lambda: + [() (unsafe-flvector-ref (flcauchy-sample c s 1) 0)] + [([n : Integer]) (flvector->list (flcauchy-sample c s n))])) + (cauchy-dist-struct pdf sample cdf inv-cdf -inf.0 +inf.0 (delay c) c s))) ) diff --git a/collects/math/private/distributions/delta-dist.rkt b/collects/math/private/distributions/delta-dist.rkt index c2a2ff2fe1..f675835224 100644 --- a/collects/math/private/distributions/delta-dist.rkt +++ b/collects/math/private/distributions/delta-dist.rkt @@ -10,13 +10,13 @@ (provide fldelta-pdf fldelta-cdf fldelta-inv-cdf - Delta-Dist delta-dist delta-dist?) + Delta-Dist delta-dist delta-dist-mean) + +(define-real-dist: delta-dist Delta-Dist + delta-dist-struct ([mean : Flonum])) (begin-encourage-inline - (define-distribution-type: Delta-Dist (Ordered-Dist Real Flonum) - delta-dist ([center : Float])) - (: delta-dist (case-> (-> Delta-Dist) (Real -> Delta-Dist))) (define (delta-dist [c 0.0]) @@ -27,7 +27,9 @@ (fldelta-cdf c (fl x) log? 1-p?))) (define inv-cdf (opt-lambda: ([p : Real] [log? : Any #f] [1-p? : Any #f]) (fldelta-inv-cdf c (fl p) log? 1-p?))) - (define (random) c) - (make-delta-dist pdf random cdf inv-cdf c c (delay c) c))) + (define sample (case-lambda: + [() c] + [([n : Integer]) (build-list n (λ (_) c))])) + (delta-dist-struct pdf sample cdf inv-cdf c c (delay c) c))) ) diff --git a/collects/math/private/distributions/discrete-dist.rkt b/collects/math/private/distributions/discrete-dist.rkt index 3f33b10d21..d4386258e5 100644 --- a/collects/math/private/distributions/discrete-dist.rkt +++ b/collects/math/private/distributions/discrete-dist.rkt @@ -1,19 +1,37 @@ #lang typed/racket/base -(require racket/promise +(require racket/performance-hint + racket/promise "../../flonum.rkt" "../statistics/statistics-utils.rkt" + "../utils.rkt" "impl/walker-table.rkt" - "dist-struct.rkt" - "utils.rkt") + "dist-struct.rkt") (provide Discrete-Dist discrete-dist discrete-dist-values discrete-dist-weights) -(define-distribution-type: (Discrete-Dist A) (Dist A A) (In Out) - discrete-dist ([values : (Listof Out)] [weights : (Listof Positive-Flonum)])) +(begin-encourage-inline + (struct: (In Out) discrete-dist-struct dist ([values : (Listof Out)] + [weights : (Listof Positive-Flonum)]) + #:property prop:custom-print-quotable 'never + #:property prop:custom-write + (λ (v port mode) + (pretty-print-constructor + 'discrete-dist (list (discrete-dist-struct-values v) (discrete-dist-struct-weights v)) + port mode))) + + (define-type (Discrete-Dist A) (discrete-dist-struct A A)) + + (: discrete-dist-values (All (A) ((Discrete-Dist A) -> (Listof A)))) + (define (discrete-dist-values d) (discrete-dist-struct-values d)) + + (: discrete-dist-weights (All (A) ((Discrete-Dist A) -> (Listof Positive-Flonum)))) + (define (discrete-dist-weights d) (discrete-dist-struct-weights d)) + + ) (: discrete-dist (All (A) (case-> ((Sequenceof A) -> (Discrete-Dist A)) @@ -30,5 +48,10 @@ (define pdf (opt-lambda: ([x : A] [log? : Any #f]) (define p (hash-ref (force hash) x (λ () 0.0))) (if log? (fllog p) p))) - (define random (λ () (walker-table-sample (force table)))) - (make-discrete-dist pdf random xs ws))) + (define sample (case-lambda: + [() (walker-table-sample (force table))] + [([n : Integer]) + (cond [(n . < . 0) (raise-argument-error 'discrete-dist-sample "Natural" n)] + [else (let ([table (force table)]) + (build-list n (λ (_) (walker-table-sample table))))])])) + (discrete-dist-struct pdf sample xs ws))) diff --git a/collects/math/private/distributions/dist-struct.rkt b/collects/math/private/distributions/dist-struct.rkt index c2ac60a676..71d6d3b9e1 100644 --- a/collects/math/private/distributions/dist-struct.rkt +++ b/collects/math/private/distributions/dist-struct.rkt @@ -1,37 +1,34 @@ #lang racket/base (provide - PDF CDF Inverse-CDF - Dist - Ordered-Dist - Real-Dist - Real-PDF - Real-CDF - Real-Inverse-CDF - dist? - ordered-dist? - dist-pdf - dist-random + (all-from-out (submod "." typed-defs)) + ;; Rename transformers dist-cdf dist-inv-cdf dist-min - dist-max - dist-median - sample - sample* - real-dist-prob) + dist-max) (module typed-defs typed/racket/base (require racket/performance-hint racket/promise "../../flonum.rkt") - (provide (all-defined-out)) + (provide + PDF Sample CDF Inverse-CDF + (struct-out dist) + (struct-out ordered-dist) + Real-Dist + dist-median sample + real-dist-prob) (define-type (PDF In) (case-> (In -> Flonum) (In Any -> Flonum))) + (define-type (Sample Out) + (case-> (-> Out) + (Integer -> (Listof Out)))) + (define-type (CDF In) (case-> (In -> Flonum) (In Any -> Flonum) @@ -42,11 +39,11 @@ (Real Any -> Out) (Real Any Any -> Out))) - (struct: (In Out) Dist ([pdf : (PDF In)] - [random : (-> Out)]) + (struct: (In Out) dist ([pdf : (PDF In)] + [sample : (Sample Out)]) #:transparent) - (struct: (In Out) Ordered-Dist Dist + (struct: (In Out) ordered-dist dist ([cdf : (CDF In)] [inv-cdf : (Inverse-CDF Out)] [min : Out] @@ -54,73 +51,99 @@ [median : (Promise Out)]) #:transparent) - (define-type Real-Dist (Ordered-Dist Real Flonum)) - (define-type Real-PDF (PDF Real)) - (define-type Real-CDF (CDF Real)) - (define-type Real-Inverse-CDF (Inverse-CDF Flonum)) + (define-type Real-Dist (ordered-dist Real Flonum)) ;; ================================================================================================= (begin-encourage-inline - (: dist-median (All (In Out) ((Ordered-Dist In Out) -> Out))) - (define (dist-median d) (force (Ordered-Dist-median d))) + (: dist-median (All (In Out) ((ordered-dist In Out) -> Out))) + (define (dist-median d) (force (ordered-dist-median d))) - (: sample (All (In Out) ((Dist In Out) -> Out))) - (define (sample d) ((Dist-random d))) - - (: sample* (All (In Out) ((Dist In Out) Integer -> (Listof Out)))) - (define (sample* d n) - (cond [(n . < . 0) (raise-argument-error 'sample* "Natural" 1 d n)] - [(index? n) (define random (Dist-random d)) - (for/list: : (Listof Out) ([_ (in-range n)]) (random))] - [else (raise-argument-error 'sample* "Index" 1 d n)])) + (: sample (All (In Out) (case-> ((dist In Out) -> Out) + ((dist In Out) Integer -> (Listof Out))))) + (define sample + (case-lambda + [(d) ((dist-sample d))] + [(d n) ((dist-sample d) n)])) ) + (: real-dist-prob* (Real-Dist Flonum Flonum Any -> Flonum)) + ;; Assumes a <= b + (define (real-dist-prob* d a b 1-p?) + (define c (dist-median d)) + (define cdf (ordered-dist-cdf d)) + (define p + (cond [(a . fl= . b) (if 1-p? 1.0 0.0)] + [1-p? (fl+ (cdf a #f #f) (cdf b #f #t))] + [(b . fl<= . c) + ;; Both less than the median; use lower tail only + (fl- (cdf b #f #f) (cdf a #f #f))] + [(a . fl>= . c) + ;; Both greater than the median; use upper tail only + (fl- (cdf a #f #t) (cdf b #f #t))] + [else + ;; Median between a and b; use lower for (a,c] and upper for (c,b] + (define P_x<=a (cdf a #f #f)) + (define P_x>b (cdf b #f #t)) + (fl+ (fl- 0.5 P_x<=a) (fl- 0.5 P_x>b))])) + (max 0.0 (min 1.0 p))) + + (: real-dist-log-prob* (Real-Dist Flonum Flonum Any -> Flonum)) + ;; Assumes a <= b + (define (real-dist-log-prob* d a b 1-p?) + (define c (dist-median d)) + (define cdf (ordered-dist-cdf d)) + (define log-p + (cond [(a . fl= . b) (if 1-p? 0.0 -inf.0)] + [1-p? (lg+ (cdf a #t #f) (cdf b #t #t))] + [(b . fl<= . c) + ;; Both less than the median; use lower tail only + (define log-P_x<=a (cdf a #t #f)) + (define log-P_x<=b (cdf b #t #f)) + (cond [(log-P_x<=b . fl< . log-P_x<=a) -inf.0] + [else (lg- log-P_x<=b log-P_x<=a)])] + [(a . fl>= . c) + ;; Both greater than the median; use upper tail only + (define log-P_x>a (cdf a #t #t)) + (define log-P_x>b (cdf b #t #t)) + (cond [(log-P_x>a . fl< . log-P_x>b) -inf.0] + [else (lg- log-P_x>a log-P_x>b)])] + [else + ;; Median between a and b; try 1-upper first + (define log-P_x<=a (cdf a #t #f)) + (define log-P_x>b (cdf b #t #t)) + (define log-p (lg1- (lg+ log-P_x<=a log-P_x>b))) + (cond [(log-p . fl> . (log 0.1)) log-p] + [else + ;; Subtracting from 1.0 (in log space) lost bits; split and add instead + (define log-P_ab) -inf.0] + [else (lg- (fllog 0.5) log-P_x>b)])) + (lg+ log-P_a (Real-Dist Real Real -> Flonum) (Real-Dist Real Real Any -> Flonum) (Real-Dist Real Real Any Any -> Flonum))) (define (real-dist-prob d a b [log? #f] [1-p? #f]) (let ([a (fl a)] [b (fl b)]) (let ([a (flmin a b)] [b (flmax a b)]) - (define c (force (Ordered-Dist-median d))) - (define cdf (Ordered-Dist-cdf d)) - (define log-p - (min (cond [1-p? (lg+ (cdf a #t #f) (cdf b #t #t))] - [(b . fl<= . c) - (define log-P_x<=a (cdf a #t #f)) - (define log-P_x<=b (cdf b #t #f)) - (cond [(log-P_x<=b . fl< . log-P_x<=a) -inf.0] - [else (lg- log-P_x<=b log-P_x<=a)])] - [(a . fl>= . c) - (define log-P_x>a (cdf a #t #t)) - (define log-P_x>b (cdf b #t #t)) - (cond [(log-P_x>a . fl< . log-P_x>b) -inf.0] - [else (lg- log-P_x>a log-P_x>b)])] - [else - (define log-P_x<=a (cdf a #t #f)) - (define log-P_x>b (cdf b #t #t)) - (define log-P_ab) -inf.0] - [else (lg- (fllog 0.5) log-P_x>b)])) - (lg+ log-P_a . +max-subnormal.0) (p . fl< . 0.9)) (fllog p)] + [else (real-dist-log-prob* d a b 1-p?)])] + [else (real-dist-prob* d a b 1-p?)])))) + ) (require (submod "." typed-defs) (for-syntax racket/base)) -(define-syntax dist? (make-rename-transformer #'Dist?)) -(define-syntax ordered-dist? (make-rename-transformer #'Ordered-Dist?)) -(define-syntax dist-pdf (make-rename-transformer #'Dist-pdf)) -(define-syntax dist-random (make-rename-transformer #'Dist-random)) -(define-syntax dist-cdf (make-rename-transformer #'Ordered-Dist-cdf)) -(define-syntax dist-inv-cdf (make-rename-transformer #'Ordered-Dist-inv-cdf)) -(define-syntax dist-min (make-rename-transformer #'Ordered-Dist-min)) -(define-syntax dist-max (make-rename-transformer #'Ordered-Dist-max)) +(define-syntax dist-cdf (make-rename-transformer #'ordered-dist-cdf)) +(define-syntax dist-inv-cdf (make-rename-transformer #'ordered-dist-inv-cdf)) +(define-syntax dist-min (make-rename-transformer #'ordered-dist-min)) +(define-syntax dist-max (make-rename-transformer #'ordered-dist-max)) diff --git a/collects/math/private/distributions/exponential-dist.rkt b/collects/math/private/distributions/exponential-dist.rkt index be768d5929..272f48e9d1 100644 --- a/collects/math/private/distributions/exponential-dist.rkt +++ b/collects/math/private/distributions/exponential-dist.rkt @@ -3,14 +3,16 @@ (require racket/performance-hint racket/promise "../../flonum.rkt" + "../../vector.rkt" + "../unsafe.rkt" "dist-struct.rkt" "utils.rkt") (provide flexponential-pdf flexponential-cdf flexponential-inv-cdf - flexponential-random - Exponential-Dist exponential-dist exponential-dist? exponential-dist-scale) + flexponential-sample + Exponential-Dist exponential-dist exponential-dist-mean) (: flexponential-pdf (Float Float Any -> Float)) (define flexponential-pdf @@ -33,18 +35,19 @@ (: flexponential-inv-cdf (Float Float Any Any -> Float)) (define flexponential-inv-cdf (make-one-sided-scale-flinv-cdf standard-flexponential-inv-cdf)) -(: flexponential-random (Float -> Float)) -(define (flexponential-random s) - (fl* s (- (fllog (random))))) +(: flexponential-sample (Flonum Integer -> FlVector)) +(define (flexponential-sample s n) + (cond [(n . < . 0) (raise-argument-error 'flexponential-sample "Natural" 1 s n)] + [else (build-flvector n (λ (_) (fl* s (- (fllog (random))))))])) ;; =================================================================================================== ;; Distribution object +(define-real-dist: exponential-dist Exponential-Dist + exponential-dist-struct ([mean : Flonum])) + (begin-encourage-inline - (define-distribution-type: Exponential-Dist (Ordered-Dist Real Flonum) - exponential-dist ([scale : Flonum])) - (: exponential-dist (case-> (-> Exponential-Dist) (Real -> Exponential-Dist))) (define (exponential-dist [s 1.0]) @@ -55,7 +58,9 @@ (flexponential-cdf s (fl x) log? 1-p?))) (define inv-cdf (opt-lambda: ([p : Real] [log? : Any #f] [1-p? : Any #f]) (flexponential-inv-cdf s (fl p) log? 1-p?))) - (define (random) (flexponential-random s)) - (make-exponential-dist pdf random cdf inv-cdf 0.0 +inf.0 (delay (fl* s (fllog 2.0))) s))) + (define sample (case-lambda: + [() (unsafe-flvector-ref (flexponential-sample s 1) 0)] + [([n : Integer]) (flvector->list (flexponential-sample s n))])) + (exponential-dist-struct pdf sample cdf inv-cdf 0.0 +inf.0 (delay (fl* s (fllog 2.0))) s))) ) diff --git a/collects/math/private/distributions/gamma-dist.rkt b/collects/math/private/distributions/gamma-dist.rkt index 228db323da..22a72d9c03 100644 --- a/collects/math/private/distributions/gamma-dist.rkt +++ b/collects/math/private/distributions/gamma-dist.rkt @@ -3,6 +3,8 @@ (require racket/performance-hint racket/promise "../../flonum.rkt" + "../../vector.rkt" + "../unsafe.rkt" "../functions/incomplete-gamma.rkt" "impl/gamma-pdf.rkt" "impl/gamma-inv-cdf.rkt" @@ -13,8 +15,8 @@ (provide flgamma-pdf flgamma-cdf flgamma-inv-cdf - flgamma-random - Gamma-Dist gamma-dist gamma-dist? gamma-dist-shape gamma-dist-scale) + flgamma-sample + Gamma-Dist gamma-dist gamma-dist-shape gamma-dist-scale) (: flgamma-pdf (Float Float Float Any -> Float)) (define (flgamma-pdf k s x log?) @@ -45,16 +47,11 @@ (standard-flgamma-inv-cdf k p log? 1-p?))) s p log? 1-p?)])) -(: flgamma-random (Float Float -> Float)) -(define (flgamma-random k s) - (cond [(k . fl< . 0.0) +nan.0] - [else (fl* s (standard-flgamma-random k))])) +(define-real-dist: gamma-dist Gamma-Dist + gamma-dist-struct ([shape : Flonum] [scale : Flonum])) (begin-encourage-inline - (define-distribution-type: Gamma-Dist (Ordered-Dist Real Flonum) - gamma-dist ([shape : Float] [scale : Float])) - (: gamma-dist (case-> (-> Gamma-Dist) (Real -> Gamma-Dist) (Real Real -> Gamma-Dist))) @@ -66,9 +63,12 @@ (flgamma-cdf k s (fl x) log? 1-p?))) (define inv-cdf (opt-lambda: ([p : Real] [log? : Any #f] [1-p? : Any #f]) (flgamma-inv-cdf k s (fl p) log? 1-p?))) - (define (random) (flgamma-random k s)) - (make-gamma-dist pdf random cdf inv-cdf - 0.0 +inf.0 (delay (flgamma-inv-cdf k s 0.5 #f #f)) - k s))) + (define sample (case-lambda: + [() (unsafe-flvector-ref (flgamma-sample k s 1) 0)] + [([n : Integer]) (flvector->list (flgamma-sample k s n))])) + (gamma-dist-struct + pdf sample cdf inv-cdf + 0.0 +inf.0 (delay (flgamma-inv-cdf k s 0.5 #f #f)) + k s))) ) diff --git a/collects/math/private/distributions/geometric-dist.rkt b/collects/math/private/distributions/geometric-dist.rkt index 73c1b0b67c..7ae905c414 100644 --- a/collects/math/private/distributions/geometric-dist.rkt +++ b/collects/math/private/distributions/geometric-dist.rkt @@ -3,14 +3,16 @@ (require racket/performance-hint racket/promise "../../flonum.rkt" + "../../vector.rkt" + "../unsafe.rkt" "dist-struct.rkt" "utils.rkt") (provide flgeometric-pdf flgeometric-cdf flgeometric-inv-cdf - flgeometric-random - Geometric-Dist geometric-dist geometric-dist? geometric-dist-prob) + flgeometric-sample + Geometric-Dist geometric-dist geometric-dist-prob) (: flgeometric-pdf (Flonum Flonum Any -> Flonum)) (define (flgeometric-pdf q k log?) @@ -49,22 +51,27 @@ [else (if log? (lg1- p) (fllog1p (- p)))])) (flmax 0.0 (fl- (flceiling (fl/ log-1-p (fllog1p (- q)))) 1.0))])) -(: flgeometric-random (Flonum -> Flonum)) -(define (flgeometric-random q) - (cond [(or (q . fl<= . 0.0) (q . fl>= . 1.0)) - (cond [(fl= q 1.0) 0.0] - [(fl= q 0.0) +inf.0] - [else +nan.0])] +(: flgeometric-sample (Flonum Integer -> FlVector)) +(define (flgeometric-sample q n) + (cond [(n . < . 0) (raise-argument-error 'flgeometric-sample "Natural" 1 q n)] + [(or (q . fl<= . 0.0) (q . fl>= . 1.0)) + (define v + (cond [(fl= q 1.0) 0.0] + [(fl= q 0.0) +inf.0] + [else +nan.0])) + (build-flvector n (λ (_) v))] [else - (define p (fl* 0.5 (random))) - (define log-1-p (if ((random) . fl> . 0.5) (fllog p) (fllog1p (- p)))) - (flmax 0.0 (fl- (flceiling (fl/ log-1-p (fllog1p (- q)))) 1.0))])) + (build-flvector + n (λ (_) + (define p (fl* 0.5 (random))) + (define log-1-p (if ((random) . fl> . 0.5) (fllog p) (fllog1p (- p)))) + (flmax 0.0 (fl- (flceiling (fl/ log-1-p (fllog1p (- q)))) 1.0))))])) + +(define-real-dist: geometric-dist Geometric-Dist + geometric-dist-struct ([prob : Flonum])) (begin-encourage-inline - (define-distribution-type: Geometric-Dist (Ordered-Dist Real Flonum) - geometric-dist ([prob : Flonum])) - (: geometric-dist (case-> (-> Geometric-Dist) (Real -> Geometric-Dist))) (define (geometric-dist [q 0.5]) @@ -75,9 +82,12 @@ (flgeometric-cdf q (fl k) log? 1-p?))) (define inv-cdf (opt-lambda: ([p : Real] [log? : Any #f] [1-p? : Any #f]) (flgeometric-inv-cdf q (fl p) log? 1-p?))) - (define (random) (flgeometric-random q)) - (make-geometric-dist pdf random cdf inv-cdf - 0.0 +inf.0 (delay (flgeometric-inv-cdf q 0.5 #f #f)) - q))) + (define sample (case-lambda: + [() (unsafe-flvector-ref (flgeometric-sample q 1) 0)] + [([n : Integer]) (flvector->list (flgeometric-sample q n))])) + (geometric-dist-struct + pdf sample cdf inv-cdf + 0.0 +inf.0 (delay (flgeometric-inv-cdf q 0.5 #f #f)) + q))) ) diff --git a/collects/math/private/distributions/impl/binomial-random.rkt b/collects/math/private/distributions/impl/binomial-random.rkt index cae63d8a42..c19dd58c56 100644 --- a/collects/math/private/distributions/impl/binomial-random.rkt +++ b/collects/math/private/distributions/impl/binomial-random.rkt @@ -4,33 +4,36 @@ Wolfgang Hormann. The Generation of Binomial Random Variates. |# -(require math/base - math/flonum - "../normal-dist.rkt" - "../dist-struct.rkt") +(require "../../../base.rkt" + "../../../flonum.rkt" + "../../../vector.rkt" + "../../unsafe.rkt" + "normal-random.rkt") -(provide flbinomial-random) +(provide flbinomial-sample) -(: flbinomial-random-small (Flonum Flonum -> Flonum)) +(: flbinomial-sample-small (Flonum Flonum Natural -> FlVector)) ;; For n*min(p,1-p) <= 30 -(define (flbinomial-random-small n p) +(define (flbinomial-sample-small n p m) (let-values ([(p q s?) (cond [(p . fl< . 0.5) (values p (fl- 1.0 p) #f)] [else (values (fl- 1.0 p) p #t)])]) (define q^n (flexpt q n)) (define r (fl/ p q)) (define g (fl* r (fl+ n 1.0))) - (define k - (let: reject : Flonum () - (let loop ([k 0.0] [f q^n] [u (random)]) - (cond [(u . fl< . f) k] - [(k . fl> . 110.0) (reject)] - [else (let ([k (fl+ k 1.0)]) - (loop k (fl* f (fl- (fl/ g k) r)) (fl- u f)))])))) - (if s? (fl- n k) k))) + (build-flvector + m (λ (_) + (define k + (let: reject : Flonum () + (let loop ([k 0.0] [f q^n] [u (random)]) + (cond [(u . fl< . f) k] + [(k . fl> . 110.0) (reject)] + [else (let ([k (fl+ k 1.0)]) + (loop k (fl* f (fl- (fl/ g k) r)) (fl- u f)))])))) + (if s? (fl- n k) k))))) -(: flbinomial-random-hormann (Flonum Flonum -> Flonum)) +(: flbinomial-sample-hormann (Flonum Flonum Natural -> FlVector)) ;; For n*min(p,1-p) >= 10 -(define (flbinomial-random-hormann n p) +(define (flbinomial-sample-hormann n p j) (let-values ([(p q s?) (cond [(p . fl< . 0.5) (values p (fl- 1.0 p) #f)] [else (values (fl- 1.0 p) p #t)])]) (define σ (flsqrt (* n p q))) @@ -41,35 +44,39 @@ Wolfgang Hormann. The Generation of Binomial Random Variates. (define c (fl+ 0.5 (fl* n p))) (define α (fl* σ (fl+ 2.83 (fl/ 5.1 b)))) (define vr (fl- 0.92 (fl/ 4.2 b))) - (define k - (let: loop : Flonum () - (define v (random)) - (define u (fl- (random) 0.5)) - (define us (fl- 0.5 (flabs u))) - (define k (flfloor (fl+ c (fl* u (fl+ b (fl/ (fl* 2.0 a) us)))))) - (cond [(or (k . fl< . 0.0) (k . fl> . n)) (loop)] - [(and (us . fl>= . 0.07) (v . fl<= . vr)) k] - [else - (let ([v (fl* v (fl/ α (fl+ b (fl/ a (fl* us us)))))]) - (define h (+ (fllog-factorial m) - (fllog-factorial (fl- n m)) - (- (fllog-factorial k)) - (- (fllog-factorial (fl- n k))) - (fl* (fl- k m) (fllog (fl/ p q))))) - (cond [((fllog v) . fl<= . h) k] - [else (loop)]))]))) - (if s? (fl- n k) k))) + (build-flvector + j (λ (_) + (define k + (let: loop : Flonum () + (define v (random)) + (define u (fl- (random) 0.5)) + (define us (fl- 0.5 (flabs u))) + (define k (flfloor (fl+ c (fl* u (fl+ b (fl/ (fl* 2.0 a) us)))))) + (cond [(or (k . fl< . 0.0) (k . fl> . n)) (loop)] + [(and (us . fl>= . 0.07) (v . fl<= . vr)) k] + [else + (let ([v (fl* v (fl/ α (fl+ b (fl/ a (fl* us us)))))]) + (define h (+ (fllog-factorial m) + (fllog-factorial (fl- n m)) + (- (fllog-factorial k)) + (- (fllog-factorial (fl- n k))) + (fl* (fl- k m) (fllog (fl/ p q))))) + (cond [((fllog v) . fl<= . h) k] + [else (loop)]))]))) + (if s? (fl- n k) k))))) -(: flbinomial-random-normal (Flonum Flonum -> Flonum)) -(define (flbinomial-random-normal n p) +(: flbinomial-sample-normal (Flonum Flonum Natural -> FlVector)) +(define (flbinomial-sample-normal n p m) (define q (fl- 1.0 p)) (define μ (fl- (fl* (fl+ n 1.0) p) 0.5)) (define σ (flsqrt (* (+ 1.0 n) p q))) (define γ (fl/ (fl- q p) σ)) - (let loop () - (define z (flnormal-random 0.0 1.0)) - (define k (flround (fl+ μ (fl* σ (fl+ z (fl/ (fl* γ (fl- (fl* z z) 1.0)) 6.0)))))) - (if (and (k . fl>= . 0.0) (k . fl<= . n)) k (loop)))) + (build-flvector + m (λ (_) + (let loop () + (define z (unsafe-flvector-ref (flnormal-sample 0.0 1.0 1) 0)) + (define k (flround (fl+ μ (fl* σ (fl+ z (fl/ (fl* γ (fl- (fl* z z) 1.0)) 6.0)))))) + (if (and (k . fl>= . 0.0) (k . fl<= . n)) k (loop)))))) (: flbinomial-normal-appx-error-bound (Flonum Flonum -> Flonum)) ;; Returns a bound on the integrated difference between the normal and binomial cdfs @@ -78,16 +85,19 @@ Wolfgang Hormann. The Generation of Binomial Random Variates. (define q (fl- 1.0 p)) (fl/ (fl* 0.4784 (fl+ (fl* p p) (fl* q q))) (flsqrt (* n p q)))) -(: flbinomial-random (Flonum Flonum -> Flonum)) -(define (flbinomial-random n p) - (cond [(not (integer? n)) +nan.0] - [(n . fl<= . 0.0) (if (fl= n 0.0) 0.0 +nan.0)] - [(p . fl<= . 0.0) (if (fl= p 0.0) 0.0 +nan.0)] - [(p . fl>= . 1.0) (if (fl= p 1.0) n +nan.0)] +(: flbinomial-sample (Flonum Flonum Integer -> FlVector)) +(define (flbinomial-sample n p m) + (cond [(m . < . 0) (raise-argument-error 'flbinomial-sample "Natural" 2 n p m)] + [(or (not (integer? n)) (n . fl< . 0.0) (p . fl< . 0.0) (p . fl> . 1.0)) + (build-flvector m (λ (_) +nan.0))] + [(or (fl= n 0.0) (fl= p 0.0)) + (build-flvector m (λ (_) 0.0))] + [(fl= p 1.0) + (build-flvector m (λ (_) n))] [(and (n . fl> . 1e8) ((flbinomial-normal-appx-error-bound n p) . fl< . (flexp -10.0))) - (flbinomial-random-normal n p)] + (flbinomial-sample-normal n p m)] [((fl* n (flmin p (fl- 1.0 p))) . fl>= . 10.0) - (flbinomial-random-hormann n p)] + (flbinomial-sample-hormann n p m)] [else - (flbinomial-random-small n p)])) + (flbinomial-sample-small n p m)])) diff --git a/collects/math/private/distributions/impl/gamma-random.rkt b/collects/math/private/distributions/impl/gamma-random.rkt index d5b6846f78..ddb4a421ec 100644 --- a/collects/math/private/distributions/impl/gamma-random.rkt +++ b/collects/math/private/distributions/impl/gamma-random.rkt @@ -19,73 +19,96 @@ For others: sum of Gamma and Exponential variables, Normal approximation. |# (require "../../../flonum.rkt" + "../../../vector.rkt" + "../../unsafe.rkt" "normal-random.rkt") -(provide standard-flgamma-random) +(provide flgamma-sample) -(: standard-flgamma-random-small (Float -> Float)) +(: flgamma-sample-small (Flonum Flonum Natural -> FlVector)) ;; Ahrens and Dieter's rejection method ;; Good for 0.0 <= k < 1.0 -(define (standard-flgamma-random-small k) - (cond [(fl= k 0.0) 0.0] - [else - (define e (fl+ 1.0 (fl* k (flexp -1.0)))) - (let loop () - (define p (fl* e (random))) - (define q (fllog (random))) - (cond [(p . fl>= . 1.0) - (define x (- (fllog (fl/ (fl- e p) k)))) - (cond [(q . fl<= . (fl* (fl- k 1.0) (fllog x))) x] - [else (loop)])] - [else - (define x (flexpt p (fl/ 1.0 k))) - (cond [(q . fl<= . (- x)) x] - [else (loop)])]))])) +(define (flgamma-sample-small k s n) + (cond + [(fl= k 0.0) (make-flvector n 0.0)] + [else + (define e (fl+ 1.0 (fl* k (flexp -1.0)))) + (define k-1 (fl- k 1.0)) + (define 1/k (fl/ 1.0 k)) + (build-flvector + n (λ (_) + (let loop () + (define p (fl* e (random))) + (define q (fllog (random))) + (cond [(p . fl>= . 1.0) + (define x (- (fllog (fl/ (fl- e p) k)))) + (cond [(q . fl<= . (fl* k-1 (fllog x))) (fl* s x)] + [else (loop)])] + [else + (define x (flexpt p 1/k)) + (cond [(q . fl<= . (- x)) (fl* s x)] + [else (loop)])]))))])) -(: standard-flgamma-random-1-2 (Float -> Float)) +(: flgamma-sample-1-2 (Flonum Flonum Natural -> FlVector)) ;; Sum of Gamma and Exponential rvs ;; Good for 1.0 <= k < 2.0 -(define (standard-flgamma-random-1-2 k) - (fl- (standard-flgamma-random-small (fl- k 1.0)) - (fllog (random)))) +(define (flgamma-sample-1-2 k s n) + (define xs (flgamma-sample-small (fl- k 1.0) s n)) + (for ([i (in-range n)]) + (define x (unsafe-flvector-ref xs i)) + (unsafe-flvector-set! xs i (fl- x (fl* s (fllog (random)))))) + xs) -(: standard-flgamma-random-2-3 (Float -> Float)) +(: flgamma-sample-2-3 (Flonum Flonum Natural -> FlVector)) ;; Sum of Gamma and two Exponential rvs ;; Good for 2.0 <= k < 3.0 -(define (standard-flgamma-random-2-3 k) - (fl- (fl- (standard-flgamma-random-small (fl- k 2.0)) - (fllog (random))) - (fllog (random)))) +(define (flgamma-sample-2-3 k s n) + (define xs (flgamma-sample-small (fl- k 2.0) s n)) + (for ([i (in-range n)]) + (define x (unsafe-flvector-ref xs i)) + (unsafe-flvector-set! xs i (fl- x (fl* s (fl+ (fllog (random)) (fllog (random))))))) + xs) -(: standard-flgamma-random-large (Float -> Float)) +(: flgamma-sample-large (Flonum Flonum Natural -> FlVector)) ;; Tadikamalla's rejection method (Laplacian candidate) ;; Good for 1.0 <= k < huge (where "huge" causes the floating-point ops to behave badly) ;; Faster than the other methods for large k when k >= 3 or so (Laplacian left tail generates too ;; many negative candidates, which are rejected, when k < 3) -(define (standard-flgamma-random-large k) +(define (flgamma-sample-large k s n) (define A (fl- k 1.0)) (define B (fl+ 0.5 (fl* 0.5 (flsqrt (fl- (fl* 4.0 k) 3.0))))) (define C (fl/ (fl* A (fl+ 1.0 B)) B)) (define D (fl/ (fl- B 1.0) (fl* A B))) - (let loop () - (define lx (flmax -max.0 (fllog (random)))) - (define x (fl+ A (fl* B (if ((random) . fl< . 0.5) (- lx) lx)))) - (cond [(x . fl< . 0.0) (loop)] - [((fllog (random)) . fl<= . (fl+ (fl- (fl- (* A (fllog (* D x))) x) lx) C)) x] - [else (loop)]))) + (build-flvector + n (λ (_) + (let loop () + (define lx (flmax -max.0 (fllog (random)))) + (define x (fl+ A (fl* B (if ((random) . fl< . 0.5) (- lx) lx)))) + (cond [(x . fl< . 0.0) + (loop)] + [((fllog (random)) . fl<= . (fl+ (fl- (fl- (fl* A (fllog (fl* D x))) x) lx) C)) + (fl* s x)] + [else + (loop)]))))) -(: standard-flgamma-random-huge (Float -> Float)) +(: flgamma-sample-huge (Flonum Flonum Natural -> FlVector)) ;; Normal approximation ;; Good for 1e10 <= k <= +inf.0 -(define (standard-flgamma-random-huge k) - (cond [(fl= k +inf.0) +inf.0] - [else (flmax 0.0 (fl+ k (fl* (flsqrt k) (standard-flnormal-random))))])) +(define (flgamma-sample-huge k s n) + (cond [(fl= k +inf.0) (build-flvector n (λ (_) +inf.0))] + [else + (define xs (flnormal-sample k (flsqrt k) n)) + (for ([i (in-range n)]) + (define x (unsafe-flvector-ref xs i)) + (unsafe-flvector-set! xs i (flmax 0.0 (fl* s x)))) + xs])) -(: standard-flgamma-random (Float -> Float)) -(define (standard-flgamma-random k) - (cond [(k . fl>= . 1e10) (standard-flgamma-random-huge k)] - [(k . fl>= . 3.0) (standard-flgamma-random-large k)] - [(k . fl>= . 2.0) (standard-flgamma-random-2-3 k)] - [(k . fl>= . 1.0) (standard-flgamma-random-1-2 k)] - [(k . fl>= . 0.0) (standard-flgamma-random-small k)] - [else +nan.0])) +(: flgamma-sample (Flonum Flonum Integer -> FlVector)) +(define (flgamma-sample k s n) + (cond [(n . < . 0) (raise-argument-error 'flgamma-sample "Natural" 2 k s n)] + [(k . fl>= . 1e10) (flgamma-sample-huge k s n)] + [(k . fl>= . 3.0) (flgamma-sample-large k s n)] + [(k . fl>= . 2.0) (flgamma-sample-2-3 k s n)] + [(k . fl>= . 1.0) (flgamma-sample-1-2 k s n)] + [(k . fl>= . 0.0) (flgamma-sample-small k s n)] + [else (build-flvector n (λ (_) +nan.0))])) diff --git a/collects/math/private/distributions/impl/normal-random.rkt b/collects/math/private/distributions/impl/normal-random.rkt index e6651c6afe..9fa37f48b9 100644 --- a/collects/math/private/distributions/impl/normal-random.rkt +++ b/collects/math/private/distributions/impl/normal-random.rkt @@ -1,9 +1,10 @@ #lang typed/racket/base (require "../../../flonum.rkt" - "../../../base.rkt") + "../../../base.rkt" + "../../../vector.rkt") -(provide box-muller-transform standard-flnormal-random) +(provide flnormal-sample) (: box-muller-transform (Float Float -> Float)) (define (box-muller-transform x y) @@ -11,14 +12,19 @@ [else (fl* (flsqrt (fl* -2.0 (fllog x))) (flsin (fl* (fl* 2.0 pi) y)))])) -(: standard-flnormal-random (-> Float)) +(: flnormal-sample (Flonum Flonum Integer -> FlVector)) ;; The Box-Muller method has an bad reputation, but undeservedly: ;; 1. There's nothing unstable about the floating-point implementation of the transform ;; 2. It has good tail behavior (i.e. it can return very unlikely numbers) ;; 3. With today's RNGs, there's no need to worry about generating two random numbers ;; 4. With today's FPUs, there's no need to worry about computing `log' and `sin' (sheesh) ;; Points in favor: it's simple and fast -(define (standard-flnormal-random) - (let loop () - (define r (box-muller-transform (random) (random))) - (if (rational? r) r (loop)))) +(define (flnormal-sample c s n) + (cond [(n . < . 0) (raise-argument-error 'flnormal-sample "Natural" 2 c s n)] + [else + (build-flvector + n (λ (_) + (let loop () + (define r (box-muller-transform (random) (random))) + (if (rational? r) (fl+ c (fl* s r)) (loop)))))])) + diff --git a/collects/math/private/distributions/impl/poisson-random.rkt b/collects/math/private/distributions/impl/poisson-random.rkt index 82775be1c1..3541393ca2 100644 --- a/collects/math/private/distributions/impl/poisson-random.rkt +++ b/collects/math/private/distributions/impl/poisson-random.rkt @@ -3,52 +3,59 @@ (require racket/fixnum "../../../flonum.rkt" "../../../base.rkt" + "../../../vector.rkt" "../../functions/log-gamma.rkt") -(provide flpoisson-random) +(provide flpoisson-sample) -(: flpoisson-random-small (Flonum -> Flonum)) +(: flpoisson-sample-small (Flonum Natural -> FlVector)) ;; Good for l < -log(+min.0); suffers from underflow otherwise ;; O(l) in time and the number of uniform random variates used -(define (flpoisson-random-small l) +(define (flpoisson-sample-small l n) (define exp-l (flexp (- l))) - (let loop ([k 0.0] [p 1.0]) - (define u (random)) - (let ([p (fl* p u)]) - (cond [(p . fl<= . exp-l) k] - [else (loop (fl+ k 1.0) p)])))) + (build-flvector + n (λ (_) + (let loop ([k 0.0] [p 1.0]) + (define u (random)) + (let ([p (fl* p u)]) + (cond [(p . fl<= . exp-l) k] + [else (loop (fl+ k 1.0) p)])))))) -(: flpoisson-random-atkinson (Flonum -> Flonum)) +(: flpoisson-sample-atkinson (Flonum Natural -> FlVector)) ;; For l < 5, converges so slowly it's not even worth considering; fast for l > 30 or so, ;; just as fast as flpoisson-random-small for l > 9 ;; For l > 9, uses 5 random variates on average, which converges to 1 as l grows -(define (flpoisson-random-atkinson l) +(define (flpoisson-sample-atkinson l n) (define c (fl- 0.767 (fl/ 3.36 l))) (define beta (fl/ pi (flsqrt (fl* 3.0 l)))) (define alpha (fl* beta l)) (define k (fl- (fl- (fllog c) l) (fllog beta))) (define log-l (fllog l)) - (let loop () - (define u (random)) - (define x (fl/ (fl- alpha (fllog (fl/ (fl- 1.0 u) u))) beta)) - (define n (flfloor (fl+ x 0.5))) - (cond [(n . fl< . 0.0) (loop)] - [else - (define v (random)) - (define y (fl- alpha (fl* beta x))) - (define 1+exp-y (fl+ 1.0 (flexp y))) - (define lhs (fl+ y (fllog (fl/ (fl/ v 1+exp-y) 1+exp-y)))) - (define rhs (fl- (fl+ k (fl* n log-l)) (fllog-gamma (fl+ n 1.0)))) - (cond [(lhs . fl<= . rhs) n] - [else (loop)])]))) + (build-flvector + n (λ (_) + (let loop () + (define u (random)) + (define x (fl/ (fl- alpha (fllog (fl/ (fl- 1.0 u) u))) beta)) + (define n (flfloor (fl+ x 0.5))) + (cond [(n . fl< . 0.0) (loop)] + [else + (define v (random)) + (define y (fl- alpha (fl* beta x))) + (define 1+exp-y (fl+ 1.0 (flexp y))) + (define lhs (fl+ y (fllog (fl/ (fl/ v 1+exp-y) 1+exp-y)))) + (define rhs (fl- (fl+ k (fl* n log-l)) (fllog-gamma (fl+ n 1.0)))) + (cond [(lhs . fl<= . rhs) n] + [else (loop)])]))))) -(: flpoisson-random (Flonum -> Flonum)) -(define (flpoisson-random l) - (cond [(l . fl<= . 0.0) (if (fl= l 0.0) 0.0 +nan.0)] - [(l . fl<= . 9.0) (flpoisson-random-small l)] - [(l . fl<= . 1e35) (flpoisson-random-atkinson l)] +(: flpoisson-sample (Flonum Integer -> FlVector)) +(define (flpoisson-sample l n) + (cond [(n . < . 0) (raise-argument-error 'flpoisson-sample "Natural" 1 l n)] + [(l . fl< . 0.0) (build-flvector n (λ (_) +nan.0))] + [(l . fl= . 0.0) (build-flvector n (λ (_) 0.0))] + [(l . fl<= . 9.0) (flpoisson-sample-small l n)] + [(l . fl<= . 1e35) (flpoisson-sample-atkinson l n)] [else ;; At this point, the flonums are so sparse that: ;; 1. The mean `l' must be an integer; it is therefore also the mode ;; 2. The only flonum integer with probability >= +min.0 is `l' - l])) + (build-flvector n (λ (_) l))])) diff --git a/collects/math/private/distributions/logistic-dist.rkt b/collects/math/private/distributions/logistic-dist.rkt index c4fb7e6684..20df61931c 100644 --- a/collects/math/private/distributions/logistic-dist.rkt +++ b/collects/math/private/distributions/logistic-dist.rkt @@ -3,14 +3,16 @@ (require racket/performance-hint racket/promise "../../flonum.rkt" + "../../vector.rkt" + "../unsafe.rkt" "dist-struct.rkt" "utils.rkt") (provide fllogistic-pdf fllogistic-cdf fllogistic-inv-cdf - fllogistic-random - Logistic-Dist logistic-dist logistic-dist? logistic-dist-center logistic-dist-scale) + fllogistic-sample + Logistic-Dist logistic-dist logistic-dist-mean logistic-dist-scale) (: fllogistic-pdf (Float Float Float Any -> Float)) (define fllogistic-pdf @@ -52,17 +54,23 @@ (: fllogistic-inv-cdf (Float Float Float Any Any -> Float)) (define fllogistic-inv-cdf (make-symmetric-location-scale-flinv-cdf standard-fllogistic-inv-cdf)) -(: fllogistic-random (Float Float -> Float)) -(define fllogistic-random (make-symmetric-location-scale-flrandom standard-fllogistic-inv-cdf)) +(: fllogistic-sample-single (Float Float -> Float)) +(define fllogistic-sample-single + (make-symmetric-location-scale-flrandom standard-fllogistic-inv-cdf)) + +(: fllogistic-sample (Flonum Flonum Integer -> FlVector)) +(define (fllogistic-sample c s n) + (cond [(n . < . 0) (raise-argument-error 'fllogistic-sample "Natural" 2 c s n)] + [else (build-flvector n (λ (_) (fllogistic-sample-single c s)))])) ;; =================================================================================================== ;; Distribution object +(define-real-dist: logistic-dist Logistic-Dist + logistic-dist-struct ([mean : Flonum] [scale : Flonum])) + (begin-encourage-inline - (define-distribution-type: Logistic-Dist (Ordered-Dist Real Flonum) - logistic-dist ([center : Float] [scale : Float])) - (: logistic-dist (case-> (-> Logistic-Dist) (Real -> Logistic-Dist) (Real Real -> Logistic-Dist))) @@ -74,7 +82,9 @@ (fllogistic-cdf c s (fl x) log? 1-p?))) (define inv-cdf (opt-lambda: ([p : Real] [log? : Any #f] [1-p? : Any #f]) (fllogistic-inv-cdf c s (fl p) log? 1-p?))) - (define (random) (fllogistic-random c s)) - (make-logistic-dist pdf random cdf inv-cdf -inf.0 +inf.0 (delay c) c s))) + (define sample (case-lambda: + [() (unsafe-flvector-ref (fllogistic-sample c s 1) 0)] + [([n : Integer]) (flvector->list (fllogistic-sample c s n))])) + (logistic-dist-struct pdf sample cdf inv-cdf -inf.0 +inf.0 (delay c) c s))) ) diff --git a/collects/math/private/distributions/normal-dist.rkt b/collects/math/private/distributions/normal-dist.rkt index c354e6a4ec..c6ae10e1d4 100644 --- a/collects/math/private/distributions/normal-dist.rkt +++ b/collects/math/private/distributions/normal-dist.rkt @@ -3,6 +3,8 @@ (require racket/performance-hint racket/promise "../../flonum.rkt" + "../../vector.rkt" + "../unsafe.rkt" "impl/normal-pdf.rkt" "impl/normal-cdf.rkt" "impl/normal-inv-cdf.rkt" @@ -13,8 +15,8 @@ (provide flnormal-pdf flnormal-cdf flnormal-inv-cdf - flnormal-random - Normal-Dist normal-dist normal-dist? normal-dist-mean normal-dist-stddev) + flnormal-sample + Normal-Dist normal-dist normal-dist-mean normal-dist-stddev) (: flnormal-pdf (Float Float Float Any -> Float)) (define flnormal-pdf @@ -37,18 +39,13 @@ (cond [log? (standard-flnormal-inv-log-cdf q)] [else (standard-flnormal-inv-cdf q)])))) -(: flnormal-random (Float Float -> Float)) -(define (flnormal-random c s) - (fl+ c (fl* s (standard-flnormal-random)))) - ;; =================================================================================================== ;; Distribution object +(define-real-dist: normal-dist Normal-Dist + normal-dist-struct ([mean : Flonum] [stddev : Flonum])) + (begin-encourage-inline - - (define-distribution-type: Normal-Dist (Ordered-Dist Real Flonum) - normal-dist ([mean : Float] [stddev : Float])) - (: normal-dist (case-> (-> Normal-Dist) (Real -> Normal-Dist) (Real Real -> Normal-Dist))) @@ -60,7 +57,9 @@ (flnormal-cdf c s (fl x) log? 1-p?))) (define inv-cdf (opt-lambda: ([p : Real] [log? : Any #f] [1-p? : Any #f]) (flnormal-inv-cdf c s (fl p) log? 1-p?))) - (define (random) (flnormal-random c s)) - (make-normal-dist pdf random cdf inv-cdf -inf.0 +inf.0 (delay c) c s))) + (define sample (case-lambda: + [() (unsafe-flvector-ref (flnormal-sample c s 1) 0)] + [([n : Integer]) (flvector->list (flnormal-sample c s n))])) + (normal-dist-struct pdf sample cdf inv-cdf -inf.0 +inf.0 (delay c) c s))) ) diff --git a/collects/math/private/distributions/poisson-dist.rkt b/collects/math/private/distributions/poisson-dist.rkt index 8f4603ae4b..fd7e6299df 100644 --- a/collects/math/private/distributions/poisson-dist.rkt +++ b/collects/math/private/distributions/poisson-dist.rkt @@ -3,6 +3,8 @@ (require racket/performance-hint racket/promise "../../flonum.rkt" + "../../vector.rkt" + "../unsafe.rkt" "../functions/incomplete-gamma.rkt" (prefix-in impl: "impl/poisson-pdf.rkt") "impl/poisson-random.rkt" @@ -13,9 +15,9 @@ (provide flpoisson-pdf flpoisson-cdf flpoisson-inv-cdf - flpoisson-random + flpoisson-sample flpoisson-median - Poisson-Dist poisson-dist poisson-dist? poisson-dist-mean) + Poisson-Dist poisson-dist poisson-dist-mean) (: flpoisson-pdf (Flonum Flonum Any -> Flonum)) (define (flpoisson-pdf l k log?) @@ -53,11 +55,11 @@ (cond [(fl= k 0.0) k] [else (if ((flpoisson-cdf l (- k 1.0) #f #f) . fl< . 0.5) k (- k 1.0))])])) +(define-real-dist: poisson-dist Poisson-Dist + poisson-dist-struct ([mean : Flonum])) + (begin-encourage-inline - (define-distribution-type: Poisson-Dist (Ordered-Dist Real Flonum) - poisson-dist ([mean : Flonum])) - (: poisson-dist (case-> (-> Poisson-Dist) (Real -> Poisson-Dist))) (define (poisson-dist [l 0.5]) @@ -68,9 +70,9 @@ (flpoisson-cdf l (fl k) log? 1-p?))) (define inv-cdf (opt-lambda: ([p : Real] [log? : Any #f] [1-p? : Any #f]) (flpoisson-inv-cdf l (fl p) log? 1-p?))) - (define (random) (flpoisson-random l)) - (make-poisson-dist pdf random cdf inv-cdf - 0.0 +inf.0 (delay (flpoisson-median l)) - l))) + (define sample (case-lambda: + [() (unsafe-flvector-ref (flpoisson-sample l 1) 0)] + [([n : Integer]) (flvector->list (flpoisson-sample l n))])) + (poisson-dist-struct pdf sample cdf inv-cdf 0.0 +inf.0 (delay (flpoisson-median l)) l))) ) diff --git a/collects/math/private/distributions/triangle-dist.rkt b/collects/math/private/distributions/triangle-dist.rkt index 28984047c6..7aa8fd9cd5 100644 --- a/collects/math/private/distributions/triangle-dist.rkt +++ b/collects/math/private/distributions/triangle-dist.rkt @@ -2,8 +2,9 @@ (require racket/performance-hint racket/promise - racket/unsafe/ops "../../flonum.rkt" + "../../vector.rkt" + "../unsafe.rkt" "../inline-sort.rkt" "dist-struct.rkt" "utils.rkt") @@ -11,13 +12,12 @@ (provide fltriangle-pdf fltriangle-cdf fltriangle-inv-cdf - fltriangle-random - Triangular-Dist triangle-dist triangle-dist? - triangle-dist-min triangle-dist-max triangle-dist-center) + fltriangle-sample + Triangle-Dist triangle-dist triangle-dist-min triangle-dist-max triangle-dist-mode) (: flsort3 (Flonum Flonum Flonum -> (Values Flonum Flonum Flonum))) (begin-encourage-inline - (define (flsort3 a b c) (inline-sort unsafe-fl< a b c))) + (define (flsort3 a b c) (inline-sort fl< a b c))) (: unsafe-fltriangle-pdf (Float Float Float Float Any -> Float)) (define (unsafe-fltriangle-pdf a b c x log?) @@ -59,8 +59,8 @@ [(q . fl= . 1.0) b] [else +nan.0]))) -(: unsafe-fltriangle-random (Float Float Float -> Float)) -(define (unsafe-fltriangle-random a b c) +(: unsafe-fltriangle-sample-single (Float Float Float -> Float)) +(define (unsafe-fltriangle-sample-single a b c) (unsafe-fltriangle-inv-cdf a b c (fl* 0.5 (random)) #f ((random) . fl> . 0.5))) (begin-encourage-inline @@ -80,25 +80,27 @@ (let-values ([(a c b) (flsort3 a b c)]) (unsafe-fltriangle-inv-cdf a b c x log? 1-p?))) - (: fltriangle-random (Float Float Float -> Float)) - (define (fltriangle-random a b c) - (let-values ([(a c b) (flsort3 a b c)]) - (unsafe-fltriangle-random a b c))) + (: fltriangle-sample (Flonum Flonum Flonum Integer -> FlVector)) + (define (fltriangle-sample a b c n) + (cond [(n . < . 0) (raise-argument-error 'fltriangle-sample "Natural" 3 a b c n)] + [else + (let-values ([(a c b) (flsort3 a b c)]) + (build-flvector n (λ (_) (unsafe-fltriangle-sample-single a b c))))])) ) ;; =================================================================================================== ;; Distribution object +(define-real-dist: triangle-dist Triangle-Dist + triangle-dist-struct ([min : Float] [max : Float] [mode : Float])) + (begin-encourage-inline - (define-distribution-type: Triangular-Dist (Ordered-Dist Real Flonum) - triangle-dist ([min : Float] [max : Float] [center : Float])) - - (: triangle-dist (case-> (-> Triangular-Dist) - (Real -> Triangular-Dist) - (Real Real -> Triangular-Dist) - (Real Real Real -> Triangular-Dist))) + (: triangle-dist (case-> (-> Triangle-Dist) + (Real -> Triangle-Dist) + (Real Real -> Triangle-Dist) + (Real Real Real -> Triangle-Dist))) (define (triangle-dist [a 0.0] [b 1.0] [c (* 0.5 (+ a b))]) (let ([a (fl a)] [b (fl b)] [c (fl c)]) (let-values ([(a c b) (flsort3 a b c)]) @@ -108,9 +110,12 @@ (unsafe-fltriangle-cdf a b c (fl x) log? 1-p?))) (define inv-cdf (opt-lambda: ([p : Real] [log? : Any #f] [1-p? : Any #f]) (unsafe-fltriangle-inv-cdf a b c (fl p) log? 1-p?))) - (define (random) (unsafe-fltriangle-random a b c)) - (make-triangle-dist pdf random cdf inv-cdf - a b (delay (unsafe-fltriangle-inv-cdf a b c 0.5 #f #f)) - a b c)))) + (define sample (case-lambda: + [() (unsafe-flvector-ref (fltriangle-sample a b c 1) 0)] + [([n : Integer]) (flvector->list (fltriangle-sample a b c n))])) + (triangle-dist-struct + pdf sample cdf inv-cdf + a b (delay (unsafe-fltriangle-inv-cdf a b c 0.5 #f #f)) + a b c)))) ) diff --git a/collects/math/private/distributions/truncated-dist.rkt b/collects/math/private/distributions/truncated-dist.rkt index 6eafd63761..3e33375f97 100644 --- a/collects/math/private/distributions/truncated-dist.rkt +++ b/collects/math/private/distributions/truncated-dist.rkt @@ -8,13 +8,12 @@ (provide Truncated-Dist truncated-dist - truncated-dist? truncated-dist-min truncated-dist-max truncated-dist-original) -(define-distribution-type: Truncated-Dist (Ordered-Dist Real Flonum) - truncated-dist ([original : Real-Dist] [min : Float] [max : Float])) +(define-real-dist: truncated-dist Truncated-Dist + truncated-dist-struct ([original : Real-Dist] [min : Flonum] [max : Flonum])) (: truncated-dist (case-> (Real-Dist -> Truncated-Dist) (Real-Dist Real -> Truncated-Dist) @@ -29,20 +28,18 @@ (min (dist-max d) (max a b)))]) (unsafe-truncated-dist d a b))])) -(define-syntax-rule (make-inv-cdf-random inv-cdf) - (λ () (inv-cdf (fl* 0.5 (random)) #f ((random) . fl> . 0.5)))) - (: unsafe-truncated-dist (Real-Dist Float Float -> Truncated-Dist)) (define (unsafe-truncated-dist d a b) (define orig-pdf (dist-pdf d)) (define orig-cdf (dist-cdf d)) (define orig-inv-cdf (dist-inv-cdf d)) - (define orig-random (dist-random d)) + (define orig-sample (dist-sample d)) (define log-P_ab (delay (orig-cdf b #t #t))) - (: pdf Real-PDF) + (: pdf (case-> (Real -> Flonum) + (Real Any -> Flonum))) (define (pdf x [log? #f]) (let ([x (fl x)]) (define log-d @@ -51,7 +48,9 @@ [else (fl- (orig-pdf x #t) log-P_a (Real -> Flonum) + (Real Any -> Flonum) + (Real Any Any -> Flonum))) (define (cdf x [log? #f] [1-p? #f]) (let ([x (fl x)]) (define log-p @@ -67,7 +66,9 @@ log-P_a (Real -> Flonum) + (Real Any -> Flonum) + (Real Any Any -> Flonum))) (define (inv-cdf p [log? #f] [1-p? #f]) (let ([log-p (if log? (fl p) (fllog (fl p)))]) (cond @@ -86,18 +87,24 @@ #t #f)])])) (min b (max a x))]))) - (: random (-> Float)) - (define random + (: sample-single (-> Flonum)) + (define sample-single (cond [(log-P_a . 0.5)))] [else - (define (random) - (define x (orig-random)) + (define (sample-single) + (define x (orig-sample)) (cond [(and (a . fl<= . x) (x . fl<= . b)) x] - [else (random)])) - random])) + [else (sample-single)])) + sample-single])) + + (: sample (Sample Flonum)) + (define sample + (case-lambda: + [() (sample-single)] + [([n : Integer]) + (cond [(n . < . 0) (raise-argument-error 'truncated-dist-sample "Natural" n)] + [else (build-list n (λ (_) (sample-single)))])])) ;; Finally put it together - (make-truncated-dist pdf random cdf inv-cdf - a b (delay (inv-cdf 0.5)) - d a b)) + (truncated-dist-struct pdf sample cdf inv-cdf a b (delay (inv-cdf 0.5)) d a b)) diff --git a/collects/math/private/distributions/uniform-dist.rkt b/collects/math/private/distributions/uniform-dist.rkt index 443d2f577b..70b5205e8d 100644 --- a/collects/math/private/distributions/uniform-dist.rkt +++ b/collects/math/private/distributions/uniform-dist.rkt @@ -3,14 +3,16 @@ (require racket/performance-hint racket/promise "../../flonum.rkt" + "../../vector.rkt" + "../unsafe.rkt" "dist-struct.rkt" "utils.rkt") (provide fluniform-pdf fluniform-cdf fluniform-inv-cdf - fluniform-random - Uniform-Dist uniform-dist uniform-dist? uniform-dist-min uniform-dist-max) + fluniform-sample + Uniform-Dist uniform-dist uniform-dist-min uniform-dist-max) (: unsafe-fluniform-pdf (Float Float Float Any -> Float)) (define (unsafe-fluniform-pdf a b x log?) @@ -44,18 +46,20 @@ [(fl= q 0.0) a] [else (fl+ (fl* (fl- b a) q) a)]))])) -(: unsafe-fluniform-random (Float Float -> Float)) +(: unsafe-fluniform-sample (Float Float Natural -> FlVector)) ;; Chooses a random flonum in [a,b] in a way that preserves the precision of the random flonum ;; returned by (random) -(define (unsafe-fluniform-random a b) +(define (unsafe-fluniform-sample a b n) (define d (fl- b a)) (cond ;; Both positive, `a' smaller in magnitude - [(a . fl>= . 0.0) (fl+ a (fl* d (random)))] + [(a . fl>= . 0.0) (build-flvector n (λ (_) (fl+ a (fl* d (random)))))] ;; Both negative, `b' smaller in magnitude - [(b . fl<= . 0.0) (fl+ b (fl* (- d) (random)))] + [(b . fl<= . 0.0) (build-flvector n (λ (_) (fl- b (fl* d (random)))))] ;; Straddle 0 case - [else (if ((fl* d (random)) . fl> . b) (fl* a (random)) (fl* b (random)))])) + [else + (build-flvector + n (λ (_) (if ((fl* d (random)) . fl> . b) (fl* a (random)) (fl* b (random)))))])) (begin-encourage-inline @@ -71,20 +75,21 @@ (define (fluniform-inv-cdf a b q log? 1-p?) (unsafe-fluniform-inv-cdf (flmin a b) (flmax a b) q log? 1-p?)) - (: fluniform-random (Float Float -> Float)) - (define (fluniform-random a b) - (unsafe-fluniform-random (flmin a b) (flmax a b))) + (: fluniform-sample (Float Float Integer -> FlVector)) + (define (fluniform-sample a b n) + (cond [(n . < . 0) (raise-argument-error 'fluniform-sample "Natural" 2 a b n)] + [else (unsafe-fluniform-sample (flmin a b) (flmax a b) n)])) ) ;; =================================================================================================== ;; Distribution object +(define-real-dist: uniform-dist Uniform-Dist + uniform-dist-struct ([min : Float] [max : Float])) + (begin-encourage-inline - (define-distribution-type: Uniform-Dist (Ordered-Dist Real Flonum) - uniform-dist ([min : Float] [max : Float])) - (: uniform-dist (case-> (-> Uniform-Dist) (Real -> Uniform-Dist) (Real Real -> Uniform-Dist))) @@ -101,7 +106,9 @@ (unsafe-fluniform-cdf a b (fl x) log? 1-p?))) (define inv-cdf (opt-lambda: ([p : Real] [log? : Any #f] [1-p? : Any #f]) (unsafe-fluniform-inv-cdf a b (fl p) log? 1-p?))) - (define (random) (fluniform-random a b)) - (make-uniform-dist pdf random cdf inv-cdf a b (delay (fl* 0.5 (fl+ a b))) a b)))])) + (define sample (case-lambda: + [() (unsafe-flvector-ref (fluniform-sample a b 1) 0)] + [([n : Integer]) (flvector->list (fluniform-sample a b n))])) + (uniform-dist-struct pdf sample cdf inv-cdf a b (delay (fl* 0.5 (fl+ a b))) a b)))])) ) diff --git a/collects/math/private/distributions/utils.rkt b/collects/math/private/distributions/utils.rkt index 97dca21538..c9a3b26a0e 100644 --- a/collects/math/private/distributions/utils.rkt +++ b/collects/math/private/distributions/utils.rkt @@ -5,6 +5,7 @@ racket/string) racket/performance-hint "../../flonum.rkt" + "../utils.rkt" "impl/delta-dist.rkt" "dist-struct.rkt") @@ -31,19 +32,40 @@ ;; =================================================================================================== -(define-syntax (define-distribution-type: stx) - (syntax-case stx () - [(_ (type-name T ...) (parent-type-name In Out) - (A B) name ([arg-names arg-opts ...] ...)) +(define-syntax (define-real-dist: stx) + (syntax-case stx (:) + [(_ name type-name struct-name ([arg-names : arg-types] ...) struct-opts ...) (let ([arg-name-lst (syntax->list #'(arg-names ...))]) - (with-syntax* ([internal-type-name (format-id #'type-name "~a-Struct" #'type-name)] - [(internal-proc-names ...) (map (λ (arg-name) + (with-syntax* ([(struct-proc-names ...) + (map (λ (arg-name) (format-id #'struct-name "~a-~a" #'struct-name arg-name)) + arg-name-lst)] + [(proc-names ...) + (map (λ (arg-name) (format-id #'name "~a-~a" #'name arg-name)) + arg-name-lst)]) + (syntax/loc stx + (begin-encourage-inline + (struct: (In Out) struct-name ordered-dist ([arg-names : arg-types] ...) struct-opts ... + #:property prop:custom-print-quotable 'never + #:property prop:custom-write + (λ (v port mode) + (pretty-print-constructor 'name (list (struct-proc-names v) ...) port mode))) + (define-type type-name (struct-name Real Flonum)) + (: proc-names (type-name -> arg-types)) ... + (define (proc-names v) (struct-proc-names v)) ...))))])) + +(define-syntax (define-distribution-type: stx) + (syntax-case stx (:) + [(_ (type-name T ...) (parent-type-name In Out) + (A B) name ([arg-names : arg-type] ...)) + (let ([arg-name-lst (syntax->list #'(arg-names ...))]) + (with-syntax* ([struct-type-name (format-id #'type-name "~a-struct" #'type-name)] + [(struct-proc-names ...) (map (λ (arg-name) (format-id #'type-name "~a-~a" - #'internal-type-name + #'struct-type-name arg-name)) arg-name-lst)] - [internal-pred-name (format-id #'type-name "~a?" #'internal-type-name)] + [struct-pred-name (format-id #'type-name "~a?" #'struct-type-name)] [make-name (format-id #'name "make-~a" #'name)] [(proc-names ...) (map (λ (arg-name) (format-id #'name "~a-~a" #'name arg-name)) @@ -55,15 +77,16 @@ ")")]) (syntax/loc stx (begin - (struct: (A B) internal-type-name parent-type-name ([arg-names arg-opts ...] ...) + (struct: (A B) struct-type-name parent-type-name ([arg-names : arg-types] ...) #:property prop:custom-print-quotable 'never #:property prop:custom-write (λ (v port write?) (fprintf port format-str 'name (proc-names v) ...))) - (define-type (type-name T ...) (internal-type-name In Out)) - (define proc-names internal-proc-names) ... - (define make-name internal-type-name) - (define pred-name internal-pred-name)))))] + (define-type (type-name T ...) (struct-type-name In Out)) + (: proc-names ((type-name T ...) -> arg-types)) ... + (define (proc-names d) (struct-proc-names d)) ... + (define make-name struct-type-name) + (define pred-name struct-pred-name)))))] [(_ type-name (parent-type-name In Out) name ([arg-names arg-opts ...] ...)) (syntax/loc stx (define-distribution-type: (type-name) (parent-type-name In Out) diff --git a/collects/math/private/flonum/flonum-log.rkt b/collects/math/private/flonum/flonum-log.rkt index aa38f08a87..7ef2abac7d 100644 --- a/collects/math/private/flonum/flonum-log.rkt +++ b/collects/math/private/flonum/flonum-log.rkt @@ -57,7 +57,8 @@ (: lg- (Float Float -> Float)) (define (lg- log-x log-y) - (cond [(log-y . fl> . log-x) +nan.0] + (cond [(log-x . fl< . log-y) +nan.0] + [(fl= log-x -inf.0) -inf.0] [else (fl+ log-x (lg1- (fl- log-y log-x)))])) ) ; begin-encourage-inline diff --git a/collects/math/private/utils.rkt b/collects/math/private/utils.rkt index 0c595506fc..29223b75d3 100644 --- a/collects/math/private/utils.rkt +++ b/collects/math/private/utils.rkt @@ -1,5 +1,61 @@ #lang typed/racket/base -(require) +(require racket/pretty) -(provide (all-defined-out)) +(provide pretty-print-constructor) + +(: port-next-column (Output-Port -> Natural)) +;; Helper to avoid the annoying #f column value +(define (port-next-column port) + (define-values (_line col _pos) (port-next-location port)) + (if col col 0)) + +(define-type Constructor-Layout (U 'one-line 'multi-line)) + +(: pretty-print-constructor (Symbol (Listof Any) Output-Port (U #t #f 0 1) -> Any)) +(define (pretty-print-constructor name args port mode) + ;; Called to print arguments; may recur (e.g. printing constructed arguments) + ;; We never have to consider the `mode' argument again after defining `recur-print' + (define recur-print + (cond [(not mode) display] + [(integer? mode) (λ: ([p : Any] [port : Output-Port]) + (print p port mode))] ; pass the quote depth through + [else write])) + + (define cols (pretty-print-columns)) + + (: print-all (Output-Port Constructor-Layout -> Any)) + (define (print-all port layout) + ;; Get current column so we can indent new lines at least that far + (define col (port-next-column port)) + ;; Print the constructor name + (write-string (format "(~a" name) port) + (for ([arg (in-list args)]) + (case layout + [(one-line) (write-string " " port)] + [else (pretty-print-newline port (assert cols integer?)) + (write-string (make-string (+ col 1) #\space) port)]) + (recur-print arg port)) + (write-string ")" port)) + + ;; See what the printer has in mind for us this time + (cond [(and (pretty-printing) (integer? cols)) + ;; Line-width-constrained pretty-printing: woo woo! + (let/ec: return : Any ; used as a return statement + ;; Wrap the port with a tentative one, in case compact layout overflows lines + (define: tport : Output-Port + (make-tentative-pretty-print-output-port + port + (max 0 (- cols 1)) ; width: make sure there's room for the closing delimiter + (λ () ; failure thunk + ;; Reset accumulated graph state + (tentative-pretty-print-port-cancel (assert tport output-port?)) + ;; Compact layout failed, so print in multi-line layout + (return (print-all port 'multi-line))))) + ;; Try printing on one line + (print-all tport 'one-line) + ;; If a line overflows, the failure thunk returns past this + (tentative-pretty-print-port-transfer tport port))] + [else + ;; No pretty printer, or printing to infinite-width lines, so print on one line + (print-all port 'one-line)])) diff --git a/collects/math/scribblings/math-array.scrbl b/collects/math/scribblings/math-array.scrbl new file mode 100644 index 0000000000..c1fbb2139e --- /dev/null +++ b/collects/math/scribblings/math-array.scrbl @@ -0,0 +1,92 @@ +#lang scribble/manual + +@(require scribble/eval + racket/sandbox + (for-label racket/base + math plot + (only-in typed/racket/base Flonum Real Boolean Any Listof Integer)) + "utils.rkt") + +@(define untyped-eval (make-untyped-math-eval)) + +@title[#:tag "arrays"]{Arrays} +@(author-neil) + +@defmodule[math/array] + +@section{Introduction} + +One of the most common ways to structure data is with an array: a grid of homogeneous, +independent elements, usually consisting of rows and columns. But an array data type +is often absent from functional languages' libraries. This is probably because arrays +are perceived as requiring users to operate on them using destructive updates, write +loops that micromanage array elements, and in general, stray far from the declarative +ideal. + +@margin-note{TODO: Cite Haskell array paper} +Normally, they do. However, experience in Python, and more recently Haskell, has shown +that providing the right data types and a rich collection of whole-array operations +allows working effectively with arrays in a functional, declarative style. As a bonus, +doing so opens the possibility of parallelizing nearly every operation. + +It requires a change in definition. The new definition is this: + +@bold{An @deftech{array} is just a function with a finite, rectangular domain.} + +Some arrays are mutable, some are lazy, some are strict, some are sparse, and some +do not even allocate space to store their elements. All are functions that can be +applied to indexes to retrieve elements. + +@subsection{Definitions} + +The domain of an array is determined by its @deftech{shape}, a vector of numbers that +describes the extent of each dimension. + +@examples[#:eval untyped-eval + (array-shape (array [0 1 2 3])) + (array-shape (array [[0 1] [2 3] [4 5]])) + (array-shape (array 0))] + +The function represented by the array is called its @deftech{procedure}, which accepts +vectors of indexes and returns elements. + +@examples[#:eval untyped-eval + (define arr (array [[0 1] [2 3]])) + (define proc (array-proc arr)) + (proc #(1 1)) + (array-ref arr #(1 1))] + +@section{Quick Start} + +@section{Array Types} + +@section{Array Constructors} + +array syntax + +(: make-array (All (A) (User-Indexes A -> (Array A)))) + +(: axis-index-array (User-Indexes Integer -> (Array Index))) + +(: index-array (User-Indexes -> (Array Index))) + +(: indexes-array (User-Indexes -> (Array Indexes))) + +(: diagonal-array (All (A) (Integer Integer A A -> (Array A)))) + + +@section{Pointwise Array Operations} + +@section{Array Folds} + +@section{Array Transformations} + +@section{Mutable Arrays} + +mutable-array syntax + +@section{Flonum Arrays} + +@section{Float-Complex Arrays} + +@(close-eval untyped-eval) diff --git a/collects/math/scribblings/math-distributions.scrbl b/collects/math/scribblings/math-distributions.scrbl new file mode 100644 index 0000000000..78d3b9d1f3 --- /dev/null +++ b/collects/math/scribblings/math-distributions.scrbl @@ -0,0 +1,517 @@ +#lang scribble/manual + +@(require scribble/eval + racket/sandbox + (for-label racket/base racket/promise + math plot + (only-in typed/racket/base + Flonum Real Boolean Any Listof Integer case-> -> Promise)) + "utils.rkt") + +@(define untyped-eval (make-untyped-math-eval)) +@interaction-eval[#:eval untyped-eval (require racket/list)] + +@title[#:tag "dist"]{Probability Distributions} +@(author-neil) + +@defmodule[math/distributions] + +@local-table-of-contents[] + +@section[#:tag "dist:intro"]{Introduction} + +The @racketmodname[math/distributions] module exports distribution objects and functions that +operate on them. + +A @deftech{distribution object} represents a probability distribution over a common domain, +such as the real numbers, integers, or a set of symbols. Their constructors correspond with +distribution families, such as the family of normal distributions. + +A distribution object has a density function (a @deftech{pdf}) and a procedure to generate random +samples. An @italic{ordered} distribution object additionally has a cumulative distribution function +(a @deftech{cdf}), and its generalized inverse (an @deftech{inverse cdf}). + +The following example creates an ordered distribution object representing a normal distribution +with mean 2 and standard deviation 5, computes an approximation of the probability of the +half-open interval (1/2,1], and computes another approximation from random samples: +@interaction[#:eval untyped-eval + (define d (normal-dist 2 5)) + (real-dist-prob d 0.5 1.0) + (define xs (sample d 10000)) + (fl (/ (count (λ (x) (and (1/2 . < . x) (x . <= . 1))) xs) + (length xs)))] + +This plots the pdf and a kernel density estimate of the pdf from random samples: +@interaction[#:eval untyped-eval + (plot (list (function (dist-pdf d) #:color 0 #:style 'dot) + (density xs)) + #:x-label "x" #:y-label "N(2,5)")] + +There are also higher-order distributions, which take other distributions as constructor +arguments. For example, the truncated distribution family returns a distribution like its +distribution argument, but sets probability outside an interval to 0 and renormalizes +the probabilities within the interval: +@interaction[#:eval untyped-eval + (define d-trunc (truncated-dist d -inf.0 5)) + (real-dist-prob d-trunc 5 6) + (real-dist-prob d-trunc 0.5 1.0) + (plot (list (function (dist-pdf d-trunc) #:color 0 #:style 'dot) + (density (sample d-trunc 1000))) + #:x-label "x" #:y-label "T(N(2,5),-∞,5)")] + +Because real distributions' cdfs represent the probability P[@italic{X} ≤ @italic{x}], they are +right-continuous: +@interaction[#:eval untyped-eval + (define d (geometric-dist 0.4)) + (define cdf (dist-cdf d)) + (plot (for/list ([i (in-range -1 7)]) + (define i+1-ε (flprev (+ i 1.0))) + (list (lines (list (vector i (cdf i)) + (vector i+1-ε (cdf i+1-ε))) + #:width 2) + (points (list (vector i (cdf i))) + #:sym 'fullcircle5 #:color 1) + (points (list (vector i+1-ε (cdf i+1-ε))) + #:sym 'fullcircle5 #:color 1 #:fill-color 0))) + #:x-min -0.5 #:x-max 6.5 #:y-min -0.05 #:y-max 1 + #:x-label "x" #:y-label "P[X ≤ x]")] +For convenience, cdfs are defined over the extended reals regardless of their distribution's +support, but their inverses return values only within the support: +@interaction[#:eval untyped-eval + (define inv-cdf (dist-inv-cdf d)) + (cdf +inf.0) + (cdf 1.5) + (cdf -inf.0) + (inv-cdf (cdf +inf.0)) + (inv-cdf (cdf 1.5)) + (inv-cdf (cdf -inf.0))] +A distribution's inverse cdf is defined on the interval [0,1] and is always left-continuous, +except possibly at 0 when its support is bounded on the left (as with @racket[geometric-dist]). + +Every pdf and cdf can return log densities and log probabilities, in case densities or probabilities +are too small to represent as flonums (i.e. are less than @racket[+min.0]): +@interaction[#:eval untyped-eval + (define d (normal-dist)) + (define pdf (dist-pdf d)) + (define cdf (dist-cdf d)) + (pdf 40.0) + (cdf -40.0) + (pdf 40.0 #t) + (cdf -40.0 #t)] +Additionally, every cdf can return upper-tail probabilities, which are always more accurate when +lower-tail probabilities are greater than @racket[0.5]: +@interaction[#:eval untyped-eval + (cdf 20.0) + (cdf 20.0 #f #t)] +Upper-tail probabilities can also be returned as logs in case probabilities are too small: +@interaction[#:eval untyped-eval + (cdf 40.0) + (cdf 40.0 #f #t) + (cdf 40.0 #t #t)] +Inverse cdfs accept log probabilities and upper-tail probabilities. + +The functions @racket[lg+] and @racket[lgsum], as well as others in @racketmodname[math/flonum], +perform arithmetic on log probabilities. + +@section{Basic Distribution Types and Operations} + +@defform[(PDF In)]{ +The type of probability density functions, or @tech{pdfs}, defined as +@racketblock[(case-> (In -> Flonum) + (In Any -> Flonum))] +The second argument defaults to @racket[#f]. When the second argument is not @racket[#f], a +pdf returns a log density. +} + +@defform[(Sample Out)]{ +The type of a distribution's sampling procedure, defined as +@racketblock[(case-> (-> Out) + (Integer -> (Listof Out)))] +When given a nonnegative integer @racket[n] as an argument, a sampling procedure returns a list of +random samples of length @racket[n]. +} + +@defstruct*[dist ([pdf (PDF In)] [sample (Sample Out)])]{ +The parent type of @tech{distribution objects}. The @racket[In] type parameter is the data type a +distribution accepts as arguments to its @tech{pdf}. The @racket[Out] type parameter is the data +type a distribution returns as random samples. + +@examples[#:eval untyped-eval + (dist? (discrete-dist '(a b c))) + (dist? (normal-dist))] +} + +@defproc*[([(sample [d (dist In Out)]) Out] + [(sample [d (dist In Out)] [n Integer]) (Listof Out)])]{ +An uncurried form of @racket[dist-sample]. +@examples[#:eval untyped-eval + (sample (exponential-dist)) + (sample (exponential-dist) 3)] +} + +@section{Ordered Distribution Types and Operations} + +@defform[(CDF In)]{ +The type of cumulative distribution functions, or @tech{cdfs}, defined as +@racketblock[(case-> (In -> Flonum) + (In Any -> Flonum) + (In Any Any -> Flonum))] +Both optional arguments default to @racket[#f]. + +Suppose @racket[cdf : (CDF Real)], and @racket[p = (cdf x log? 1-p?)]. The flonum @racket[p] +represents a probability. If @racket[log?] is a true value, @racket[p] is a log probability. +If @racket[1-p?] is a true value, @racket[p] is an upper-tail probability or upper-tail +log probability. +} + +@defform[(Inverse-CDF Out)]{ +The type of inverse cumulative distribution functions, or @tech{inverse cdfs}, defined as +@racketblock[(case-> (Real -> Out) + (Real Any -> Out) + (Real Any Any -> Out))] +Both optional arguments default to @racket[#f]. + +Suppose @racket[inv-cdf : (Inverse-CDF Flonum)], and @racket[x = (inv-cdf p log? 1-p?)]. +If @racket[log?] is a true value, @racket[inv-cdf] interprets @racket[p] as a log probability. +If @racket[1-p?] is a true value, @racket[inv-cdf] interprets @racket[p] as an upper-tail probability +or upper-tail log probability. +} + +@defstruct*[(ordered-dist dist) ([cdf (CDF In)] + [inv-cdf (Inverse-CDF Out)] + [min Out] + [max Out] + [median (Promise Out)])]{ +The parent type of @italic{ordered} @tech{distribution objects}. + +Similarly to @racket[dist], the @racket[In] type parameter is the data type an ordered distribution +accepts as arguments to its @tech{pdf}; it is also the data type its @tech{cdf} accepts. + +Also similarly to @racket[dist], the @racket[Out] type parameter is the data type an ordered +distribution returns as random samples; it is also the data type its @tech{inverse cdf} returns. + +@examples[#:eval untyped-eval + (ordered-dist? (discrete-dist '(a b c))) + (ordered-dist? (normal-dist))] + +The median is stored in an @racket[ordered-dist] to allow interval probabilities to be calculated +as accurately as possible. For example, for @racket[d = (normal-dist)], whose median is @racket[0.0], +@racket[(real-dist-prob d -2.0 -1.0)] is calculated using lower-tail probabilities, and +@racket[(real-dist-prob d 1.0 2.0)] is calculated using upper-tail probabilities. +} + +@defidform[Real-Dist]{ +The parent type of real-valued distributions, such as any distribution returned by +@racket[normal-dist]. Equivalent to the type @racket[(ordered-dist Real Flonum)]. +} + +@defproc[(real-dist-prob [d Real-Dist] [a Real] [b Real] [log? Any #f] [1-p? Any #f]) Flonum]{ +Computes the probability of the half-open interval (@racket[a], @racket[b]]. (If @racket[b < a], +the two endpoints are swapped first.) The @racket[log?] and @racket[1-p?] arguments determine how the +return value should be interpreted, just as the arguments to a function of type +@racket[(CDF Real Flonum)]. +} + +@deftogether[(@defproc[(dist-cdf [d (ordered-dist In Out)]) (CDF In)] + @defproc[(dist-inv-cdf [d (ordered-dist In Out)]) (Inverse-CDF Out)] + @defproc[(dist-min [d (ordered-dist In Out)]) Out] + @defproc[(dist-max [d (ordered-dist In Out)]) Out] + @defproc[(dist-median [d (ordered-dist In Out)]) Out])]{ +The first four are synonyms for @racket[ordered-dist-cdf], @racket[ordered-dist-inv-cdf], +@racket[ordered-dist-min] and @racket[ordered-dist-max]. The last is equivalent to +@racket[(force (ordered-dist-median d))]. + +@examples[#:eval untyped-eval + (dist-min (normal-dist)) + (dist-min (geometric-dist 1/3)) + (dist-max (uniform-dist 1 2)) + (dist-median (gamma-dist 4 2))] +} + +@section{Discrete Distributions} + +@subsection{Discrete Distribution} + +@subsection{Categorical Distribution} + +@section{Integer Distributions} + +@subsection{Bernoulli Distribution} + +@subsection{Binomial Distribution} + +@subsection{Geometric Distribution} + +@subsection{Poisson Distribution} + +@section{Real Distributions} + +@subsection{Beta Distribution} + +@margin-note{Wikipedia: + @hyperlink["http://wikipedia.org/wiki/Beta_distribution"]{Beta Distribution}.} +@deftogether[(@defidform[Beta-Dist] + @defproc[(beta-dist [alpha Real] [beta Real]) Beta-Dist] + @defproc[(beta-dist-alpha [d Beta-Dist]) Flonum] + @defproc[(beta-dist-beta [d Beta-Dist]) Flonum])]{ +Represents the beta distribution family parameterized by two shape parameters, or pseudocounts. + +@examples[#:eval untyped-eval + (plot (for/list ([α (in-list '(1 2 3 1/2))] + [β (in-list '(1 3 1 1/2))] + [i (in-naturals)]) + (function (dist-pdf (beta-dist α β)) + #:color i #:label (format "Beta(~a,~a)" α β))) + #:x-min 0 #:x-max 1 #:y-max 4 #:y-label "density") + + (plot (for/list ([α (in-list '(1 2 3 1/2))] + [β (in-list '(1 3 1 1/2))] + [i (in-naturals)]) + (function (dist-cdf (beta-dist α β)) + #:color i #:label (format "Beta(~a,~a)" α β))) + #:x-min 0 #:x-max 1 #:y-label "probability")] +} + +@subsection{Cauchy Distribution} + +@margin-note{Wikipedia: + @hyperlink["http://wikipedia.org/wiki/Cauchy_distribution"]{Cauchy Distribution}.} +@deftogether[(@defidform[Cauchy-Dist] + @defproc[(cauchy-dist [mode Real 0] [scale Real 1]) Cauchy-Dist] + @defproc[(cauchy-dist-mode [d Cauchy-Dist]) Flonum] + @defproc[(cauchy-dist-scale [d Cauchy-Dist]) Flonum])]{ +Represents the Cauchy distribution family parameterized by mode and scale. + +@examples[#:eval untyped-eval + (plot (for/list ([m (in-list '(0 -1 0 2))] + [s (in-list '(1 1/2 2.25 0.7))] + [i (in-naturals)]) + (function (dist-pdf (cauchy-dist m s)) + #:color i #:label (format "Cauchy(~a,~a)" m s))) + #:x-min -8 #:x-max 8 #:y-label "density" + #:legend-anchor 'top-right) + + (plot (for/list ([m (in-list '(0 -1 0 2))] + [s (in-list '(1 1/2 2.25 0.7))] + [i (in-naturals)]) + (function (dist-cdf (cauchy-dist m s)) + #:color i #:label (format "Cauchy(~a,~a)" m s))) + #:x-min -8 #:x-max 8 #:y-label "probability")] +} + +@subsection{Delta Distribution} + +@margin-note{Wikipedia: + @hyperlink["http://wikipedia.org/wiki/Dirac_delta_function"]{Dirac Delta Function}.} +@deftogether[(@defidform[Delta-Dist] + @defproc[(delta-dist [mean Real 0]) Delta-Dist] + @defproc[(delta-dist-mean [d Delta-Dist]) Flonum])]{ +Represents the family of distributions whose densities are Dirac delta functions. + +@examples[#:eval untyped-eval + ((dist-pdf (delta-dist)) 0) + ((dist-pdf (delta-dist)) 1) + (plot (for/list ([μ (in-list '(-1 0 1))] + [i (in-naturals)]) + (function (dist-cdf (delta-dist μ)) + #:color i #:style i #:label (format "δ(~a)" μ))) + #:x-min -2 #:x-max 2 #:y-label "probability")] +When a distribution with a scale parameter has scale zero, it behaves like a delta distribution: +@interaction[#:eval untyped-eval + ((dist-pdf (normal-dist 0 0)) 0) + ((dist-pdf (normal-dist 0 0)) 1) + (plot (function (dist-cdf (normal-dist 0 0)) -1e-300 1e-300))] +} + +@subsection{Exponential Distribution} + +@margin-note{Wikipedia: + @hyperlink["http://wikipedia.org/wiki/Exponential_distribution"]{Exponential + Distribution}.} +@deftogether[(@defidform[Exponential-Dist] + @defproc[(exponential-dist [mean Real 1]) Exponential-Dist] + @defproc[(exponential-dist-mean [d Exponential-Dist]) Flonum])]{ +Represents the exponential distribution family parameterized by mean, or scale. + +@bold{Warning:} The exponential distribution family is often parameterized by @italic{rate}, which +is the reciprocal of mean or scale. If you have rates, construct exponential distributions using +@racketblock[(exponential-dist (/ 1.0 rate))] + +@examples[#:eval untyped-eval + (plot (for/list ([μ (in-list '(2/3 1 2))] + [i (in-naturals)]) + (function (dist-pdf (exponential-dist μ)) + #:color i #:label (format "Exponential(~a)" μ))) + #:x-min 0 #:x-max 5 #:y-label "density" + #:legend-anchor 'top-right) + + (plot (for/list ([μ (in-list '(2/3 1 2))] + [i (in-naturals)]) + (function (dist-cdf (exponential-dist μ)) + #:color i #:label (format "Exponential(~a)" μ))) + #:x-min 0 #:x-max 5 #:y-label "probability" + #:legend-anchor 'bottom-right)] +} + +@subsection{Gamma Distribution} + +@margin-note{Wikipedia: + @hyperlink["http://wikipedia.org/wiki/Gamma_distribution"]{Gamma Distribution}.} +@deftogether[(@defidform[Gamma-Dist] + @defproc[(gamma-dist [shape Real 1] [scale Real 1]) Gamma-Dist] + @defproc[(gamma-dist-shape [d Gamma-Dist]) Flonum] + @defproc[(gamma-dist-scale [d Gamma-Dist]) Flonum])]{ +Represents the gamma distribution family parameterized by shape and scale. + +@bold{Warning:} The gamma distribution family is often parameterized by shape and @italic{rate}, +which is the reciprocal of scale. If you have rates, construct gamma distributions using +@racketblock[(gamma-dist shape (/ 1.0 rate))] + +@examples[#:eval untyped-eval + (plot (for/list ([k (in-list '(1 2 3 9))] + [s (in-list '(2 2 3 1/2))] + [i (in-naturals)]) + (function (dist-pdf (gamma-dist k s)) + #:color i #:label (format "Gamma(~a,~a)" k s))) + #:x-min 0 #:x-max 15 #:y-label "density" + #:legend-anchor 'top-right) + + (plot (for/list ([k (in-list '(1 2 3 9))] + [s (in-list '(2 2 3 1/2))] + [i (in-naturals)]) + (function (dist-cdf (gamma-dist k s)) + #:color i #:label (format "Gamma(~a,~a)" k s))) + #:x-min 0 #:x-max 15 #:y-label "probability" + #:legend-anchor 'bottom-right)] +} + +@subsection{Logistic Distribution} + +@margin-note{Wikipedia: + @hyperlink["http://wikipedia.org/wiki/Logistic_distribution"]{Logistic Distribution}.} +@deftogether[(@defidform[Logistic-Dist] + @defproc[(logistic-dist [mean Real 0] [scale Real 1]) Logistic-Dist] + @defproc[(logistic-dist-mean [d Logistic-Dist]) Flonum] + @defproc[(logistic-dist-scale [d Logistic-Dist]) Flonum])]{ +Represents the logistic distribution family parameterized by mean (also called ``location'') +and scale. In this parameterization, the variance is @racket[(* 1/3 (sqr (* pi scale)))]. + +@examples[#:eval untyped-eval + (plot (for/list ([μ (in-list '(0 -1 0 2))] + [s (in-list '(1 1/2 2.25 0.7))] + [i (in-naturals)]) + (function (dist-pdf (logistic-dist μ s)) + #:color i #:label (format "Logistic(~a,~a)" μ s))) + #:x-min -8 #:x-max 8 #:y-label "density" + #:legend-anchor 'top-right) + + (plot (for/list ([μ (in-list '(0 -1 0 2))] + [s (in-list '(1 1/2 2.25 0.7))] + [i (in-naturals)]) + (function (dist-cdf (logistic-dist μ s)) + #:color i #:label (format "Logistic(~a,~a)" μ s))) + #:x-min -8 #:x-max 8 #:y-label "probability")] +} + +@subsection{Normal Distribution} + +@margin-note{Wikipedia: + @hyperlink["http://wikipedia.org/wiki/Normal_distribution"]{Normal Distribution}.} +@deftogether[(@defidform[Normal-Dist] + @defproc[(normal-dist [mean Real 0] [stddev Real 1]) Normal-Dist] + @defproc[(normal-dist-mean [d Normal-Dist]) Flonum] + @defproc[(normal-dist-stddev [d Normal-Dist]) Flonum])]{ +Represents the normal distribution family parameterized by mean and standard deviation. + +@bold{Warning:} The normal distribution family is often parameterized by mean and @italic{variance}, +which is the square of standard deviation. If you have variances, construct normal distributions +using +@racketblock[(normal-dist mean (sqrt var))] + +@examples[#:eval untyped-eval + (plot (for/list ([μ (in-list '(0 -1 0 2))] + [σ (in-list '(1 1/2 2.25 0.7))] + [i (in-naturals)]) + (function (dist-pdf (normal-dist μ σ)) + #:color i #:label (format "N(~a,~a)" μ σ))) + #:x-min -5 #:x-max 5 #:y-label "density") + + (plot (for/list ([μ (in-list '(0 -1 0 2))] + [σ (in-list '(1 1/2 2.25 0.7))] + [i (in-naturals)]) + (function (dist-cdf (normal-dist μ σ)) + #:color i #:label (format "N(~a,~a)" μ σ))) + #:x-min -5 #:x-max 5 #:y-label "probability")] +} + +@subsection{Triangular Distribution} + +@margin-note{Wikipedia: + @hyperlink["http://wikipedia.org/wiki/Triangular_distribution"]{Triangular + Distribution}.} +@deftogether[(@defidform[Triangle-Dist] + @defproc[(triangle-dist [min Real 0] [max Real 1] [mode Real (* 0.5 (+ min max))]) + Triangle-Dist] + @defproc[(triangle-dist-min [d Triangle-Dist]) Flonum] + @defproc[(triangle-dist-max [d Triangle-Dist]) Flonum] + @defproc[(triangle-dist-mode [d Triangle-Dist]) Flonum])]{ +Represents the triangular distribution family parameterized by minimum, maximum and mode. + +If @racket[min], @racket[mode] and @racket[max] are not in ascending order, they are sorted +before constructing the distribution object. + +@examples[#:eval untyped-eval + (plot (for/list ([a (in-list '(-3 -1 -2))] + [b (in-list '(0 1 3))] + [m (in-list '(-2 0 2))] + [i (in-naturals)]) + (function (dist-pdf (triangle-dist a b m)) #:color i + #:label (format "Triangle(~a,~a,~a)" a b m))) + #:x-min -3 #:x-max 3 #:y-label "density") + + (plot (for/list ([a (in-list '(-3 -1 -2))] + [b (in-list '(0 1 3))] + [m (in-list '(-2 0 2))] + [i (in-naturals)]) + (function (dist-cdf (triangle-dist a b m)) #:color i + #:label (format "Triangle(~a,~a,~a)" a b m))) + #:x-min -3 #:x-max 3 #:y-label "probability")] +} + +@subsection{Truncated Distribution} + +@subsection{Uniform Distribution} + +@margin-note{ +Wikipedia: +@hyperlink["http://wikipedia.org/wiki/Uniform_distribution_%28continuous%29"]{Uniform Distribution}.} + +@deftogether[(@defidform[Uniform-Dist] + @defproc*[([(uniform-dist) Uniform-Dist] + [(uniform-dist [max Real]) Uniform-Dist] + [(uniform-dist [min Real] [max Real]) Uniform-Dist])] + @defproc[(uniform-dist-min [d Uniform-Dist]) Flonum] + @defproc[(uniform-dist-max [d Uniform-Dist]) Flonum])]{ +Represents the uniform distribution family parameterized by minimum and maximum. + +@racket[(uniform-dist)] is equivalent to @racket[(uniform-dist 0 1)]. +@racket[(uniform-dist max)] is equivalent to @racket[(uniform-dist 0 max)]. If @racket[max < min], +they are swapped before constructing the distribution object. +@;{ +@examples[#:eval untyped-eval + (plot (for/list ([a (in-list '(-3 -1 -2))] + [b (in-list '(0 1 3))] + [i (in-naturals)]) + (function (dist-pdf (uniform-dist a b)) #:color i + #:label (format "Uniform(~a,~a)" a b))) + #:x-min -3 #:x-max 3 #:y-label "density") + + (plot (for/list ([a (in-list '(-3 -1 -2))] + [b (in-list '(0 1 3))] + [i (in-naturals)]) + (function (dist-cdf (uniform-dist a b)) #:color i + #:label (format "Uniform(~a,~a)" a b))) + #:x-min -3 #:x-max 3 #:y-label "probability")] +} +} +@(close-eval untyped-eval) diff --git a/collects/math/scribblings/math-number-theory.scrbl b/collects/math/scribblings/math-number-theory.scrbl index b14d2e1ab4..09ee2245aa 100644 --- a/collects/math/scribblings/math-number-theory.scrbl +++ b/collects/math/scribblings/math-number-theory.scrbl @@ -180,7 +180,7 @@ modulus. For example, the code corresponds with the mathematical statement @italic{a}@superscript{@italic{b}} = @italic{c} (mod @italic{n}). -The current modulus is stored in a @tech{parameter} that, for performance reasons, can only be +The current modulus is stored in a parameter that, for performance reasons, can only be set using @racket[with-modulus]. (The basic modular operators cache parameter reads, and this restriction guarantees that the cached values are current.) diff --git a/collects/math/scribblings/math.scrbl b/collects/math/scribblings/math.scrbl index ada16029c3..e2a7f01eac 100644 --- a/collects/math/scribblings/math.scrbl +++ b/collects/math/scribblings/math.scrbl @@ -39,3 +39,5 @@ be used in untyped Racket. Exceptions and performance warnings are in @bold{bold @include-section["math-special-functions.scrbl"] @include-section["math-number-theory.scrbl"] @include-section["math-bigfloat.scrbl"] +@include-section["math-array.scrbl"] +@include-section["math-distributions.scrbl"] diff --git a/collects/math/scribblings/utils.rkt b/collects/math/scribblings/utils.rkt index c3fba841f9..a4f7f77c17 100644 --- a/collects/math/scribblings/utils.rkt +++ b/collects/math/scribblings/utils.rkt @@ -5,6 +5,7 @@ (provide author-neil author-jens-axel + make-plain-math-eval make-math-eval make-untyped-math-eval) @@ -14,10 +15,14 @@ (define (author-jens-axel) @author{@(author+email "Jens Axel Søgaard" "jensaxel@soegaard.net")}) -(define (make-math-eval) +(define (make-plain-math-eval) (define eval (make-base-eval)) (eval '(require typed/racket/base)) (eval '(require math)) + eval) + +(define (make-math-eval) + (define eval (make-plain-math-eval)) (eval '(require math/scribblings/rename-defines)) (λ (v) (cond [(syntax? v) (eval #`(rename-defines #,v))]