juliaでよくつかうパッケージまとめ

Author: @takuizum
Date: 5th Dec 2019

※このドキュメントはjuliaのWeave.jlで作成しています。

どうも“@takuizum”です。

5日目のカレンダーが空白だったので埋めます。この記事では個人的によく使うパッケージを紹介します。あまり深く突っ込んだことは書いてないです。「ほーん，こんなもんか」くらいに読んでください。

パッケージの追加方法

コードを評価する場合はこちら。 using Pkg;Pkg.add("PkgName")

REPLに直接打ち込むのであれば，PkgREPLが便利です。

  1. REPLに]を打ち込んでPkgREPLに移動
  2. add hoge

パッケージ名のあとに@数字を入力するとバージョンを指定できます。#でブランチ名やコミットの数字（ハッシュ？よくわからん～）を指定することもできます。

パッケージの使い方

julia REPLでusing hogeとすれば，必要に応じてprecompileが走ってすぐ（物によっては1分くらいかかる）使えるようになります。

using Plots
y = randn(100)
plot(y)

 

複数のパッケージで関数の名前が衝突して，指示が曖昧になるのを避けるには，明示的にhoge.fuga()と実行すればよいです。

WARNING: using Gadfly.plot in module Main conflicts with an existing identifier.

using Gadfly, DataFrames
# 2019/12/07 Weaveの実行関数でエラーが出るためコメントアウト
# SystemError("opening file ~
Gadfly.plot(DataFrame(x = 1:100, y = y), x = :x, y = :y , Geom.line)

 

import hoge.fugaとすればhogeパッケージのhuga関数だけをインポートできます。

import Statistics.std
std([1:1:10;])

3.0276503540974917

var() # 同じStatisticsにあっても，こちらは使えない。

ERROR: UndefVarError: var not defined

よく使うパッケージ一覧

いち心理学系ユーザーが個人的によく使う程度なので，範囲は十分ではないです。 野良ではない公式のパッケージはここで調べることができます。

例えばおなじベイズ推論を扱うパッケージでも，Turingは公式にありますがGenはgitから引っ張ってくる必要がある野良パッケージに当たります。

Statistics

記述統計用の関数が揃っています。平均，中央値，quantile，標準偏差，分散，相関，共分散が含まれます。

Rと同じで分散・標準偏差は，デフォルトだと不偏のほうを求めます。

using Statistics
Statistics.mean([1 2 3])
Statistics.std(1:10, corrected = false) # デフォルトだと不偏標準偏差を計算する。
[Statistics.var(1:10, corrected = true)
 Statistics.var(1:10, corrected = false)]

2-element Array{Float64,1}: 9.166666666666666 8.25

NewFeature julia1.3からmean関数に，任意の関数fを適用してから平均をとる機能が追加されました（誰得？）

mean(√,[1 2 3;4 5 6], dims = 2)

2×1 Array{Float64,2}: 1.3820881233139908 2.2285192400943226

StatsFuns

ロジスティック関数logistic(x::Real)やソフトマックス関数sodtma(x::Array)などが含まれています。心理系だとよくお世話になる関数なのでたまに使います。

そのほかオーバーフロー回避で使われるlogsumexp(x)もあります。

using StatsFuns
Gadfly.plot(DataFrame(x = -10:.1:10, y = logistic.(-10:.1:10)), x = :x, y = :y , Geom.line)

 

StatsBase

積率や度数，共分散などを計算することができるパッケージです。

using StatsBase
# 任意のレンジで度数をカウント
counts(sample([1:1:100;], 100), 1:100)

100-element Array{Int64,1}: 2 1 0 1 0 2 0 1 1 2 ⋮ 4 0 2 0 0 0 1 2 1

おすすめはPlotsとのあわせ技で，かんたんに2変量のヒストグラムをかけるfit(Histogram, (X, Y)) |> Plots.plotです

SpecialFunstions

ベータ分布系のモデルを微分したりするとお世話になる，digammaやtrigammaなどの関数を扱うパッケージです。それ以外の関数は私にはよくわかりません。

\(\ln{Γ}\)を計算したいときにlgammaよりもlogabsgammaを使うとより高速に計算できます。ちなみにlgammaを使うと以下のような警告が出ます。 "lgamma(x::Real)is deprecated, use(logabsgamma(x))[1]instead."

でもloggammaを使ったほうがより速いみたいです。logabsgammaはベクトルで計算するとTuplesのベクトルが出てきて，ちょっと厄介です。

using SpecialFunctions
log(1);gamma(1);lgamma(1);loggamma(1);logabsgamma(1);# precompile
@time log.(gamma.(0.001:0.001:10))

0.071152 seconds (189.39 k allocations: 9.678 MiB)

@time lgamma.(0.001:0.001:10)

3.418047 seconds (1.55 M allocations: 139.278 MiB, 0.70% gc time)

@time loggamma.(0.001:0.001:10)

0.034178 seconds (140.87 k allocations: 7.096 MiB)

@time (logabsgamma.(0.001:0.001:10));

0.035795 seconds (147.86 k allocations: 7.603 MiB)

Distributions

様々な確率分布を扱うパッケージです。統計モデルを扱う上では必須です。 乱数を発生させたいときはRandomも一緒にインポートしておき，randを使います。確率密度はpdf.です（.をつけることが推奨される）。

using Distributions
using Random
Random.seed!(1205)
d = rand(Beta(1, 3), 1000)
histogram(d, normalize = true)
plot!(0:.01:1, pdf.(Beta(1, 3), 0:.01:1))

 

