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

Y_src = readdlm("../../Exercises/azdiabetes.dat", header=true)
Y = Y_src[1]
header = Y_src[2] |> vec
df = DataFrame(Y, header)

col_int = [:npreg, :glu, :bp, :skin, :age]
col_float = [:bmi, :ped]
col_str = :diabetes
df[!,col_int] = Int.(df[!,col_int])
df[!,col_float] = Float32.(df[!,col_float])
df[!,col_str] = String.(df[!,col_str])

function ResidualSE(y,X)
    fit = lm(X, y)
    return sum(map(x -> x^2, residuals(fit))) / (size(X, 1) - size(X, 2))
end

function lmGprior(y, X; g=length(y), ν₀=1, σ₀²=ResidualSE(y,X), S=1000)
    n, p = size(X)

    Hg = (g/(g+1)) * X * inv(X'X) * X'
    SSRg = y' * (diagm(ones(n)) - Hg) * y
    σ² = 1 ./ rand(Gamma((ν₀+n)/2, 2/(ν₀*σ₀²+SSRg)), S)

    Vb = g * inv(X'X) /(g+1)
    Eb = Vb * X' * y

    β = Matrix{Float64}(undef, 0, p)
    for i in 1:S
        s = σ²[i]
        β = vcat(β, rand(MvNormal(Eb, Symmetric(s * Vb)))')
    end
    return β, σ²
end

# set prior parameters
g = size(df, 1)
ν₀ = 2
σ₀² = 1

y = df[:, :glu]
X = select(df, Not([:diabetes, :glu])) |> Matrix
X = [ones(size(X, 1)) X] # add intercept
S=5000
BETA, SIGMA² = lmGprior(y, X; g=g, ν₀=ν₀, σ₀²=σ₀², S=S)

# posterior confidence intervals for σ²
ci_σ² = quantile(SIGMA², [0.025, 0.975])

# %%
# posterior confidence intervals for β
ci_β = map(x -> quantile(x, [0.025, 0.975]), eachcol(BETA))

function lpyX(y, X; g=length(y), ν₀=1)

    n = size(X, 1)
    p = size(X, 2)

    if p ==0
        σ̂² = mean(y.^2)
        H0 = zeros(n, n)
    else
        fit = lm(X, y)
        σ̂² = sum(map(x -> x^2, residuals(fit))) / (n - p)
        H0 = (g/(g+1)) * X * inv(X'X) * X'
    end
    SS0 = y' * (diagm(ones(n)) - H0) * y
    return(
        -0.5*n*log(2*π) + logabsgamma(0.5*(ν₀+n))[1] - logabsgamma(0.5*ν₀)[1] - 0.5*p*log(1+g) +
            0.5*ν₀*log(0.5*ν₀*σ̂²) - 0.5*(ν₀+n)*log(0.5*(ν₀*σ̂²+SS0))
    )
end

function BayesianModelAveraging(y, X; S=5000)
    n, p = size(X)
    BETA = Matrix{Float64}(undef, S, p)
    Z = Matrix{Int64}(undef, S, p)
    z = ones(p)
    lpy_c = lpyX(y, X[:, findall(z .== 1)])

    for s in 1:S
        for j in shuffle(1:p)
            zp = copy(z)
            zp[j] = 1 .- zp[j]
            lpy_p = lpyX(y, X[:, findall(zp .== 1)])
            r = (lpy_p - lpy_c)*(-1)^(zp[j]==0)
            z[j] = rand(Bernoulli(1/(1+exp(-r)))) |> Int
            if z[j] == zp[j]
                lpy_c = lpy_p
            end
        end

        β = copy(z)
        if sum(z) > 0
            β[findall(z .== 1)] = lmGprior(y, X[:, findall(z .== 1)], S=1)[1]
        end
        Z[s, :] = z
        BETA[s, :] = β

        if s % 10 == 0
            bpm = mean(BETA[1:s, :], dims=1) |> vec
            mse = mean((y - X * bpm).^2)
            println("s = ", s, ", mse = ", mse)
        end
    end
    return BETA, Z
end

BETA2, Z = BayesianModelAveraging(y, X; S=5000)
z_prob = mean(Z, dims=1) |> vec
ci_β_bma = map(x -> quantile(x, [0.025, 0.975]), eachcol(BETA2))
df_ci_bma = DataFrame(
    param = coef_name,
    prob = z_prob,
    ci_lower = map(x -> x[1], ci_β_bma),
    ci_upper = map(x -> x[2], ci_β_bma)
)