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

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

# 作りたいもの
みどりぼん（データ解析のための統計モデリング入門)にあるグラフ。
![スクリーンショット 2018-10-19 13.17.45.png](https://qiita-image-store.s3.amazonaws.com/0/196857/7ccc2749-c075-8e9e-87f9-441d0938b9dd.png)


(本当の目的)
サンプルコードは既にあるが、よくわからないので一つ一つ理解する。
### 補足
ヒストグラムは生データのようなものであり、生データがポアソン分布に従うと仮定して、平均がいくらなのか、視覚的にみせたい。
そこで、平均を変えながらプロットしようと考えた。
また各平均のポアソン分布がどれぐらい適しているのか、対数尤度を導入して、この値もプロットする。

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

データは[みどりぼん作者のサイト](http://hosho.ees.hokudai.ac.jp/~kubo/ce/IwanamiBook.html)を参照(第2章のdata.RData)。

```r:前準備
# 生データにあたるものをロード
data <-load("data.RData")
```


```r:グラフ作成
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)
```

# 仕組み
```r:グラフ作成(一行目)
logL <- function(m) sum(dpois(data, m, log = TRUE))
```
logLに対数尤度用の関数を定義。対数尤度なので、和を取る。
ポアソン分布

```math
P(X=k)=\frac{\lambda^k e^{-\lambda}}{k!}
```
[wikipedia](https://ja.wikipedia.org/wiki/%E3%83%9D%E3%82%A2%E3%82%BD%E3%83%B3%E5%88%86%E5%B8%83)によると
> P(X = k) は、「所与の時間中に平均で λ 回発生する事象がちょうど k 回（k は非負の整数）発生する確率

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

```math
L(\lambda)\equiv\prod_{i = 0}^{n}P(data[i]\ |\ \lambda) \ (\textrm{R}ならi=1から)
```
と表すことができる。

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

```r:グラフ作成(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 を参照)

<a href="refMatrix"></a>

```r:グラフ作成(layout)
layout(matrix(1:9, byrow = T, ncol = 3))
```

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

```r: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であれば、そこはスキップされるようだ。

```r: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という変数に代入する必要はないんじゃないか、と思うが、代入しなければ

```r:変数に代入しない場合
> sapply(seq(2, 5.2, 0.4), plot.poisson)
[[1]]
NULL

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

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

```r:もとに戻す
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