他にもサンプルデータに確率分布をフィットさせて，パラメタを最尤推定することも簡単にできます（任意の制約でMAP推定もできるみたいです）。

fit(Beta, d)

Distributions.Beta{Float64}(α=1.0566418122764945, β=3.0216405369600863)

2019/12/07 アップデートしました。

juliaプログラミングクックブックの9章を参考にし，非負のパラメタ推定を行う際にlogをとった値を推定することで，関数で計算できないパラメタが探索されることを避けるように変更しました。

# optimで最尤推定
using Optim
fn(x, logm, logn) = -sum(logpdf.(Beta(exp(logm), exp(logn)), x))
opt = optimize(x -> fn(d, x[1], x[2]), [-1.0, 2.0])# 結構初期値に敏感....
exp.(opt.minimizer)

2-element Array{Float64,1}: 1.0426010627969013 2.9947359227621457

ForwardDiff & LinearAlgebra

juliaで高速な数値微分を提供してくれるパッケージです。gradientやhessianの評価を近似的に，お気楽にできるのでとても重宝しています。

一緒に使っているLinearAlgebraも，固有値や行列の操作をするときによく使います。心理学徒は大好きですよね，固有値。

試しに先程最適化した値の標準誤差を計算してみます。

using ForwardDiff, LinearAlgebra
hess = ForwardDiff.hessian(x -> fn(d, x), opt.minimizer)

ERROR: MethodError: no method matching fn(::Array{Float64,1}, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{Main.WeaveSandBox0.var"#3#4",Float64},ForwardDiff.Dual{ForwardDiff.Tag{Main.WeaveSandBox0.var"#3#4",Float64},Float64,2},2},1})
Closest candidates are:
  fn(::Any, ::Any, !Matched::Any) at none:1

sqrt.(LinearAlgebra.diag(inv(hess)))

ERROR: UndefVarError: hess not defined

Gadfly

LINK

R言語にあるggplot2ライクなプロットを書くことができます。

GadflyはThe Grammar Of Graphicsというグラフの文法を元に設計されており，グラフの様々な要素を重ねていくスタイルです。

下のプロットの例ではpointにlmで回帰直線を引いています。グループのラベルがあるときに，グループごとに可視化…というのも楽にできます。

using Gadfly, RDatasets
iris = RDatasets.dataset("datasets", "iris")
Gadfly.plot(iris, x=:SepalLength, y=:PetalLength, Geom.point, color=:Species, Geom.smooth(method=:lm))

 

StatsPlots

LINK

統計分析用のプロットパッケージです。

周辺分布をかんたんに書けたりします。他にもインタラクティブなプロットもできるみたいです。あまり使ったことがないですが，個人的に注目しているパッケージです。

using StatsPlots, Distributions
StatsPlots.marginalhist(sort(rand(Normal(), 1000)), sort(rand(Normal(), 1000)))

 

データフレーム操作パッケージQueryと連携することで，データをざっと要約することもできます。

using Query
# データのざっとした要約
@df iris corrplot([:SepalLength :SepalWidth :PetalLength :PetalWidth], grid = false)

 

Turing

LINK

先日のjulia Tokyo 10でも紹介があった，juliaによるjuliaのためのjuliaなPPLです。特徴は生成モデルの記述によるかんたんな確率モデルの表現だけで，比較的複雑な（階層的な）モデルを表現できるところです。

下のコードは項目反応理論の2パラメタロジスティックモデルです。項目反応モデルはテストやアンケートに受検者が反応する確率をロジスティックモデルで表現しています。

モデルは書きましたが，このモデルをMCMCすると，パラメタ数が多すぎて推定にめっちゃ時間取られるので，ものすごくサンプル数とチェイン数を少なくしています。なので推定精度はダメダメです。

なお，本当はこのモデルを分析するためには，2000人くらいのデータが必要です。

using Turing, StatsFuns, Distributions

# simulation data
function respGen(N, J)
    θ = rand(Normal(), N)
    a = rand(LogNormal(), J)
    b = rand(Uniform(-2, 2), J)
    resp = zeros(Int64, N, J)
    for i in 1:N
        for j in 1:J
            resp[i,j] = rand(Uniform(0.0,1.0),1)[1] < StatsFuns.logistic(a[j]*(θ[i]-b[j])) ? 1 : 0
        end
    end
    return resp, θ, a, b
end

@model irt2pl(data,  N, J) = begin
   θ = Vector{Float64}(undef, N)
   α = β = Vector{Float64}(undef, J)
   # assign distributon to each element
   θ ~ [Normal(0, 1)]
   α ~ [LogNormal(0, 1)]
   β ~ [Normal(0, 2)]
   for i = 1:N
       for j = 1:J
           p = logistic(α[j]*(θ[i]-β[j]))
           data[i,j] ~ Bernoulli(p)
       end
   end
end

data, θ, a, b = respGen(200, 20);
N, J = size(data)
iterations = 100
num_chains = 1
Turing.setadbackend(:forward_diff)
chain_HMCDA = sample(irt2pl(data, N, J), HMCDA(0.15, 0.65), iterations) # Fast!!

# visualize
using Plots
est = mean(chain_HMCDA[:θ][:,:,1], dims = 3)
Plots.scatter(θ, est.df.mean, ylims = [-4,4], xlims = [-4,4], xlabel = "true", ylabel = "estimate")

 

Published from summaryPkg.jmd using Weave.jl on 2019-12-16.