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)
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)
が、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()
サブプロットの要素には リストと同じようにアクセスできる。[
でのアクセスは ggmultiplot
インスタンスに、[[
での アクセスは ggplot
インスタンスになる。代入もできるため、一部のプロットにだけ処理を適用したい場合は以下のように書ける。
p[2:3] <- p[2:3] + stat_smooth(method = 'lm') p
+
演算子の右項に ggplot
インスタンスを指定するとサブプロットが追加できる。
p + (ggplot(iris, aes(Sepal.Width, fill = Species)) + geom_density())
ggmultiplot
作成時に ncol
もしくは nrow
キーワードを指定することで、サブプロットの列数 / 行数が指定できる。
p <- new('ggmultiplot', plots = list(p1, p2, p3, p4), ncol = 3) p
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")
一行目は今なら {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によるグラフ作成のレシピ集
- 作者: Winston Chang,石井弓美子,河内崇,瀬戸山雅人,古畠敦
- 出版社/メーカー: オライリージャパン
- 発売日: 2013/11/30
- メディア: 大型本
- この商品を含むブログ (2件) を見る
{purrr} でリストデータを操作する <2>
前の記事に続けて、{purrr}
で {rlist}
相当の処理を行う。今回はレコードの選択とソート。
サンプルデータは前回と同じものを利用する。リストの表示も同じく 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}
には対応する関数はないが、組み込みの head
と tail
でよいのでは...。
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.cases
で logical
に変換する。
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.clean
は recursive = 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.map
や purrr::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}
を使ったほうがよさそう。
さらに続きます。
- 作者: Hadley Wickham,石田基広,市川太祐,高柳慎一,福島真太朗
- 出版社/メーカー: 共立出版
- 発売日: 2016/01/23
- メディア: 単行本
- この商品を含むブログ (25件) を見る
{purrr} でリストデータを操作する <1>
R で関数型プログラミングを行うためのパッケージである {purrr}
、すこし使い方がわかってきたので整理をしたい。RStudio のブログの記載をみると、とくにデータ処理フローを関数型のように記述することが目的のようだ。
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 標準の lapply
他 apply
系の関数に対応するのは rlist::list.map
と purrr::map
。purrr::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.iter
と purrr::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::map
と purrr::keep
だけでもパイプ演算子 + ラムダ式と組み合わせて幅広い処理ができそうだ。
11/28追記 続きです。
- 作者: Hadley Wickham,石田基広,市川太祐,高柳慎一,福島真太朗
- 出版社/メーカー: 共立出版
- 発売日: 2016/01/23
- メディア: 単行本
- この商品を含むブログ (25件) を見る
Python Jupyter + pandas で DataFrame 表示をカスタマイズする
先日 pandas
v0.17.1 がリリースされた。v0.17.0 に対するバグフィックスがメインだが、以下の追加機能もあるため その内容をまとめたい。
HTML 表示のカスタマイズ
Jupyer
上では pandas
の DataFrame
は自動的に 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
v0.17.0 にて DataFrame
に DataFrame.style
アクセサが追加され、これを経由して種々のカスタマイズができる。例えば、各セルの値に応じてカラーマップを適用するには Styler.background_gradient()
。
type(df.style) # pandas.core.style.Styler df.style.background_gradient(cmap='winter')
カラーマップの一部範囲の色のみ利用する場合は low
, high
キーワードで範囲を指定する。
df.style.background_gradient(cmap='winter', low=2.)
特定の列にのみ Styler
を適用する場合は subset
キーワード。
df.style.background_gradient(cmap='winter', low=2., subset=['values1'])
さらに特定のセルにのみ適用する場合は 対象セルを指定する pd.IndexSlice
インスタンスを subset
として指定する。
pd.IndexSlice[2:5, ['values1']] # (slice(2, 5, None), ['values1']) # .loc に渡した場合は対象セルのみをスライス df.loc[pd.IndexSlice[2:5, ['values1']]]
df.style.background_gradient(cmap='winter', low=2., subset=pd.IndexSlice[2:5, ['values1']])
Styler
は background_gradient
以外にも 以下のようなメソッドをもつ。
Styler.highlight_max()
: 列もしくは行の最大値を指定して色分け。Styler.highlight_min()
: 列もしくは行の最小値を指定して色分け。Styler.highlight_null()
:NaN
を指定して色分け。Styler.background_gradient()
: 上述。Styler.bar()
: 各セルの値に応じて棒グラフのように背景色表示。
df.iloc[1, 1] = np.nan df.style.highlight_null()
さらに、Styler.set_properties
を利用すると任意の CSS を全セルに対して適用できる。数値以外のセルを色分けするにはこれを使えばよい。
df.style.set_properties(**{'background-color': '#00ff7f'})
また、Styler
はチェインできる。
(df.style. set_properties(**{'background-color': '#00ff7f'}). background_gradient(cmap='winter', low=2.). highlight_null())
より柔軟な条件分けのため、 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)
各行について、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)
最後に、Jupyter
の widget と組み合わせることで動的なスタイル変更ができる。以下ではスライドバーの位置によって フォントの大きさを変更している。
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)}))
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を使ったデータ処理
- 作者: Wes McKinney,小林儀匡,鈴木宏尚,瀬戸山雅人,滝口開資,野上大介
- 出版社/メーカー: オライリージャパン
- 発売日: 2013/12/26
- メディア: 大型本
- この商品を含むブログ (11件) を見る
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 を利用するツールとして、以下のようなものもある。
- r-appveyor: Windows 用の CI 環境である AppVeyor が利用できる。
- r-builder: Travis CI を R から利用する。R-devel が利用できる。
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 で書き直した。
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
時のセンタリングを予測データから行なってしまっているが、これもおかしいので修正。
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()
まずは単純な折れ線グラフを描く。
df.plot()
次に、適当な単位 (ここでは月次) でグルーピングして棒グラフを書きたい。集計には pd.Grouper
を利用する。
参考 Python pandas アクセサ / Grouperで少し高度なグルーピング/集計 - StatsFragments
monthly = df.groupby(pd.Grouper(level=0, freq='M')).mean() monthly
このとき、index
が 時系列 ( DatetimeIndex
) の場合、そのままでは x 軸がタイムスタンプとして表示されてしまう。
monthly.plot.bar(color=['#348ABD', '#7A68A6', '#A60628'])
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'])
また、各棒の塗り分けには colormap
を指定することもできる。
monthly.plot.bar(cmap='rainbow')
折れ線グラフ / 棒グラフを一つのプロットとして描画する場合は以下のようにする。.plot
メソッドは matplotlib.axes.Axes
インスタンスを返すため、続くプロットの描画先として その Axes
を指定すればよい。
ax = monthly[u'札幌'].plot(legend=True) monthly[[u'東京', u'福岡']].plot.bar(ax=ax, rot=30)
また、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));
分布描画系
各都市の気温の分布が見たい、といった場合にはヒストグラムを描画する。
df.plot.hist(bins=100, alpha=0.5)
より細かい単位でグループ分けしてグラフを描きたい場合は、まず 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)
同様に箱ヒゲ図も描ける。
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)
散布図系
通常の散布図は以下のようにして描ける。
df.plot(kind='scatter', x=u'札幌', y=u'福岡')
各点を適当に色分けしたい場合、c
キーワードで各点の色を指定する。また、colormap
も合わせて指定できる。ここで指定した値は連続値として扱われるため colorbar
が表示されている (非表示にすることもできる)。
df.plot(kind='scatter', x=u'札幌', y=u'福岡', c=df.index.month, cmap='winter')
各点をグループ別に塗り分ける ( 離散値として扱う ) 場合は、単一の 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)
また、ダミーのカラムを作成することで月次の気温の分布を散布図としてみることもできる。
df['month'] = df.index.month df.plot(kind='scatter', x='month', y=u'札幌')
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'札幌')
バブルチャートも描ける。前処理として札幌の気温を離散化し、月次のデータ数を集計する。
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()
バブルの大きさは s
キーワードで指定する。
df2.plot.scatter(x='month', y=u'札幌2', s=df2['count'] * 10) df.plot(kind='scatter', x='month2', y=u'札幌')
さらに変形して plt.imshow
でヒートマップを描画する。imshow
は pandas
の API にはないため、DataFrame.pipe
でデータを渡す。
piv = pd.pivot_table(df, index=u'札幌2', columns=df.index.month, values=u'札幌', aggfunc='count') piv = piv.fillna(0) piv
piv.pipe(plt.imshow, cmap='winter') ax = plt.gca() ax.invert_yaxis() ax.set_yticks(np.arange(4)) ax.set_yticklabels(piv.index);
まとめ
pandas
での前処理 + 可視化機能の組み合わせを利用して、より柔軟にプロットを行う方法を記載した。pandas
の裏側は ndarray
のため、最後の例のように pandas
側に API がないプロットも簡単に描ける。
Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理
- 作者: Wes McKinney,小林儀匡,鈴木宏尚,瀬戸山雅人,滝口開資,野上大介
- 出版社/メーカー: オライリージャパン
- 発売日: 2013/12/26
- メディア: 大型本
- この商品を含むブログ (11件) を見る
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 と比較すると、型のために Vec
や slice
が透過的に扱えず戸惑う。また、
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追記 上では逆行列を使ったが、利用するパッケージに線形方程式を解く関数 / メソッドが用意されていればそちらを使った方が数値的に安定 / 高速。
岩波データサイエンス Vol.1
ご恵贈いただきました。 ありがとうございます! あわせてタスクもいただきました (下部)。
書籍のコンテンツ
各章ごとの内容は id:sfchaos さんが詳しく紹介されています。
まだ すべて読めていないのですが、以下 3 点がよいポイントだと思います。
- 理論 と サンプルプログラム 両方の記載がある
BUGS
,Stan
,PyMC3
と主要なパッケージが網羅されている- サンプルは単純な回帰だけでなく 時系列 / 空間ベイズを含む
補足 書籍には コラム "Pythonとは" という データ分析視点での Python 紹介があるのですが、中身は結構な pandas
推しでした。著者の方、いったい何者なんだ...。
Stan 入門
依頼により、著者の松浦さんが作成した RStan
サンプルの PyStan
版を作成させていただきました。
以下リポジトリの "Example..." フォルダ中に含まれる Jupyter Notebook ( .ipynb
) を開くと、 GitHub 上でプログラムと出力を確認することができます。
RStan
と PyStan
の API 差異については公式ドキュメントに記載がありますがアップデートがされていません。下表にサンプル中で利用したものを簡単にまとめます。
RStan |
PyStan |
概要 |
---|---|---|
stan(...) |
pystan.stan(...) |
モデルのコンパイルとサンプリングの実行 |
stan_model(...) |
pystan.StanModel(...) |
モデルのコンパイル |
sampling(...) |
StanModel.sampling(...) |
サンプリングの実行 |
extract(...) |
StanFit4model.extract() |
サンプリング結果を取得 |
traceplot(...) |
StanFit4model.plot(...) |
サンプリングの経過をプロット |
書籍や Rstan
のサンプルと合わせてご参照ください。
- 作者: 岩波データサイエンス刊行委員会
- 出版社/メーカー: 岩波書店
- 発売日: 2015/10/08
- メディア: 単行本(ソフトカバー)
- この商品を含むブログ (4件) を見る