StatsFragments

Python, R, Rust, 統計, 機械学習とか

R {ggplot2} で 少し複雑なサブプロットが描きたい

この記事は R Advent Calendar 2015 4 日目の記事です。


{ggplot2} でのサブプロット

{ggplot2} でサブプロットを描画したいことがある。同じ種類のプロットを水準別に描画するなど、単純なものであれば facet で描ける。例えば 適当な散布図を Species 別に描きたい場合、

library(dplyr)
library(ggplot2)

ggplot(iris, aes(Sepal.Width, Sepal.Length)) + 
  geom_point() + 
  facet_wrap(~ Species)

f:id:sinhrks:20151204000924p:plain

facet では異なる種類のプロットをサブプロットとすることはできない。例えば散布図の縦軸/横軸をそれぞれ変えたい、といった場合には gridExtra::grid.arrange を使う。

library(gridExtra)
p1 <- ggplot(iris, aes(Sepal.Width, Sepal.Length)) + geom_point()
p2 <- ggplot(iris, aes(Petal.Width, Petal.Length)) + geom_point()
p3 <- ggplot(iris, aes(Sepal.Length, Petal.Length)) + geom_point()
p4 <- ggplot(iris, aes(Sepal.Width, Petal.Width)) + geom_point()
grid.arrange(p1, p2, p3, p4)

f:id:sinhrks:20151204000934p:plain

が、grid.arrange は引数を ドット ... で受け取ること、また呼び出すと直ちに描画を行うことから、以下のような場合は少し面倒だ。

  • 描画が直ちに行われるため、サブプロットの描画設定をまとめて変更できない。単一のテーマを使う場合でも、個々のプロットそれぞれについてあらかじめ処理を行っておく必要がある。
  • リストが渡せない。サブプロットの数を動的に変更したい場合は、リストに入れて do.call が必要である。

というわけで、自作パッケージ {ggfortify} でこの辺の操作を少し便利にできるようにした。サブプロット周りの処理は ggmultiplot クラスで行う。

コンストラクタの引数にはリストが渡せるため、サブプロット数が変わる場合でも楽だと思う。また、インスタンスに対して + 演算子を用いることで ggplot と同じ処理を各プロットに対して適用できる。

# install.packages('ggfortify')
library(ggfortify)
p <- new('ggmultiplot', plots = list(p1, p2, p3, p4))
p
# 略。上記 grid.arrange と同じ出力

# theme_bw を全サブプロットに適用
p + theme_bw()

f:id:sinhrks:20151204000945p:plain

