Julia に翻訳--148 重回帰分析の標準化残差

このページのランキングは？？

 

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

重回帰分析の標準化残差
http://aoki2.si.gunma-u.ac.jp/R/stdres.html

ファイル名: stdres.jl　　関数名: stdres

翻訳するときに書いたメモ

h の計算はもう少し簡単になるかもしれない

==========#

using Statistics, LinearAlgebra

function stdres(x, y)
    nr, nc = size(x)
    fittedvalues, residuals = mreg(hcat(x, y))
    x = hcat(ones(size(x, 1)), x)
    s2 = sum(residuals .^ 2) / (nr - nc - 1)
    h = vec(sqrt.(s2 .* (1 .- sum((x'*x \ x') .* x', dims=1))))
    hcat(y, fittedvalues, residuals, residuals ./ h)
end

function mreg(dat)  # 回帰分析を行い，予測値，残差を返す関数
    n, nc = size(dat)
    r = cor(dat)
    betas = r[1:end-1, 1:end-1] \ r[1:end-1, end]
    prop = var(dat, dims=1, corrected=false) .* n
    b = betas .* sqrt.(prop[end] ./ prop[1:end-1])
    means = mean(dat, dims=1)
    b0 = means[end] - sum(b .* means[1:end-1])
    fittedvalues = dat[:, 1:end-1] * b .+ b0
    residuals = dat[:, end] .- fittedvalues
    fittedvalues, residuals
end

x = [
    1 6 4
    3 2 5
    4 3 3
    2 1 3
    5 4 7
    3 4 2];

y = [4, 2, 6, 5, 4, 3];

stdres(x, y)
# 6×4 Array{Float64,2}:
#  4.0  3.46951   0.530488   1.17902
#  2.0  3.69512  -1.69512   -1.06417
#  6.0  4.56707   1.43293    0.943228
#  5.0  3.9939    1.0061     0.870388
#  4.0  3.68293   0.317073   0.393939
#  3.0  4.59146  -1.59146   -1.17902