元ページ

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

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

準備

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

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

$ 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>

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

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をインストール？中に以下のエラーで中断してしまいました。

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がうまくいってないんじゃないかな〜？！

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してみたところ以下のようなエラーメッセージが。

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します。

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から作業します。
前回の記事で「IJulia」はインストール済みなので jupyter notebookを起動し、「New」->「Julia 1.1.0」で新規のJuliaノートブックを作成します。
実際のグラフなどはこちらのnbviewerからご覧ください。

Julia 1.1 + Jupyter で回帰分析サンプル

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

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

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

ここは問題なし。

using RDatasets
ds = RDatasets.datasets()

データセットの一覧も表示されました。

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

実際のデータも取得おっけい。

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

散布図もプロットできます。

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

GLMの準備はOKです。

# 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
配列の"ブロードキャスティング"ってやつを明示的にした、というのが新しい点でしょうか。

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

ここの@formulaの書き方もだいぶ探しました。GLMのドキュメントにたどり着けたのでなんとか動かせました。とりあえず最新のGLMのAPIは以下です。

http://juliastats.github.io/GLM.jl/latest/manual/

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

そしてここの.+もブロードキャストしてます。

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パッケージを使用するように修正しています。

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

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

図のプロットは問題なし。

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

ここは先ほどの@formulaと同じですね。また、LogLink()はデフォルトで省略可能でしたね。

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 で回帰分析サンプル

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