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

  • 26
    いいね
  • 1
    コメント
この記事は最終更新日から1年以上が経過しています。

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

この投稿は Julia Advent Calendar 201519日目の記事です。