Skip to content

Commit

Permalink
Fix residual calculation for default priors (#78)
Browse files Browse the repository at this point in the history 
* Update turing_model.jl

* fix residual calculation

* updating tests after change in residual definition

Co-authored-by: Jose Storopoli <[email protected]>
  • Loading branch information
@burtonjosh @storopoli
burtonjosh and storopoli authored Jan 20, 2023
1 parent cbcb0d0 commit 3869546
Show file tree
Hide file tree
Showing 2 changed files with 24 additions and 24 deletions.
8 changes: 4 additions & 4 deletions src/turing_model.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -177,7 +177,7 @@ function _model(μ_X, σ_X, prior, intercept_ranef, idx, ::Type{Normal})
μ_X=μ_X,
σ_X=σ_X,
prior=prior,
residual=1 / std(y),
residual=std(y),
mad_y=mad(y; normalize=true),
)
α ~ prior.intercept
Expand All @@ -197,7 +197,7 @@ function _model(μ_X, σ_X, prior, intercept_ranef, idx, ::Type{Normal})
end
function _model(μ_X, σ_X, prior, ::Type{Normal})
@model function normal_model(
y, X; predictors=size(X, 2), μ_X=μ_X, σ_X=σ_X, prior=prior, residual=1 / std(y)
y, X; predictors=size(X, 2), μ_X=μ_X, σ_X=σ_X, prior=prior, residual=std(y)
)
α ~ prior.intercept
β ~ filldist(prior.predictors, predictors)
Expand All @@ -219,7 +219,7 @@ function _model(μ_X, σ_X, prior, intercept_ranef, idx, ::Type{TDist})
μ_X=μ_X,
σ_X=σ_X,
prior=prior,
residual=1 / std(y),
residual=std(y),
mad_y=mad(y; normalize=true),
)
α ~ prior.intercept
Expand All @@ -240,7 +240,7 @@ function _model(μ_X, σ_X, prior, intercept_ranef, idx, ::Type{TDist})
end
function _model(μ_X, σ_X, prior, ::Type{TDist})
@model function student_model(
y, X; predictors=size(X, 2), μ_X=μ_X, σ_X=σ_X, prior=prior, residual=1 / std(y)
y, X; predictors=size(X, 2), μ_X=μ_X, σ_X=σ_X, prior=prior, residual=std(y)
)
α ~ prior.intercept
β ~ filldist(prior.predictors, predictors)
Expand Down
40 changes: 20 additions & 20 deletions test/turing_model.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,9 +9,9 @@
@testset "standardize=false" begin
m = turing_model(f, kidiq)
chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
@test summarystats(chn)[:α, :mean] ≈ 29.30 atol = 2.0
@test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.533 atol = 0.2
@test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ 0.593 atol = 0.2
@test summarystats(chn)[:α, :mean] ≈ 31.80 atol = 2.0
@test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.507 atol = 0.2
@test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ 0.22 atol = 0.2
end

@testset "standardize=true" begin
Expand All @@ -33,30 +33,30 @@
@testset "explicit calling Normal" begin
m = turing_model(f, kidiq; model=Normal)
chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
@test summarystats(chn)[:α, :mean] ≈ 29.30 atol = 2.0
@test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.533 atol = 0.2
@test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ 0.593 atol = 0.2
@test summarystats(chn)[:α, :mean] ≈ 31.80 atol = 2.0
@test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.507 atol = 0.2
@test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ 0.22 atol = 0.2
end
end
@timed_testset "TDist Model" begin
f = @formula(kid_score ~ mom_iq * mom_hs)
@testset "standardize=false" begin
m = turing_model(f, kidiq; model=TDist)
chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
@test summarystats(chn)[:α, :mean] ≈ 40.380 atol = 2.0
@test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.478 atol = 0.2
@test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ 0.736 atol = 0.2
@test quantile(chn)[:ν, Symbol("50.0%")] ≈ 1.039 atol = 0.5
@test summarystats(chn)[:α, :mean] ≈ 33.31 atol = 2.0
@test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.519 atol = 0.2
@test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ 0.340 atol = 0.2
@test quantile(chn)[:ν, Symbol("50.0%")] ≈ 2.787 atol = 0.5
end

@testset "custom_priors" begin
priors = CustomPrior(Normal(), Normal(28, 5), Exponential(2))
m = turing_model(f, kidiq; model=TDist, priors)
chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
@test summarystats(chn)[:α, :mean] ≈ 35.506 atol = 2.0
@test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.522 atol = 0.2
@test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ 0.628 atol = 0.2
@test quantile(chn)[:ν, Symbol("50.0%")] ≈ 1.178 atol = 0.5
@test summarystats(chn)[:α, :mean] ≈ 28.565 atol = 2.0
@test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.551 atol = 0.2
@test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ 0.255 atol = 0.2
@test quantile(chn)[:ν, Symbol("50.0%")] ≈ 10.339 atol = 0.5
end
end
@timed_testset "Bernoulli Model" begin
Expand Down Expand Up @@ -112,20 +112,20 @@
priors = CustomPrior(Normal(), Normal(2, 5), Exponential(0.5))
m = turing_model(f, roaches; model=NegativeBinomial, priors)
chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
@test summarystats(chn)[:α, :mean] ≈ 2.422 atol = 0.5
@test summarystats(chn)[:α, :mean] ≈ 2.401 atol = 0.5
@test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.013 atol = 0.2
@test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ -0.732 atol = 0.2
@test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ -0.723 atol = 0.2
@test quantile(chn)[:ϕ⁻, Symbol("50.0%")] ≈ 3.56 atol = 0.2
end
end
@timed_testset "Hierarchical Model" begin
f = @formula(y ~ (1 | cheese) + background)
m = turing_model(f, cheese)
chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
@test summarystats(chn)[:α, :mean] ≈ 68.33 atol = 2.0
@test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 6.928 atol = 0.2
@test summarystats(chn)[Symbol("zⱼ[1]"), :mean] ≈ 0.306 atol = 0.2
@test quantile(chn)[Symbol("zⱼ[2]"), Symbol("50.0%")] ≈ -1.422 atol = 0.5
@test summarystats(chn)[:α, :mean] ≈ 68.07 atol = 2.0
@test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 6.60 atol = 0.2
@test summarystats(chn)[Symbol("zⱼ[1]"), :mean] ≈ 0.348 atol = 0.2
@test quantile(chn)[Symbol("zⱼ[2]"), Symbol("50.0%")] ≈ -1.376 atol = 0.5
end
@testset "Unsupported Model Likelihoods" begin
@test_throws ArgumentError turing_model(@formula(y ~ x), nt_str; model=Binomial)
Expand Down

0 comments on commit 3869546

Please sign in to comment.