Help us understand the problem. What is going on with this article?

Juliaで学ぶ統計学入門(第3章:2次元のデータ)

More than 5 years have passed since last update.

はじめに

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関数でエラーが起こらなくなった。

nezuq
ITエンジニア。仕事や趣味で学んだ事を共有しています。Enjoy Tech! Update World!
https://www.slideshare.net/nezuQ/presentations
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away