Rで条件を変えながら、3×3のレイアウトに9個のグラフを一括でプロットする。

注

3×3のレイアウトになっているグラフ方法が知りたい場合は記事下部、layoutについて説明している部分以降のみを見て下さい。

作りたいもの

みどりぼん（データ解析のための統計モデリング入門)にあるグラフ。

(本当の目的)
サンプルコードは既にあるが、よくわからないので一つ一つ理解する。

補足

ヒストグラムは生データのようなものであり、生データがポアソン分布に従うと仮定して、平均がいくらなのか、視覚的にみせたい。
そこで、平均を変えながらプロットしようと考えた。
また各平均のポアソン分布がどれぐらい適しているのか、対数尤度を導入して、この値もプロットする。

　コード

http://rpubs.com/kaz_yos/kubobook2 を引用（一部改変）。

データはみどりぼん作者のサイトを参照(第2章のdata.RData)。

前準備

# 生データにあたるものをロード
data <-load("data.RData")

グラフ作成

logL <- function(m) sum(dpois(data, m, log = TRUE))
plot.poisson <- function(lambda) {
    y <- 0:9
    prob <- dpois(y, lambda = lambda)

    hist(data, breaks = seq(-0.5, 9.5, 1), ylim = c(0, 15),
         main = "", xlab = "", ylab = "")
    points(y, prob * 50)
    lines(y,  prob * 50, lty = 2)

    title(sprintf("lambda= %.1f\n logL= %.1f", lambda, logL(lambda)))
}

layout(matrix(1:9, byrow = T, ncol = 3))
junk <- sapply(seq(2, 5.2, 0.4), plot.poisson)

仕組み

グラフ作成(一行目)

logL <- function(m) sum(dpois(data, m, log = TRUE))

logLに対数尤度用の関数を定義。対数尤度なので、和を取る。
ポアソン分布

P(X=k)=\frac{\lambda^k e^{-\lambda}}{k!}

wikipediaによると

P(X = k) は、「所与の時間中に平均で λ 回発生する事象がちょうど k 回（k は非負の整数）発生する確率

今回は$P(X=k)$は、平均で$\lambda$である事象に対して、$k$が観測される確率、ぐらいの解釈で良いと思う。
$\lambda$も考慮して$P(k\ |\ \lambda)$と表記することにする。
尤度はその事象が実現する確率の積なので、$P(k\ |\ \lambda)$を使って

L(\lambda)\equiv\prod_{i = 0}^{n}P(data[i]\ |\ \lambda) \ (\textrm{R}ならi=1から)

と表すことができる。

前置きが非常に長くなったが、この$L(\lambda)$に対して対数尤度を考えていく、というわけだ。
Rのコードに戻る。
logL(1)とすれば、平均1のポアソン分布に対する対数尤度の値が代入される。

グラフ作成(plot.poisson)

plot.poisson <- function(lambda) {
    y <- 0:9
    prob <- dpois(y, lambda = lambda)

    hist(data, breaks = seq(-0.5, 9.5, 1), ylim = c(0, 15),
         main = "", xlab = "", ylab = "")
    points(y, prob * 50)
    lines(y,  prob * 50, lty = 2)

    title(sprintf("lambda= %.1f\n logL= %.1f", lambda, logL(lambda)))
}

この部分は、関数である、ということを除けば(関数については[1])、

• 単にヒストグラムと折れ線をプロットする (ylimに関しては https://qiita.com/7of9/items/21654c15ceee84bd6801 を参照)
• タイトルをつける

の二点に分解される。
プロット部分については、記事(まだ未完成、追記予定)を参照。

タイトルについて[2]。
titleはグラフにタイトルのラベルをつけるための関数。
sprintf()は他の言語にもあるので…。(Rのsprintfについてはhttp://cse.naro.affrc.go.jp/takezawa/r-tips/r/10.html を参照)

グラフ作成(layout)

layout(matrix(1:9, byrow = T, ncol = 3))

layout関数の中身にmatrixを指定すれば、そのmatrixの値の順にプロットされる[2]。
今回のmatrixの中身は以下のようになっている。
コメントで対応する位置の図を書いておく。

matrixの中身(layout)

> print(matrix(1:9, byrow = T, ncol = 3))
     [,1] [,2] [,3]
[1,]    1    2    3 #(1番目の図) (2番目の図) (3番目の図)
[2,]    4    5    6 #(4番目の図) (5番目の図) (6番目の図)
[3,]    7    8    9 #(7番目の図) (8番目の図) (9番目の図)

ちなみに行列の要素が0であれば、そこはスキップされるようだ。

sapply

junk <- sapply(seq(2, 5.2, 0.4), plot.poisson) 

sapplyはapply関数の仲間[3]。
リストに対して関数を実行する。
今回の例では、リストはseq(2, 5.2, 0.4) = (2, 2.4, 2.8, ..., 5.2)、関数はplot.poisson (=　平均値を引数にとり、ポアソン分布とかをプロットする関数)。
つまり、平均値を2, 2.4, 2.8, ..., 5.2と変えながら、ポアソン分布等をプロットすることになる。
これがまさにほしかったもの。
別にjunkという変数に代入する必要はないんじゃないか、と思うが、代入しなければ

変数に代入しない場合

> sapply(seq(2, 5.2, 0.4), plot.poisson)
[[1]]
NULL

[[2]]
NULL
…

となってしまう。取り敢えず代入しなければいけないものだと思っておく。
原因がわかれば追記する。（すみません）

もとに戻す

layoutを上記のように変えたまま別のグラフをプロットすると、一つのグラフが9分割されたうちの一つのスペースにプロットされてしまう。
もとに戻すためには

もとに戻す

layout(1)

とすれば良い。

参考文献

[1]https://stats.biopapyrus.jp/r/basic/function.html
[2]http://cse.naro.affrc.go.jp/takezawa/r-tips/r/55.html
[3]http://cse.naro.affrc.go.jp/takezawa/r-tips/r/24.html