-
Notifications
You must be signed in to change notification settings - Fork 0
/
rasch.jl
50 lines (40 loc) · 1.54 KB
/
rasch.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
# Import Packages
using Random, Distributions, Turing, Memoization, LazyArrays, ReverseDiff
using CSV, DataFrames, CategoricalArrays
# helper functions
lazyarray(f, x) = LazyArray(Base.broadcasted(f, x))
# Settings for Turing
Random.seed!(1234)
Turing.setrdcache(true)
Turing.setadbackend(:reversediff)
# Data Handling stuff
url = "https://raw.githubusercontent.com/t-alfers/Rasch-Turing/master/dichotom.txt"
data = CSV.read(download(url), DataFrame, delim = "\t")
data.person = collect(1:nrow(data));
data_responses = data[:, Not([:Anger, :Sex])];
data_long = DataFrames.stack(data_responses, Not(:person), :person);
data_long.variable = CategoricalArray(data_long.variable)
levels!(data_long.variable, unique(data_long.variable))
data_long.variable = levelcode.(data_long.variable)
# statistical modeling
@model rasch(y, ii, jj) = begin
I = maximum(ii);
J = maximum(jj);
N = I * J;
# prior item parameter
β_σ ~ truncated(Cauchy(0, 5), 0, Inf);
β ~ filldist(Normal(0, β_σ), I);
# prior person parameter
θ_σ ~ truncated(Cauchy(0, 5), 0, Inf);
θ ~ filldist(Normal(0, θ_σ), J);
# target distribution
y ~ arraydist(lazyarray(x -> BernoulliLogit(x), θ[jj] .- β[ii]))
end
@time chn = sample(rasch(data_long.value, data_long.variable, data_long.person),
NUTS(1000, 0.8), 2000, progress = true)
# LazyArrays 0.16.16
# NUTS(1000, 0.8), 2000 :: 43.783335 seconds (15.58 M allocations: 3.886 GiB, 0.78% gc time)
#
# LaryArrays 0.22.2
# NUTS(1000, 0.8), 2000 :: 354.621508 seconds (16.82 M allocations: 3.857 GiB, 0.09% gc time)
#