Rでsummary()を使わずに四分位数を計算する

### 概要

Rで四分位数を計算するには、summary関数やquantile関数を使えばよいです。ここでは、勉強のために四分位数をRの組み込み関数を使わずに計算してみます。


### はじめに

あるデータの四分位数を知りたい場合、Rではsummary関数やquantile関数を使います。例えば、

```r:R
> x <- rpois(10, 2.5)  # 個数10で平均2.5のポアソン乱数を発生してxに代入せよ
> x
 [1] 2 2 2 3 2 6 1 1 1 3
```

のように、まず10個の整数データをポアソン乱数をrpois関数で作ります。この生成した乱数データxの四分位数は、

```r:R
> summary(x)  # 要約統計量を計算せよ
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1.00    1.25    2.00    2.30    2.75    6.00 
```

のように、summary関数の出力する要約統計量から第一四分位数(1st Quantile)などが分かります。quantile関数を使っても四分位数が計算できます。

```r:R
> quantile(x)  # 四分位数を計算せよ
  0%  25%  50%  75% 100% 
1.00 1.25 2.00 2.75 6.00 
```

このようにRではコマンド1つで四分位数を計算できるのですが、勉強のためにsummary関数やquantile関数を使わずに四分位数を計算してみようと思います。


### 四分位数を順を追って計算する

まず、サンプルデータを用意します。サンプルデータは、標本数50、平均2.8のポアソン乱数の整数データとします。サンプルデータの設定に特に意味はありません。平均などの値は何でもよいです。

```r:R
> x <- rpois(50, 2.8)
> options(width=24)  # 一行の表示文字数を24にせよ
> x
 [1] 4 3 2 1 1 1 4 3 3 0
[11] 0 4 1 4 3 3 4 0 3 2
[21] 2 2 4 0 1 1 2 1 4 1
[31] 2 2 1 6 5 3 4 3 2 5
[41] 1 2 4 3 0 2 4 2 2 0
```

サンプルデータxの要素は0から6までの整数になりました。xをプロットすると下図になります。

```
> plot(x)  # xをプロットせよ
```

