| PasteRack.org | ||
| Paste # 29444 | ||
| 2019-07-23 02:20:39 | ||
Fork as a new paste. | ||
Paste viewed 452 times. | ||
Tweet | ||
Embed: | ||
#lang typed/racket/base
;; These functions accept non-integral `k' because they're used for the gamma distribution as well
;; as the Poisson distribution
(require math/flonum
math/base
math/private/functions/gamma
math/private/functions/log-gamma
math/private/functions/stirling-error
math/private/distributions/normal-dist)
(provide flpoisson-pdf flpoisson-log-pdf)
(: flpoisson-pdf-limits (Flonum Flonum Any -> Flonum))
(define (flpoisson-pdf-limits l k log?)
(cond [(l . fl<= . 0.0)
(cond [(fl= l 0.0) (cond [(fl= k 0.0) (if log? 0.0 1.0)]
[else (if log? -inf.0 0.0)])]
[else +nan.0])]
[(k . fl<= . 0.0)
(cond [(fl= k 0.0) (if log? (- l) (flexp (- l)))]
[else (if log? -inf.0 0.0)])]
[else +nan.0]))
(: flpoisson-in-bounds? (Flonum Flonum -> Boolean))
(define (flpoisson-in-bounds? l k)
(and (l . fl> . 0.0) (k . fl> . 0.0) (l . fl< . +inf.0) (k . fl< . +inf.0)))
(define log-c (fllog (flsqrt +max-subnormal.0)))
(: get-div (Flonum Flonum -> Flonum))
;; Returns a divisor d that ensures l^(k/d)/k^(k/d)*exp((k-l)/d) won't underflow; assumes l and
;; k are positive
(define (get-div l k)
(define d0 (fl/ (fl* k (fllog (fl/ l k))) log-c))
(define d1 (fl/ (fl- k l) log-c))
(define d (max d0 d1))
(flmax 1.0 (flexpt 2.0 (flfloor (fl/ (fllog d) (fllog 2.0))))))
(: flpoisson-pdf (Flonum Flonum -> Flonum))
;; Maximum relative error depends on l and k:
;; l,k < 4000 5 ulps
;; l,k < 1e4 18 ulps
;; l,k < 1e5 65 ulps
;; l,k < 1e6 170 ulps
;; l,k < 1e7 600 ulps
;; This is due to the argument reduction used in these areas to keep exponentials from underflowing
;; (search for `get-div' below)
(define (flpoisson-pdf l k)
(cond [(not (flpoisson-in-bounds? l k))
(flpoisson-pdf-limits l k #f)]
;; When l is large enough, there's no density > +min.0/2 within 48 stddevs
[(and (l . fl> . 500.0) ((flabs (fl- k l)) . fl> . (fl* 48.0 (flsqrt l)))) 0.0]
;; Use the direct definition when it won't over/underflow
[(and (l . fl< . 700.0) (k . fl> . 1.0) (k . fl< . 100.0))
(fl/ (fl/ (fl* (flexpt l k) (flexp (- l))) (flgamma k)) k)]
[(k . fl<= . 1.0)
;; Using Stirling error works great here
(define-values (l-k l-k-lo) (fast-fl-/error l k))
(fl/ (* (flexp (- l-k))
(flexpt l k)
(flexpt k (- k))
(flexp (- l-k-lo))
(/ 1.0 (flsqrt (fl* 2.0 pi))))
(fl* (flexp-stirling k) (flsqrt k)))]
;; When l > 1e21, the normal approximation has *lower* relative error than the algorithm
;; below. It *does not* have *low* relative error! It's about 2e9...
;; The normal approximation is perfect starting around 4.2e34, which is incidentally the
;; same time that 48 standard deviations from `l' becomes its floating-point neighbor.
[(l . fl>= . 1e21) (flnormal-pdf l (flsqrt l) k #f)]
[else
(define-values (l-k l-k-lo) (fast-fl-/error l k))
(define-values (l/k l/k-lo) (fast-fl//error l k))
(define div (get-div l k))
(define y
(fl* (flexpt (fl* (flexpt+ l/k l/k-lo (fl/ k div))
(flexp (fl/ (- l-k) div)))
div)
(flexp (- l-k-lo))))
(fl/ (fl* y (fl/ 1.0 (flsqrt (fl* 2.0 pi))))
(fl* (flexp-stirling k) (flsqrt k)))]))
(: fllog-gamma1p-taylor-0 (Flonum -> Flonum))
;; Good for k < 0.005
(define fllog-gamma1p-taylor-0
(make-flpolyfun
(0.0
-5.772156649015328606065120900824024310422e-1
8.224670334241132182362075833230125946095e-1
-4.00685634386531428466579387170483330255e-1
2.705808084277845478790009241352919756937e-1
-2.073855510286739852662730972914068336114e-1
1.69557176997408189952419654965153421317e-1)))
(: flpoisson-log-pdf (Flonum Flonum -> Flonum))
(define (flpoisson-log-pdf l k)
(cond [(not (flpoisson-in-bounds? l k))
(flpoisson-pdf-limits l k #t)]
[(k . fl<= . 0.005)
;; Don't subtract log(k)+log(gamma(k)) (which suffers from cancellation here), subtract
;; terms from its Taylor series
(fl- (fl+ (fl* k (fllog l)) (- l))
(fllog-gamma1p-taylor-0 k))]
[(and (k . fl< . 1.0) (l . fl< . 1.0))
(cond [(or (k . fl<= . 0.96) (l . fl> . 1e-300))
(fl- (fllog (fl* (flexpt l k) (flexp (- l))))
(fl+ (fllog k) (fllog-gamma k)))]
[else
(- (fl* k (fllog l)) (fllog k) (fllog-gamma k) l)])]
[(or (k . fl< . 2.0) (l . fl< . 2.0))
(- (fl* k (fllog l)) (fllog k) (fllog-gamma k) l)]
[(and (k . fl> . 2.0) (l . fl> . 2.0)
((fl- k l) . fl< . (fl* 40.0 (flsqrt l)))
((fl- l k) . fl< . (fl* 29.0 (flsqrt l))))
;; Error <= 1 ulp when flpoisson-pdf error is just a few ulps
(fllog (flpoisson-pdf l k))]
[(and (or (l . fl> . 1e18) (k . fl> . 1e18))
(or ((fl- k l) . fl> . (fl* 0.5 k))
((fl- l k) . fl> . (fl* 0.5 l))))
;; Asymptotic expansion had by replacing log(gamma(k)) with k*log(k)-k
(- (fl* k (fl+ 1.0 (fllog (fl/ l k)))) (fllog k) l)]
[else
;; Log version of the above; relative error here also climbs with distance from 0
(define-values (l-k l-k-lo) (fast-fl-/error l k))
(define-values (l/k l/k-lo) (fast-fl//error l k))
(define div (get-div l k))
(- (fl* div (fllog (fl* (flexpt+ l/k l/k-lo (fl/ k div))
(flexp (fl/ (- l-k) div)))))
(flstirling k)
(fllog (* (flsqrt k) (flsqrt (fl* 2.0 pi)) (flexp l-k-lo))))]))
(: flpoisson-log-pdf+ (Flonum Flonum -> Flonum))
(define (flpoisson-log-pdf+ l k)
(cond [(not (flpoisson-in-bounds? l k))
(flpoisson-pdf-limits l k #t)]
[(k . fl<= . 0.005)
;; Don't subtract log(k)+log(gamma(k)) (which suffers from cancellation here), subtract
;; terms from its Taylor series
(fl- (fl+ (fl* k (fllog l)) (- l))
(fllog-gamma1p-taylor-0 k))]
[(and (k . fl< . 1.0) (l . fl< . 1.0))
(cond [(or (k . fl<= . 0.96) (l . fl> . 1e-300))
(fl- (fllog (fl* (flexpt l k) (flexp (- l))))
(fl+ (fllog k) (fllog-gamma k)))]
[else
(- (fl* k (fllog l)) (fllog k) (fllog-gamma k) l)])]
[(or (k . fl< . 2.0) (l . fl< . 2.0))
(- (fl* k (fllog l)) (fllog k) (fllog-gamma k) l)]
[(and (k . fl> . 2.0) (l . fl> . 2.0)
((fl- k l) . fl< . (fl* 40.0 (flsqrt l)))
((fl- l k) . fl< . (fl* 38.0 (flsqrt l))))
;; Error <= 1 ulp when flpoisson-pdf error is just a few ulps
(define ans (fllog (flpoisson-pdf l k)))
(cond
[(fl< ans -700.0)
(define bns (flpoisson-log-pdf-logassymp l k))
;(println (list ans bns))
;bail out -> rounding errors to big in this region
bns]
[else ans])]
[(and (or (l . fl> . 1e18) (k . fl> . 1e18))
(or ((fl- k l) . fl> . (fl* 0.5 k))
((fl- l k) . fl> . (fl* 0.5 l))))
;; Asymptotic expansion had by replacing log(gamma(k)) with k*log(k)-k
(- (fl* k (fl+ 1.0 (fllog (fl/ l k)))) (fllog k) l)]
[else
;; Log version of the above; relative error here also climbs with distance from 0
(flpoisson-log-pdf-logassymp l k)]))
(: flpoisson-log-pdf++ (Flonum Flonum -> Flonum))
(define (flpoisson-log-pdf++ l k)
(cond [(not (flpoisson-in-bounds? l k))
(flpoisson-pdf-limits l k #t)]
[(k . fl<= . 0.005)
;; Don't subtract log(k)+log(gamma(k)) (which suffers from cancellation here), subtract
;; terms from its Taylor series
(fl- (fl+ (fl* k (fllog l)) (- l))
(fllog-gamma1p-taylor-0 k))]
[(and (k . fl< . 1.0) (l . fl< . 1.0))
(cond [(or (k . fl<= . 0.96) (l . fl> . 1e-300))
(fl- (fllog (fl* (flexpt l k) (flexp (- l))))
(fl+ (fllog k) (fllog-gamma k)))]
[else
(- (fl* k (fllog l)) (fllog k) (fllog-gamma k) l)])]
[(or (k . fl< . 2.0) (l . fl< . 2.0))
(- (fl* k (fllog l)) (fllog k) (fllog-gamma k) l)]
[(and (k . fl> . 2.0) (l . fl> . 2.0)
(or ((fl- k l) . fl< . (fl* 38.0 (flexpt l 0.499)))
(and (l . fl< . 12500.)
((fl- k l) . fl< . (fl* 40.0 (flexpt l 0.5)))))
(or (l . fl< . 735.)
(and (l . fl> . 2000.0)
((fl- l k) . fl< . (fl- (fl* 38.0 (flexpt l 0.499)) 400.0)))
(k . fl> . (fl+ (fl* (flexpt (fl- l 735.) 1.2) 0.12221088561886605) 2.0))))
;; Error <= 1 ulp when flpoisson-pdf error is just a few ulps
(fllog (flpoisson-pdf l k))]
[(and (or (l . fl> . 1e18) (k . fl> . 1e18))
(or ((fl- k l) . fl> . (fl* 0.5 k))
((fl- l k) . fl> . (fl* 0.5 l))))
;; Asymptotic expansion had by replacing log(gamma(k)) with k*log(k)-k
(- (fl* k (fl+ 1.0 (fllog (fl/ l k)))) (fllog k) l)]
[else
;; Log version of the above; relative error here also climbs with distance from 0
(define-values (l-k l-k-lo) (fast-fl-/error l k))
(define-values (l/k l/k-lo) (fast-fl//error l k))
(define div (get-div l k))
(- (fl* div (fllog (fl* (flexpt+ l/k l/k-lo (fl/ k div))
(flexp (fl/ (- l-k) div)))))
(flstirling k)
(fllog (* (flsqrt k) (flsqrt (fl* 2.0 pi)) (flexp l-k-lo))))]))
(: flpoisson-log-pdf1 (Flonum Flonum -> Flonum))
(define (flpoisson-log-pdf1 l k)
(cond [(not (flpoisson-in-bounds? l k))
(flpoisson-pdf-limits l k #t)]
[else
;; Error <= 1 ulp when flpoisson-pdf error is just a few ulps
(fllog (flpoisson-pdf l k))]))
(: flpoisson-log-pdf2 (Flonum Flonum -> Flonum))
(define (flpoisson-log-pdf2 l k)
(cond [(not (flpoisson-in-bounds? l k))
(flpoisson-pdf-limits l k #t)]
[else
(flpoisson-log-pdf-logassymp l k)]))
(: flpoisson-log-pdf-logassymp (Flonum Flonum -> Flonum))
(define (flpoisson-log-pdf-logassymp l k)
;; Log version of the above; relative error here also climbs with distance from 0
(define-values (l-k l-k-lo) (fast-fl-/error l k))
(define-values (l/k l/k-lo) (fast-fl//error l k))
(define div (get-div l k))
(- (fl* div (fllog (fl* (flexpt+ l/k l/k-lo (fl/ k div))
(flexp (fl/ (- l-k) div)))))
(flstirling k)
(fllog (* (flsqrt k) (flsqrt (fl* 2.0 pi)) (flexp l-k-lo)))))
(module+ test0
(require plot)
(define lowerborder
;goes from k=2 along line l-k=29√l
(let* ([n 29.]
[k0 2.]
[l0 (flexpt (fl/ (fl+ n (flsqrt (fl+ (flexpt n 2.) (fl* 4. k0))))
2.)
2.)])
(λ ([l : Flonum])
(cond
[(< l l0) (list l 2.0)]
[else (list l (fl- l (fl* n (flsqrt l))))]))))
(define highrborder
;goes from l=2 along line k-l=40√l
(let* ([n 40.0]
[l0 2.0]
[k0 (fl+ l0 (fl* n (flsqrt l0)))])
(λ ([l : Flonum])
(cond
[(< l l0) (list 2.0 (fl* k0 (fl/ l 2.0)))]
[else
(list l (fl+ l (fl* n (flsqrt l))))]))))
(define altr0border
;goes from k=2 along line l-k=29√l
(let* ([n 38.]
[k0 2.]
[l0 (flexpt (fl/ (fl+ n (flsqrt (fl+ (flexpt n 2.) (fl* 4. k0))))
2.)
2.)])
(λ ([l : Flonum])
(cond
[(< l l0) (list l 2.0)]
[else (list l (fl- l (fl* n (flsqrt l))))]))))
(define altr1border
;goes from l=2 along line k-l=40√l
(let* ([n 38.0]
[a 0.499]
[b 000.0];300
[c 400.0];500
[d 0.0]
[l0 735.0]
[k0 2.0]
[l1 4050.0];2000
[k1 (fl- l1
(fl- (fl* n (flexpt (fl- l1 b) a))
(fl* c (flexpt l1 d))))]
[e 1.2]
[f (fl/ (fl- k1 k0) (flexpt (fl- l1 l0) e))]
)
(println f)
(λ ([l : Flonum])
(cond
[(< l l0) (list l k0)]
[(< l l1) (list l (fl+ (fl* (flexpt (fl- l l0) e) f) k0))]
[else
(list l (fl- l
(fl- (fl* n (flexpt (fl- l b) a))
(fl* c (flexpt l d)))))]))))
(define altrLborder
;goes from l=2 along line k-l=40√l
(let* ([n 38.0]
[a 0.499]
[b 300.0];300
[c 500.0];500
[d -0.1]
[l0 745.0]
[k0 2.0]
[l1 2000.0];2000
[k1 (fl- l1
(fl- (fl* n (flexpt (fl- l1 b) a))
(fl* c (flexpt l1 d))))]
[e 1.2]
[f (fl/ (fl- k1 k0) (flexpt (fl- l1 l0) e))]
)
(λ ([l : Flonum])
(cond
[(< l l0) (list l k0)]
[(< l l1) (list l (fl+ (fl* (flexpt (fl- l l0) e) f) k0))]
[else
(list l (fl- l
(fl- (fl* n (flexpt (fl- l b) a))
(fl* c (flexpt l d)))))]))))
(define altr2border
;goes from l=2 along line k-l=40√l
(let* ([n 39.0]
[a 0.499]
[l0 2.0]
[k0 (fl+ l0 (fl* n (flexpt l0 a)))])
(λ ([l : Flonum])
(cond
[(< l l0) (list l0 (fl* k0 (fl/ l l0)))]
[else
(list l (fl+ l (fl* n (flexpt l a))))]))))
#|
(plot
(list
(parametric (λ (t) (lowerborder (fl t))) 0 1000)
(parametric (λ (t) (highrborder (fl t))) 0 1800)
(parametric (λ (t) (altr1border (fl t))) 0 1000)
(parametric (λ (t) (altr2border (fl t))) 0 1000)))
(plot
(list
(parametric (λ ([t : Real]) (lowerborder (fl t))) 0 160 #:color 'green)
(parametric (λ ([t : Real]) (highrborder (fl t))) 0 260 #:color 'green)
(contours (λ ([l : Real] [k : Real]) (flpoisson-log-pdf (fl l) (fl k)))
0 1000
0 1000
#:label "p-pdf"
#:levels '(0 -100 -200 -500 -700 -740 -1000 -2000))))
(plot
(list
(parametric (λ ([t : Real]) (lowerborder (fl t))) 0 160 #:color 'green)
(parametric (λ ([t : Real]) (highrborder (fl t))) 0 260 #:color 'green)
(contours (λ ([l : Real] [k : Real]) (flpoisson-log-pdf+ (fl l) (fl k)))
0 1000
0 1000
#:label "p-pdf"
#:levels '(0 -100 -200 -500 -700 -740 -1000 -2000))))
(let ([t0 0]
[t1 100])
(plot
(list
(parametric (λ ([t : Real]) (lowerborder (fl t))) t0 t1 #:color 'green)
;(parametric (λ ([t : Real]) (highrborder (fl t))) t0 t1 #:color 'green)
(contours (λ ([l : Real] [k : Real]) (flpoisson-log-pdf+ (fl l) (fl k)))
(car (lowerborder (fl t0)))(car (lowerborder (fl t1)))
(cadr (lowerborder (fl t0)))(cadr (lowerborder (fl t1)))
#:label "p-pdf"
#:levels '(-720 -730 -740 -750 -760 -770 -780 -790 -800)))))
(let ([t0 23000]
[t1 25000])
(plot
(list
;(parametric (λ ([t : Real]) (lowerborder (fl t))) t0 t1 #:color 'green)
(parametric (λ ([t : Real]) (highrborder (fl t))) t0 t1 #:color 'green)
(contours (λ ([l : Real] [k : Real]) (flpoisson-log-pdf (fl l) (fl k)))
(car (highrborder (fl t0)))(car (highrborder (fl t1)))
(cadr (highrborder (fl t0)))(cadr (highrborder (fl t1)))
#:label "p-pdf"
#:levels '(-740 -741 -742 -743 -744 -745 -746 -747 -748)))))
;|#
(let
#;([t0 .99999999999e27]
[t1 1.0000000000e27])
#;([t0 .99999999999e26]
[t1 1.0000000000e26])
([t0 .999999999e22]
[t1 1.00000000e22])
#;([t0 .9999999e18]
[t1 1.000000e18])
#;([t0 .999e10]
[t1 1.0e10])
#;([t0 .9e5]
[t1 1.0e5])
#;([t0 0]
[t1 2e4])
(define l- (max 0.0 (min (car (lowerborder (fl t0)))(car (highrborder (fl t0))))))
(define l+ (max (car (lowerborder (fl t1)))(car (highrborder (fl t1)))))
(define k- (max 0.0 (min (cadr (lowerborder (fl t0)))(cadr (highrborder (fl t0))))))
(define k+ (max (cadr (lowerborder (fl t1)))(cadr (highrborder (fl t1)))))
(parameterize ()
#;([plot-x-transform log-transform]
[plot-x-ticks (log-ticks)]
[plot-y-transform log-transform]
[plot-y-ticks (log-ticks)])
(plot
(list
(parametric (λ ([t : Real]) (lowerborder (fl t))) t0 t1 #:color 'green)
(parametric (λ ([t : Real]) (highrborder (fl t))) t0 t1 #:color 'green)
(parametric (λ ([t : Real]) (altr0border (fl t))) t0 t1 #:color 'Red)
(parametric (λ ([t : Real]) (altr1border (fl t))) t0 t1 #:color 'black)
(parametric (λ ([t : Real]) (altrLborder (fl t))) t0 t1 #:color 'blue)
(parametric (λ ([t : Real]) (altr2border (fl t))) t0 t1 #:color 'gray)
(isoline (λ ([l : Real] [k : Real]) (flpoisson-log-pdf+ (fl l) (fl k)))
-740
l- l+
k- k+))
#:width 800
#:height 800
#:x-min l- #:x-max l+
#:y-min k- #:y-max k+)))
;#|
(define (plotit [fct : (-> Flonum (List Flonum Flonum))] [t0 : Real] [t1 : Real])
(list
(plot
(list
(function (λ ([t : Real]) (apply flpoisson-log-pdf1 (fct (fl t)))) t0 t1 #:color 'blue)
(function (λ ([t : Real]) (apply flpoisson-log-pdf2 (fct (fl t)))) t0 t1 #:color 'green)
(function (λ ([t : Real]) (apply flpoisson-log-pdf (fct (fl t)))) t0 t1 #:color 'gray)))
(plot
(function (λ ([t : Real]) (abs (- (apply flpoisson-log-pdf1 (fct (fl t)))
(apply flpoisson-log-pdf2 (fct (fl t)))))) t0 t1))
(format "from ~a to ~a\n" (fct (fl t0)) (fct (fl t1))))
)
(define t0 10)
(define t1 40000)
(plotit lowerborder t0 t1)
(plotit highrborder t0 t1)
(plotit altr1border t0 t1)
(plotit altrLborder t0 t1)
(plotit altr2border t0 t1)
;|#
)
(module+ test1
(require plot
math/bigfloat
racket/list)
(define (bfpois [l : Bigfloat] [k : Bigfloat]) : Bigfloat
(bf/ (bf* (bfexpt l k)
(bfexp (bf- l)))
(bfgamma (bf+ k (bf 1)))))
(define (bflogpois [l : Bigfloat] [k : Bigfloat]) : Bigfloat
(bf- (bf- (bf* k (bflog l))
l)
(bflog-gamma (bf+ k (bf 1)))))
(define nr 1000000)
(define pts1
(for*/list : (Listof (List Real Real Real))
([i : Integer (in-range nr)]
[k : Real (in-value (* 1e6 (random)))]
[l : Real (in-value (* 1e6 (random)))]
[diff : Real
(in-value (abs
(- ;(bigfloat->flonum (bflogpois (bf l)(bf k)))
(bigfloat->flonum (bflog (bfpois (bf l)(bf k))))
(flpoisson-log-pdf (fl l)(fl k)))))]
#:unless (nan? diff))
(list l k diff)))
(define NR (length pts1))
NR
(define pts2
((inst sort (List Real Real Real))
pts1 <
#:key (λ ([x : (List Real Real Real)]) (list-ref x 2))))
(define L000% (last (list-ref pts2 (floor (* 000/100 NR)))))
(define L010% (last (list-ref pts2 (floor (* 010/100 NR)))))
(define L050% (last (list-ref pts2 (floor (* 050/100 NR)))))
(define L090% (last (list-ref pts2 (floor (* 090/100 NR)))))
(define L100% (last (list-ref pts2 (floor (- (* 100/100 NR) 1)))))
(list L000% L010% L050% L090% L100%)
(last pts2)
(define-values (pts-ok pts3)
(partition (λ ([x : (List Real Real Real)])
(< (list-ref x 2) 1e-15))
pts2))
(define-values (pts-b1 pts4)
(partition (λ ([x : (List Real Real Real)])
(< (list-ref x 2) 1e-10))
pts3))
(define-values (pts-b2 pts-no)
(partition (λ ([x : (List Real Real Real)])
(< (list-ref x 2) 1e-7))
pts4))
(plot (list (points pts-ok #:size 1 #:color 'Green)
(points pts-b1 #:size 1 #:color 'OliveDrab)
(points pts-b2 #:size 1 #:color 'Blue)
(points pts-no #:size 2 #:color 'Red))
#:width 800 #:height 800)
)
(module+ test2
(require plot
math/bigfloat
racket/list)
(define (bfpois [l : Bigfloat] [k : Bigfloat]) : Bigfloat
(bf/ (bf* (bfexpt l k)
(bfexp (bf- l)))
(bfgamma (bf+ k (bf 1)))))
(define (bflogpois [l : Bigfloat] [k : Bigfloat]) : Bigfloat
(bf- (bf- (bf* k (bflog l))
l)
(bflog-gamma (bf+ k (bf 1)))))
(define nr 100000)
(define l- 600)(define l+ 1100)(define k- 0)(define k+ 80)
;(define l- 0)(define l+ 5e3)(define k- 0)(define k+ 5e3)
;(define l- 0)(define l+ 5e3)(define k- 0)(define k+ 5e3)
;(define l- 1e4)(define l+ 2.5e4)(define k- 1e4)(define k+ 2.5e4)
;(define l- 0)(define l+ 1.1e6)(define k- 0)(define k+ 1e6)
;(define l- 0)(define l+ 1.1e7)(define k- 0)(define k+ 1e7)
;(define l- 9e9)(define l+ 1e10)(define k- 9e9)(define k+ 1e10)(bf-precision 2048)
(struct ans ([kl : (List Real Real)]
[bf : (List Real Real)]
[lbf : (List Real Real Real Real Real)]
[diffs : (List Real Real Real Real Real Real Real)])
#:transparent)
(define (adif [x : ans] [y : Integer]) (list-ref (ans-diffs x) y))
(define pts
(for*/list : (Listof ans)
([i : Integer (in-range nr)]
[k : Real (in-value (+ k- (* k+ (random))))]
[l : Real (in-value (+ l- (* l+ (random))))]
[vb : Real (in-value (bigfloat->flonum (bfpois (bf l)(bf k))))]
[vf : Real (in-value (flpoisson-pdf (fl l)(fl k)))]
[vlb1 : Real (in-value (bigfloat->flonum (bflog (bfpois (bf l)(bf k)))))]
[vlb2 : Real (in-value (bigfloat->flonum (bflogpois (bf l)(bf k))))]
[vlf1 : Real (in-value (flpoisson-log-pdf (fl l)(fl k)))]
[vlf2 : Real (in-value (flpoisson-log-pdf+ (fl l)(fl k)))]
[vlf3 : Real (in-value (flpoisson-log-pdf++ (fl l)(fl k)))]
[diff0 : Real (in-value (abs (- vb vf)))]
[diff1 : Real (in-value (abs (- vlb1 vlb2)))]
[diff2 : Real (in-value (abs (- vlb1 vlf1)))]
[diff3 : Real (in-value (abs (- vlb1 vlf2)))]
[diff4 : Real (in-value (abs (- vlb2 vlf1)))]
[diff5 : Real (in-value (abs (- vlb2 vlf2)))]
[diff6 : Real (in-value (abs (- vlb2 vlf3)))])
(ans (list l k)
(list vb vf)
(list vlb1 vlb2 vlf1 vlf2 vlf3)
(list diff0 diff1 diff2 diff3 diff4 diff5 diff6))))
(define (->pts [x : (Listof ans)]) (map ans-kl x))
(define (sort&print [pts0 : (Listof ans)] [slct : (-> ans Real)])
(define-values (prnt1 pts1)
(partition (λ ([x : ans])
(< (slct x) 1e-15))
(filter (λ ([x : ans])(not (nan? (slct x)))) pts0)))
(define-values (prnt2 pts2)
(partition (λ ([x : ans])
(< (slct x) 1e-13))
pts1))
(define-values (prnt3 pts3)
(partition (λ ([x : ans])
(< (slct x) 1e-10))
pts2))
(define-values (prnt4 pts4)
(partition (λ ([x : ans])
(< (slct x) 1e-7))
pts3))
(parameterize ()
#;([plot-x-transform log-transform]
[plot-x-ticks (log-ticks)]
[plot-y-transform log-transform]
[plot-y-ticks (log-ticks)])
(plot (list (points (->pts prnt1) #:size 1 #:color 'Green)
(points (->pts prnt2) #:size 1 #:color 'Orange)
(points (->pts prnt3) #:size 1 #:color 'DarkGreen)
(points (->pts prnt4) #:size 1 #:color 'Blue)
(points (->pts pts4) #:size 2 #:color 'Red)
;borders of norm calculation
(lines `((500 ,k+)(500 1573.312629)) #:color 'Gray)
(function (λ ([l : Real]) (+ l (* 48 (sqrt (abs l))))) 500.0 l+ #:color 'Gray)
(function (λ ([l : Real]) (- l (* 48 (sqrt (abs l))))) 500.0 l+ #:color 'Gray)
(lines `((0 1)(700 1)(700 100)(0 100)) #:color 'Gray)
(lines `((0 1)(,l+ 1)) #:color 'Gray)
;borders of log calculation
(lines `((,l- 0.05)(,l+ 0.05)) #:color 'Black)
(lines `((1 0)(1 1)(0 1)) #:color 'Black)
(lines `((,l+ 2)(2 2)(2 ,k+)) #:color 'Black)
(function (λ ([l : Real]) (- l (* 29 (sqrt (abs l))))) 2.0 l+ #:color 'Black)
(function (λ ([l : Real]) (+ l (* 40 (sqrt (abs l))))) 2.0 l+ #:color 'Black)
)
;#:width 800 #:height 800
#:x-min l- #:x-max l+
#:y-min k- #:y-max k+
)))
;(sort&print pts (λ ([x : ans]) (adif x 0)))
;(sort&print pts (λ ([x : ans]) (adif x 1)))
;(sort&print pts (λ ([x : ans]) (adif x 2)))
;(sort&print pts (λ ([x : ans]) (adif x 3)))
(sort&print pts (λ ([x : ans]) (adif x 4)))
(sort&print pts (λ ([x : ans]) (adif x 5)))
(sort&print pts (λ ([x : ans]) (adif x 6)))
)
(module+ test3
(define b : Positive-Integer (+ 1 (random 100000)))
(define l 2.5e5)
(define k 2.5e5)
(random-seed b)(collect-garbage)(collect-garbage)
(time (for ([i (in-range 1000000)])
(flpoisson-log-pdf (fl* l (random))(fl* k (random)))))
(random-seed b)(collect-garbage)(collect-garbage)
(time (for ([i (in-range 1000000)])
(flpoisson-log-pdf+ (fl* l (random))(fl* k (random)))))
(random-seed b)(collect-garbage)(collect-garbage)
(time (for ([i (in-range 1000000)])
(flpoisson-log-pdf++ (fl* l (random))(fl* k (random)))))
)