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]);

50-element Array{Float64,1}:
 -0.6062895051354262 
 -0.6866992538271283 
 -0.20424076167692148
 -1.4103869920524357 
  0.5998567252400879 
 -0.28465051036862354
  1.2431347147736962 
  0.43903722785668664
  2.9317394372994143 
  0.27821773047328247
  ⋮                  
 -0.6866992538271283 
 -0.6866992538271283 
 -2.214484478969445  
  0.6802664739317872 
  0.27821773047328247
 -0.12383101298522224
 -0.8475187512105296 
  0.19780798178158324
 -1.4907967407441378

Show the first six rows of the dataset.

first(df, 6)

6 rows × 14 columns

	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ℝ₊))

problem_transformation (generic function with 1 method)

Wrap the problem with a transformation, then use Flux for the gradient.

P = TransformedLogDensity(problem_transformation(p), p)
∇P = ADgradient(:ForwardDiff, P);

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.TransformTuple{NamedTuple{(:β, :σ),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)

Object of type Chains, with data of type 1000×3×4 Array{Float64,3}

Iterations        = 1:1000
Thinning interval = 1
Chains            = 1, 2, 3, 4
Samples per chain = 1000
parameters        = a, bA, sigma

2-element Array{ChainDataFrame,1}

Summary Statistics
. Omitted printing of 1 columns
│ Row │ parameters │ mean     │ std      │ naive_se   │ mcse       │ ess     │
│     │ Symbol     │ Float64  │ Float64  │ Float64    │ Float64    │ Any     │
├─────┼────────────┼──────────┼──────────┼────────────┼────────────┼─────────┤
│ 1   │ a          │ 9.68736  │ 0.209196 │ 0.00330768 │ 0.00371194 │ 3721.28 │
│ 2   │ bA         │ -1.08234 │ 0.213982 │ 0.00338335 │ 0.00325931 │ 3592.15 │
│ 3   │ sigma      │ 1.49522  │ 0.154237 │ 0.00243869 │ 0.00324173 │ 3045.84 │

Quantiles

│ Row │ parameters │ 2.5%     │ 25.0%    │ 50.0%   │ 75.0%     │ 97.5%     │
│     │ Symbol     │ Float64  │ Float64  │ Float64 │ Float64   │ Float64   │
├─────┼────────────┼──────────┼──────────┼─────────┼───────────┼───────────┤
│ 1   │ a          │ 9.27672  │ 9.54732  │ 9.68918 │ 9.82899   │ 10.0901   │
│ 2   │ bA         │ -1.50173 │ -1.22384 │ -1.0799 │ -0.940167 │ -0.672096 │
│ 3   │ sigma      │ 1.23234  │ 1.38745  │ 1.483   │ 1.5882    │ 1.83725   │

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)

2-element Array{ChainDataFrame,1}

Summary Statistics
. Omitted printing of 1 columns
│ Row │ parameters │ mean     │ std      │ naive_se   │ mcse       │ ess     │
│     │ Symbol     │ Float64  │ Float64  │ Float64    │ Float64    │ Any     │
├─────┼────────────┼──────────┼──────────┼────────────┼────────────┼─────────┤
│ 1   │ a          │ 9.68736  │ 0.209196 │ 0.00330768 │ 0.00371194 │ 3721.28 │
│ 2   │ bA         │ -1.08234 │ 0.213982 │ 0.00338335 │ 0.00325931 │ 3592.15 │
│ 3   │ sigma      │ 1.49522  │ 0.154237 │ 0.00243869 │ 0.00324173 │ 3045.84 │

Quantiles

│ Row │ parameters │ 2.5%     │ 25.0%    │ 50.0%   │ 75.0%     │ 97.5%     │
│     │ Symbol     │ Float64  │ Float64  │ Float64 │ Float64   │ Float64   │
├─────┼────────────┼──────────┼──────────┼─────────┼───────────┼───────────┤
│ 1   │ a          │ 9.27672  │ 9.54732  │ 9.68918 │ 9.82899   │ 10.0901   │
│ 2   │ bA         │ -1.50173 │ -1.22384 │ -1.0799 │ -0.940167 │ -0.672096 │
│ 3   │ sigma      │ 1.23234  │ 1.38745  │ 1.483   │ 1.5882    │ 1.83725   │

Plot the chains

plot(chns)

end of m4.5d.jl#- This page was generated using Literate.jl.