Skip to content

Commit e60a581

Browse files
jw3126devmotion
andauthored
speed up Semicirle sampling (#1580)
* speed up Semicirle sampling * fix rand for Semicircle * add semicircle rand tests * Update src/univariate/continuous/semicircle.jl Co-authored-by: David Widmann <devmotion@users.noreply.github.com> * improve semicirle rand tests Co-authored-by: David Widmann <devmotion@users.noreply.github.com>
1 parent 254f27f commit e60a581

File tree

2 files changed

+40
-1
lines changed

2 files changed

+40
-1
lines changed

src/univariate/continuous/semicircle.jl

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -72,4 +72,14 @@ function cdf(d::Semicircle, x::Real)
7272
end
7373
end
7474

75+
function rand(rng::AbstractRNG, d::Semicircle)
76+
# Idea:
77+
# sample polar coodinates r,θ
78+
# of point uniformly distributed on radius d.r half disk
79+
# project onto x axis
80+
θ = rand(rng) # multiple of π
81+
r = d.r * sqrt(rand(rng))
82+
return cospi(θ) * r
83+
end
84+
7585
@quantile_newton Semicircle

test/semicircle.jl

Lines changed: 30 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
using Distributions
2-
using Test
2+
using Random: MersenneTwister
3+
using Test
34

45
d = Semicircle(2.0)
56

@@ -37,3 +38,31 @@ d = Semicircle(2.0)
3738
@test quantile(d, .0) == -2.0
3839
@test quantile(d, .5) == .0
3940
@test quantile(d, 1.0) == +2.0
41+
42+
rng = MersenneTwister(0)
43+
for r in rand(rng, Uniform(0,10), 5)
44+
N = 10^4
45+
semi = Semicircle(r)
46+
sample = rand(rng, semi, N)
47+
mi, ma = extrema(sample)
48+
@test -r <= mi < ma <= r
49+
50+
# test order statistic of sample min is sane
51+
d_min = Beta(1, N)
52+
lo = quantile(d_min, 0.01)
53+
hi = quantile(d_min, 0.99)
54+
@test lo < cdf(semi, mi) < hi
55+
56+
# test order statistic of sample max is sane
57+
d_max = Beta(N, 1)
58+
lo = quantile(d_max, 0.01)
59+
hi = quantile(d_max, 0.99)
60+
@test lo < cdf(semi, ma) < hi
61+
62+
# central limit theorem
63+
dmean = Normal(mean(semi), std(semi)/√(N))
64+
@test quantile(dmean, 0.01) < mean(sample) < quantile(dmean, 0.99)
65+
66+
pvalue = pvalue_kolmogorovsmirnoff(sample, semi)
67+
@test pvalue > 1e-2
68+
end

0 commit comments

Comments
 (0)