![図1 サンプルデータx](https://qiita-image-store.s3.amazonaws.com/0/197109/c8a71248-b882-357f-8fe8-8baefe958e67.png)

横軸はxのIndex、縦軸は各Indexの要素の値です。Index 0がありますが、xのIndexは1から50までです。50個の○が散らばっています。

四分位数は昇順に並べ替えたデータを4等分する値なので、サンプルデータxを昇順にソートします。ソートしたデータをxsとし、xと同様にプロットします。

```r:R
> xs <- sort(x)  # xを昇順にソートせよ
> xs
 [1] 0 0 0 0 0 0 1 1 1 1
[11] 1 1 1 1 1 1 2 2 2 2
[21] 2 2 2 2 2 2 2 2 3 3
[31] 3 3 3 3 3 3 3 4 4 4
[41] 4 4 4 4 4 4 4 5 5 6
> plot(xs)
```
![図2 ソートしたサンプルデータxs](https://qiita-image-store.s3.amazonaws.com/0/197109/97fb6319-a481-ecb1-461e-7d6fdf85670a.png)

0から6までの整数データを昇順にソートしたので、Index 1から50まで階段状に増加しています。

このソートしたデータxsのIndexを4等分する点を求めます。下図は、Index 1からnを4等分した点を示しています。

![図3 標本数nの4等分](https://qiita-image-store.s3.amazonaws.com/0/197109/ced49096-5b69-05dc-38d8-50e114efc173.png)

Index 1を0%、Index nを100%として、Indexを4等分する25%、50%、75%の点をそれぞれ$i_{0.25}$、$i_{0.5}$、$i_{0.75}$とします。4等分は (n - 1)/4 ですので[^10]、$i_{0.25}$、$i_{0.5}$、$i_{0.75}$はそれぞれ以下で表されます[^20]。

[^10]: 0からnを4等分するならばn/4ですが、1からnを4等分するので(n-1)/4となります。
[^20]: 内分点の公式を用いても同じことです。

$$ i_{0.25} = 1 + \frac{n - 1}{4} = \frac{n + 3}{4} $$

$$ i_{0.5} = 1 + 2 \cdot \frac{n - 1}{4} = \frac{n + 1}{2} $$

$$ i_{0.75} = 1 + 3 \cdot \frac{n - 1}{4} = \frac{3 n + 1}{4} $$

これらの式を用いて$i_{0.25}$、$i_{0.5}$、$i_{0.75}$をそれぞれRで計算します。

```r:R
> n <- length(x)  # xの長さ(サンプルデータxの要素数)

> i25 <- (n+3)/4  # 25%点
> i25
[1] 13.25

> i50 <- (n+1)/2  # 50%点
> i50
[1] 25.5

> i75 <- (3*n+1)/4  # 75%点
> i75
[1] 37.75
```

この4等分点からIndexを4等分する線をソートしたxの図に加えると下図になります。

```r:R
> lines(c(1,1), c(-1,7), type="l", col="red")  # (1,-1)から(1,7)に赤い線を引け
> lines(c(i25,i25), c(-1,7), type="l", col="red")
> lines(c(i50,i50), c(-1,7), type="l", col="red")
> lines(c(i75,i75), c(-1,7), type="l", col="red")
> lines(c(n,n), c(-1,7), type="l", col="red")
> text(3, 5.5, "0%", col="red")  # (3, 5.5)に赤い字で25%と書け
> text(i25+2, 5.5, "25%", col="red")
> text(i50+2, 5.5, "50%", col="red")
> text(i75+2, 5.5, "75%", col="red")
> text(48, 5.5, "100%", col="red")
```

![図4 Indexを4等分する線を加えたxsのグラフ](https://qiita-image-store.s3.amazonaws.com/0/197109/67c12379-dbc2-23d6-0a74-52b9b4207acb.png)

サンプルデータxの第一四分位数は、図4の25%線とデータxsが交差する点から求まります。同様に第二四分位数（中央値）は50%線との交点、第三四分位数は75%線との交点から求まります。

第一四分位数を計算します。

```r:R
> i25f <- floor(i25)  # i25 = 13.25 以下の最大の整数を求めよ
> i25f
[1] 13
> qu1 <- xs[i25f] + (xs[i25f + 1] - xs[i25f]) * (i25 - i25f)  # 第一四分位数を計算せよ
> qu1
[1] 1
```

Rコードは変数で計算していますが、今回のサンプルデータのケースでは次の計算をしています。

$$ Q_1 = xs_{13} + (xs_{14} - xs_{13}) \cdot 0.25 $$

計算ををわかりやすくするために拡大した図を示します。

```r:R
> plot(xs, xlim=c(i25-2,i25+2), ylim=c(0,2))
> lines(c(i25,i25),c(-1,7),type="l",col="red")
> lines(c(13,14),c(xs[13],xs[14]),type="l",col="blue")
> points(i25, qu1, pch=19)

> par(xpd=T)  # プロット領域外への描画を許可せよ
> text(locator(1), "1st Qu.")
> text(locator(1), "25%", col="red")
> text(locator(1), i25)
```

![図5 第一四分位数付近のグラフ"](https://qiita-image-store.s3.amazonaws.com/0/197109/07ea26fd-31e2-24e3-148c-949bbdec77bd.png)

Index 1から50の25％は先に求めた$i_{0.25} = 13.25$ですので、Index 13と14の間にあります。赤線とIndex 13と14の要素を線形補間した青線との交点が第一四分位数です。このケースではxs[13]とxs[14]がどちらも1ですので、計算せずとも第一四分位数は$Q_1 = 1$です。

Indexの4等分線を加えた図から第二四分位数(中央値)は2になることが明らかですので、計算を飛ばします。

第一四分位数と同様に、第三四分位数を計算します。

```r:R
> i75f <- floor(i75)  # i75 = 37.75
> i75f
[1] 37

> qu3 <- xs[i75f] + (xs[i75f + 1] - xs[i75f]) * (i75 - i75f)  # 第三四分位数を計算せよ
> qu3
[1] 3.75
```

具体的には次の計算をしています。

$$ Q_3 = xs_{37} + (xs_{38} - xs_{37}) \cdot 0.75 $$

第三四分位数あたりを拡大したグラフを下図に示します。

![図6 第三四分位数付近のグラフ](https://qiita-image-store.s3.amazonaws.com/0/197109/abcdb98e-c8e2-d672-13fb-0c7d04ac6912.png)

Index 1から50の75％は先に求めた$i_{0.75} = 37.75$ですので、Index 37と38の間にあります。赤線とIndex 37と38の要素を線形補間した青線との交点が第三四分位数です。このケースでは$xs[37] = 3$で$xs[38] = 4$で両者に差があります。赤線と青線の交点を計算して、第三四分位数は$Q_3 = 3.75$ です。


### summary関数で答え合わせ

第一四分位数（1st Qu.）と第三四分位数（3rd Qu.）を計算したので、summary関数の出力と比べてみます。

```r:R
> qu1;qu3  # 自分で計算した 1st Qu. と 3rd Qu.
[1] 1
[1] 3.75
> summary(x)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.00    1.00    2.00    2.34    3.75    6.00 
```

両者は一致しています。


### まとめ

- Rでsummary関数やquantile関数などの組み込み関数を使わずに四分位数を計算しました。
- 随時グラフを書きながら、順を追って四分位数を計算しました。
- 最後に計算結果とsummary関数を比較して両者が一致していることを確かめました。