はじめに
Facebookグループ『東大の統計学教科書を読み進める会』で、
統計学入門（通称：赤本）の読み進めを行っています。
赤本の統計処理をJuliaでどう実装するかが気になったので、コードを書き連ねています。
統計学入門（通称：赤本）とは？
「東京大学教養学部統計学教室 編」の統計学入門書です。古典です。
参考：統計学入門
Juliaとは？
科学計算向けのプログラミング言語です。ポストR・ベターMATLABと呼ばれています。
参考：Juliaの公式サイト
調査対象の章
第3章「2次元のデータ」
コード
sec03.jl
blood_mtrx = [35 114; 45 124; 55 143; 65 158; 75 166;]
# [散布図 scattergram]
Pkg.add("Gadfly")
using Gadfly
plot(x=blood_mtrx[:, 1], y=blood_mtrx[:, 2], Guide.XLabel("年齢"), Guide.YLabel("血圧"), Geom.point)
# [分割表 contigency table]
## カウント
count(x -> 150 < x, blood_mtrx[:, 2])
## 合計
sum(blood_mtrx, 1) #縦計
sum(blood_mtrx, 2) #横計
# [共分散 covariance]
cov(blood_mtrx[:, 1], blood_mtrx[:, 2], corrected = false)
#sum([(blood_mtrx[i, 1] - mean(blood_mtrx[:, 1])) * (blood_mtrx[i, 2] - mean(blood_mtrx[:, 2])) for i = 1:size(blood_mtrx, 1)]) / size(blood_mtrx,1)
# [ピアソンの積率相関係数 product-moment correlation coefficient]
cor(blood_mtrx[:, 1], blood_mtrx[:, 2])
#sum([(blood_mtrx[i, 1] - mean(blood_mtrx[:, 1])) * (blood_mtrx[i, 2] - mean(blood_mtrx[:, 2])) for i = 1:size(blood_mtrx, 1)]) / size(blood_mtrx,1) / (sqrt(sum([(x - mean(blood_mtrx[:, 1])) ^ 2 for x = blood_mtrx[:, 1]]) / size(blood_mtrx,1)) * sqrt(sum([(y - mean(blood_mtrx[:, 2])) ^ 2 for y = blood_mtrx[:, 2]]) / size(blood_mtrx,1)))
# [偏相関係数 partial correlation coefficient]
Pkg.add("RDatasets")
using RDatasets
iris = dataset("datasets", "iris")
pcor(x, y, z) = (cor(x, y) - (cor(x, z) * cor(y, z))) / (sqrt(1 - (cor(x, z) ^ 2)) * sqrt(1 - (cor(y, z) ^ 2)))
pcor(iris[:, :SepalLength], iris[:, :PetalWidth], iris[:, :PetalLength])
# [並び替え sort]
sort(blood_mtrx[:, 2])
blood_mtrx[:, 2][sortperm(blood_mtrx[:, 2])]
# [順位相関係数 rank correlation coefficient]
Pkg.add("DataFrames")
using DataFrames
flower_mtrx = [1 3; 2 1; 3 2; 4 5; 5 4; 6 7; 7 6; 8 8;]
## スピアマン Spearman
Pkg.add("StatsBase")
using StatsBase
corspearman(flower_mtrx[:, 1], flower_mtrx[:, 2])
#1 - ((6 / ((size(flower_mtrx, 1) ^ 3) - size(flower_mtrx, 1))) * sum([(flower_mtrx[i, 1] - flower_mtrx[i, 2]) ^ 2 for i = 1:size(flower_mtrx, 1)]))
## ケンドール Kendall
using StatsBase
corkendall(flower_mtrx[:, 1], flower_mtrx[:, 2])
#flower_mtrx |>
# x -> [(x[i, 1] < x[i + 1, 1]) == (x[i, 2] < x[i + 1, 2]) for i = 1:(size(x, 1) - 1)] |>
# x -> (countnz(x) + count(x1 -> !x1, x)) / (((length(x) + 1) * length(x)) / 2)
# [時系列データのプロット]
using RDatasets
arbuthnot = dataset("HistData", "Arbuthnot")
using Gadfly
plot(arbuthnot, x = :Year, y = :Mortality, Geom.line)
# [自己相関係数 auto correlation coefficient]
using StatsBase
autocor(arbuthnot[:, :Mortality], 1)
#arbuthnot[:, :Mortality] |>
# ts -> sum([(ts[i] - mean(ts)) * (ts[i + 1] - mean(ts)) / (length(ts) - 1) for i = 1:(size(ts, 1) - 1)]) / sum([(x - mean(ts)) ^ 2 / length(ts) for x = ts])
# [コレログラム correlogram]
using Gadfly
plot(x = 1:(size(arbuthnot, 1) - 1), y = [autocor(arbuthnot[:, :Mortality], i) for i = 1:(size(arbuthnot, 1) - 1)], Geom.line)
# [線形回帰 liner regression]
using RDatasets
iris = dataset("datasets", "iris")
sepal_length = convert(Array{Float64,1}, iris[:, :SepalLength])
petal_length = convert(Array{Float64,1}, iris[:, :PetalLength])
a, b = linreg(sepal_length, petal_length) #y = a + b * x
# [線形回帰のプロット]
plot(
layer(x=sepal_length, y=petal_length, Geom.point),
layer(x=ifloor(minimum(sepal_length)):iceil(maximum(sepal_length)), y=[a + b * x for x = ifloor(minimum(sepal_length)):iceil(maximum(sepal_length))], Geom.line)
)
# [多項式回帰 polynomial regression]
Pkg.add("GLM")
using GLM
using RDatasets
LifeCycleSavings = dataset("datasets", "LifeCycleSavings")
fm2 = fit(LinearModel, SR ~ Pop15 + Pop75 + DPI + DDPI, LifeCycleSavings)
# [練習問題]
# 3.1 社会経済指標と投票行動
vhrate_mtrx = {
41.4 52.8;
76.3 71.2;
59.2 72.6;
51.8 63.7;
52.5 81.3;
53.2 81.8;
62.4 70.9;
55.0 74.0;
57.7 73.2;
63.2 72.9;
37.5 66.7;
48.5 65.7;
32.4 43.7;
20.5 55.5;
47.9 79.6;
68.9 85.7;
68.5 75.3;
52.5 80.5;
63.3 73.0;
58.8 77.0;
59.7 77.5;
48.4 69.2;
40.7 60.0;
51.0 78.2;
50.9 79.5;
34.3 61.8;
25.8 49.6;
32.1 59.6;
34.4 72.1;
55.1 71.0;
60.3 76.3;
57.0 72.8;
45.6 71.8;
54.2 60.7;
55.1 67.0;
55.7 71.8;
70.3 71.2;
61.8 68.3;
47.6 68.5;
42.5 54.8;
71.3 76.0;
55.2 65.8;
65.2 69.4;
42.9 66.9;
54.7 69.7;
62.0 71.2;
48.2 59.6;
}
plot(x = vhrate_mtrx[:, 1], y = vhrate_mtrx[:, 2], Guide.XLabel("自民得票率"), Guide.YLabel("持ち家比率"), Geom.point)
cor(vhrate_mtrx[:, 1], vhrate_mtrx[:, 2])
# 3.2 統計的な関連
# (略)
# 3.3 社会的リスクの順位づけ
riskeval_mtrx = [
1 1 8 20;
2 5 3 1;
3 2 1 4;
4 3 4 2;
5 6 2 6;
6 7 5 3;
7 15 11 12;
8 8 7 17;
9 4 15 8;
10 11 9 5;
11 10 6 18;
12 14 13 13;
13 18 10 23;
14 13 22 26;
15 22 12 29;
16 24 14 15;
17 16 18 16;
18 19 19 9;
19 30 17 10;
20 9 22 11;
21 25 16 30;
22 17 24 7;
23 26 21 27;
24 23 20 19;
25 12 28 14;
26 20 30 21;
27 28 25 28;
28 21 26 24;
29 27 27 22;
30 29 29 25;
]
## スピアマン Spearman
using StatsBase
corspearman(riskeval_mtrx[:, 1], riskeval_mtrx[:, 2])
## ケンドール Kendall
using StatsBase
corkendall(riskeval_mtrx[:, 1], riskeval_mtrx[:, 2])
## 3.4 ブートストラップ
# ⅰ)
iceil(rand()*11)
## ⅱ)
height_mtrx = [
71 69;
68 64;
66 65;
67 63;
70 65;
71 62;
70 65;
73 64;
72 66;
65 59;
66 62;
]
bootstrap(mtrx, cnt) =
reduce(vcat, [mtrx[iceil(rand()*11), :] for i = 1:cnt]) |>
boot_mtrx -> cor(boot_mtrx[:, 1], boot_mtrx[:, 2])
bootstrap(height_mtrx, 11)
## ⅲ)
plot(x=[bootstrap(height_mtrx, i) for i = 1:200], Geom.histogram)
感想
GLM.jlのfit関数でエラーが発生する。
GithubのReadme.mdのサンプルコードを実行してエラーになっているので、
v0.3.0に伴う言語の仕様変更が原因だと思える。
v1.0.0が待たれる。
[2014/9/12]
Glm.jlをv0.4.2に更新した所、fit関数でエラーが起こらなくなった。