chisq.jl

######
# Homework on Chi² Tests / Loglinear Models on 2x2 tables
# 
# instead of computing it by hand or in R, I wrote
# it in Julia, an awesome new language for technical
# computing. I discovered it yesterday, and it is
# truly intriguing. Do check it out! http://julialang.org/
######

# imports; functions in those modules are global now
using Distributions
using Compat, Mamba, Jags # very cool; Python still lacks a JAGS interface!


####
# Frequentist Statistics
######
function expected(m::Array) # optional type checking
    r1, r2 = mapslices(sum, m, 2)
    c1, c2 = mapslices(sum, m, 1)
    exp = [r1 * c1 r1 * c2;
           r2 * c1 r2 * c2] / sum(m)
    return exp
end


# concise function definitions & matrix operations
pearson(emp, exp) = sum((emp .- exp).^2 ./ exp)
deviance(emp, exp) = 2 * sum(emp .* log(emp ./ exp))
df(m) = map(-, size(m), (1, 1)) |> x -> x[1] * x[2] # pipes! awesome


function test(stat::Function, m::Array)
    d = Chisq(df(m))
    chi = stat(m, expected(m))
    p = pdf(d, chi) * 2
    return Dict(["stat" => chi, "df" => df(m), "p" => p])
end

# Bayesian Statistics
# Instead of doing a Chi²-Test, we can also run a binomial proportion test.
# Given this frequency table below
#
#                correct   incorrect
# right-handed     13         7
#  left-handed     11         9
#
# we can ask if the proportion of people who solved some task.
# is the same for right- and left-handed people
# e.g. ask if 13/20 ~ 11/20

function bayes_prop(m::Array)
    x = m[:, 1]
    n = collect(mapslices(sum, m, 2))

    line = "
    model {
       for (i in 1:length(n)) {
            theta[i] ~ dbeta(1, 1)
            x[i] ~ dbinom(theta[i], n[i])
            x_pred[i] ~ dbinom(theta[i], n[i])
        }
    }"

    len = length(n)
    # comprehensions! ohh yes
    inits = [@Compat.Dict("theta" => rand(Beta(1, 1), len)) for i in 1:len]

    # needed, see https://github.com/goedman/Jags.jl/issues/13
    inits = map((x) -> convert(Dict{ASCIIString, Any}, x), inits)

    data = @Compat.Dict("x" => x, "n" => n)
    monitors = (ASCIIString => Bool)["theta" => true] # typed dictionary
    
    model = Jagsmodel(
        name="prop", model=line,
        monitor=monitors, thin=10,
        update=10000, adapt=1000, nchains=4)

    sim = jags(model, data, inits)
    return sim
end


function diff(sim::Chains)
    n = length(sim.chains)
    res = zeros(n)
    for i in 1:n
        val = sim.value[:, :, i]
        res[i] = mean(val[:, 1] - val[:, 2])
    end
    return mean(res)
end

# example data
m1 = [13 7; 11 9]
m2 = [13 5; 7 10]
m3 = [16 4; 4 16]
ms = (m1, m2, m3)

res_bayes = map((m) -> diff(bayes_prop(m)), ms)
res_pearson = map((m) -> test(pearson, m)["p"], ms)
res_deviance = map((m) -> test(deviance, m)["p"], ms)