Moved flvector functions into math/flonum
Sped up normal distribution sampling procedure (2x for large samples)
This commit is contained in:
Neil Toronto
parent
cd002d5830
commit
6009eed8d2
|
|
@ -15,7 +15,8 @@
|
|||
"private/flonum/flonum-syntax.rkt"
|
||||
"private/flonum/expansion/expansion-base.rkt"
|
||||
"private/flonum/expansion/expansion-exp.rkt"
|
||||
"private/flonum/expansion/expansion-log.rkt")
|
||||
"private/flonum/expansion/expansion-log.rkt"
|
||||
"private/flonum/flvector.rkt")
|
||||
|
||||
(provide (all-from-out
|
||||
"private/flonum/flonum-bits.rkt"
|
||||
|
|
@ -32,7 +33,8 @@
|
|||
"private/flonum/flonum-syntax.rkt"
|
||||
"private/flonum/expansion/expansion-base.rkt"
|
||||
"private/flonum/expansion/expansion-exp.rkt"
|
||||
"private/flonum/expansion/expansion-log.rkt")
|
||||
"private/flonum/expansion/expansion-log.rkt"
|
||||
"private/flonum/flvector.rkt")
|
||||
lg* lg/ lgprod)
|
||||
|
||||
(define-syntax lg* (make-rename-transformer #'fl+))
|
||||
|
|
|
|||
|
|
@ -10,7 +10,6 @@
|
|||
(module defs typed/racket/base
|
||||
|
||||
(require "../../flonum.rkt"
|
||||
"../vector/flvector.rkt"
|
||||
"../unsafe.rkt"
|
||||
"array-struct.rkt"
|
||||
"utils.rkt"
|
||||
|
|
|
|||
|
|
@ -2,7 +2,7 @@
|
|||
|
||||
(require racket/flonum
|
||||
(for-syntax racket/base)
|
||||
"../vector/flvector.rkt"
|
||||
"../../flonum.rkt"
|
||||
"array-struct.rkt"
|
||||
"array-broadcast.rkt"
|
||||
"array-pointwise.rkt"
|
||||
|
|
|
|||
|
|
@ -10,7 +10,6 @@
|
|||
|
||||
(require "../../flonum.rkt"
|
||||
"../unsafe.rkt"
|
||||
"../vector/flvector.rkt"
|
||||
"array-struct.rkt"
|
||||
"utils.rkt"
|
||||
"for-each.rkt")
|
||||
|
|
|
|||
|
|
@ -3,7 +3,6 @@
|
|||
(require racket/performance-hint
|
||||
racket/promise
|
||||
"../../flonum.rkt"
|
||||
"../../vector.rkt"
|
||||
"../unsafe.rkt"
|
||||
"dist-struct.rkt"
|
||||
"utils.rkt")
|
||||
|
|
|
|||
|
|
@ -4,7 +4,6 @@
|
|||
racket/performance-hint
|
||||
racket/promise
|
||||
"../../flonum.rkt"
|
||||
"../../vector.rkt"
|
||||
"../unsafe.rkt"
|
||||
"../functions/beta.rkt"
|
||||
"../functions/incomplete-beta.rkt"
|
||||
|
|
|
|||
|
|
@ -3,7 +3,6 @@
|
|||
(require racket/performance-hint
|
||||
racket/promise
|
||||
"../../flonum.rkt"
|
||||
"../../vector.rkt"
|
||||
"../unsafe.rkt"
|
||||
"../functions/incomplete-beta.rkt"
|
||||
(prefix-in impl: "impl/binomial-pdf.rkt")
|
||||
|
|
|
|||
|
|
@ -4,7 +4,6 @@
|
|||
racket/promise
|
||||
"../../flonum.rkt"
|
||||
"../../base.rkt"
|
||||
"../../vector.rkt"
|
||||
"../unsafe.rkt"
|
||||
"dist-struct.rkt"
|
||||
"utils.rkt")
|
||||
|
|
|
|||
|
|
@ -3,7 +3,6 @@
|
|||
(require racket/performance-hint
|
||||
racket/promise
|
||||
"../../flonum.rkt"
|
||||
"../../vector.rkt"
|
||||
"../unsafe.rkt"
|
||||
"dist-struct.rkt"
|
||||
"utils.rkt")
|
||||
|
|
|
|||
|
|
@ -3,7 +3,6 @@
|
|||
(require racket/performance-hint
|
||||
racket/promise
|
||||
"../../flonum.rkt"
|
||||
"../../vector.rkt"
|
||||
"../unsafe.rkt"
|
||||
"../functions/incomplete-gamma.rkt"
|
||||
"impl/gamma-pdf.rkt"
|
||||
|
|
|
|||
|
|
@ -3,7 +3,6 @@
|
|||
(require racket/performance-hint
|
||||
racket/promise
|
||||
"../../flonum.rkt"
|
||||
"../../vector.rkt"
|
||||
"../unsafe.rkt"
|
||||
"dist-struct.rkt"
|
||||
"utils.rkt")
|
||||
|
|
|
|||
|
|
@ -6,7 +6,6 @@ Wolfgang Hormann. The Generation of Binomial Random Variates.
|
|||
|
||||
(require "../../../base.rkt"
|
||||
"../../../flonum.rkt"
|
||||
"../../../vector.rkt"
|
||||
"../../unsafe.rkt"
|
||||
"normal-random.rkt")
|
||||
|
||||
|
|
|
|||
|
|
@ -19,7 +19,6 @@ For others: sum of Gamma and Exponential variables, Normal approximation.
|
|||
|#
|
||||
|
||||
(require "../../../flonum.rkt"
|
||||
"../../../vector.rkt"
|
||||
"../../unsafe.rkt"
|
||||
"normal-random.rkt")
|
||||
|
||||
|
|
|
|||
|
|
@ -1,30 +1,68 @@
|
|||
#lang typed/racket/base
|
||||
|
||||
(require "../../../flonum.rkt"
|
||||
;; Return random samples from a normal distribution using the Box-Muller transform
|
||||
|
||||
(require racket/fixnum
|
||||
"../../../flonum.rkt"
|
||||
"../../../base.rkt"
|
||||
"../../../vector.rkt")
|
||||
"../../unsafe.rkt")
|
||||
|
||||
(provide flnormal-sample)
|
||||
|
||||
(: box-muller-transform (Float Float -> Float))
|
||||
(define (box-muller-transform x y)
|
||||
(cond [(and (fl= x 0.0) (fl= y 0.0)) 0.0]
|
||||
[else (fl* (flsqrt (fl* -2.0 (fllog x)))
|
||||
(flsin (fl* (fl* 2.0 pi) y)))]))
|
||||
;; Leaving these in, in case we discover in the future that it's actually important for them
|
||||
;; to be accurate
|
||||
#|
|
||||
(: flsin2pix (Flonum -> Flonum))
|
||||
;; Computes sin(2*pi*x) accurately in the range [0,1]
|
||||
(define (flsin2pix x)
|
||||
(let*-values ([(x s) (if (x . fl> . 0.5) (values (fl- x 0.5) -1.0) (values x 1.0))]
|
||||
[(x) (if (x . fl> . 0.25) (fl- 0.5 x) x)])
|
||||
(fl* s (flsin (fl* (fl* 2.0 pi) x)))))
|
||||
|
||||
(: flcos2pix (Flonum -> Flonum))
|
||||
;; Computes cos(2*pi*x) accurately in the range [0,1]
|
||||
(define (flcos2pix x)
|
||||
(let*-values ([(x) (if (x . fl> . 0.5) (fl- 1.0 x) x)]
|
||||
[(x s) (if (x . fl> . 0.25) (values (fl- 0.5 x) 1.0) (values x -1.0))])
|
||||
(fl* s (flsin (fl* (fl* 2.0 pi) (fl- x 0.25))))))
|
||||
|#
|
||||
|
||||
(: nonzero-random (-> Flonum))
|
||||
(define (nonzero-random)
|
||||
(let ([u (random)])
|
||||
(if (fl= u 0.0) (nonzero-random) u)))
|
||||
|
||||
(: 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)
|
||||
;; 3. 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 (flnormal-sample c s n)
|
||||
(cond [(n . < . 0) (raise-argument-error 'flnormal-sample "Natural" 2 c s n)]
|
||||
(cond [(not (index? n)) (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)))))]))
|
||||
|
||||
(define xs (make-flvector n))
|
||||
(cond
|
||||
[(fx= n 0) xs]
|
||||
[else
|
||||
(define n-1 (fx- n 1))
|
||||
(let loop ([#{i : Nonnegative-Fixnum} 0])
|
||||
(cond [(i . fx< . n-1)
|
||||
(define u1 (nonzero-random))
|
||||
(define u2 (random))
|
||||
(define t (flsqrt (fl* -2.0 (fllog u1))))
|
||||
(define z (fl* (fl* 2.0 pi) u2))
|
||||
(define x (fl* t (flcos z)))
|
||||
(define y (fl* t (flsin z)))
|
||||
(unsafe-flvector-set! xs i (fl+ c (fl* s x)))
|
||||
(unsafe-flvector-set! xs (fx+ i 1) (fl+ c (fl* s y)))
|
||||
(loop (fx+ i 2))]
|
||||
[(i . fx= . n-1)
|
||||
(define u1 (nonzero-random))
|
||||
(define u2 (random))
|
||||
(define x (fl* (flsqrt (fl* -2.0 (fllog u1)))
|
||||
(flsin (fl* (fl* 2.0 pi) u2))))
|
||||
(unsafe-flvector-set! xs i (fl+ c (fl* s x)))
|
||||
xs]
|
||||
[else
|
||||
xs]))])]))
|
||||
|
|
|
|||
|
|
@ -3,7 +3,6 @@
|
|||
(require racket/fixnum
|
||||
"../../../flonum.rkt"
|
||||
"../../../base.rkt"
|
||||
"../../../vector.rkt"
|
||||
"../../functions/log-gamma.rkt")
|
||||
|
||||
(provide flpoisson-sample)
|
||||
|
|
|
|||
|
|
@ -3,7 +3,6 @@
|
|||
(require racket/performance-hint
|
||||
racket/promise
|
||||
"../../flonum.rkt"
|
||||
"../../vector.rkt"
|
||||
"../unsafe.rkt"
|
||||
"dist-struct.rkt"
|
||||
"utils.rkt")
|
||||
|
|
|
|||
|
|
@ -3,7 +3,6 @@
|
|||
(require racket/performance-hint
|
||||
racket/promise
|
||||
"../../flonum.rkt"
|
||||
"../../vector.rkt"
|
||||
"../unsafe.rkt"
|
||||
"impl/normal-pdf.rkt"
|
||||
"impl/normal-cdf.rkt"
|
||||
|
|
|
|||
|
|
@ -3,7 +3,6 @@
|
|||
(require racket/performance-hint
|
||||
racket/promise
|
||||
"../../flonum.rkt"
|
||||
"../../vector.rkt"
|
||||
"../unsafe.rkt"
|
||||
"../inline-sort.rkt"
|
||||
"dist-struct.rkt"
|
||||
|
|
|
|||
|
|
@ -3,7 +3,6 @@
|
|||
(require racket/performance-hint
|
||||
racket/promise
|
||||
"../../flonum.rkt"
|
||||
"../../vector.rkt"
|
||||
"../unsafe.rkt"
|
||||
"dist-struct.rkt"
|
||||
"utils.rkt")
|
||||
|
|
|
|||
|
|
@ -3,7 +3,7 @@
|
|||
(require racket/fixnum
|
||||
racket/string
|
||||
(for-syntax racket/base syntax/parse)
|
||||
"../../flonum.rkt"
|
||||
"flonum-functions.rkt"
|
||||
"../unsafe.rkt"
|
||||
"flvector-syntax.rkt")
|
||||
|
||||
|
|
@ -12,9 +12,6 @@
|
|||
;; Construction
|
||||
unsafe-flvector-copy!
|
||||
flvector-copy!
|
||||
;; Loops
|
||||
for/flvector:
|
||||
for*/flvector:
|
||||
;; Conversion
|
||||
list->flvector
|
||||
flvector->list
|
||||
|
|
@ -90,47 +87,6 @@
|
|||
[else
|
||||
(unsafe-flvector-copy! dest dest-start src src-start src-end)])]))
|
||||
|
||||
;; ===================================================================================================
|
||||
;; Loops
|
||||
|
||||
(define-syntax (base-for/flvector: stx)
|
||||
(syntax-parse stx
|
||||
[(_ for: #:length n-expr:expr (clauses ...) body ...+)
|
||||
(syntax/loc stx
|
||||
(let: ([n : Integer n-expr])
|
||||
(cond [(n . > . 0)
|
||||
(define xs (make-flvector n))
|
||||
(define: i : Nonnegative-Fixnum 0)
|
||||
(let/ec: break : Void
|
||||
(for: (clauses ...)
|
||||
(unsafe-flvector-set! xs i (let () body ...))
|
||||
(set! i (unsafe-fx+ i 1))
|
||||
(when (i . unsafe-fx>= . n) (break (void)))))
|
||||
xs]
|
||||
[else (flvector)])))]
|
||||
[(_ for: (clauses ...) body ...+)
|
||||
(syntax/loc stx
|
||||
(let ()
|
||||
(define n 4)
|
||||
(define xs (make-flvector 4))
|
||||
(define i 0)
|
||||
(for: (clauses ...)
|
||||
(let: ([x : Float (let () body ...)])
|
||||
(cond [(unsafe-fx= i n) (define new-n (unsafe-fx* 2 n))
|
||||
(define new-xs (make-flvector new-n x))
|
||||
(unsafe-flvector-copy! new-xs 0 xs 0 n)
|
||||
(set! n new-n)
|
||||
(set! xs new-xs)]
|
||||
[else (unsafe-flvector-set! xs i x)]))
|
||||
(set! i (unsafe-fx+ i 1)))
|
||||
(flvector-copy xs 0 i)))]))
|
||||
|
||||
(define-syntax-rule (for/flvector: e ...)
|
||||
(base-for/flvector: for: e ...))
|
||||
|
||||
(define-syntax-rule (for*/flvector: e ...)
|
||||
(base-for/flvector: for*: e ...))
|
||||
|
||||
;; ===================================================================================================
|
||||
;; Conversion
|
||||
|
||||
|
|
@ -25,7 +25,6 @@ Approximations:
|
|||
"../../flonum.rkt"
|
||||
"../../base.rkt"
|
||||
"../number-theory/factorial.rkt"
|
||||
"../vector/flvector.rkt"
|
||||
"lanczos.rkt")
|
||||
|
||||
(provide gamma flgamma)
|
||||
|
|
|
|||
|
|
@ -2,7 +2,6 @@
|
|||
|
||||
(require "../../../flonum.rkt"
|
||||
"../../../base.rkt"
|
||||
"../../../vector.rkt"
|
||||
"../../unsafe.rkt"
|
||||
"../../distributions/impl/normal-cdf.rkt"
|
||||
"../stirling-error.rkt"
|
||||
|
|
|
|||
|
|
@ -2,8 +2,7 @@
|
|||
|
||||
(require racket/fixnum
|
||||
"../../flonum.rkt"
|
||||
"../vector/vector.rkt"
|
||||
"../vector/flvector.rkt")
|
||||
"../vector/vector.rkt")
|
||||
|
||||
(provide fltan-diff/y flcot-diff/y)
|
||||
|
||||
|
|
|
|||
|
|
@ -5,10 +5,9 @@
|
|||
racket/list
|
||||
racket/future
|
||||
"../../base.rkt"
|
||||
"../../flonum.rkt"
|
||||
"../../parameters.rkt"
|
||||
"../unsafe.rkt"
|
||||
"flvector.rkt"
|
||||
"vector.rkt")
|
||||
"../unsafe.rkt")
|
||||
|
||||
(provide vector-fft flvector-fft!
|
||||
vector-inverse-fft flvector-inverse-fft!)
|
||||
|
|
|
|||
|
|
@ -1,10 +1,8 @@
|
|||
#lang typed/racket/base
|
||||
|
||||
(require "private/vector/flvector.rkt"
|
||||
"private/vector/vector.rkt"
|
||||
(require "private/vector/vector.rkt"
|
||||
"private/vector/vector-fft.rkt")
|
||||
|
||||
(provide (all-from-out
|
||||
"private/vector/flvector.rkt"
|
||||
"private/vector/vector.rkt"
|
||||
"private/vector/vector-fft.rkt"))
|
||||
|
|
|
|||
Loading…
Reference in New Issue
Block a user