JupyterでJuliaを動かして回帰分析をやってみる。

Julia上でデータセットを取得して回帰分析(線形回帰、ポアソン回帰)までやってみるという記事です。

MacにJuliaを導入して、Jupyter Notebookで動かすという最初の最初から書いていきます。なぜなら僕がJulia初心者だからです Jupyterの"Ju"はJuliaの"Ju"らしいので、ぜひJupyterを使いましょう

ちなみに、Juliaの本を探そうとAmazonで"Julia"を検索すると痛い目にあうので（特に職場では）気をつけましょうw (参考： https://twitter.com/Kenmatsu4/status/663990102478028800)

##環境

• OSX Yosemite 10.10.5

コードの全文

本記事で扱ったコードの全文はGitHub
https://github.com/matsuken92/Qiita_Contents/blob/master/Julia/Regression_with_Julia.ipynb
にアップしています。(一部の画像は壊れてしまっているので、そちらについては本記事で画像をご確認ください. D3.jsのグラフをGitHub上で見れるよう画像にコンバートしている部分がうまくいっていないようです。）

Juliaのインストール

まずは公式ダウンロードページへ。
http://julialang.org/downloads/

Current Release (v0.4.2)が最新のバージョンですね(2015/12/16現在)
Mac用の
　　Mac OS X Package (.dmg)
からイメージをダウンロードします。

julia-0.4.2-osx10.7+.dmg をダブルクリックして、開かれるディレクトリにある Julia-0.4.1.app をアプリケーションフォルダにコピーして完了です。これを起動すると、ターミナルが開いて、Juliaが実行できるようになりました。簡単ですね

Jupyter (iPython Notebook)での実行

※ Jupyter (iPython Notebook)はすでに導入済みと想定します。Anaconda等で事前にインストールしてください。

先ほど立ち上がったJuliaのターミナルにてインストールを実施します。

Pkg.add("IJulia")

ありゃ、エラーが出ました。

out

===============================[ ERROR: IJulia ]================================

LoadError: __precompile__(true) but require failed to create a precompiled cache file
while loading /Users/XXXXXXXX/.julia/v0.4/IJulia/deps/build.jl, in expression starting on line 2

================================================================================

=============================================[ BUILD ERRORS ]==============================================

WARNING: IJulia had build errors.

 - packages with build errors remain installed in /Users/XXXXXXXX/.julia/v0.4
 - build the package(s) and all dependencies with `Pkg.build("IJulia")`
 - build a single package by running its `deps/build.jl` script

===========================================================================================================
INFO: Package database updated :grin: 

依存関係含めビルドするには、

Pkg.build("IJulia")

を実行せよということらしいので、やってみる。

out

INFO: Building Homebrew
HEAD is now at d872c99 tap: warn already tapped
HEAD is now at b178725 Merge pull request #83 from staticfloat/staging
INFO: Building Nettle
INFO: Building ZMQ
INFO: Building IJulia
INFO: Found Jupyter version 3.2.0: ipython
Writing IJulia kernelspec to /Users/XXXXXXXX/.julia/v0.4/IJulia/deps/julia-0.4/kernel.json ...
Installing julia kernelspec julia-0.4

できた！

では、JupyterでJuliaを起動してみましょう。detached=trueにしておくと、ブラウザで自動起動します。

using IJulia
notebook(detached=true)

できた！！！

ちょっと動かしてみる

========
# この辺のパッケージインストールはREPLからやった方がいいかも。
Pkg.update()
Pkg.status()

Pkg.add("Distributions")    # 乱数生成や、確率分布を使う時用
Pkg.add("Gadfly")           # グラフ描画用
Pkg.add("Cairo")            # グラフ描画用
Pkg.add("RDatasets")        # Rのデータセットを使用するため
Pkg.add("GLM")              # 回帰の実行用
Pkg.add("Requests")         # HTTP通信用

=======
# パッケージのインポート
# 初回のみコンパイルが実行されるので時間がかかる。
using Distributions
using Gadfly
using DataFrames

======

# 正規分布オブジェクトの作成
d = Normal()

# 標準正規分布に従う100個の乱数を2セット生成する
x = rand(d, 100)
y = rand(d, 100)

set_default_plot_size(7inch, 7inch/golden);

# 散布図の描画
plot(x=x, y=y)

# 折れ線グラフ
plot(x=rand(50), y=rand(50), Geom.point, Geom.line)

できた！！！！！
しかも、実際のグラフはD3.jsなので、インタラクティブに動かせます！

Rのデータセットを使ってみる

まずRにあるデータセットを取得できるRDatasetsというパッケージがあるのでこれを導入します。

Pkg.add("RDatasets")        # RのデータをJuliaで使いたい時用
using RDatasets

それで、下記のdatasetsというメソッドを呼ぶと使えるデータ一式を見ることができます。

ds = RDatasets.datasets()

	Package	Dataset	Title	Rows	Columns
1	COUNT	affairs	affairs	601	18
2	COUNT	azdrg112	azdrg112	1798	4
3	COUNT	azpro	azpro	3589	6
4	COUNT	badhealth	badhealth	1127	3
5	COUNT	fasttrakg	fasttrakg	15	9
6	COUNT	lbw	lbw	189	10
7	COUNT	lbwgrp	lbwgrp	6	7
…	…	…	…	…	…

こんな感じでデータの名称とデータ数を見ることができます。ただ、これだと途中で切れてしまうので、下記のように全てファイルに書き出してみるとどのようなデータがあるかを確認できます。

open( "datalist.txt", "w" ) do fp
    write(fp, "Package\tDataset\tTitle\tRows\tColumns\n")
    for i in 1:length(ds[2]);
        for j in 1:5;
            tmp = ds[j][i]
            write(fp, "$tmp\t");
        end
        write(fp, "\n");
    end
end

Irisデータを試す

iris = dataset("datasets", "iris")
println(iris)

おなじみ、みんな大好きirisのデータが使えます

Row	SepalLength	SepalWidth	PetalLength	PetalWidth	Species
1	5.1	3.5	1.4	0.2	"setosa"
2	4.9	3.0	1.4	0.2	"setosa"
3	4.7	3.2	1.3	0.2	"setosa"
4	4.6	3.1	1.5	0.2	"setosa"
5	5.0	3.6	1.4	0.2	"setosa"
…	…	…	…	…	…
145	6.7	3.3	5.7	2.5	"virginica"
146	6.7	3.0	5.2	2.3	"virginica"
147	6.3	2.5	5.0	1.9	"virginica"
148	6.5	3.0	5.2	2.0	"virginica"
149	6.2	3.4	5.4	2.3	"virginica"
150	5.9	3.0	5.1	1.8	"virginica"

plot(dataset("datasets", "iris"), x="SepalLength", y="SepalWidth", color="Species",Geom.point)

なかなか綺麗なグラフです。

回帰分析とGLM(ポアソン回帰)

GLMパッケージを使って、直線回帰とポアソン回帰を行ってみます。まずはパッケージのインポートから。

using GLM

直線回帰

housing = dataset("Ecdat", "Housing")

列はいくつもあるデータですが、今回はシンプルに部屋のサイズだけから価格を推定します。

	Price	LotSize
1	42000	5850
2	38500	4000
3	49500	3060
4	60500	6650
…	…	…

# まずは散布図のプロット
plot(housing, x="LotSize", y="Price", Geom.point)

散布図をみるとlogを取ってみた方がよさそうなのでコンバートする。

# logをとってみた方がよさそう
log_housing = DataFrame(LotSize=log(housing[:,2]), Price=log(housing[:,1]))

#プロットしてみる
plot(log_housing, x="LotSize", y="Price", Geom.point,Guide.xlabel("LotSize(log)"), Guide.ylabel("Price(log)"))

このデータに対して線形回帰を行います。

# 線形回帰の実行
lm = fit(LinearModel, Price ~ LotSize, log_housing)

out

DataFrames.DataFrameRegressionModel{GLM.LinearModel{GLM.DensePredQR{Float64}},Float64}:

Coefficients:
             Estimate Std.Error t value Pr(>|t|)
(Intercept)   6.46853  0.276741  23.374   <1e-83
LotSize      0.542179 0.0326501 16.6057   <1e-49

# グラフに直線を描画
intercept = coef(lm)[1]
slope     = coef(lm)[2]

xx = collect(7:0.1:10)
yy = intercept + xx*slope 

plot(
     layer(log_housing, x="LotSize", y="Price", Geom.point,
           order=1,Theme(default_color=color(colorant"#7799ff"))),
     layer(x=xx, y=yy, Geom.line, order=2,Theme(default_color=color(colorant"green"), 
           line_width=3px)),
     Guide.Title("Regression of Housing price"),
     Guide.xlabel("LotSize(log)"), Guide.ylabel("Price(log)")
)

いい感じで線が引けました

ポアソン回帰

ポアソン回帰用のデータをweb上から取ってきて加工、DataFrameにしてポアソン回帰を実行します。

using Requests

# Children Ever Born Data
# Lecture Note: http://data.princeton.edu/wws509/notes/c4.pdf
# Dataset: http://data.princeton.edu/wws509/datasets/#ceb
url = "http://data.princeton.edu/wws509/datasets/ceb.dat"
res = get(url)

# DataFrame化できるように加工する
ceb = readcsv(IOBuffer(res.data));
n_col = length(split(ceb[2]))
n_row = length(ceb)

ceb_list = []
header = true
for i in 1:n_col;
    if i == 7 || i == 8;
        # ちゃんとここでArrayの型をAnyではなく、Float64にしておかないとglmの実行時にエラーになる
        push!(ceb_list, Float64[]) 
    else;
        push!(ceb_list, []) 
    end
end

col_name = Symbol[]
for data in ceb;
    if header == true
        header = false 
        
        # ヘッダーリストを作成
        names = split("id " * data)
        for n in 1:length(names);
            push!(col_name, symbol(names[n]))
        end
        
        continue
    end
    
    tmp = split(data)
    tmp[2] = split(tmp[2], "-")[2]
    for i in 1:n_col;
        if i == 7 || i == 8;
            push!(ceb_list[i],round(parse(Float64, tmp[i])))
        else;
            push!(ceb_list[i], tmp[i])        
        end
    end
end

# DataFrame化
df_ceb = DataFrame(ceb_list, Vector(col_name))[:,2:8]

さて、データを無事DataFrameに出来たので、グラフ描画とポアソン回帰を実行します。

# まずは散布図のプロット
plot(df_ceb, x="n", y="y", Geom.point)

# ポアソン回帰の実行
gm1 = fit(GeneralizedLinearModel, y ~ n, df_ceb, Poisson(), LogLink())

out

DataFrames.DataFrameRegressionModel{GLM.GeneralizedLinearModel{GLM.GlmResp{Array{Float64,1},Distributions.Poisson,GLM.LogLink},GLM.DensePredChol{Float64,Base.LinAlg.Cholesky{Float64,Array{Float64,2}}}},Float64}:

Coefficients:
              Estimate   Std.Error z value Pr(>|z|)
(Intercept)    4.03177   0.0160541 251.137   <1e-99
n            0.0175205 0.000155328 112.797   <1e-99

結果も良好で有意になっています。これを使って曲線を当てはめてみます。

# 散布図の上のポアソン回帰の結果の曲線を載せてみる
xx = collect(0:1:200)
yy = exp(4.03177 + 0.0175205 * xx)
plot(layer(df_ceb, x="n", y="y", Geom.point, order=1,Theme(default_color=color(colorant"#7799ff"))),
     layer(x=xx, y=yy, Geom.line, order=2,Theme(default_color=color(colorant"green"), line_width=3px) ),
     Guide.Title("Regression of Housing price"),
     Guide.xlabel("LotSize(log)"), Guide.ylabel("Price(log)")
)

参考 (いろいろ調べた履歴です。ここが一番有用かも？)

X分で学ぶJulia -りんごがでている-
　http://bicycle1885.hatenablog.com/entry/2014/12/01/050522

Mathematical Operations and Elementary Functions
　http://docs.julialang.org/en/release-0.4/manual/mathematical-operations/

DataFrames.jl Overview
　https://dataframesjl.readthedocs.org/en/latest/

RDatasetsパッケージ
　https://github.com/johnmyleswhite/RDatasets.jl#data-sets

Gadflyパッケージ
　http://dcjones.github.io/Gadfly.jl/

Gadfly Examples
　https://github.com/dcjones/Gadfly.jl/tree/master/doc

Introducing Julia/Plotting
　https://en.wikibooks.org/wiki/Introducing_Julia/Plotting#Gadfly

文字列のSplit
　https://en.wikibooks.org/wiki/Introducing_Julia/Strings_and_characters#Splitting_and_joining_strings

配列Arrayの使い方
　http://docs.julialang.org/en/release-0.4/manual/arrays/

各種統計パッケージ
　https://github.com/JuliaStats

Julia for R Programmers
　http://www.stat.wisc.edu/~bates/JuliaForRProgrammers.pdf