Edited at
JuliaDay 19

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

More than 3 years have passed since last update.

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

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

ちなみに、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が実行できるようになりました。簡単ですね :laughing:

julia-REPL-compressor.png


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


できた! :grin:

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

using IJulia

notebook(detached=true)

スクリーンショット 2015-11-10 20.21.36.png

できた!!! :grin::grin::grin:


ちょっと動かしてみる

========

# この辺のパッケージインストールは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)

スクリーンショット 2015-11-10 21.09.28.png

# 折れ線グラフ

plot(x=rand(50), y=rand(50), Geom.point, Geom.line)

スクリーンショット 2015-11-12 01.13.46.png

できた!!!!! :grin::grin::grin::grin::grin:

しかも、実際のグラフは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のデータが使えます :kissing_closed_eyes:

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)

スクリーンショット 2015-11-12 01.17.13.png

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


回帰分析と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)

reg_plot.png

散布図をみると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)"))

reg_plot_log.png

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

# 線形回帰の実行

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)")
)

reg_plot_log_line.png

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


ポアソン回帰

ポアソン回帰用のデータを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)

pois_plot.png

# ポアソン回帰の実行

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)")
)

pois_plot_curve.png


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

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