#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。
Excel にある二変量統計関数
http://aoki2.si.gunma-u.ac.jp/R/bivariate.html
ファイル名: bivariate.jl 関数名: correl, covar, forcast, growth,
intercept, logest, pearson,
rsq, slope, steyx, trend
翻訳するときに書いたメモ
Excel に用意する関数は，あまり網羅的でもないんだなあ。
==========#
using Statistics
x = [3, 5, 4, 8, 2, 1];
y = [2, 4, 5, 7, 1, 6];
correl(x, y) = cor(x, y)
correl(x, y) # 0.492515750317586
covar(x, y) = cov(x, y, corrected=false)
covar(x, y) # 2.3611111111111107
forecast(data, y, x) = mean(y) .+ cov(x, y)/var(x) .* (data .- mean(x))
forecast(6, y, x) # 5.162162162162162
forecast(8, y, x) # 6.081081081081081
forecast([6, 8, 10], y, x)
# 3-element Array{Float64,1}:
# 5.162162162162162
# 6.081081081081081
# 7.0
function growth(y, x, data, one = false)
all(y .> 0) || error("any y is less than or equal to zero.")
y = log.(y)
if one
b = sum(x .* y) / sum(x .^ 2)
constant = 0
else
b = cov(x, y) / var(x)
constant = mean(y) - b * mean(x)
end
exp.(constant .+ b .* data)
end
growth(y, x, 6) # 4.677106652484691
growth(y, x, [1, 2, 3])
# 3-element Array{Float64,1}:
# 2.3140857250280833
# 2.6637840201891647
# 3.0663277636912243
growth(y, x, 6, true) # 5.228738441530447
growth(y, x, [1, 2, 3], true)
# 3-element Array{Float64,1}:
# 1.3174459890112413
# 1.7356639339618076
# 2.2866434880694557
intercept(y, x) = mean(y) - cov(x, y) / var(x) * mean(x)
intercept(y, x) # 2.4054054054054057
function logest(y, x, one = false)
all(y .> 0) || error("any y is less than or equal to zero.")
y = log.(y)
if one
b = sum(x .* y) / sum(x .^ 2)
constant = 0
else
b = cov(x, y) / var(x)
constant = mean(y) - b * mean(x)
end
Dict(:model => "a*b^x", :a => exp(constant), :b => exp(b))
end
logest(y, x)
# Dict{Symbol,Any} with 3 entries:
# :a => 2.0103
# :b => 1.15112
# :model => "a*b^x"
logest(y, x, true)
# Dict{Symbol,Any} with 3 entries:
# :a => 1.0
# :b => 1.31745
# :model => "a*b^x"
pearson(x, y) = cor(x, y)
pearson(x, y) # 0.492515750317586
rsq(y, x) = cor(y, x)^2
rsq(y, x) # 0.2425717643108947
slope(y, x, zero = false) = zero ? sum(x .* y) / sum(x .^ 2) : cov(x, y) / var(x)
slope(y, x) # 0.45945945945945943
slope(y, x, true) # 0.9243697478991597
function steyx(y, x)
n = length(x)
sqrt((var(y) - cov(x, y) ^ 2 / var(x)) * (n - 1) / (n - 2))
end
steyx(y, x) # 2.254125347242491
trend(y, x, data, zero = false) = zero ? sum(x .* y) / sum(x .^ 2) .* data :
intercept(y, x) .+ (cov(x, y) / var(x)) .* data
trend(y, x, 6) # 5.162162162162162
trend(y, x, [4, 5, 6])
# 3-element Array{Float64,1}:
# 4.243243243243244
# 4.7027027027027035
# 5.162162162162162
trend(y, x, 6, true) # 5.546218487394958
trend(y, x, [4, 5, 6], true)
# 3-element Array{Float64,1}:
# 3.697478991596639
# 4.621848739495799
# 5.546218487394958
