TuringModels

Admit or reject

Data

import Random
import TuringModels

using CSV
using DataFrames
using StatsFuns
using Turing

Random.seed!(1)

file_path = joinpath(TuringModels.project_root, "data", "UCBadmit.csv")
df = CSV.read(file_path, DataFrame; delim=';')
12×5 DataFrame
 Row │ dept       gender     admit  reject  applications
     │ InlineSt…  InlineSt…  Int64  Int64   Int64
─────┼───────────────────────────────────────────────────
   1 │ A          male         512     313           825
   2 │ A          female        89      19           108
   3 │ B          male         353     207           560
   4 │ B          female        17       8            25
   5 │ C          male         120     205           325
   6 │ C          female       202     391           593
   7 │ D          male         138     279           417
   8 │ D          female       131     244           375
   9 │ E          male          53     138           191
  10 │ E          female        94     299           393
  11 │ F          male          22     351           373
  12 │ F          female        24     317           341

Model

@model m_pois(admit, reject) = begin
   α₁ ~ Normal(0,100)
   α₂ ~ Normal(0,100)

   for i ∈ 1:length(admit)
       λₐ = exp(α₁)
       λᵣ = exp(α₂)
       admit[i] ~ Poisson(λₐ)
       reject[i] ~ Poisson(λᵣ)
   end
end;

Output

chns = sample(m_pois(df.admit, df.reject), NUTS(), 1000)
Chains MCMC chain (1000×14×1 Array{Float64, 3}):

Iterations        = 501:1:1500
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 3.97 seconds
Compute duration  = 3.97 seconds
parameters        = α₂, α₁
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size

Summary Statistics
  parameters      mean       std   naive_se      mcse        ess      rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64    Float64   Float64       Float64

          α₁    4.9851    0.0252     0.0008    0.0008   818.3047    0.9991      206.0183
          α₂    5.4417    0.0185     0.0006    0.0006   979.1879    1.0008      246.5226

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5%
      Symbol   Float64   Float64   Float64   Float64   Float64

          α₁    4.9370    4.9683    4.9848    5.0028    5.0317
          α₂    5.4041    5.4307    5.4421    5.4531    5.4782

using StatsPlots

StatsPlots.plot(chns)

Original output

m_10_yyt_result = "
    mean   sd 5.5% 94.5% n_eff Rhat
 a1 4.99 0.02 4.95  5.02  2201    1
 a2 5.44 0.02 5.41  5.47  2468    1
";