Exercise 7.5 Solution Example - Hoff, A First Course in Bayesian Statistical Methods
標準ベイズ統計学 演習問題 7.5 解答例

Y_src = readdlm("../../Exercises/interexp.dat", skipstart=1)
function parse_matrix(input_matrix)
    nrow, ncol = size(input_matrix)
    output_matrix = Matrix{Union{Missing, Float64}}(missing, nrow, ncol)

    for i in 1:nrow
        for j in 1:ncol
            if input_matrix[i, j] == "NA"
                output_matrix[i, j] = missing
            else
                output_matrix[i, j] = Float64(input_matrix[i, j])
            end
        end
    end
    return output_matrix
end
Y = parse_matrix(Y_src)
θ̂_A = mean(skipmissing(Y[:, 1]))
θ̂_B = mean(skipmissing(Y[:, 2]))
σ̂²_A = var(skipmissing(Y[:, 1]))
σ̂²_B = var(skipmissing(Y[:, 2]))

# get index of not missing
ind_A = ismissing.(Y[:, 1])
ind_B = ismissing.(Y[:, 2])
Y_comp = Y[ind_A.+ind_B.==0, :]
ρ̂ = cor(Y_comp[:, 1], Y_comp[:, 2])

julia> θ̂_A

24.200487804878044

julia> θ̂_B
24.805348837209298

julia> ρ̂
0.6164509013184667

julia> σ̂²_A
4.092799756097562

julia> σ̂²_B
4.69157785160576

function impute(Y)
    nrow = size(Y, 1)
    Y_imp = copy(Y)

    for i in 1:nrow
        if ismissing(Y[i, 1])
            Y_imp[i, 1] = θ̂_A + (Y[i, 2] - θ̂_B) * ρ̂ * sqrt(σ̂²_A / σ̂²_B)
        end
        if ismissing(Y[i, 2])
            Y_imp[i, 2] = θ̂_B + (Y[i, 1] - θ̂_A) * ρ̂ * sqrt(σ̂²_B / σ̂²_A)
        end
    end
    return Y_imp
end

Y_imp = impute(Y)

# Do a paired sample t-test
using HypothesisTests
OneSampleTTest(Y_imp[:, 1] .- Y_imp[:, 2])