[2019年2月版] Julia + Jupyterで回帰分析してみる

# 元ページ
@kenmatsu4 さんの以下の記事を参考に、最新版のJulia 1.1で回帰分析&グラフ描画の動作を再現してみました。

- [JupyterでJuliaを動かして回帰分析をやってみる。](https://qiita.com/kenmatsu4/items/b77ecec75418cbaa2212)

# 準備
REPLで必要なパッケージのインストールを行います。

まず例によってREPLコンソールを立ち上げます。

```shell
$ julia
               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.1.0 (2019-01-21)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia> 
```

1. ```using```で```Pkg```を有効にし、以下のように各種パッケージをインストールします。
2. インストールしたパッケージはすぐさま```using```して初回のコンパイルを実施しておきます。

```julia
julia> using Pkg

julia> Pkg.update()

julia> Pkg.add("Distributions")

julia> using Distributions

julia> Pkg.add("Gadfly")

julia> using Gadfly

julia> Pkg.add("RDatasets")

julia> using RDatasets

julia> Pkg.add("GLM")

julia> using GLM

julia> Pkg.add("DataFrames")

julia> using DataFrames

julia> Pkg.add("HTTP")

julia> using HTTP

julia> Pkg.add("CSV")

julia> using CSV

```
## インストールに失敗したら
実は```Cairo```のインストール中にネットワークの調子が悪く、```homebrew```をインストール？中に以下のエラーで中断してしまいました。

```shell
julia> Pkg.add("Cairo")
 Resolving package versions...
 Installed Graphics ──────── v0.4.0
 Installed WinRPM ────────── v0.4.2
 Installed Homebrew ──────── v0.7.1
 Installed LibExpat ──────── v0.5.0
 Installed Libz ──────────── v1.0.0
 Installed HTTPClient ────── v0.2.1
 Installed BufferedStreams ─ v1.0.0
 Installed LibCURL ───────── v0.4.1
 Installed Cairo ─────────── v0.5.6
  Updating `~/.julia/environments/v1.1/Project.toml`
  [159f3aea] + Cairo v0.5.6
  Updating `~/.julia/environments/v1.1/Manifest.toml`
  [e1450e63] + BufferedStreams v1.0.0
  [159f3aea] + Cairo v0.5.6
  [a2bd30eb] + Graphics v0.4.0
  [0862f596] + HTTPClient v0.2.1
  [d9be37ee] + Homebrew v0.7.1
  [b27032c2] + LibCURL v0.4.1
  [522f3ed2] + LibExpat v0.5.0
  [2ec943e9] + Libz v1.0.0
  [c17dfb99] + WinRPM v0.4.2
  Building LibCURL ─→ `~/.julia/packages/LibCURL/OoXMv/deps/build.log`
  Building WinRPM ──→ `~/.julia/packages/WinRPM/Y9QdZ/deps/build.log`
  Building Homebrew → `~/.julia/packages/Homebrew/s09IX/deps/build.log`
^[


┌ Error: Error building `Homebrew`: 
         -=O=-                                                           # # #  
│ Initialized empty Git repository in /Users/yoshinori/.julia/packages/Homebrew/s09IX/deps/usr/.git/
│ Updating Homebrew...
│ From https://github.com/Homebrew/brew
│  * [new branch]      dependabot/bundler/Library/Homebrew/test/rubocop-0.64.0 -> origin/dependabot/bundler/Library/Homebrew/test/rubocop-0.64.0
│  * [new branch]      dependabot/bundler/Library/Homebrew/vendor/rubocop-0.64.0 -> origin/dependabot/bundler/Library/Homebrew/vendor/rubocop-0.64.0
│  * [new branch]      master     -> origin/master
│  * [new tag]         1.9.0      -> 1.9.0
│  * [new tag]         1.9.1      -> 1.9.1
│  * [new tag]         1.9.2      -> 1.9.2
│  * [new tag]         1.9.3      -> 1.9.3
│  * [new tag]         2.0.0      -> 2.0.0
│  * [new tag]         2.0.1      -> 2.0.1
│ HEAD is now at 6ae600b Merge pull request #5689 from Homebrew/dependabot/bundler/Library/Homebrew/test/parallel_tests-2.28.0
│ ==> Tapping homebrew/core
│ Cloning into '/Users/yoshinori/.julia/packages/Homebrew/s09IX/deps/usr/Library/Taps/homebrew/homebrew-core'...
│ error: RPC failed; curl 56 LibreSSL SSL_read: SSL_ERROR_SYSCALL, errno 60
│ fatal: The remote end hung up unexpectedly
│ fatal: early EOF
│ fatal: index-pack failed
│ Error: Failure while executing; `git clone https://github.com/Homebrew/homebrew-core /Users/yoshinori/.julia/packages/Homebrew/s09IX/deps/usr/Library/Taps/homebrew/homebrew-core` exited with 128.
│ [ Info: Downloading brew...
│ ERROR: LoadError: InitError: failed process: Process(`/Users/yoshinori/.julia/packages/Homebrew/s09IX/deps/usr/bin/brew tap --full homebrew/core`, ProcessExited(1)) [1]
│ Stacktrace:
│  [1] error(::String, ::Base.Process, ::String, ::Int64, ::String) at ./error.jl:42
│  [2] pipeline_error at ./process.jl:785 [inlined]
│  [3] #run#515(::Bool, ::Function, ::Cmd) at ./process.jl:726
│  [4] run at ./process.jl:724 [inlined]
│  [5] #brew#4(::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Function, ::Cmd) at /Users/yoshinori/.julia/packages/Homebrew/s09IX/src/API.jl:19
│  [6] #brew at ./none:0 [inlined]
│  [7] #tap#53(::Bool, ::Bool, ::Function, ::String) at /Users/yoshinori/.julia/packages/Homebrew/s09IX/src/API.jl:619
│  [8] tap at /Users/yoshinori/.julia/packages/Homebrew/s09IX/src/API.jl:599 [inlined]
│  [9] install_brew() at /Users/yoshinori/.julia/packages/Homebrew/s09IX/src/private_API.jl:52
│  [10] __init__() at /Users/yoshinori/.julia/packages/Homebrew/s09IX/src/Homebrew.jl:37
│  [11] _include_from_serialized(::String, ::Array{Any,1}) at ./loading.jl:633
│  [12] _require_from_serialized(::String) at ./loading.jl:684
│  [13] _require(::Base.PkgId) at ./loading.jl:967
│  [14] require(::Base.PkgId) at ./loading.jl:858
│  [15] require(::Module, ::Symbol) at ./loading.jl:853
│  [16] top-level scope at /Users/yoshinori/.julia/packages/Homebrew/s09IX/deps/build.jl:2
│  [17] include at ./boot.jl:326 [inlined]
│  [18] include_relative(::Module, ::String) at ./loading.jl:1038
│  [19] include(::Module, ::String) at ./sysimg.jl:29
│  [20] include(::String) at ./client.jl:403
│  [21] top-level scope at none:0
│ during initialization of module Homebrew
│ in expression starting at /Users/yoshinori/.julia/packages/Homebrew/s09IX/deps/build.jl:1
└ @ Pkg.Operations /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/Pkg/src/Operations.jl:1075
  Building Cairo ───→ `~/.julia/packages/Cairo/CXPG1/deps/build.log`
┌ Error: Error building `Cairo`: 
│ ==> Tapping homebrew/core
│ Cloning into '/Users/yoshinori/.julia/packages/Homebrew/s09IX/deps/usr/Library/Taps/homebrew/homebrew-core'...
│ fatal: unable to access 'https://github.com/Homebrew/homebrew-core/': Could not resolve host: github.com
│ Error: Failure while executing; `git clone https://github.com/Homebrew/homebrew-core /Users/yoshinori/.julia/packages/Homebrew/s09IX/deps/usr/Library/Taps/homebrew/homebrew-core` exited with 128.
│ ERROR: LoadError: InitError: failed process: Process(`/Users/yoshinori/.julia/packages/Homebrew/s09IX/deps/usr/bin/brew tap --full homebrew/core`, ProcessExited(1)) [1]
│ Stacktrace:
│  [1] error(::String, ::Base.Process, ::String, ::Int64, ::String) at ./error.jl:42
│  [2] pipeline_error at ./process.jl:785 [inlined]
│  [3] #run#515(::Bool, ::Function, ::Cmd) at ./process.jl:726
│  [4] run at ./process.jl:724 [inlined]
│  [5] #brew#4(::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Function, ::Cmd) at /Users/yoshinori/.julia/packages/Homebrew/s09IX/src/API.jl:19
│  [6] #brew at ./none:0 [inlined]
│  [7] #tap#53(::Bool, ::Bool, ::Function, ::String) at /Users/yoshinori/.julia/packages/Homebrew/s09IX/src/API.jl:619
│  [8] tap at /Users/yoshinori/.julia/packages/Homebrew/s09IX/src/API.jl:599 [inlined]
│  [9] install_brew() at /Users/yoshinori/.julia/packages/Homebrew/s09IX/src/private_API.jl:52
│  [10] __init__() at /Users/yoshinori/.julia/packages/Homebrew/s09IX/src/Homebrew.jl:37
│  [11] _include_from_serialized(::String, ::Array{Any,1}) at ./loading.jl:633
│  [12] _require_search_from_serialized(::Base.PkgId, ::String) at ./loading.jl:713
│  [13] _require(::Base.PkgId) at ./loading.jl:937
│  [14] require(::Base.PkgId) at ./loading.jl:858
│  [15] require(::Module, ::Symbol) at ./loading.jl:853
│  [16] top-level scope at /Users/yoshinori/.julia/packages/Cairo/CXPG1/deps/build.jl:41
│  [17] include at ./boot.jl:326 [inlined]
│  [18] include_relative(::Module, ::String) at ./loading.jl:1038
│  [19] include(::Module, ::String) at ./sysimg.jl:29
│  [20] include(::String) at ./client.jl:403
│  [21] top-level scope at none:0
│ during initialization of module Homebrew
│ in expression starting at /Users/yoshinori/.julia/packages/Cairo/CXPG1/deps/build.jl:40
└ @ Pkg.Operations /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/Pkg/src/Operations.jl:1075
```
再度、```Pkg.add```してみたところ、あっさり終了。なんか怪しいので```Pkg.rm```で一旦削除してから```Pkg.add```してみたところ、あっさり終了。というより```Pkg.rm```がうまくいってないんじゃないかな〜？！

```shell
julia> Pkg.add("Cairo")
 Resolving package versions...
  Updating `~/.julia/environments/v1.1/Project.toml`
 [no changes]
  Updating `~/.julia/environments/v1.1/Manifest.toml`
 [no changes]

julia> Pkg.rm("Cairo")
  Updating `~/.julia/environments/v1.1/Project.toml`
  [159f3aea] - Cairo v0.5.6
  Updating `~/.julia/environments/v1.1/Manifest.toml`
  [e1450e63] - BufferedStreams v1.0.0
  [159f3aea] - Cairo v0.5.6
  [a2bd30eb] - Graphics v0.4.0
  [0862f596] - HTTPClient v0.2.1
  [d9be37ee] - Homebrew v0.7.1
  [b27032c2] - LibCURL v0.4.1
  [522f3ed2] - LibExpat v0.5.0
  [2ec943e9] - Libz v1.0.0
  [c17dfb99] - WinRPM v0.4.2

julia> Pkg.add("Cairo")
 Resolving package versions...
  Updating `~/.julia/environments/v1.1/Project.toml`
  [159f3aea] + Cairo v0.5.6
  Updating `~/.julia/environments/v1.1/Manifest.toml`
  [e1450e63] + BufferedStreams v1.0.0
  [159f3aea] + Cairo v0.5.6
  [a2bd30eb] + Graphics v0.4.0
  [0862f596] + HTTPClient v0.2.1
  [d9be37ee] + Homebrew v0.7.1
  [b27032c2] + LibCURL v0.4.1
  [522f3ed2] + LibExpat v0.5.0
  [2ec943e9] + Libz v1.0.0
  [c17dfb99] + WinRPM v0.4.2
```
試しに ```using```してみたところ以下のようなエラーメッセージが。

```shell
julia> using Cairo
[ Info: Precompiling Cairo [159f3aea-2a34-519c-b102-8c37f9878175]
ERROR: LoadError: Cairo not properly installed. Please run
Pkg.build("Cairo")
Stacktrace:
 [1] error(::String, ::String) at ./error.jl:42
 [2] top-level scope at none:0
 [3] include at ./boot.jl:326 [inlined]
 [4] include_relative(::Module, ::String) at ./loading.jl:1038
 [5] include(::Module, ::String) at ./sysimg.jl:29
 [6] top-level scope at none:2
 [7] eval at ./boot.jl:328 [inlined]
 [8] eval(::Expr) at ./client.jl:404
 [9] top-level scope at ./none:3
in expression starting at /Users/yoshinori/.julia/packages/Cairo/CXPG1/src/Cairo.jl:9
ERROR: Failed to precompile Cairo [159f3aea-2a34-519c-b102-8c37f9878175] to /Users/yoshinori/.julia/compiled/v1.1/Cairo/l6vnT.ji.
Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] compilecache(::Base.PkgId, ::String) at ./loading.jl:1197
 [3] _require(::Base.PkgId) at ./loading.jl:960
 [4] require(::Base.PkgId) at ./loading.jl:858
 [5] require(::Module, ::Symbol) at ./loading.jl:853
```
- ```Pkg.build```してみろとのことなので、ご指示に従って```Pkg.build```します。

```shell
julia> Pkg.build("Cairo")
  Building LibCURL ─→ `~/.julia/packages/LibCURL/OoXMv/deps/build.log`
  Building WinRPM ──→ `~/.julia/packages/WinRPM/Y9QdZ/deps/build.log`
  Building Homebrew → `~/.julia/packages/Homebrew/s09IX/deps/build.log`
  Building Cairo ───→ `~/.julia/packages/Cairo/CXPG1/deps/build.log`

julia> using Cairo
[ Info: Precompiling Cairo [159f3aea-2a34-519c-b102-8c37f9878175]
[ Info: Loading Cairo backend into Compose.jl
```
Homebrewのビルドにがっつり時間がかかったので今回はうまくいったようです。```using```したらめでたくプリコンパイルが走りました。

# 本編
さて、事前の準備はここまでで、これ以降はJupyter notebookから作業します。
[前回の記事](https://qiita.com/koinori/items/bdb01ec789d5e49aed0f)で「IJulia」はインストール済みなので jupyter notebookを起動し、「New」->「Julia 1.1.0」で新規のJuliaノートブックを作成します。
実際のグラフなどはこちらのnbviewerからご覧ください。
- [Julia 1.1 + Jupyter で回帰分析サンプル](https://nbviewer.jupyter.org/gist/Yoshinori-Koide-PRO/15da84e0d21c29f0c63ab4ee57face82)

それではさっそく、進めていきます。

```Julia
using Distributions
using Gadfly
using DataFrames

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

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

set_default_plot_size(7inch, 7inch/MathConstants.golden);

# 散布図の描画
plot(x=x, y=y)
```
ここで```Base.golden```は```Base.MathConstants.golden```に移動していました。
- https://docs.julialang.org/en/v1.1/base/numbers/#Base.MathConstants.golden

```Julia
# 折れ線グラフ
plot(x=rand(50), y=rand(50), Geom.point, Geom.line)
```
ここは問題なし。

```Julia
using RDatasets
ds = RDatasets.datasets()
```
データセットの一覧も表示されました。

```Julia
iris = dataset("datasets", "iris")
println(iris)
```
実際のデータも取得おっけい。

```Julia
plot(dataset("datasets", "iris"), x="SepalLength", y="SepalWidth", color="Species",Geom.point)
```
散布図もプロットできます。

```Julia
using GLM
housing = dataset("Ecdat", "Housing")
# まずは散布図のプロット
plot(housing, x="LotSize", y="Price", Geom.point)
```
GLMの準備はOKです。

```Julia
# 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)"))
```
上記の```log.(housing[:,2])```で例の```.(dot)```演算子の登場です。詳しくは以下のページを見つけました。
- https://docs.julialang.org/en/v1.1/manual/mathematical-operations/#man-dot-operators-1
- https://docs.julialang.org/en/v1.1/manual/arrays/#Broadcasting-1
- https://docs.julialang.org/en/v1.1/base/arrays/#Base.Broadcast.broadcast
配列の"ブロードキャスティング"ってやつを明示的にした、というのが新しい点でしょうか。

```Julia
# 線形回帰の実行
olm = fit(LinearModel, @formula(Price ~ LotSize), log_housing)
```
ここの```@formula```の書き方もだいぶ探しました。GLMのドキュメントにたどり着けたのでなんとか動かせました。とりあえず最新のGLMのAPIは以下です。
- http://juliastats.github.io/GLM.jl/latest/manual/

```Julia
# グラフに直線を描画
intercept = coef(olm)[1]
slope     = coef(olm)[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)")
)
```
そしてここの```.+```もブロードキャストしてます。

```Julia
using HTTP
using CSV

# Children Ever Born Data
# Lecture Note: http://data.princeton.edu/wws509/notes/c4.pdf
# Dataset: http://data.princeton.edu/wws509/datasets/#ceb
url = "https://data.princeton.edu/wws509/datasets/ceb.dat"
res = HTTP.request("GET", url)
body = String(res.body)
```
元の記事では ```Requests```パッケージを使っていましたがこれは非推奨になっていたので```HTTP```パッケージを使用するように修正しています。

```Julia
df_ceb = (
    CSV.File(
        IOBuffer(
            # 行頭の空白+数字を数字だけにしておく
            replace(String(body), r"^ (\d)"m => s"\1")
            # スペース区切りで連続したスペースは無視
            ), delim=' ', ignorerepeated=true, types=Dict(:n=>Float64,:y=>Float64))
    |> DataFrame
)[:,2:8]
# n と y はFloat64型のまま整数値にする
df_ceb[:n] = round.(df_ceb[:n])
df_ceb[:y] = round.(df_ceb[:y])
df_ceb
```
Http通信で受け取ったBodyをCSVにパースしてDataFrameに変換すること自体は簡単になりました。```CSV```パッケージのおかげです。
- http://juliadata.github.io/CSV.jl/stable/index.html

また、元データの1〜10レコード目まではインデックスの数字が空白+数字で2文字に埋めているので、空白区切でパースさせるとカラムがずれちゃうんですよね。。。そんなわけで正規表現で行頭の空白+数字を数字だけに変換させるようにいたしました。
- https://docs.julialang.org/en/v1.1/manual/strings/#Regular-Expressions-1

そしてDataFrameの列ごと取得＆代入はpandasとだいたい同じような操作ができるので助かります。
- https://juliadata.github.io/DataFrames.jl/stable/index.html

```Julia
# まずは散布図のプロット
plot(df_ceb, x="n", y="y", Geom.point)
```
図のプロットは問題なし。

```Julia
# ポアソン回帰の実行
gm1 = fit(GeneralizedLinearModel, @formula(y ~ n), df_ceb, Poisson())
```
ここは先ほどの```@formula```と同じですね。また、```LogLink()```はデフォルトで省略可能でしたね。

```Julia
intercept = coef(gm1)[1]
slope     = coef(gm1)[2]

# 散布図の上のポアソン回帰の結果の曲線を載せてみる
xx = collect(0:1:200)
yy = exp.(intercept .+ slope * 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))
)
```
めでたくポアソン回帰の曲線が重ねられました。

# 動く Notebookはこちらから

再掲いたしますがGistにあげた notebookをnbviewer形式でちゃんと動くようにいたしました。以下のリンクからどうぞ。

- [Julia 1.1 + Jupyter で回帰分析サンプル](https://nbviewer.jupyter.org/gist/Yoshinori-Koide-PRO/15da84e0d21c29f0c63ab4ee57face82)

グラフのパンニングなどもヌルヌル動くのでいいですね〜！