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

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](https://qiita-image-store.s3.amazonaws.com/0/50670/56d34aad-f885-f654-93ef-cdc7b82d0ba8.png)


# Jupyter (iPython Notebook)での実行

※ Jupyter (iPython Notebook)はすでに導入済みと想定します。[Anaconda](https://www.continuum.io/downloads)等で事前にインストールしてください。

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

```jl
Pkg.add("IJulia")

```

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

```log: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: 
```

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

```jl
Pkg.build("IJulia")
```

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

```log: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`にしておくと、ブラウザで自動起動します。

```jl
using IJulia
notebook(detached=true)
```

<img width="799" alt="スクリーンショット 2015-11-10 20.21.36.png" src="https://qiita-image-store.s3.amazonaws.com/0/50670/72caa573-6bd1-8278-ea15-57e091677b6c.png">


できた！！！ :grin::grin::grin: 

### ちょっと動かしてみる


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


<img width="793" alt="スクリーンショット 2015-11-10 21.09.28.png" src="https://qiita-image-store.s3.amazonaws.com/0/50670/6c8a9494-690b-7dd9-2a4f-dea291a156a7.png">

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

<img width="775" alt="スクリーンショット 2015-11-12 01.13.46.png" src="https://qiita-image-store.s3.amazonaws.com/0/50670/7bcc18ec-7812-85ac-8af5-f10fa00ca567.png">


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


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

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

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

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

```jl
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|
|…|…|…|…|…|…|

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

```jl
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データを試す

```jl
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" |

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

<img width="919" alt="スクリーンショット 2015-11-12 01.17.13.png" src="https://qiita-image-store.s3.amazonaws.com/0/50670/2b3d77c8-281e-db10-da15-1a4aa5db6662.png">

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


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

```jl
using GLM
```

## 直線回帰

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

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

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


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

<img width="675" alt="reg_plot.png" src="https://qiita-image-store.s3.amazonaws.com/0/50670/faaa5d37-ce14-326f-9d52-b14cbf6e041d.png">

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

```jl
# 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)"))
```

<img width="670" alt="reg_plot_log.png" src="https://qiita-image-store.s3.amazonaws.com/0/50670/fb688ff4-5966-6d4d-9ca6-e7053c5fe9d5.png">


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

```jl
# 線形回帰の実行
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

```

```jl
# グラフに直線を描画
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)")
)
```

<img width="692" alt="reg_plot_log_line.png" src="https://qiita-image-store.s3.amazonaws.com/0/50670/99ef4dfd-2380-48d9-1cd4-7131aabfed19.png">

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

## ポアソン回帰

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

```jl
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に出来たので、グラフ描画とポアソン回帰を実行します。

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

<img width="669" alt="pois_plot.png" src="https://qiita-image-store.s3.amazonaws.com/0/50670/d4c43e18-3cd2-2efb-f551-cb43ff5fa755.png">

```jl
# ポアソン回帰の実行
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
```

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


```jl
# 散布図の上のポアソン回帰の結果の曲線を載せてみる
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)")
)
```

<img width="666" alt="pois_plot_curve.png" src="https://qiita-image-store.s3.amazonaws.com/0/50670/99d87486-536e-c25d-6faa-b7558648c56d.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