Linear regression
using DynamicHMCModels, MCMCChains
ProjDir = rel_path_d("..", "scripts", "05")
cd(ProjDir)
Import the dataset.
snippet 5.1
wd = CSV.read(rel_path("..", "data", "WaffleDivorce.csv"), delim=';')
df = convert(DataFrame, wd);
mean_ma = mean(df[:MedianAgeMarriage])
df[:MedianAgeMarriage_s] = convert(Vector{Float64},
(df[:MedianAgeMarriage]) .- mean_ma)/std(df[:MedianAgeMarriage]);
Show the first six rows of the dataset.
first(df, 6)
Location | Loc | Population | MedianAgeMarriage | Marriage | Marriage SE | Divorce | Divorce SE | WaffleHouses | South | Slaves1860 | Population1860 | PropSlaves1860 | MedianAgeMarriage_s | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
String⍰ | String⍰ | Float64⍰ | Float64⍰ | Float64⍰ | Float64⍰ | Float64⍰ | Float64⍰ | Int64⍰ | Int64⍰ | Int64⍰ | Int64⍰ | Float64⍰ | Float64 | |
1 | Alabama | AL | 4.78 | 25.3 | 20.2 | 1.27 | 12.7 | 0.79 | 128 | 1 | 435080 | 964201 | 0.45 | -0.60629 |
2 | Alaska | AK | 0.71 | 25.2 | 26.0 | 2.93 | 12.5 | 2.05 | 0 | 0 | 0 | 0 | 0.0 | -0.686699 |
3 | Arizona | AZ | 6.33 | 25.8 | 20.3 | 0.98 | 10.8 | 0.74 | 18 | 0 | 0 | 0 | 0.0 | -0.204241 |
4 | Arkansas | AR | 2.92 | 24.3 | 26.4 | 1.7 | 13.5 | 1.22 | 41 | 1 | 111115 | 435450 | 0.26 | -1.41039 |
5 | California | CA | 37.25 | 26.8 | 19.1 | 0.39 | 8.0 | 0.24 | 0 | 0 | 0 | 379994 | 0.0 | 0.599857 |
6 | Colorado | CO | 5.03 | 25.7 | 23.5 | 1.24 | 11.6 | 0.94 | 11 | 0 | 0 | 34277 | 0.0 | -0.284651 |
Model $y ∼ Normal(y - Xβ, σ)$. Flat prior for β
, half-T for σ
.
struct WaffleDivorceProblem{TY <: AbstractVector, TX <: AbstractMatrix}
"Observations."
y::TY
"Covariates"
X::TX
end
Make the type callable with the parameters as a single argument.
function (problem::WaffleDivorceProblem)(θ)
@unpack y, X, = problem # extract the data
@unpack β, σ = θ # works on the named tuple too
ll = 0.0
ll += logpdf(Normal(10, 10), X[1]) # a = X[1]
ll += logpdf(Normal(0, 1), X[2]) # b1 = X[2]
ll += logpdf(TDist(1.0), σ)
ll += loglikelihood(Normal(0, σ), y .- X*β)
ll
end
Instantiate the model with data and inits.
N = size(df, 1)
X = hcat(ones(N), df[:MedianAgeMarriage_s]);
y = convert(Vector{Float64}, df[:Divorce])
p = WaffleDivorceProblem(y, X);
p((β = [1.0, 2.0], σ = 1.0))
-2225.6614871340917
Write a function to return properly dimensioned transformation.
problem_transformation(p::WaffleDivorceProblem) =
as((β = as(Array, size(p.X, 2)), σ = asℝ₊))
Wrap the problem with a transformation, then use Flux for the gradient.
P = TransformedLogDensity(problem_transformation(p), p)
∇P = LogDensityRejectErrors(ADgradient(:ForwardDiff, P));
LogDensityRejectErrors{InvalidLogDensityException,LogDensityProblems.ForwardDiffLogDensity{TransformedLogDensity{TransformVariables.TransformNamedTuple{(:β, :σ),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ShiftedExp{true,Float64}}},Main.ex-m5.1d1.WaffleDivorceProblem{Array{Float64,1},Array{Float64,2}}},ForwardDiff.GradientConfig{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##1#2")){TransformedLogDensity{TransformVariables.TransformNamedTuple{(:β, :σ),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ShiftedExp{true,Float64}}},Main.ex-m5.1d1.WaffleDivorceProblem{Array{Float64,1},Array{Float64,2}}}},Float64},Float64,3,Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##1#2")){TransformedLogDensity{TransformVariables.TransformNamedTuple{(:β, :σ),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ShiftedExp{true,Float64}}},Main.ex-m5.1d1.WaffleDivorceProblem{Array{Float64,1},Array{Float64,2}}}},Float64},Float64,3},1}}}}(ForwardDiff AD wrapper for TransformedLogDensity of dimension 3, w/ chunk size 3)
Create an array to hold 1000 samples of 3 parameters in 4 chains
a3d = create_a3d(1000, 3, 4);
trans = as( (β = as(Array, 2), σ = asℝ));
TransformVariables.TransformNamedTuple{(:β, :σ),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.Identity}}((TransformVariables.ArrayTransform{TransformVariables.Identity,1}(TransformVariables.Identity(), (2,)), TransformVariables.Identity()), 3)
Sample from 4 chains and store the draws in the a3d array
for j in 1:4
chain, NUTS_tuned = NUTS_init_tune_mcmc(∇P, 3000);
posterior = TransformVariables.transform.(Ref(problem_transformation(p)),
get_position.(chain));
insert_chain!(a3d, j, posterior, trans);
end;
cnames = ["a", "bA", "sigma"]
chns = Chains(a3d, cnames)
MCMC, adapting ϵ (75 steps)
5.7e-5 s/step ...done
MCMC, adapting ϵ (25 steps)
5.8e-5 s/step ...done
MCMC, adapting ϵ (50 steps)
5.7e-5 s/step ...done
MCMC, adapting ϵ (100 steps)
3.9e-5 s/step ...done
MCMC, adapting ϵ (200 steps)
0.00013 s/step ...done
MCMC, adapting ϵ (400 steps)
3.5e-5 s/step ...done
MCMC, adapting ϵ (50 steps)
6.8e-5 s/step ...done
MCMC (3000 steps)
7.9e-5 s/step ...done
MCMC, adapting ϵ (75 steps)
0.0003 s/step ...done
MCMC, adapting ϵ (25 steps)
8.8e-5 s/step ...done
MCMC, adapting ϵ (50 steps)
5.0e-5 s/step ...done
MCMC, adapting ϵ (100 steps)
4.0e-5 s/step ...done
MCMC, adapting ϵ (200 steps)
3.2e-5 s/step ...done
MCMC, adapting ϵ (400 steps)
6.3e-5 s/step ...done
MCMC, adapting ϵ (50 steps)
3.6e-5 s/step ...done
MCMC (3000 steps)
5.2e-5 s/step ...done
MCMC, adapting ϵ (75 steps)
4.4e-5 s/step ...done
MCMC, adapting ϵ (25 steps)
4.9e-5 s/step ...done
MCMC, adapting ϵ (50 steps)
5.2e-5 s/step ...done
MCMC, adapting ϵ (100 steps)
4.0e-5 s/step ...done
MCMC, adapting ϵ (200 steps)
8.6e-5 s/step ...done
MCMC, adapting ϵ (400 steps)
3.4e-5 s/step ...done
MCMC, adapting ϵ (50 steps)
4.0e-5 s/step ...done
MCMC (3000 steps)
7.6e-5 s/step ...done
MCMC, adapting ϵ (75 steps)
4.9e-5 s/step ...done
MCMC, adapting ϵ (25 steps)
4.9e-5 s/step ...done
MCMC, adapting ϵ (50 steps)
4.7e-5 s/step ...done
MCMC, adapting ϵ (100 steps)
0.00014 s/step ...done
MCMC, adapting ϵ (200 steps)
3.3e-5 s/step ...done
MCMC, adapting ϵ (400 steps)
3.2e-5 s/step ...done
MCMC, adapting ϵ (50 steps)
0.00025 s/step ...done
MCMC (3000 steps)
6.7e-5 s/step ...done
Object of type Chains, with data of type 1000×3×4 Array{Float64,3}
Log evidence = 0.0
Iterations = 1:1000
Thinning interval = 1
Chains = 1, 2, 3, 4
Samples per chain = 1000
parameters = a, bA, sigma
parameters
Mean SD Naive SE MCSE ESS
a 9.6839 0.2128 0.0034 0.0035 1000
bA -1.0859 0.2202 0.0035 0.0033 1000
sigma 1.4934 0.1634 0.0026 0.0034 1000
cmdstan result
cmdstan_result = "
Iterations = 1:1000
Thinning interval = 1
Chains = 1,2,3,4
Samples per chain = 1000
Empirical Posterior Estimates:
Mean SD Naive SE MCSE ESS
a 9.6882466 0.22179190 0.0035068378 0.0031243061 1000
bA -1.0361742 0.21650514 0.0034232469 0.0034433245 1000
sigma 1.5180337 0.15992781 0.0025286807 0.0026279593 1000
Quantiles:
2.5% 25.0% 50.0% 75.0% 97.5%
a 9.253141 9.5393175 9.689585 9.84221500 10.11121000
bA -1.454571 -1.1821025 -1.033065 -0.89366925 -0.61711705
sigma 1.241496 1.4079225 1.504790 1.61630750 1.86642750
";
"\nIterations = 1:1000\nThinning interval = 1\nChains = 1,2,3,4\nSamples per chain = 1000\n\nEmpirical Posterior Estimates:\n Mean SD Naive SE MCSE ESS\n a 9.6882466 0.22179190 0.0035068378 0.0031243061 1000\n bA -1.0361742 0.21650514 0.0034232469 0.0034433245 1000\nsigma 1.5180337 0.15992781 0.0025286807 0.0026279593 1000\n\nQuantiles:\n 2.5% 25.0% 50.0% 75.0% 97.5%\n a 9.253141 9.5393175 9.689585 9.84221500 10.11121000\n bA -1.454571 -1.1821025 -1.033065 -0.89366925 -0.61711705\nsigma 1.241496 1.4079225 1.504790 1.61630750 1.86642750\n"
Extract the parameter posterior means: β
,
describe(chns)
Log evidence = 0.0
Iterations = 1:1000
Thinning interval = 1
Chains = 1, 2, 3, 4
Samples per chain = 1000
parameters = a, bA, sigma
Empirical Posterior Estimates
─────────────────────────────────────────────
parameters
Mean SD Naive SE MCSE ESS
a 9.6839 0.2128 0.0034 0.0035 1000
bA -1.0859 0.2202 0.0035 0.0033 1000
sigma 1.4934 0.1634 0.0026 0.0034 1000
Quantiles
─────────────────────────────────────────────
parameters
2.5% 25.0% 50.0% 75.0% 97.5%
a 8.9484 9.5475 9.6843 9.8237 10.6310
bA -1.9209 -1.2337 -1.0845 -0.9385 -0.3098
sigma 1.0548 1.3744 1.4811 1.5921 2.2590
Plot the chains
plot(chns)
end of m4.5d.jl#- This page was generated using Literate.jl.