|
| 1 | +using StableRNGs: StableRNG |
| 2 | + |
1 | 3 | @timed_testset "turing_model" begin |
2 | 4 | DATA_DIR = joinpath("..", "data") |
3 | 5 | kidiq = CSV.read(joinpath(DATA_DIR, "kidiq.csv"), DataFrame) |
|
8 | 10 | f = @formula(kid_score ~ mom_iq * mom_hs) |
9 | 11 | @testset "standardize=false" begin |
10 | 12 | m = turing_model(f, kidiq) |
11 | | - chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2) |
| 13 | + chn = sample(StableRNG(123), m, NUTS(), MCMCThreads(), 2_000, 2) |
12 | 14 | @test summarystats(chn)[:α, :mean] ≈ 31.80 atol = 2.0 |
13 | 15 | @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.507 atol = 0.2 |
14 | 16 | @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ 0.22 atol = 0.2 |
15 | 17 | end |
16 | 18 |
|
17 | 19 | @testset "standardize=true" begin |
18 | 20 | m = turing_model(f, kidiq; standardize=true) |
19 | | - chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2) |
| 21 | + chn = sample(StableRNG(123), m, NUTS(), MCMCThreads(), 2_000, 2) |
20 | 22 | @test summarystats(chn)[:α, :mean] ≈ 0.000 atol = 0.2 |
21 | 23 | @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.648 atol = 0.2 |
22 | 24 | @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ 0.849 atol = 0.2 |
|
25 | 27 | @testset "custom_priors" begin |
26 | 28 | priors = CustomPrior(Normal(), Normal(28, 5), nothing) |
27 | 29 | m = turing_model(f, kidiq; priors) |
28 | | - chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2) |
| 30 | + chn = sample(StableRNG(123), m, NUTS(), MCMCThreads(), 2_000, 2) |
29 | 31 | @test summarystats(chn)[:α, :mean] ≈ 28.758 atol = 2.0 |
30 | 32 | @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.539 atol = 0.2 |
31 | 33 | @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ 0.3863 atol = 0.2 |
32 | 34 | end |
33 | 35 | @testset "explicit calling Normal" begin |
34 | 36 | m = turing_model(f, kidiq; model=Normal) |
35 | | - chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2) |
| 37 | + chn = sample(StableRNG(123), m, NUTS(), MCMCThreads(), 2_000, 2) |
36 | 38 | @test summarystats(chn)[:α, :mean] ≈ 31.80 atol = 2.0 |
37 | 39 | @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.507 atol = 0.2 |
38 | 40 | @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ 0.22 atol = 0.2 |
|
42 | 44 | f = @formula(kid_score ~ mom_iq * mom_hs) |
43 | 45 | @testset "standardize=false" begin |
44 | 46 | m = turing_model(f, kidiq; model=TDist) |
45 | | - chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2) |
| 47 | + chn = sample(StableRNG(123), m, NUTS(), MCMCThreads(), 2_000, 2) |
46 | 48 | @test summarystats(chn)[:α, :mean] ≈ 33.31 atol = 2.0 |
47 | 49 | @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.519 atol = 0.2 |
48 | 50 | @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ 0.340 atol = 0.2 |
|
52 | 54 | @testset "custom_priors" begin |
53 | 55 | priors = CustomPrior(Normal(), Normal(28, 5), Exponential(2)) |
54 | 56 | m = turing_model(f, kidiq; model=TDist, priors) |
55 | | - chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2) |
| 57 | + chn = sample(StableRNG(123), m, NUTS(), MCMCThreads(), 2_000, 2) |
56 | 58 | @test summarystats(chn)[:α, :mean] ≈ 28.565 atol = 2.0 |
57 | 59 | @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.551 atol = 0.2 |
58 | 60 | @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ 0.255 atol = 0.2 |
|
63 | 65 | f = @formula(switch ~ arsenic + dist + assoc + educ) |
64 | 66 | @testset "standardize=false" begin |
65 | 67 | m = turing_model(f, wells; model=Bernoulli) |
66 | | - chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2) |
| 68 | + chn = sample(StableRNG(123), m, NUTS(), MCMCThreads(), 2_000, 2) |
67 | 69 | @test summarystats(chn)[:α, :mean] ≈ -0.153 atol = 0.2 |
68 | 70 | @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.467 atol = 0.2 |
69 | 71 | @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ -0.009 atol = 0.2 |
|
72 | 74 | @testset "custom_priors" begin |
73 | 75 | priors = CustomPrior(Normal(), Normal(), nothing) |
74 | 76 | m = turing_model(f, wells; model=Bernoulli, priors) |
75 | | - chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2) |
| 77 | + chn = sample(StableRNG(123), m, NUTS(), MCMCThreads(), 2_000, 2) |
76 | 78 | @test summarystats(chn)[:α, :mean] ≈ -0.155 atol = 0.2 |
77 | 79 | @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.468 atol = 0.2 |
78 | 80 | @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ -0.009 atol = 0.2 |
|
82 | 84 | f = @formula(y ~ roach1 + treatment + senior + exposure2) |
83 | 85 | @testset "standardize=false" begin |
84 | 86 | m = turing_model(f, roaches; model=Poisson) |
85 | | - chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2) |
| 87 | + chn = sample(StableRNG(123), m, NUTS(), MCMCThreads(), 2_000, 2) |
86 | 88 | @test summarystats(chn)[:α, :mean] ≈ 2.969 atol = 0.5 |
87 | 89 | @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.006 atol = 0.2 |
88 | 90 | @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ -0.5145 atol = 0.2 |
|
91 | 93 | @testset "custom_priors" begin |
92 | 94 | priors = CustomPrior(Normal(2, 5), Normal(), nothing) |
93 | 95 | m = turing_model(f, roaches; model=Poisson, priors) |
94 | | - chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2) |
| 96 | + chn = sample(StableRNG(123), m, NUTS(), MCMCThreads(), 2_000, 2) |
95 | 97 | @test summarystats(chn)[:α, :mean] ≈ 2.963 atol = 0.5 |
96 | 98 | @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.006 atol = 0.2 |
97 | 99 | @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ -0.5145 atol = 0.2 |
|
101 | 103 | f = @formula(y ~ roach1 + treatment + senior + exposure2) |
102 | 104 | @testset "standardize=false" begin |
103 | 105 | m = turing_model(f, roaches; model=NegativeBinomial) |
104 | | - chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2) |
| 106 | + chn = sample(StableRNG(123), m, NUTS(), MCMCThreads(), 2_000, 2) |
105 | 107 | @test summarystats(chn)[:α, :mean] ≈ 2.448 atol = 0.5 |
106 | 108 | @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.013 atol = 0.2 |
107 | 109 | @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ -0.734 atol = 0.2 |
|
111 | 113 | @testset "custom_priors" begin |
112 | 114 | priors = CustomPrior(Normal(), Normal(2, 5), Exponential(0.5)) |
113 | 115 | m = turing_model(f, roaches; model=NegativeBinomial, priors) |
114 | | - chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2) |
| 116 | + chn = sample(StableRNG(123), m, NUTS(), MCMCThreads(), 2_000, 2) |
115 | 117 | @test summarystats(chn)[:α, :mean] ≈ 2.401 atol = 0.5 |
116 | 118 | @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.013 atol = 0.2 |
117 | 119 | @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ -0.723 atol = 0.2 |
|
121 | 123 | @timed_testset "Hierarchical Model" begin |
122 | 124 | f = @formula(y ~ (1 | cheese) + background) |
123 | 125 | m = turing_model(f, cheese) |
124 | | - chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2) |
| 126 | + chn = sample(StableRNG(123), m, NUTS(), MCMCThreads(), 2_000, 2) |
125 | 127 | @test summarystats(chn)[:α, :mean] ≈ 68.07 atol = 2.0 |
126 | 128 | @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 6.60 atol = 0.2 |
127 | 129 | @test summarystats(chn)[Symbol("zⱼ[1]"), :mean] ≈ 0.348 atol = 0.2 |
|
0 commit comments