サブプロットの要素には リストと同じようにアクセスできる。[ でのアクセスは ggmultiplot インスタンスに、[[ での アクセスは ggplot インスタンスになる。代入もできるため、一部のプロットにだけ処理を適用したい場合は以下のように書ける。

p[2:3] <- p[2:3] + stat_smooth(method = 'lm')
p

f:id:sinhrks:20151204000953p:plain

+ 演算子の右項に ggplot インスタンスを指定するとサブプロットが追加できる。

p + (ggplot(iris, aes(Sepal.Width, fill = Species)) + geom_density())

f:id:sinhrks:20151204002618p:plain

ggmultiplot 作成時に ncol もしくは nrow キーワードを指定することで、サブプロットの列数 / 行数が指定できる。

p <- new('ggmultiplot', plots = list(p1, p2, p3, p4), ncol = 3)
p

f:id:sinhrks:20151204002529p:plain

autoplot との組み合わせ

{ggfortify} を使うと主要なクラスを ggplot2::autoplot 関数で描画できるようになる。autoplot にリストを渡した場合は、リストの各要素が autoplot 可能であればそれらをサブプロットとして描画する。

サンプルとして、{broom} パッケージの vignettes にある例をやってみたい。ここでは、適当なデータに対して クラスタ数を変えながら kmeans クラスタリングを行った結果をプロットしている。

これは autoplot を使って以下のように書ける。

kclusts <- data.frame(k=1:9) %>% 
  dplyr::group_by(k) %>%
  dplyr::do(kclust = kmeans(iris[-5], .$k))
autoplot(kclusts$kclust, data = iris[-5], ncol = 3) + theme(legend.position = "none")

f:id:sinhrks:20151204001005p:plain

一行目は今なら {purrr} を使ったほうが簡潔でわかりやすくなると思う。

library(purrr)
kclusts <- purrr::map(1:9, ~ kmeans(iris[-5], .))
autoplot(kclusts, data = iris[-5], ncol = 3) + theme(legend.position = "none")
# 略

まとめ

{ggfortify} では ggmultiplot クラス、ならびに autoplot を利用してサブプロットを簡単に扱うことができる。

Rグラフィックスクックブック ―ggplot2によるグラフ作成のレシピ集

Rグラフィックスクックブック ―ggplot2によるグラフ作成のレシピ集

{purrr} でリストデータを操作する <2>

前の記事に続けて、{purrr}{rlist} 相当の処理を行う。今回はレコードの選択とソート。

sinhrks.hatenablog.com

サンプルデータは前回と同じものを利用する。リストの表示も同じく Hmisc::list.tree を使う。

library(rlist)
library(pipeR)
library(purrr)

packages <- list(
  list(name = 'dplyr', star = 979L, maintainer = 'hadley' , authors = c('hadley', 'romain')),
  list(name = 'ggplot2', star = 1546L, maintainer = 'hadley' , authors = c('hadley')),
  list(name = 'knitr', star = 1047L, maintainer = 'yihui' , authors = c('yihui', 'hadley', '...and all'))
)
packages
##  x = list 3 (2632 bytes)
## .  [[1]] = list 4
## . .  name = character 1= dplyr 
## . .  star = integer 1= 979
## . .  maintainer = character 1= hadley 
## . .  authors = character 2= hadley romain 
## .  [[2]] = list 4
## . .  name = character 1= ggplot2 
## . .  star = integer 1= 1546
## . .  maintainer = character 1= hadley 
## . .  authors = character 1= hadley 
## .  [[3]] = list 4
## . .  name = character 1= knitr 
## . .  star = integer 1= 1047
## . .  maintainer = character 1= yihui 
## . .  authors = character 3= yihui hadley  ...

レコードの選択

rlist Tutorial: Filtering 後半の内容。

ある条件にあてはまるレコードのうち最初の N 件だけを取得したいとき、{rlist} には専用の関数 list.find がある。{purrr} には対応する関数はないが、keep %>% head で同じ結果が得られる。

packages %>>% rlist::list.find(star >= 1000, 1)
##  x = list 1 (840 bytes)
## .  [[1]] = list 4
## . .  name = character 1= ggplot2 
## . .  star = integer 1= 1546
## . .  maintainer = character 1= hadley 
## . .  authors = character 1= hadley

packages %>%
  purrr::keep(~ .$star >= 1000) %>%
  head(n = 1)
# 略

また、{rlist} には list.findi という インデックスを返す関数がある。purrr::keep では、引数に logical 型のベクトルを渡せば TRUE に対応する要素のみを残すことができる。

packages %>>% rlist::list.findi(star >= 1000, 2)
## [1] 2 3

purrr::keep(seq_along(packages),
            flatmap(packages, ~ .$star >= 1000)) %>%
  head(n = 2)
# 略

条件にあてはまる最初のレコードだけを取得したいときは rlist::first もしくは purrr::detect

rlist::list.first(packages, star >= 1000)
##  x = list 4 (792 bytes)
## .  name = character 1= ggplot2 
## .  star = integer 1= 1546
## .  maintainer = character 1= hadley 
## .  authors = character 1= hadley

purrr::detect(packages, ~ .$star >= 1000)
# 略

レコード数がそこまで多くない場合は keep %>% head(n = 1) でよいと思うが、結果に 1 レベル目が入れ子として残ってしまう。上と同じ結果を得るには purrr::flatten ( unlist(recursive = FALSE) と同じ ) を通す。

purrr::keep(packages, ~ .$star >= 1000) %>% head(n = 1)
##  x = list 1 (840 bytes)
## .  [[1]] = list 4
## . .  name = character 1= ggplot2 
## . .  star = integer 1= 1546
## . .  maintainer = character 1= hadley 
## . .  authors = character 1= hadley

purrr::keep(packages, ~ .$star >= 1000) %>%
  head(n = 1) %>%
  purrr::flatten()
##  x = list 4 (792 bytes)
## .  name = character 1= ggplot2 
## .  star = integer 1= 1546
## .  maintainer = character 1= hadley 
## .  authors = character 1= hadley

同じく、最後のレコードを取得したいときは rlist::last もしくは purrr::detect(.right = TRUE)。これも自分としては keep %>% tail(n = 1) のほうが覚えやすい。

list.last(packages, star >= 1000)
##  x = list 4 (920 bytes)
## .  name = character 1= knitr 
## .  star = integer 1= 1047
## .  maintainer = character 1= yihui 
## .  authors = character 3= yihui hadley  ...

purrr::detect(packages, ~ .$star >= 1000, .right = TRUE)
# 略

purrr::keep(packages, ~ .$star >= 1000) %>%
  tail(n = 1) %>%
  purrr::flatten()
# 略

類似の処理として、{rlist} には list.take という最初の N レコードを取得する関数と list.skip という最初の N レコードを除外する関数がある。{purrr} には対応する関数はないが、組み込みの headtail でよいのでは...。

packages %>>% rlist::list.take(2)
##  x = list 2 (1696 bytes)
## .  [[1]] = list 4
## . .  name = character 1= dplyr 
## . .  star = integer 1= 979
## . .  maintainer = character 1= hadley 
## . .  authors = character 2= hadley romain 
## .  [[2]] = list 4
## . .  name = character 1= ggplot2 
## . .  star = integer 1= 1546
## . .  maintainer = character 1= hadley 
## . .  authors = character 1= hadley

packages %>% head(n = 2)
# 略

packages %>% rlist::list.skip(1)
##  x = list 2 (1768 bytes)
## .  [[1]] = list 4
## . .  name = character 1= ggplot2 
## . .  star = integer 1= 1546
## . .  maintainer = character 1= hadley 
## . .  authors = character 1= hadley 
## .  [[2]] = list 4
## . .  name = character 1= knitr 
## . .  star = integer 1= 1047
## . .  maintainer = character 1= yihui 
## . .  authors = character 3= yihui hadley  ...

packages %>% tail(length(.) - 1)
# 略

ある条件が初めて偽となるまでレコードを取得する場合は rlist::list.takeWhile もしくは purrr::head_while

packages %>>% rlist::list.takeWhile(star <= 1500)
##  x = list 1 (896 bytes)
## .  [[1]] = list 4
## . .  name = character 1= dplyr 
## . .  star = integer 1= 979
## . .  maintainer = character 1= hadley 
## . .  authors = character 2= hadley romain

packages %>% purrr::head_while(~ .$star <= 1500)
# 略

ある条件が初めて偽となったレコード以降を取得したい場合は rlist::list.skipWhile{purrr} には対応する処理はないが、purrr::detect_index を使って以下のように書ける。

packages %>>% rlist::list.skipWhile(star <= 1500)
##  x = list 2 (1768 bytes)
## .  [[1]] = list 4
## . .  name = character 1= ggplot2 
## . .  star = integer 1= 1546
## . .  maintainer = character 1= hadley 
## . .  authors = character 1= hadley 
## .  [[2]] = list 4
## . .  name = character 1= knitr 
## . .  star = integer 1= 1047
## . .  maintainer = character 1= yihui 
## . .  authors = character 3= yihui hadley  ...

purrr::detect_index(packages, ~ .$star > 1500)
## [1] 2
packages[purrr::detect_index(packages, ~ .$star > 1500):length(packages)]
# 略 

一方、{purrr} には tail_while という末尾からみて条件にあてはまるレコードだけを取得する関数がある。

packages %>% purrr::tail_while(~ .$star <= 1500)
##  x = list 1 (968 bytes)
## .  [[1]] = list 4
## . .  name = character 1= knitr 
## . .  star = integer 1= 1047
## . .  maintainer = character 1= yihui 
## . .  authors = character 3= yihui hadley  ...

ある条件を満たす要素がいくつあるかを調べるには rlist::list.count{purrr} なら keep %>% length でよい。

list.count(packages, star > 1000)
## [1] 2

purrr::keep(packages, ~ .$star > 1000) %>% length()
# 略

{rlist} には リストの要素の名前を正規表現で抽出する list.match がある。{purrr} なら keep でできる。stringr::str_match はマッチした結果を matrix で返すため、complete.caseslogical に変換する。

data <- list(p1 = 1, p2 = 2, a1 = 3, a2 = 4)
list.match(data, "p[12]")
##  x = list 2 (416 bytes)
## .  p1 = double 1= 1
## .  p2 = double 1= 2

complete.cases(stringr::str_match(names(data), 'p[12]'))
## [1]  TRUE  TRUE FALSE FALSE
keep(data, complete.cases(stringr::str_match(names(data), 'p[12]')))
# 略

logical の処理

条件に当てはまるかどうかを logical ベクトルとして返すには rlist::list.is。これは rlist::list.mapv と同じなのでは...? {purrr} なら flatmap もしくは map_lgl

packages %>>% rlist::list.is(star < 1500)
## [1]  TRUE FALSE  TRUE

packages %>>% rlist::list.mapv(star < 1500)
# 略

packages %>% purrr::flatmap(~ .$star < 1500)
# 略

packages %>% purrr::map_lgl(~ .$star < 1500)
# 略

また、logical のリスト全体の真偽値を以下のように取得できる

  • 要素全てが TRUE かどうか: rlist::list.all もしくは purrr::every
  • 要素いずれかが TRUE かどうか: rlist::list.any もしくは purrr::some

自分は purrr::flatmap %>% all もしくは purrr::flatmap %>% any のほうが覚えやすい。

ただし、{purrr} v0.1.0 の every, some にはバグがあり正しく動作しない。下のスクリプトを実行するためには GitHub から最新の master をインストールする必要がある。

rlist::list.all(packages, star >= 500)
## [1] TRUE

# install_github('hadley/purrr') # が必要
purrr::every(packages, ~ .$star >= 500) 
# 略

packages %>% purrr::flatmap(~ .$star >= 500) %>% all()
# 略

rlist::list.any(packages, star >= 1500)
## [1] TRUE

# install_github('hadley/purrr') # が必要
purrr::some(packages, ~ .$star >= 1500)

packages %>% purrr::flatmap(~ .$star >= 1500) %>% any()
# 略

要素の削除

{rlist} では list.remove で指定した名前をリストから除外できる。{purrr} には discard という条件にマッチする要素をリストから除外する関数がある ( purrr::keep とは逆 )。条件の否定をとれば keep でも書ける。

rlist::list.remove(data, c("p1", "p2"))
##  x = list 2 (416 bytes)
## .  a1 = double 1= 3
## .  a2 = double 1= 4

purrr::discard(data, names(data) %in% c("p1", "p2"))
# 略

purrr::keep(data, !(names(data) %in% c("p1", "p2")))
# 略

{rlist} で、名前ではなく 条件式を使ってリストから除外したい場合は list.exclude{purrr} の場合は 上と同じく discard もしくは keep

packages %>>% rlist::list.exclude("yihui" %in% authors)
##  x = list 2 (1696 bytes)
## .  [[1]] = list 4
## . .  name = character 1= dplyr 
## . .  star = integer 1= 979
## . .  maintainer = character 1= hadley 
## . .  authors = character 2= hadley romain 
## .  [[2]] = list 4
## . .  name = character 1= ggplot2 
## . .  star = integer 1= 1546
## . .  maintainer = character 1= hadley 
## . .  authors = character 1= hadley

packages %>% purrr::discard(~ "yihui" %in% .$authors)
# 略

また、{rlist} には 欠損などの異常値を リストから除去する list.clean という関数がある。下の例では リストの 1 レベル目にある欠損 "b" を除去している。{purrr} では同じく discard

x <- list(a = 1, b = NULL, c = list(x = 1, y = NULL, z = logical(0L), w = c(NA, 1)))
x
##  x = list 3 (1040 bytes)
## .  a = double 1= 1
## .  b = NULL 0
## .  c = list 4
## . .  x = double 1= 1
## . .  y = NULL 0
## . .  z = logical 0
## . .  w = double 2= NA 1

rlist::list.clean(x)
##  x = list 2 (960 bytes)
## .  a = double 1= 1
## .  c = list 4
## . .  x = double 1= 1
## . .  y = NULL 0
## . .  z = logical 0
## . .  w = double 2= NA 1

purrr::discard(x, ~ is.null(.))
# 略

rlist::list.cleanrecursive = TRUE を指定することで 処理を再帰的に適用できる。{purrr} ではこういった再帰的な処理は (自分で関数を書かない限り) できないと思う。

rlist::list.clean(x, recursive = TRUE)
##  x = list 2 (912 bytes)
## .  a = double 1= 1
## .  c = list 3
## . .  x = double 1= 1
## . .  z = logical 0
## . .  w = double 2= NA 1

要素の更新

rlist Tutorial: Updating に相当する内容。

rlist::list.mappurrr::map を用いたリストの選択では、結果に含める属性を明示的に指定する必要があった。 これは対象のリストに選択したい属性が多い場合はすこし面倒だ。

rlist::list.update を使うと、もともとのリストの属性を残した上で 指定した要素を追加 / 更新できる。{purrr} で同じことをするには map + update_list

packages %>>% rlist::list.update(lang = 'R', star = star + 1)
##  x = list 3 (3160 bytes)
## .  [[1]] = list 5
## . .  name = character 1= dplyr 
## . .  star = double 1= 980
## . .  maintainer = character 1= hadley 
## . .  authors = character 2= hadley romain 
## . .  lang = character 1= R 
## .  [[2]] = list 5
## . .  name = character 1= ggplot2 
## . .  star = double 1= 1547
## . .  maintainer = character 1= hadley 
## . .  authors = character 1= hadley 
## . .  lang = character 1= R 
## .  [[3]] = list 5
## . .  name = character 1= knitr 
## . .  star = double 1= 1048
## . .  maintainer = character 1= yihui 
## . .  authors = character 3= yihui hadley  ... 
## . .  lang = character 1= R

packages %>% purrr::map(~ purrr::update_list(., lang = 'R', star = ~ .$star + 1))
# 略

一部の要素を取り除きたい場合は NULL を指定する。これは {purrr} も同じ。

packages %>>% rlist::list.update(maintainer = NULL, authors = NULL)
##  x = list 3 (1464 bytes)
## .  [[1]] = list 2
## . .  name = character 1= dplyr 
## . .  star = integer 1= 979
## .  [[2]] = list 2
## . .  name = character 1= ggplot2 
## . .  star = integer 1= 1546
## .  [[3]] = list 2
## . .  name = character 1= knitr 
## . .  star = integer 1= 1047

packages %>% purrr::map(~ purrr::update_list(., maintainer = NULL, authors = NULL))
# 略

ソート

rlist Tutorial: Sorting に相当する内容。

{rlist}, {purrr} とも 標準の sort ( 値のソート ) と order ( インデックスのソート ) それぞれに対応した関数を持つ。

x = c('a', 'c', 'd', 'b')

sort(x)
## [1] "a" "b" "c" "d"

order(x)
## [1] 1 4 2 3

sort のようにソートした値を得るためには rlist::list.sort もしくは purrr::sort_by

packages %>>% rlist::list.sort(star)
##  x = list 3 (2632 bytes)
## .  [[1]] = list 4
## . .  name = character 1= dplyr 
## . .  star = integer 1= 979
## . .  maintainer = character 1= hadley 
## . .  authors = character 2= hadley romain 
## .  [[2]] = list 4
## . .  name = character 1= knitr 
## . .  star = integer 1= 1047
## . .  maintainer = character 1= yihui 
## . .  authors = character 3= yihui hadley  ... 
## .  [[3]] = list 4
## . .  name = character 1= ggplot2 
## . .  star = integer 1= 1546
## . .  maintainer = character 1= hadley 
## . .  authors = character 1= hadley

packages %>% purrr::sort_by('star')
# 略

{rlist} では 列名を括弧でくくると逆順となる。{purrr} には対応する記法がないが、結果を rev で逆順にすればよい。

packages %>>% rlist::list.sort((star))
##  x = list 3 (2632 bytes)
## .  [[1]] = list 4
## . .  name = character 1= ggplot2 
## . .  star = integer 1= 1546
## . .  maintainer = character 1= hadley 
## . .  authors = character 1= hadley 
## .  [[2]] = list 4
## . .  name = character 1= knitr 
## . .  star = integer 1= 1047
## . .  maintainer = character 1= yihui 
## . .  authors = character 3= yihui hadley  ... 
## .  [[3]] = list 4
## . .  name = character 1= dplyr 
## . .  star = integer 1= 979
## . .  maintainer = character 1= hadley 
## . .  authors = character 2= hadley romain

packages %>% purrr::sort_by('star') %>% rev()
# 略

order のように結果をインデックスで得るためには rlist::list.order もしくは purrr::order_by

rlist::list.order(packages, star)
## [1] 1 3 2

purrr::order_by(packages, 'star')
## [1] 1 3 2

まとめ

前回と同じく、{purrr} では map (flatmap)keep (discard) でだいたいのことはできる。処理を少し変えたい場合は 組み込み関数を適当に組み合わせればよいため、追加での学習コストは低いと思う。

一方、{rlist} はパッケージの中だけで処理が完結するので、そちらにメリットを感じる方は {rlist} を使ったほうがよさそう。

さらに続きます。

R言語徹底解説

R言語徹底解説

{purrr} でリストデータを操作する <1>

R で関数型プログラミングを行うためのパッケージである {purrr}、すこし使い方がわかってきたので整理をしたい。RStudio のブログの記載をみると、とくにデータ処理フローを関数型のように記述することが目的のようだ。

purrr 0.1.0 | RStudio Blog

The core of purrr is a set of functions for manipulating vectors (atomic vectors, lists, and data frames). The goal is similar to dplyr: help you tackle the most common 90% of data manipulation challenges.

ここでいう"関数型プログラミング言語"とは Haskell のような静的型の言語を想定しており、型チェック、ガード、リフトなど断片的に影響を受けたと思われる関数群をもつ。ただし、同時に Haskell を目指すものではない、とも明言されている。

そもそも R は Scheme の影響をうけているためか、関数型らしい言語仕様やビルトイン関数群 Map, Reduce, Filter などをすでに持っている。ただ、これらの関数群は引数の順序からパイプ演算子と組み合わせて使いにくい。{purrr} でこのあたりが使いやすくなっているとうれしい。

R: Common Higher-Order Functions in Functional Programming...

※ 上記の Map など、S のリファレンスに見当たらなかったため R 独自だと思っていますが、間違っていたら教えてください。

{purrr} によるリストの操作

{purrr} を使うとリストやベクトルの処理が書けるのかをまとめたい。リスト操作のためのパッケージである {rlist} とできること、記法をくらべてみる。

{purrr} を利用する目的は 処理を簡単に、見やすく書くことであるため、{rlist} と比べて可読性が悪いもの、複雑になるものは "できない" 処理として扱う。サンプルデータとして、R の著名パッケージに関連した情報のリストを用意した。

packages <- list(
  list(name = 'dplyr', star = 979L, maintainer = 'hadley' , authors = c('hadley', 'romain')),
  list(name = 'ggplot2', star = 1546L, maintainer = 'hadley' , authors = c('hadley')),
  list(name = 'knitr', star = 1047L, maintainer = 'yihui' , authors = c('yihui', 'hadley', '...and all'))
)
packages
##  x = list 3 (2632 bytes)
## .  [[1]] = list 4
## . .  name = character 1= dplyr 
## . .  star = integer 1= 979
## . .  maintainer = character 1= hadley 
## . .  authors = character 2= hadley romain 
## .  [[2]] = list 4
## . .  name = character 1= ggplot2 
## . .  star = integer 1= 1546
## . .  maintainer = character 1= hadley 
## . .  authors = character 1= hadley 
## .  [[3]] = list 4
## . .  name = character 1= knitr 
## . .  star = integer 1= 1047
## . .  maintainer = character 1= yihui 
## . .  authors = character 3= yihui hadley  ...

1レベル目がレコード、以降のレベルが各レコードの要素となっている。この記事では主に以下二つの操作に関連した機能について記載する。

  • 要素に対する操作: 特定の要素を選択する。特定の要素を元に新しい要素をつくる (マッピングする)。
  • レコードに対する操作: 特定のレコードを選択する。

準備: リストの表示

既定の print ではリストの階層構造がわかりにくく、出力も長い。わかりやすくするため Hmisc::list.tree を利用した出力を記載する。また、出力が同一のものは省略する。

library(Hmisc)
l1 <- list(a = 1, b = c(3L, 4L), c = list(x = c(5, 6), y = 'AAA'))

# 既定での表示
l1
## $a
## [1] 1
## 
## $b
## [1] 3 4
## 
## $c
## $c$x
## [1] 5 6
## 
## $c$y
## [1] "AAA"

# list.tree での表示 (名前 = 型 要素数 = 値 と読む)
Hmisc::list.tree(l1)
## l1 = list 3 (968 bytes)
## .  a = double 1= 1
## .  b = integer 2= 3 4
## .  c = list 2
## . .  x = double 2= 5 6
## . .  y = character 1= AAA 

パッケージのロード

library(rlist)
library(pipeR)
library(purrr)
library(dplyr)

要素の選択とマッピング

rlist Tutorial: Mapping に相当する内容。

要素の選択

R 標準の lapplyapply 系の関数に対応するのは rlist::list.mappurrr::mappurrr::map では第二引数に文字列を渡すとその要素の抽出、ラムダ式 ( ~ .$name )を渡すとその式の適用になる。

lapply(packages, function(x) { x$name })
##  x = list 3 (360 bytes)
## .  [[1]] = character 1= dplyr 
## .  [[2]] = character 1= ggplot2 
## .  [[3]] = character 1= knitr

rlist::list.map(packages, name)
# 略

purrr::map(packages, 'name')
# 略

purrr::map(packages, ~ .$name)
# 略

複数の要素を一度に選択したい場合、{rlist} には rlist::list.select という選択専用の関数がある。{purrr} では渡せる式が一つのため、ラムダ式を利用する。

rlist::list.select(packages, maintainer, star)
##  x = list 3 (1488 bytes)
## .  [[1]] = list 2
## . .  maintainer = character 1= hadley 
## . .  star = integer 1= 979
## .  [[2]] = list 2
## . .  maintainer = character 1= hadley 
## . .  star = integer 1= 1546
## .  [[3]] = list 2
## . .  maintainer = character 1= yihui 
## . .  star = integer 1= 1047

purrr::map(packages, ~ .[c('maintainer', 'star')])
# 略

補足 第二引数にベクトルを渡した場合の処理は、以下のように階層的な選択になる ( l[[c('a', 'b')]] と一緒)。

l2 <- list(list(a=list(b=1)), list(a=list(b=2)))
l2
##  x = list 2 (1176 bytes)
## .  [[1]] = list 1
## . .  a = list 1
## . . .  b = double 1= 1
## .  [[2]] = list 1
## . .  a = list 1
## . . .  b = double 1= 2

purrr::map(l2, c('a', 'b'))
##  x = list 2 (152 bytes)
## .  [[1]] = double 1= 1
## .  [[2]] = double 1= 2

{rlist}, {purrr} でのラムダ式

f <- function(.) {. + 1} に対応する無名関数を {rlist}, {purrr} それぞれの記法で書く。形式はチルダの有無以外は同じ。ドットが引数に対応する。

nums <- c(a = 3, b = 2, c = 1)
rlist::list.map(nums, . + 1)
##  x = list 3 (544 bytes)
## .  a = double 1= 4
## .  b = double 1= 3
## .  c = double 1= 2

purrr::map(nums, ~ . + 1)
# 略

{rilst} では ラムダ式中の .i で元のリストの位置 (インデックス) を、.name で名前 (names) をそれぞれ参照できる。purrrラムダ式には対応する記法がないため、近いことをやるためには直接 map の引数として渡す必要がある。

rlist::list.map(nums, .i) 
##  x = list 3 (544 bytes)
## .  a = integer 1= 1
## .  b = integer 1= 2
## .  c = integer 1= 3

purrr::map(seq_along(nums), ~ .)  # namesは変わってしまう
##  x = list 3 (216 bytes)
## .  [[1]] = integer 1= 1
## .  [[2]] = integer 1= 2
## .  [[3]] = integer 1= 3
rlist::list.map(nums, paste0("Name: ", .name))
##  x = list 3 (688 bytes)
## .  a = character 1= Name: a 
## .  b = character 1= Name: b 
## .  c = character 1= Name: c

purrr::map(names(nums), ~ paste0("Name: ", .)) # namesは変わってしまう
##  x = list 3 (360 bytes)
## .  [[1]] = character 1= Name: a 
## .  [[2]] = character 1= Name: b 
## .  [[3]] = character 1= Name: c

また、内部的にはラムダ式eval で評価をしているため、変数として扱われない位置にあるドットは評価されない。

purrr::map(nums, ~ list(. = 1))
##  x = list 3 (1312 bytes)
## .  a = list 1
## . .  . = double 1= 1
## .  b = list 1
## . .  . = double 1= 1
## .  c = list 1
## . .  . = double 1= 1

要素の追加、変更

元となるリストの値から新しいリストを作りたい場合はラムダ式でリストを返す。

rlist::list.map(packages, list(star = star, had = 'hadley' %in% authors))
##  x = list 3 (1320 bytes)
## .  [[1]] = list 2
## . .  star = integer 1= 979
## . .  had = logical 1= TRUE
## .  [[2]] = list 2
## . .  star = integer 1= 1546
## . .  had = logical 1= TRUE
## .  [[3]] = list 2
## . .  star = integer 1= 1047
## . .  had = logical 1= TRUE

purrr::map(packages, ~ list(star = .$star, had = 'hadley' %in% .$authors))
# 略

結果の型の変更

結果をベクトルで取得したい場合、{rilst} では list.mapv{purrr} では flatmap もしくは map_int ( map_xxx のように結果の型を指定できる関数群がある )。

rlist::list.mapv(packages, star)
## [1]  979 1546 1047

purrr::flatmap(packages, 'star')
# 略

purrr::map_int(packages, 'star')
# 略

data.frame としたい場合は rlist::list.stack{purrr} では dplyr::bind_rows に渡す。

packages %>>%
  rlist::list.select(name, star) %>>%
  rlist::list.stack()
##      name star
## 1   dplyr  979
## 2 ggplot2 1546
## 3   knitr 1047

packages %>%
  purrr::map(~ .[c('name', 'star')]) %>%
  dplyr::bind_rows())
# 略

関数の適用

返り値を変更せずに関数を適用するには rlist::list.iterpurrr::walk。パイプ処理の間に出力やプロットなど意図する返り値を持たない処理を挟み込む場合に利用する。

r <- rlist::list.iter(packages, cat(name, ":", star, "\n"))
## dplyr : 979 
## ggplot2 : 1546 
## knitr : 1047

r <- purrr::walk(packages, ~ cat(.$name, ":", .$star, "\n"))
# 略

レコードの選択

rlist Tutorial: Filtering 前半の単純な条件での選択。

{rilst} では list.filter{purrr} では keepラムダ式が使えるのは map などと同じ。

packages %>>%
  list.filter(star >= 1500)
##  x = list 1 (840 bytes)
## .  [[1]] = list 4
## . .  name = character 1= ggplot2 
## . .  star = integer 1= 1546
## . .  maintainer = character 1= hadley 
## . .  authors = character 1= hadley

packages %>% 
   purrr::keep(~ .$star >= 1500)
# 略
packages %>>%
  list.filter("yihui" %in% authors)
##  x = list 1 (968 bytes)
## .  [[1]] = list 4
## . .  name = character 1= knitr 
## . .  star = integer 1= 1047
## . .  maintainer = character 1= yihui 
## . .  authors = character 3= yihui hadley  ...

packages %>%
   purrr::keep(~ "yihui" %in% .$authors)
# 略

まとめ

{purrr} によるリストデータの属性、レコードに対する操作を記載した。purrr::mappurrr::keep だけでもパイプ演算子 + ラムダ式と組み合わせて幅広い処理ができそうだ。

11/28追記 続きです。

sinhrks.hatenablog.com

R言語徹底解説

R言語徹底解説

Python Jupyter + pandas で DataFrame 表示をカスタマイズする

先日 pandas v0.17.1 がリリースされた。v0.17.0 に対するバグフィックスがメインだが、以下の追加機能もあるため その内容をまとめたい。

HTML 表示のカスタマイズ

Jupyer 上では pandasDataFrame は自動的に HTML として描画される。この HTML に対して、さまざまな CSS を柔軟に設定できるようになった。

このエントリでは、添付した公式ドキュメントとは少し違う例を記載する。

重要 公式ドキュメントにも記載がされているが v0.17.1 時点で開発中 / Experimental な追加のため、今後 破壊的な変更が発生する可能性がある。ご要望やお気づきの点があれば GitHub issue か @ ください。

以降、Jupyter にて実行する。まず DataFrame を作成して表示する。

import pandas as pd
import numpy as np

np.__version__
# '1.10.1'

pd.__version__
# u'0.17.1'

np.random.seed(1)

df = pd.DataFrame({'name': list('abcdefg'),
                   'values1': np.random.randn(7), 
                   'values2': np.random.randn(7)})
df

f:id:sinhrks:20151122200445p:plain

v0.17.0 にて DataFrameDataFrame.style アクセサが追加され、これを経由して種々のカスタマイズができる。例えば、各セルの値に応じてカラーマップを適用するには Styler.background_gradient()

type(df.style)
# pandas.core.style.Styler

df.style.background_gradient(cmap='winter')

f:id:sinhrks:20151122200539p:plain

カラーマップの一部範囲の色のみ利用する場合は low, high キーワードで範囲を指定する。

df.style.background_gradient(cmap='winter', low=2.)

f:id:sinhrks:20151122200546p:plain

特定の列にのみ Styler を適用する場合は subset キーワード。

df.style.background_gradient(cmap='winter', low=2., subset=['values1'])

f:id:sinhrks:20151122200646p:plain

さらに特定のセルにのみ適用する場合は 対象セルを指定する pd.IndexSlice インスタンスsubset として指定する。

pd.IndexSlice[2:5, ['values1']]
# (slice(2, 5, None), ['values1'])

# .loc に渡した場合は対象セルのみをスライス
df.loc[pd.IndexSlice[2:5, ['values1']]]

f:id:sinhrks:20151122201309p:plain

df.style.background_gradient(cmap='winter', low=2., subset=pd.IndexSlice[2:5, ['values1']])

f:id:sinhrks:20151122200655p:plain

Stylerbackground_gradient 以外にも 以下のようなメソッドをもつ。

df.iloc[1, 1] = np.nan
df.style.highlight_null()

f:id:sinhrks:20151122200712p:plain

さらに、Styler.set_properties を利用すると任意の CSS を全セルに対して適用できる。数値以外のセルを色分けするにはこれを使えばよい。

df.style.set_properties(**{'background-color': '#00ff7f'})

f:id:sinhrks:20151122200741p:plain

また、Styler はチェインできる。

(df.style.
 set_properties(**{'background-color': '#00ff7f'}).
 background_gradient(cmap='winter', low=2.).
 highlight_null())

f:id:sinhrks:20151122200724p:plain

より柔軟な条件分けのため、 Styler.apply 系のメソッドも持つ。使い方は通常の DataFrame.apply.applymap と同じだが、値そのものではなく CSS を返す関数を適用する。

  • Styler.applymap: 各セルに対して CSS を返す関数を適用 (関数への入力はスカラー)。
  • Styler.apply: 各列もしくは各行に対して CSS を返す関数を適用 (関数への入力は 各列/各行の値からなる Series )。

参考 Python pandas データのイテレーションと関数適用、pipe - StatsFragments

セルの値が 0 より小さい場合は 赤 / 太字で表示する場合は Styler.applymap

def highlight_negative(val):
    if val < 0:
        return 'color: {0}; font-weight: bold'.format('red')
    else:
        return 'color: {0}'.format('black')

df.style.applymap(highlight_negative)

f:id:sinhrks:20151122200800p:plain

各行について、values1 列の値が values2 列の値より大きいときに values1 列のみ 赤背景で表示したい。行ごとに関数を適用するため、axis=1 を指定する。highlight_values1 関数が返す リストの各要素が、各列に適用される CSS に対応している。

def highlight_values1(s):
    if s['values1'] > s['values2']:
        # 
        return ['', 'background-color: red', '']
    else:
        # スタイル変更しない場合は空の文字列のリストを返す
        return [''] * len(s)

df.style.apply(highlight_values1, axis=1)

f:id:sinhrks:20151122200808p:plain

最後に、Jupyterwidget と組み合わせることで動的なスタイル変更ができる。以下ではスライドバーの位置によって フォントの大きさを変更している。

from IPython.html import widgets
@widgets.interact
def f(s=(5, 30)):
    return (df.style.background_gradient('winter', low=2.).set_properties(**{'font-size': '{0}px'.format(s)}))

f:id:sinhrks:20151122200816p:plain

f:id:sinhrks:20151122200825p:plain

Jupyter 以外に 同様の出力を HTML として埋め込みたい場合は、Styler.render() を利用する。

style = df.style.background_gradient('winter', low=2.)
style.render()
# 略 (HTML が文字列として出力される)

ピボットテーブル + グラフ表示

また Jupyter 上での DataFrame の扱いに関連して、データそのものを GUIインタラクティブに操作 / グラフ表示ができるpivottablejs というパッケージがある。

まとめ

pandas v0.17.1 で追加された DataFrame の HTML 表示のカスタマイズ方法について記載した。

Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理

Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理

R パッケージを CRAN で公開する

少し前に 自作パッケージを CRAN で公開したのだが ブログに書くのを忘れていた。CRAN 公開時の注意点に関して、日本語の説明があまりない / 情報が古いので簡単にまとめたい。

パッケージの作成

この資料を読みましょう。

www.slideshare.net

継続的インテグレーション (CI)

Travis CI は R をサポート (community supportだが) しているため、.travis.yml に2行記載するだけで利用できる。CI 上でパッケージのチェック (R CMD check) も走るので利用したほうが楽。

複数の環境でテストを実行したい場合、Travis CI では Build Matrix という仕組みがある。R では実行したいテストの数だけ env: タグ中に環境変数を記載するのがよいと思う。自分のパッケージでは {ggplot2} の公式版 / GitHub 上の開発版で 2 つのテストを実行している。

そのほかの CI ツール

ほか、R から CI を利用するツールとして、以下のようなものもある。

vignettes の作成

11/30追記 {knitr} で作るのがカンタン。vignettes 化したい {knitr}ドキュメント中に以下のようなタグを含めればよい。VignetteIndexEntry は vignettes へのリンクテキストとして利用される。ASCII文字以外を含める場合、DESCRIPTION ファイル中で Encoding: UTF-8 を指定する。

<!--
%\VignetteEngine{knitr::knitr}
%\VignetteIndexEntry{Introduction to ggfortify package}
-->

参考 Package vignettes: How to build package vignettes with knitr | knitr

パッケージのチェック

開発版 (R-devel) の日次ビルド を利用して、R CMD Check --as-cran package_name.tar.gz で行う。開発版のほうが警告が厳しく出るようだ。

自分は 開発版環境を AWS 上に構築している。 Amazon Linux であれば以下のスクリプトでインストールできると思う。

R CMD check では ERROR, WARNING はもちろん、NOTE も "checking CRAN incoming feasibility" 1 つを除いて全て表示されないようにするのが望ましい。

特に、{dplyr} 等で NSE を使っている場合は xxx: no visible binding for global variable 'xxx' といった NOTE が表示されるが、これらについても CRAN 管理者から指摘される場合がある。 SE を使うか、関数中に (ダミーとして) 同名の変数を作るなどで回避しておくのがよい。

パッケージのサブミット

FTP で直接あげてはダメ。Submit package to CRAN に必要事項を記載して Upload package をする。その後、確認のメールが来たら Confirmation する。

公開待ちのパッケージは 以下の FTP サーバにて確認できる。

パッケージによっては 拡張子の後に以下のようなステータス情報が付けられる場合がある。これは CRAN 公式の情報ではなく、いくつかのフォーラム文書をみて自分が推測したものなので間違っているかもしれない。

  • .save: Upload されたが Confirmation されていないもの。
  • .pending: CRAN 管理者にて何らかの理由でペンディングされているもの (放っておくしかない)。
  • .noemail: Confirmation はされているが確認時のメールが CRAN に飛んでいないもの。

参考 What are the steps in submitting an R-package to CRAN and how long does each step take? - Stack Overflow (ほぼ情報ない)

その後の対応

CRAN 管理者から指摘がきたら適切に対応してください。

パッケージの公開

11/30追記 CRAN 上でのテストとレビューに問題がなければ CRAN 公開後にメールがくる。その後 1日程度すると Windows 用のバイナリのビルドが自動的に行われ、改めて通知がくる。

Rust で階層的クラスタリング

こちらの続き。

階層的クラスタリングについてはこちらを。サンプルデータもリンク先のものを利用させていただく。

再帰的なデータ構造の作成

階層的クラスタリングでは もっとも距離の小さいクラスタ同士を順に併合していく。そのため併合元となるクラスタ (子) と 作成される新しいクラスタ (親) に親子関係をもたせたい。

こういったとき、親から子に参照を持たせればいいかな? と思うのだが、Rust でこれをやろうとすると "親のインスタンスが存在している間は子への参照が有効となっている" ということを lifetime として記述しなければならず、すこし面倒である。

そのため、参照ではなく所有権を移してしまうのが一般的なようだ。

参考 The idea behind a recursive struct - help - The Rust Programming Language Forum

階層的クラスタリング

場合分けをシンプルにするため、ほぼ同一のロジックである最近隣法、最遠隣法、群平均法を書いた。

extern crate csv;
extern crate nalgebra;
extern crate num;
extern crate rand;

use csv::Reader;
use nalgebra::{DVec, DMat, RowSlice, Iterable};
use num::{Signed, Zero, Float};
use std::f64;
use std::io;
use std::ops::Index;
use std::str;
use std::str::FromStr;

fn main() {

    let data = "名前,算数,理科,国語,英語,社会
田中,89,90,67,46,50
佐藤,57,70,80,85,90
鈴木,80,90,35,40,50
本田,40,60,50,45,55
川端,78,85,45,55,60
吉野,55,65,80,75,85
斉藤,90,85,88,92,95";

    // http://burntsushi.net/rustdoc/csv/
    let mut reader = csv::Reader::from_string(data).has_headers(true);
    // http://nalgebra.org/doc/nalgebra/struct.DMat.html
    let dx = read_csv_f64(&mut reader);

    println!("最近隣法");
    let mut hclust = HClust::new(ClusterDistance::Single);
    hclust.fit(&dx);

    println!("最遠隣法");
    let mut hclust = HClust::new(ClusterDistance::Complete);
    hclust.fit(&dx);

    println!("群平均法");
    let mut hclust = HClust::new(ClusterDistance::Average);
    hclust.fit(&dx);
}


enum ClusterDistance {
    Single,                 // 最近隣法
    Complete,               // 最遠隣法
    Average,                // 群平均法
}

struct HClust {
    dist_mat: DMat<f64>,        // 距離行列
    method: ClusterDistance,    // クラスタ間距離の計算方法
    clusters: Vec<Cluster>
}

impl HClust {

    fn new(method: ClusterDistance) -> HClust {
        HClust {
            // dummy
            dist_mat: DMat::from_elem(1, 1, 1.),
            method: method,
            clusters: vec![]
        }
    }

    fn fit(&mut self, data: &DMat<f64>) {
        self.dist_mat = self.get_dist_matrix(&data);

        // クラスタを初期化
        for i in 0..data.nrows() {
            let c = Cluster::from_nodes(vec![i]);
            self.clusters.push(c);
        }
        while self.clusters.len() > 1 {
            self.fit_step();
        }
    }

    /// もっとも距離の小さいクラスタを併合
    fn fit_step(&mut self) {
        let mut tmp_i = 0;
        let mut tmp_j = 0;
        let mut current_dist = f64::MAX;

        for i in 0..self.clusters.len() {
            for j in (i + 1)..self.clusters.len() {
                let d = self.get_cluster_dist(&self.clusters[i], &self.clusters[j]);
                if d < current_dist {
                    current_dist = d;
                    tmp_i = i;
                    tmp_j = j;
                }
            }
        }
        // R と比較するための出力を表示
        println!("Mergeing clusters: {:?} and {:?}: Distance {:?}",
                 &self.clusters[tmp_i].nodes,
                 &self.clusters[tmp_j].nodes,
                 (current_dist * 1e+6).round() / 1e+6);

        let mut new_clusters: Vec<Cluster> = vec![];
        for (i, n) in self.clusters.iter().enumerate() {
            if i != tmp_i && i != tmp_j {
                let n2 = Cluster::from_nodes(n.nodes.clone());
                new_clusters.push(n2);
            }
        }
        // Vec から要素を取り出し & 所有権を取得
        // コンストラクタに所有権ごと渡す
        let new = Cluster::from_clusters(self.clusters.swap_remove(tmp_j),
                                         self.clusters.swap_remove(tmp_i));

        new_clusters.push(new);
        self.clusters = new_clusters;
    }

    /// クラスタに含まれるノード間の距離を計算
    fn get_cluster_dist(&self, c1: &Cluster, c2: &Cluster) -> f64 {
        let mut dists: Vec<f64> = vec![];
        for i in &c1.nodes {
            for j in &c2.nodes {
                dists.push(self.get_node_dist(*i, *j));
            }
        }
        return self.fold_dist_vec(dists);
    }

    /// クラスタに含まれるノード間の距離から、クラスタ間の距離を計算
    fn fold_dist_vec(&self, dists: Vec<f64>) -> f64 {
        match self.method {
            ClusterDistance::Single => {
                return dists.iter().fold(f64::MAX, |a, b| a.min(*b));
            },
            ClusterDistance::Complete => {
                return dists.iter().fold(f64::MIN, |a, b| a.max(*b));
            },
            ClusterDistance::Average => {
                return dists.iter().fold(0., |a, b| a + b) / (dists.len() as f64);
            },
        }
    }

    /// 距離行列を作成
    fn get_dist_matrix(&self, data: &DMat<f64>) -> DMat<f64> {
        // 各列は 0番目 から ノード数-1 番目までのノードに対応 -
        // 各行は 1番目 から ノード数番目 までのノードに対応
        return DMat::from_fn(data.nrows() - 1, data.nrows() - 1,
                             |i, j| if i >= j {
                                euc_dist(&data.row_slice(i + 1, 0, data.ncols()),
                                         &data.row_slice(j, 0, data.ncols()))
                                } else { 0. });
    }

    /// 距離行列を利用して ノード間の距離を求める
    fn get_node_dist(&self, i: usize, j: usize) -> f64 {
        match i > j {
            true => *self.dist_mat.index((i - 1, j)),
            false => *self.dist_mat.index((j - 1, i))
        }
    }
}

struct Cluster {
    nodes: Vec<usize>,
    children: Vec<Cluster>
}

impl Cluster {

    fn from_nodes(nodes: Vec<usize>) -> Cluster {
        Cluster {
            nodes: nodes,
            children: vec![]
        }
    }

    /// 二つのクラスタを併合
    fn from_clusters(left: Cluster, right: Cluster) -> Cluster {
        let mut id = vec![];
        for i in &left.nodes {
            id.push(*i);
        }
        for j in &right.nodes {
            id.push(*j);
        }
        Cluster {
            nodes: id,
            children: vec![left, right]
        }
    }
}

// ユークリッド距離
fn euc_dist<T: Float + Signed>(vec1: &DVec<T>, vec2: &DVec<T>) -> T {
    let mut val: T = Zero::zero();
    for (v1, v2) in vec1.iter().zip(vec2.iter()) {
        val = val + num::pow(num::abs(*v1 - *v2), 2);
    }
    return val.sqrt();
}

/// csv::Reder から f64 に変換できるカラムのみ読み込み
pub fn read_csv_f64<R: io::Read>(reader: &mut Reader<R>) -> DMat<f64> {
    

    let mut x:Vec<f64> = vec![];
    let mut nrows: usize = 0;

    for record in reader.byte_records().map(|r| r.unwrap()) {
        // f64 に変換できる列のみ読み込み
        for item in record.iter().map(|i| str::from_utf8(i).unwrap()) {
            match f64::from_str(item) {
                Ok(v) => x.push(v),
                Err(_) => {}
            };
        }
        nrows += 1;
    }
    let ncols = x.len() / nrows;

    // http://nalgebra.org/doc/nalgebra/struct.DMat.html
    return DMat::from_row_vec(nrows, ncols, &x);
}

出力。

# 最近隣法
# Mergeing clusters: [1] and [5]: Distance 12.409674
# Mergeing clusters: [2] and [4]: Distance 21.307276
# Mergeing clusters: [0] and [4, 2]: Distance 28.478062
# Mergeing clusters: [6] and [5, 1]: Distance 38.105118
# Mergeing clusters: [3] and [4, 2, 0]: Distance 47.106263
# Mergeing clusters: [5, 1, 6] and [4, 2, 0, 3]: Distance 54.313902

# 最遠隣法
# Mergeing clusters: [1] and [5]: Distance 12.409674
# Mergeing clusters: [2] and [4]: Distance 21.307276
# Mergeing clusters: [0] and [4, 2]: Distance 33.778692
# Mergeing clusters: [6] and [5, 1]: Distance 45.585085
# Mergeing clusters: [3] and [4, 2, 0]: Distance 60.133186
# Mergeing clusters: [5, 1, 6] and [4, 2, 0, 3]: Distance 91.531415

# 群平均法
# Mergeing clusters: [1] and [5]: Distance 12.409674
# Mergeing clusters: [2] and [4]: Distance 21.307276
# Mergeing clusters: [0] and [4, 2]: Distance 31.128377
# Mergeing clusters: [6] and [5, 1]: Distance 41.845102
# Mergeing clusters: [3] and [4, 2, 0]: Distance 53.305906
# Mergeing clusters: [5, 1, 6] and [4, 2, 0, 3]: Distance 69.922956

R の結果。コメントは ノードが 0 から始まるとして (上の出力にあわせて ) 記載。作成されるクラスタはどの方法でも同じだが、クラスタ間の距離 (3 列目) は異なる。

# 最近隣法
hc <- hclust(seiseki.d, method = "single")
cbind(hc$merge, hc$height)
#      [,1] [,2]     [,3]
# [1,]   -2   -6 12.40967   # ノード 1 と ノード 5 を併合
# [2,]   -3   -5 21.30728   # ノード 2 と ノード 4 を併合
# [3,]   -1    2 28.47806   # ノード 0 と ノード [2, 4] からなるクラスタを併合 
# [4,]   -7    1 38.10512   # ノード 6 と ノード [1, 5] からなるクラスタを併合
# [5,]   -4    3 47.10626   # ノード 3 と ノード [0, 2, 4] からなるクラスタを併合
# [6,]    4    5 54.31390   # ノード [1, 5, 6] と ノード [0, 2, 3, 4] からなるクラスタを併合

# 最遠隣法
hc <- hclust(seiseki.d, method = "complete")
cbind(hc$merge, hc$height)
#      [,1] [,2]     [,3]
# [1,]   -2   -6 12.40967   # ノード 1 と ノード 5 を併合
# [2,]   -3   -5 21.30728   # ノード 2 と ノード 4 を併合
# [3,]   -1    2 33.77869   # ノード 0 と ノード [2, 4] からなるクラスタを併合 
# [4,]   -7    1 45.58509   # ノード 6 と ノード [1, 5] からなるクラスタを併合
# [5,]   -4    3 60.13319   # ノード 3 と ノード [0, 2, 4] からなるクラスタを併合
# [6,]    4    5 91.53142   # ノード [1, 5, 6] と ノード [0, 2, 3, 4] からなるクラスタを併合

# 群平均法
hc <- hclust(seiseki.d, method = "average")
cbind(hc$merge, hc$height)
#      [,1] [,2]     [,3]
# [1,]   -2   -6 12.40967   # ノード 1 と ノード 5 を併合
# [2,]   -3   -5 21.30728   # ノード 2 と ノード 4 を併合
# [3,]   -1    2 31.12838   # ノード 0 と ノード [2, 4] からなるクラスタを併合 
# [4,]   -7    1 41.84510   # ノード 6 と ノード [1, 5] からなるクラスタを併合
# [5,]   -4    3 53.30591   # ノード 3 と ノード [0, 2, 4] からなるクラスタを併合
# [6,]    4    5 69.92296   # ノード [1, 5, 6] と ノード [0, 2, 3, 4] からなるクラスタを併合

2016/12/20 追記 上記以外の手法を追加 / Lance-Williams formulation で書き直した。

github.com

Rust で主成分分析

Rust で重回帰に続き、今日は 主成分分析をやりたい。

  • Rust(Nightly) rustc 1.6.0-nightly (d5fde83ae 2015-11-12)
  • rust-csv 0.14.3 : CSV の読み込みに利用
  • nalgebra 0.3.2 : 行列/ベクトルの処理に利用

rust-csv での CSV ファイルの読み込み

今回は ローカルのCSV ファイル (iris.csv) を読みとる処理を加える。CSV が単一の型しか含まない場合、 前の記事 のように読み取った値をそのまま Vec に追加していけばよい。

iris のように複数の型を含むデータから特定の型のみを抜き出すには以下のような処理を書く必要がある。

まず Reader.byte_records() で各行のエントリをバイト列として読み取る。

let mut reader = csv::Reader::from_file(path).unwrap().has_headers(false);
for record in reader.byte_records().map(|r| r.unwrap()) {
    println!("{:?}", record);
}
// [[53, 46, 49], [51, 46, 53], [49, 46, 52], [48, 46, 50], [73, 114, 105, 115, 45, 115, 101, 116, 111, 115, 97]]
// [[52, 46, 57], [51, 46, 48], [49, 46, 52], [48, 46, 50], [73, 114, 105, 115, 45, 115, 101, 116, 111, 115, 97]]
// ...

読み取ったバイト列を文字列に変換。

for record in reader.byte_records().map(|r| r.unwrap()) {
    for item in record.iter().map(|i| str::from_utf8(i).unwrap()) {
        println!("{:?}", item);
    }
}
// "5.1"
// "3.5"
// "1.4"
// "0.2"
// "Iris-setosa"
// "4.9"
// ...

文字列を f64 に変換し、成功したもののみ残す。

for record in reader.byte_records().map(|r| r.unwrap()) {
    // f64 に変換できる列のみ読み込み
    for item in record.iter().map(|i| str::from_utf8(i).unwrap()) {
        match f64::from_str(item) {
            Ok(v) => println!("{}", v),
            Err(e) => {}
        };
    }
}
// 5.1
// 3.5
// 1.4
// 0.2
// 4.9
// ...

主成分分析

nalgebra では 行列の固有ベクトル/固有値nalgebra::eigen_qr で計算できるが、この関数は 動的サイズの行列 DMat では利用できない。そのため、データの入力次元にあわせた行列 Mat4 を使った。

現状の nalgebra::DMat で任意の入力次元に対応した処理を書くためには固有値計算を自力でやる必要がある...と思う。

extern crate csv;
extern crate nalgebra;
extern crate num;

use std::f64;
use std::str;
use std::str::FromStr;
use std::vec::Vec;
use nalgebra::{DVec, DMat, Mat4, Mean, ColSlice, Iterable, Transpose};

fn main() {
    // "iris" データを使用
    // http://aima.cs.berkeley.edu/data/iris.csv
    let path = "./data/iris.csv";
    // http://burntsushi.net/rustdoc/csv/
    let mut reader = csv::Reader::from_file(path).unwrap().has_headers(false);

    let mut x:Vec<f64> = vec![];
    let mut nrows: usize = 0;

    for record in reader.byte_records().map(|r| r.unwrap()) {
        // f64 に変換できる列のみ読み込み
        for item in record.iter().map(|i| str::from_utf8(i).unwrap()) {
            match f64::from_str(item) {
                Ok(v) => x.push(v),
                Err(e) => {}
            };
        }
        nrows += 1;
    }

    let ncols = x.len() / nrows;

    // http://nalgebra.org/doc/nalgebra/struct.DMat.html
    let dx = DMat::from_row_vec(nrows, ncols, &x);

    // 主成分分析
    let mut pca = PCA::new(ncols, true);
    pca.fit(&dx);

    // 結果は小数点以下 5 桁に丸めて表示
    println!("主成分 (center=true)\n{:?}", &round(&mut pca.rotation, 5));
    println!("主成分得点\n{:?}", &round(&mut pca.transform(&dx), 5));
}

struct PCA {
    center: bool,          // センタリングするかどうか
    nfeatures: usize,      // 入力の次元
    rotation: DMat<f64>    // 主成分
}

impl PCA {

    fn new(nfeatures: usize, center: bool) -> PCA {
        PCA {
            center: center,
            // 固有値計算に Mat4 を使うため、次元は 4 に固定
            nfeatures: match nfeatures {
                4 => nfeatures,
                _ => panic!("Number of features must be 4")
            },
            rotation: DMat::new_ones(nfeatures, nfeatures)
        }
    }

    fn get_centers(&self, data: &DMat<f64>) -> DVec<f64> {
        // センタリングに用いるベクトルを計算
        return match self.center {
            true => data.mean(),
            false => DVec::from_elem(self.nfeatures, 0.)
        };
    }

    fn fit(&mut self, data: &DMat<f64>) {
        let nrows = data.nrows();
        let centers = self.get_centers(&data);

        // 偏差平方和積和行列
        // Dmat::from_fn で、行番号 i, 列番号 j を引数とする関数から各要素の値を生成できる
        let smx = DMat::from_fn(self.nfeatures, self.nfeatures,
                                |i, j| sum_square(&data.col_slice(i, 0, nrows),
                                                  &data.col_slice(j, 0, nrows),
                                                  centers[i], centers[j]));

        // DMat では eigen_qr が使えないため Mat4 に変換
        let smxv = smx.to_vec();
        let smx4 = Mat4::new(smxv[0], smxv[1], smxv[2], smxv[3],
                             smxv[4], smxv[5], smxv[6], smxv[7],
                             smxv[8], smxv[9], smxv[10], smxv[11],
                             smxv[12], smxv[13], smxv[14], smxv[15]);
        // 固有ベクトル、固有値を計算
        let (evec, eval) = nalgebra::eigen_qr(&smx4, &1e-8, 10000);

        let mut vals: Vec<f64> = vec![];
        for row in evec.transpose().as_array() {
            vals.extend(row);
        }
        self.rotation = DMat::from_row_vec(self.nfeatures, self.nfeatures, &vals);
    }

    fn transform(&mut self, data: &DMat<f64>) -> DMat<f64> {
        // 主成分得点を計算
        let centers = self.get_centers(&data);
        let cdata = DMat::from_fn(data.nrows(), data.ncols(),
                                |i, j| data[(i, j)] - centers[j]);
        return cdata * &self.rotation;
    }
}

fn sum_square(vec1: &DVec<f64>, vec2: &DVec<f64>, m1: f64, m2: f64) -> f64 {
    // 平方和
    let mut val: f64 = 0.;
    for (v1, v2) in vec1.iter().zip(vec2.iter()) {
        val += (v1 - m1) * (v2 - m2);
    }
    return val;
}

fn round(data: &mut DMat<f64>, decimals: usize) -> DMat<f64> {
    // DMat の各要素を指定された有効数字で丸め
    let nrows = data.nrows();
    let ncols = data.ncols();
    let d: f64 = num::pow(10., decimals);
    // round() は常に整数で丸めるため、有効数字の処理を別途行う
    let vals: Vec<f64> = data.as_mut_vec().iter().map(|x| (x * d).round() / d).collect();

    return DMat::from_col_vec(nrows, ncols, &vals);
}

実行結果。

# 主成分 (center=true)
#  0.36159 -0.65654  0.581    0.31725 
# -0.08227 -0.72971 -0.59642 -0.32409 
#  0.85657  0.17577 -0.07252 -0.47972 
#  0.35884  0.07471 -0.54906  0.75112 

# 主成分得点
# -2.68421 -0.32661  0.02151  0.00101 
# -2.71539  0.16956  0.20352  0.0996 
# -2.88982  0.13735 -0.02471  0.0193 
# -2.74644  0.31112 -0.03767 -0.07596 
# -2.72859 -0.33392 -0.09623 -0.06313 
# ...

R の結果。

biris = read.csv("iris.csv", header=FALSE)

# 主成分
prcomp(biris[-5])
# Standard deviations:
# [1] 2.0554417 0.4921825 0.2802212 0.1538929
# 
# Rotation:
#            PC1         PC2         PC3        PC4
# V1  0.36158968 -0.65653988  0.58099728  0.3172545
# V2 -0.08226889 -0.72971237 -0.59641809 -0.3240944
# V3  0.85657211  0.17576740 -0.07252408 -0.4797190
# V4  0.35884393  0.07470647 -0.54906091  0.7511206

# 主成分得点
head(prcomp(biris[-5])$x, n = 5)
#            PC1        PC2         PC3          PC4
# [1,] -2.684207 -0.3266073  0.02151184  0.001006157
# [2,] -2.715391  0.1695568  0.20352143  0.099602424
# [3,] -2.889820  0.1373456 -0.02470924  0.019304543
# [4,] -2.746437  0.3111243 -0.03767198 -0.075955274
# [5,] -2.728593 -0.3339246 -0.09622970 -0.063128733

補足 R 組み込みの iris データと berkeley.edu からダウンロードできる iris データは一部のレコードが異なる。そのため、R の iris データを利用した場合は上と同じ結果にはならない。

iris[apply(iris[-5] != biris[-5], 1, any), ]
#    Sepal.Length Sepal.Width Petal.Length Petal.Width Species
# 35          4.9         3.1          1.5         0.2  setosa
# 38          4.9         3.6          1.4         0.1  setosa

biris[apply(iris[-5] != biris[-5], 1, any), ]
#     V1  V2  V3  V4          V5
# 35 4.9 3.1 1.5 0.1 Iris-setosa
# 38 4.9 3.1 1.5 0.1 Iris-setosa

2016/12/20追記 rulinalg の SVD を使って書き直した。また、上のプログラムでは .transform 時のセンタリングを予測データから行なってしまっているが、これもおかしいので修正。

github.com

Python pandas プロット機能を使いこなす

pandas は可視化のための API を提供しており、折れ線グラフ、棒グラフといった基本的なプロットを簡易な API で利用することができる。一般的な使い方は公式ドキュメントに記載がある。

これらの機能は matplotlib に対する 薄い wrapper によって提供されている。ここでは pandas 側で一処理を加えることによって、ドキュメントに記載されているプロットより少し凝った出力を得る方法を書きたい。

補足 サンプルデータに対する見せ方として不適切なものがあるが、プロットの例ということでご容赦ください。

パッケージのインポート

import matplotlib.pyplot as plt
plt.style.use('ggplot')

import matplotlib as mpl
mpl.__version__
# '1.5.0'

import numpy as np
np.__version__
# '1.10.1'

import pandas as pd
pd.__version__
# u'0.17.0'

折れ線グラフ、棒グラフ系

サンプルデータとして 気象庁から 東京、札幌、福岡の 2014 年の日別平均気温データをダウンロードし、pd.read_csv で読み込んだ。

df = pd.read_csv('data.csv', index_col=0, header=[0, 1, 2], skiprows=[0], encoding='shift-jis')
df = df.iloc[:, [0, 3, 6]]
df.columns = [u'東京', u'札幌', u'福岡']
df.index = pd.to_datetime(df.index)
df.head()

f:id:sinhrks:20151115214421p:plain

まずは単純な折れ線グラフを描く。

df.plot()

f:id:sinhrks:20151115214525p:plain

次に、適当な単位 (ここでは月次) でグルーピングして棒グラフを書きたい。集計には pd.Grouper を利用する。

参考 Python pandas アクセサ / Grouperで少し高度なグルーピング/集計 - StatsFragments

monthly = df.groupby(pd.Grouper(level=0, freq='M')).mean()
monthly

f:id:sinhrks:20151115215001p:plain

このとき、index が 時系列 ( DatetimeIndex ) の場合、そのままでは x 軸がタイムスタンプとして表示されてしまう。

monthly.plot.bar(color=['#348ABD', '#7A68A6', '#A60628'])

f:id:sinhrks:20151115215123p:plain

x軸 を簡単にフォーマットしたい場合、Index.format() メソッドを利用して index を文字列表現に変換するとよい。

monthly.index.format()
# ['2014-01-31', '2014-02-28', '2014-03-31', '2014-04-30', '2014-05-31', '2014-06-30',
#  '2014-07-31', '2014-08-31', '2014-09-30', '2014-10-31', '2014-11-30', '2014-12-31']
monthly.index = monthly.index.format()
monthly.plot.bar(color=['#348ABD', '#7A68A6', '#A60628'])

f:id:sinhrks:20151115214908p:plain

また、各棒の塗り分けには colormap を指定することもできる。

monthly.plot.bar(cmap='rainbow')

f:id:sinhrks:20151115215206p:plain

折れ線グラフ / 棒グラフを一つのプロットとして描画する場合は以下のようにする。.plot メソッドmatplotlib.axes.Axes インスタンスを返すため、続くプロットの描画先として その Axes を指定すればよい。

ax = monthly[u'札幌'].plot(legend=True)
monthly[[u'東京', u'福岡']].plot.bar(ax=ax, rot=30)

f:id:sinhrks:20151115215325p:plain

また、matplotlib 側で 極座標Axes を作成しておけば、レーダーチャートのような描画もできる。

補足 一部 線が切れているのは負値のため。

fig = plt.figure()
ax = fig.add_axes([0.1, 0.1, 0.9, 0.9], polar=True)
indexer = np.arange(1, 13) * 2 * np.pi / 12
monthly.index = indexer
monthly.append(monthly.iloc[0]).plot(ax=ax)
ax.set_xticks(indexer)
ax.set_xticklabels(np.arange(1, 13));

f:id:sinhrks:20151115215334p:plain

分布描画系

各都市の気温の分布が見たい、といった場合にはヒストグラムを描画する。

df.plot.hist(bins=100, alpha=0.5)

f:id:sinhrks:20151115215608p:plain

より細かい単位でグループ分けしてグラフを描きたい場合は、まず plt.subplots で必要な Axes を生成する。続けて、pandas でグループ化したデータを各 Axes に対して描画すればよい。

ここでは月ごとに サブプロットにして カーネル密度推定したグラフを描画してみる。

fig, axes = plt.subplots(4, 3, figsize=(8, 6))
plt.subplots_adjust(wspace=0.5, hspace=0.5)
for (ax, (key, group)) in zip(axes.flatten(), df.groupby(pd.Grouper(level=0, freq='M'))):
    ax = group.plot.kde(ax=ax, legend=False, fontsize=8)
    ax.set_ylabel('')
    ax.set_title(key, fontsize=8)

f:id:sinhrks:20151115215728p:plain

同様に箱ヒゲ図も描ける。

fig, axes = plt.subplots(4, 3, figsize=(8, 6))
plt.subplots_adjust(wspace=0.5, hspace=0.5)
for (ax, (key, group)) in zip(axes.flatten(), df.groupby(pd.Grouper(level=0, freq='M'))):
    ax = group.plot.box(ax=ax)
    ax.set_ylabel('')
    ax.set_title(key, fontsize=8)

f:id:sinhrks:20151115215752p:plain

散布図系

通常の散布図は以下のようにして描ける。

df.plot(kind='scatter', x=u'札幌', y=u'福岡')

f:id:sinhrks:20151115220121p:plain

各点を適当に色分けしたい場合、c キーワードで各点の色を指定する。また、colormap も合わせて指定できる。ここで指定した値は連続値として扱われるため colorbar が表示されている (非表示にすることもできる)。

df.plot(kind='scatter', x=u'札幌', y=u'福岡', c=df.index.month, cmap='winter')

f:id:sinhrks:20151115220154p:plain

各点をグループ別に塗り分ける ( 離散値として扱う ) 場合は、単一の Axes に対して グループ化したデータを順番に描画していけばよい。

cmap = plt.get_cmap('rainbow')
colors = [cmap(c / 12.0) for c in np.arange(1, 13)]
colors
# [(0.33529411764705885, 0.25584277759443558, 0.99164469551074275, 1.0),
#  (0.17058823529411765, 0.49465584339977881, 0.96671840426918743, 1.0),
#  ...
#  (1.0, 0.25584277759443586, 0.12899921653020341, 1.0),
#  (1.0, 1.2246467991473532e-16, 6.123233995736766e-17, 1.0)]

fig, ax = plt.subplots(1, 1)
for i, (key, group) in enumerate(df.groupby(pd.Grouper(level=0, freq='M')), start=1):
    group.plot(kind='scatter', x=u'札幌', y=u'福岡', color=cmap(i / 12.0), ax=ax, label=i)

f:id:sinhrks:20151115220859p:plain

また、ダミーのカラムを作成することで月次の気温の分布を散布図としてみることもできる。

df['month'] = df.index.month
df.plot(kind='scatter', x='month', y=u'札幌')

f:id:sinhrks:20151115220521p:plain

x 軸方向に散らす。

df['month2'] = df.index.month + np.random.normal(loc=0, scale=0.05, size=len(df.index.month))
df.plot(kind='scatter', x='month2', y=u'札幌')

f:id:sinhrks:20151115220531p:plain

バブルチャートも描ける。前処理として札幌の気温を離散化し、月次のデータ数を集計する。

df2 = df.copy()
df2[u'札幌2'] = df2[u'札幌'] // 10 * 10 
df2['count'] = 1
df2 = df2.groupby(['month', u'札幌2'])['count'].count().reset_index()
df2.head()

f:id:sinhrks:20151115220541p:plain

バブルの大きさは s キーワードで指定する。

df2.plot.scatter(x='month', y=u'札幌2', s=df2['count'] * 10)
df.plot(kind='scatter', x='month2', y=u'札幌')

f:id:sinhrks:20151115220549p:plain

さらに変形して plt.imshow でヒートマップを描画する。imshowpandasAPI にはないため、DataFrame.pipe でデータを渡す。

piv = pd.pivot_table(df, index=u'札幌2', columns=df.index.month, values=u'札幌', aggfunc='count')
piv = piv.fillna(0)
piv

f:id:sinhrks:20151115220601p:plain

piv.pipe(plt.imshow, cmap='winter')
ax = plt.gca()
ax.invert_yaxis()
ax.set_yticks(np.arange(4))
ax.set_yticklabels(piv.index);

f:id:sinhrks:20151115220608p:plain

まとめ

pandas での前処理 + 可視化機能の組み合わせを利用して、より柔軟にプロットを行う方法を記載した。pandas の裏側は ndarray のため、最後の例のように pandas 側に API がないプロットも簡単に描ける。

Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理

Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理

Rust で重回帰

簡単なアルゴリズムを実装しつつ Rust を学びたい。環境はこちら。

  • Rust(Nightly) rustc 1.6.0-nightly (d5fde83ae 2015-11-12)
  • rust-csv 0.14.3 : CSV の読み込みに利用
  • nalgebra 0.3.1 : 行列/ベクトルの処理に利用

Python や R と比較すると、型のために Vecslice が透過的に扱えず戸惑う。また、 Generics を使った場合に所有権まわりのエラーが出て解決できなかったため、浮動小数f64 のみの定義とした。もう少し理解できたら書き直したい。

extern crate csv;
extern crate nalgebra;

use std::vec::Vec;
use nalgebra::{DVec, DMat, Inv, Mean, ColSlice, Iterable};

fn main() {

    // R の "trees" データを使用
    let data = "Girth,Height,Volume
8.3,70,10.3
8.6,65,10.3
8.8,63,10.2
10.5,72,16.4
10.7,81,18.8
10.8,83,19.7
11.0,66,15.6
11.0,75,18.2
11.1,80,22.6
11.2,75,19.9
11.3,79,24.2
11.4,76,21.0
11.4,76,21.4
11.7,69,21.3
12.0,75,19.1
12.9,74,22.2
12.9,85,33.8
13.3,86,27.4
13.7,71,25.7
13.8,64,24.9
14.0,78,34.5
14.2,80,31.7
14.5,74,36.3
16.0,72,38.3
16.3,77,42.6
17.3,81,55.4
17.5,82,55.7
17.9,80,58.3
18.0,80,51.5
18.0,80,51.0
20.6,87,77.0";

    let mut ncols = 0;

    // http://burntsushi.net/rustdoc/csv/
    let mut reader = csv::Reader::from_string(data).has_headers(true);

    let mut x: Vec<f64> = vec![];
    for row in reader.decode() {
        let vals: Vec<f64> = row.unwrap();
        // 1 列目を目的変数, 残りは説明変数とする
        for i in 0 .. vals.len() {
            x.push(vals[i]);
        }
        // 説明変数の数 (各行は必ず同数のデータを含むと仮定)
        ncols = vals.len();
    }

    // レコード数
    let nrows: usize = x.len() / ncols;

    // ベクトルから nrows x ncols の DMat (動的サイズの行列) を作成
    // http://nalgebra.org/doc/nalgebra/struct.DMat.html
    let dx = DMat::from_row_vec(nrows, ncols, &x);

    // 重回帰
    let coefs = lm_fit(&dx);
    println!("結果 {:?}", &coefs);
}

fn lm_fit(data: &DMat<f64>) -> Vec<f64> {
    // 重回帰

    // 列ごとの平均値を計算
    let means = data.mean();

    let nrows = data.nrows();
    let nfeatures = data.ncols() - 1;

    // 偏差平方和積和行列
    // DMat::from_fn で、行番号 i, 列番号 j を引数とする関数から各要素の値を生成できる
    let smx = DMat::from_fn(nfeatures, nfeatures,
                            |i, j| sum_square(&data.col_slice(i + 1, 0, nrows),
                                              &data.col_slice(j + 1, 0, nrows),
                                              means[i+1], means[j+1]));
    // 偏差積和行列
    // DMat と Vec では演算ができないため、こちらも一列の DMat として生成
    let smy = DMat::from_fn(nfeatures, 1,
                            |i, j| sum_square(&data.col_slice(i + 1, 0, nrows),
                                              &data.col_slice(0, 0, nrows),
                                              means[i], means[0]));
    // 偏回帰係数を計算し、Vec に変換
    let mut res = (smx.inv().unwrap() * smy).to_vec();

    // 切片を計算し、0 番目の要素として挿入
    let intercept = (1..means.len()).fold(means[0], |m, i| m - res[i-1]*means[i]);
    res.insert(0, intercept);
    return res
}

fn sum_square(vec1: &DVec<f64>, vec2: &DVec<f64>, m1: f64, m2: f64) -> f64 {
    // 平方和
    let mut val: f64 = 0.;
    for (v1, v2) in vec1.iter().zip(vec2.iter()) {
        val += (v1 - m1) * (v2 - m2);
    }
    return val;
}

実行結果。左から切片と各説明変数の偏回帰係数。

# 結果 [10.81637050444831, -0.04548349100904309, 0.19517974893553702]

R の結果。

lm(Girth ~ Height + Volume, data = trees)
# 
# Call:
# lm(formula = Girth ~ Height + Volume, data = trees)
# 
# Coefficients:
# (Intercept)       Height       Volume  
#    10.81637     -0.04548      0.19518

2016/12/26追記 上では逆行列を使ったが、利用するパッケージに線形方程式を解く関数 / メソッドが用意されていればそちらを使った方が数値的に安定 / 高速。

github.com

岩波データサイエンス Vol.1

ご恵贈いただきました。 ありがとうございます! あわせてタスクもいただきました (下部)。

書籍のコンテンツ

各章ごとの内容は id:sfchaos さんが詳しく紹介されています。

d.hatena.ne.jp

まだ すべて読めていないのですが、以下 3 点がよいポイントだと思います。

  • 理論 と サンプルプログラム 両方の記載がある
  • BUGS, Stan, PyMC3 と主要なパッケージが網羅されている
  • サンプルは単純な回帰だけでなく 時系列 / 空間ベイズを含む

補足 書籍には コラム "Pythonとは" という データ分析視点での Python 紹介があるのですが、中身は結構な pandas 推しでした。著者の方、いったい何者なんだ...。

Stan 入門

依頼により、著者の松浦さんが作成した RStan サンプルの PyStan 版を作成させていただきました。

以下リポジトリの "Example..." フォルダ中に含まれる Jupyter Notebook ( .ipynb ) を開くと、 GitHub 上でプログラムと出力を確認することができます。

github.com

RStanPyStanAPI 差異については公式ドキュメントに記載がありますがアップデートがされていません。下表にサンプル中で利用したものを簡単にまとめます。

RStan PyStan 概要
stan(...) pystan.stan(...) モデルのコンパイルとサンプリングの実行
stan_model(...) pystan.StanModel(...) モデルのコンパイル
sampling(...) StanModel.sampling(...) サンプリングの実行
extract(...) StanFit4model.extract() サンプリング結果を取得
traceplot(...) StanFit4model.plot(...) サンプリングの経過をプロット

書籍や Rstan のサンプルと合わせてご参照ください。

岩波データサイエンス Vol.1

岩波データサイエンス Vol.1