Skip to content

Commit f69552c

Browse files
bdeketsoegaard
bdeket
authored and
soegaard
committed
Add student-t distribution
1 parent fa46ced commit f69552c

File tree

8 files changed

+1227
-2
lines changed

8 files changed

+1227
-2
lines changed

‎math-doc/math/scribblings/math-distributions.scrbl

Lines changed: 33 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -677,6 +677,39 @@ which is the square of standard deviation. Construct normal distributions from v
677677
#:x-min -5 #:x-max 5 #:y-label "probability")]
678678
}
679679

680+
@subsection{Student-t Distributions}
681+
682+
@margin-note{Wikipedia:
683+
@hyperlink["https://en.wikipedia.org/wiki/Student's_t-distribution"]{Student-t Distribution}.}
684+
@deftogether[(@defidform[Student-t-Dist]
685+
@defproc[(student-t-dist [freedom Positive-Real] [mean Real 0] [scale Real 1]) Normal-Dist]
686+
@defproc[(student-t-dist-freedom [d Normal-Dist]) Positive-Flonum]
687+
@defproc[(student-t-dist-mean [d Normal-Dist]) Flonum]
688+
@defproc[(student-t-dist-scale [d Normal-Dist]) Flonum])]{
689+
Represents the student-t distribution family parameterized by degree of freedom, mean (also called ``location'') and scale.
690+
691+
@examples[#:eval untyped-eval
692+
(plot (for/list ([ν (in-list '(1 2 3 5.5))]
693+
[μ (in-list '(0 0 0 2.3))]
694+
[σ (in-list '(1 1 1 0.4))]
695+
[i (in-naturals)])
696+
(function (distribution-pdf (student-t-dist ν μ σ))
697+
#:color i #:label (if (and (= μ 0) (= σ 1))
698+
(format "Stud-t(~a)" ν)
699+
(format "Stud-t(~a,~a,~a)" ν μ σ))))
700+
#:x-min -5 #:x-max 10 #:y-label "density")
701+
702+
(plot (for/list ([ν (in-list '(1 2 3 5.5))]
703+
[μ (in-list '(0 0 0 2.3))]
704+
[σ (in-list '(1 1 1 0.4))]
705+
[i (in-naturals)])
706+
(function (ordered-dist-cdf (student-t-dist ν μ σ))
707+
#:color i #:label (if (and (= μ 0) (= σ 1))
708+
(format "Stud-t(~a)" ν)
709+
(format "Stud-t(~a,~a,~a)" ν μ σ))))
710+
#:x-min -5 #:x-max 10 #:y-label "probability")]
711+
}
712+
680713
@subsection{Triangular Distributions}
681714

682715
@margin-note{Wikipedia:

‎math-lib/math/distributions.rkt

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,8 @@
1616
"private/distributions/binomial-dist.rkt"
1717
"private/distributions/bernoulli-dist.rkt"
1818
"private/distributions/poisson-dist.rkt"
19-
"private/distributions/discrete-dist.rkt")
19+
"private/distributions/discrete-dist.rkt"
20+
"private/distributions/student-t-dist.rkt")
2021

2122
(provide (all-from-out
2223
"private/distributions/dist-struct.rkt"
@@ -35,4 +36,5 @@
3536
"private/distributions/binomial-dist.rkt"
3637
"private/distributions/bernoulli-dist.rkt"
3738
"private/distributions/poisson-dist.rkt"
38-
"private/distributions/discrete-dist.rkt"))
39+
"private/distributions/discrete-dist.rkt"
40+
"private/distributions/student-t-dist.rkt"))

‎math-lib/math/private/distributions/impl/student-t-cdf.rkt

Lines changed: 247 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,247 @@
1+
#lang typed/racket
2+
(require "student-t-utils.rkt"
3+
"normal-cdf.rkt"
4+
"../dist-struct.rkt"
5+
"../../functions/incomplete-beta.rkt"
6+
"../../../flonum.rkt")
7+
8+
(provide make-student-t-cdf)
9+
10+
;;; Student t distribution
11+
12+
;; Parameters
13+
; ν - degrees of freedom
14+
; μ - location parameter
15+
; σ - scale parameter
16+
17+
;; Domains
18+
; ν - any positive real number
19+
; μ - any real number
20+
; σ - any real number
21+
22+
23+
;; The Cumulative distribution function (CDF)
24+
25+
; For t>0
26+
; F(t) = 1 - 0.5 I_x(t) (ν/2, 0.5),
27+
; where x(t) = ν/(t²+ν).
28+
29+
; Here I is the regularized incomplete beta function.
30+
31+
32+
33+
;; *************************************************************************************************
34+
;; ** CDF IMPLEMENTAION
35+
;; *************************************************************************************************
36+
(: make-student-t-cdf : (case-> (Real -> (CDF Real))
37+
(Real Real Real -> (CDF Real))))
38+
(define make-student-t-cdf
39+
(case-lambda
40+
; Y ~ σX+μ
41+
[(ν μ σ)
42+
(define F (make-student-t-cdf ν))
43+
(let ([μ (fl μ)] [σ (fl σ)])
44+
(λ (y [log? #f] [1-p? #f])
45+
(define x (/ (- (fl y) μ) σ))
46+
(F x log? 1-p?)))]
47+
48+
; X ~ t(ν)
49+
[(ν)
50+
(let ([ν (fl ν)])
51+
(define ν/2 (fl/ ν 2.))
52+
53+
;; *************************************************************************************************
54+
;; ** CDF FUNCTIONS
55+
;; *************************************************************************************************
56+
(: cdf4 (Flonum -> Flonum))
57+
(define (cdf4 x)
58+
(flbeta1/2-regularized ν x (/ (+ 1. (/ (/ ν x) x))) (/ (+ 1. (* x (/ x ν)))) #f #t #t))
59+
60+
(: cdf4-hypergeom (Flonum -> Flonum))
61+
(define (cdf4-hypergeom x)
62+
(define X (/ (+ 1. (/ (/ ν x) x))))
63+
(define Y (/ (+ 1. (* x (/ x ν)))))
64+
(define log-X (if (< X 0.5) (fllog X) (fllog1p (- Y))))
65+
(define log-Y (if (< Y 0.5) (fllog Y) (fllog1p (- X))))
66+
(define log? #f) (define 1-? #t) (define r? #t)
67+
68+
(if (or (< 1e4 ν) (< -0.5 x)) +nan.0
69+
(let ([A ν][x Y][y X][log-x log-Y][log-y log-X])
70+
(define a+2/2 (fl+ A 2.0))
71+
(define: z : Flonum
72+
(let loop ([z 0.0] [dz 1.0] [2n 0.0] [i -1.0])
73+
;(printf "z = ~v dz = ~v i = ~v~n" z dz i)
74+
(define new-dz (fl* (fl- dz (fl/ dz (fl+ a+2/2 2n))) x))
75+
(cond [(zero? i) (fl+ z new-dz)]
76+
[else
77+
(let ([i (if (and (i . fl< . 0.0)
78+
((flabs new-dz) . fl<= . (flabs dz))
79+
((flabs new-dz) . fl<= . (fl* (fl* 0.5 epsilon.0) (flabs z))))
80+
3.0
81+
(fl- i 1.0))])
82+
(loop (fl+ z new-dz) new-dz (fl+ 2n 2.0) i))])))
83+
(define c (/ (* (flexpt x (* 0.5 A)) (flsqrt y)) (beta1/2 A)))
84+
;(println (list a+b/2 a+2/2 (/ a+b/2 a+2/2) z c))
85+
(fl/ (fl+ c (fl* z c)) A))))
86+
87+
(: cdf-largex (Flonum -> Flonum))
88+
(define (cdf-largex x)
89+
;;(B (ν/ν+x²) ; ν/2 1/2)
90+
;; 0th therm of continuous fraction when ν+x²≈x²
91+
(define _x_ (flabs x))
92+
(fl/ (flexpt (fl/ (flsqrt ν) _x_) ν)
93+
(fl* ν (beta1/2 ν))))
94+
95+
(: cdf-scalednorm (Flonum -> Flonum))
96+
(define (cdf-scalednorm x)
97+
(define V (/ 1. 4. ν))
98+
(standard-flnormal-cdf (/ (* x (- 1. V))
99+
(flsqrt (+ 1. (* x x 2. V))))))
100+
101+
(: invert ((Flonum -> Flonum) -> (Flonum -> Flonum)))
102+
(define (invert f)
103+
(λ (x)
104+
(cond
105+
[(flnan? x) +nan.0]
106+
[else
107+
((if (< x 0.) values (λ ([x : Flonum]) (fl- 1. x)))
108+
(let ([x (if (< x 0.) x (- x))])
109+
(cond
110+
[(= x -inf.0) 0.]
111+
[(< -1e-17 x) 0.5]
112+
[else (f x)])))])))
113+
114+
(: cdf : (Flonum -> Flonum))
115+
(define cdf
116+
(cond
117+
[(flnan? ν) (λ (x) +nan.0)]
118+
[(< ν 1e-19) (λ (x) 0.5)]
119+
[(< ν 1e3)
120+
(define large-lim (if (< ν 1.) (max -1e21 (* -1e8 ν)) -1e21))
121+
(define geo-lim (if (< ν 2e2) -1e0 -1e1))
122+
(invert
123+
(λ (x)
124+
(cond
125+
[(< x large-lim) (cdf-largex x)]
126+
[(< x geo-lim) (cdf4-hypergeom x)]
127+
[else (cdf4 x)])))]
128+
[(< ν 1e7)
129+
(invert
130+
(λ (x)
131+
(cond
132+
[(< x -1e2) 0.]
133+
[else (cdf4 x)])))]
134+
[(< ν 1e20)
135+
(define norm-lim (* -1e-3 (flexpt ν #i1/3)))
136+
(invert
137+
(λ (x)
138+
(cond
139+
[(< norm-lim x) (cdf-scalednorm x)]
140+
[(< x -1e2) 0.]
141+
[else (cdf4 x)])))]
142+
[else
143+
(λ (x) (standard-flnormal-cdf x))]
144+
))
145+
146+
;; *************************************************************************************************
147+
;; ** CDF LOG FUNCTIONS
148+
;; *************************************************************************************************
149+
(: lcdf1 (Flonum -> Flonum))
150+
(define (lcdf1 x)
151+
(fllog (cdf x)))
152+
153+
(: lcdf1.1 (Flonum -> Flonum))
154+
(define (lcdf1.1 x)
155+
(fllog1p (- (cdf (- x)))))
156+
157+
(: lcdf2 (Flonum -> Flonum))
158+
(define (lcdf2 x)
159+
(define z (fl/ ν (fl+ (fl* x x) ν)))
160+
(define ν/2 (/ ν 2.))
161+
(if (< x 0.)
162+
(+ (fllog 0.5) (log-beta-inc ν/2 0.5 z #f #t))
163+
(fllog1p (* -.5 (beta-inc ν/2 0.5 z #f #t)))))
164+
165+
(: lcdf5 (Flonum -> Flonum))
166+
(define (lcdf5 x)
167+
(define z (fl/ (fl+ 1. (fl/ (fl/ ν x) x))))
168+
(define ν/2 (/ ν 2.))
169+
(if (< x 0.)
170+
(+ (fllog 0.5) (log-beta-inc 0.5 ν/2 z #t #t))
171+
(fllog1p (* -.5 (beta-inc 0.5 ν/2 z #t #t)))))
172+
173+
(: lcdf9 (Flonum -> Flonum))
174+
(define (lcdf9 x)
175+
;;(B (ν/ν+x²) ; ν/2 1/2)
176+
;; 0th therm of continuous fraction when ν+x²≈x²
177+
(define _x_ (flabs x))
178+
(fl- (fl* ν (fl- (fl* 0.5 (fllog ν)) (fllog _x_)))
179+
(fl+ (fllog ν) (logbeta1/2 ν))))
180+
181+
(define lg1/2 (fllog 0.5))
182+
(define 0.5sqrtmax (* .5 (flsqrt +max.0)))
183+
184+
(: logcheck ((Flonum -> Flonum) -> (Flonum -> Flonum)))
185+
(define (logcheck f)
186+
(λ (x)
187+
(cond
188+
[(flnan? x) +nan.0]
189+
[(< (flabs x) 1e-17) lg1/2]
190+
[else (f x)])))
191+
192+
(: log-cdf : (Flonum -> Flonum))
193+
(define log-cdf
194+
(cond
195+
[(nan? ν) (λ (x) +nan.0)]
196+
[(< ν 1e-20) (λ (x) lg1/2)]
197+
[(< ν 1.)
198+
(logcheck
199+
(λ (x)
200+
(cond
201+
[(< x 0.) (lcdf1 x)]
202+
[else (lcdf1.1 x)])))]
203+
[(< ν 1e18)
204+
(define large-lim (* -1e7 (flexpt ν 0.5)))
205+
(define log-lim (* -1e0 (flexpt ν 0.5)))
206+
(logcheck
207+
(λ (x)
208+
(cond
209+
[(< x large-lim) (lcdf9 x)]
210+
[(< x log-lim) (lcdf2 x)]
211+
[(< x -10.) (lcdf5 x)]
212+
[(< x 0.) (lcdf1 x)]
213+
[else (lcdf1.1 x)])))]
214+
[(< ν 0.5sqrtmax)
215+
(define large-lim (if (< ν 1e20) (* -1e7 (flexpt ν 0.5)) (- 0.5sqrtmax)))
216+
(define log-lim (* -1e0 (flexpt ν 0.5)))
217+
(define mid-lim (* -1e-8 (flexpt ν 0.5)))
218+
(logcheck
219+
(λ (x)
220+
(cond
221+
[(< x large-lim) (lcdf9 x)]
222+
[(< x log-lim) (lcdf2 x)]
223+
[(< x mid-lim) (lcdf5 x)]
224+
[else
225+
(standard-flnormal-log-cdf x)])))]
226+
[(< ν +inf.0)
227+
(define large-lim (- 0.5sqrtmax))
228+
(define log-lim (* -1e0 (flexpt ν 0.5)))
229+
(logcheck
230+
(λ (x)
231+
(cond
232+
[(< x large-lim) (lcdf9 x)]
233+
[(< x log-lim) (lcdf2 x)]
234+
[else
235+
(standard-flnormal-log-cdf x)])))]
236+
[else
237+
(λ (x) (standard-flnormal-log-cdf x))]))
238+
239+
(: result-cdf : (CDF Real))
240+
(define (result-cdf x [log? #f] [1-p? #f])
241+
(let* ([x (fl (if 1-p? (- x) x))])
242+
(cond
243+
[log? (log-cdf x)]
244+
[else (cdf x)])))
245+
result-cdf)]
246+
))
247+

0 commit comments

Comments
 (0)