講義（山本裕樹）

この文書は読取専用です。文書のソースを閲覧することは可能ですが、変更はできません。もし変更したい場合は管理者に連絡してください。
====== 【社会情報処理】第4回 基本統計量 ======
__**基本統計量**__とは分布の特徴を表すのに使われる基本的な統計量である。

  * 分布の中心を表す統計量：平均、中央値、最頻値
  * 分布の広がりを表す統計量：分散、標準偏差
  * その他の統計量：最大値、最小値、四分位数

以下の数式では $X_{1},X_{2},\cdots,X_{n}$ の $n$ 個の観測値のデータがあるとする。

R の例として用いる標本データは、データセット airquality の Ozone から $n=21$ 個の標本を抽出したものとする。
<code rsplus>
> n <- 21
> set.seed(12); (data <- sample(na.omit(airquality$Ozone), n))
 [1]  35  73   9  91  27  89   7  47  12   8  30 168  13  32  18  48  18  20 122  45  52
</code>

''set.seed()'' は乱数の種（seed）を設定するためのもので、毎回同じ標本抽出を再現するために使用している。


===== 度数分布表とヒストグラム =====
それぞれの観測値がある範囲（階級）で何回観測されたかをカウントした数を__**度数**__といい、それぞれの階級で上限値と下限値を足して2で割ったものを__**階級値**__という。
度数を全データ数で割ったものを__**相対度数**__という。

階級ごとの度数のデータを__**級別データ**__といい、級別データを表にしたものが__**度数分布表**__である。

度数分布表の例（階級幅は 20）

^  階級番号  ^  階級[ppb]  ^  階級値[ppb]  ^  度数  ^
|  1         |  0-20       |  10           |  7     |
|  2         |  20-40      |  30           |  5     |
|  3         |  40-60      |  50           |  4     |
|  4         |  60-80      |  70           |  1     |
|  5         |  80-100     |  90           |  2     |
|  6         |  100-120    |  110          |  0     |
|  7         |  120-140    |  130          |  1     |
|  8         |  140-160    |  150          |  0     |
|  9         |  160-180    |  170          |  1     |
|  合計      |             |               |  21    |


横軸に階級、縦軸に度数をとって描いた柱状グラフを__**ヒストグラム**__という。

R でヒストグラムを描くには関数 ''hist()'' を使う。
<code rsplus>
> hg <- hist(data)
</code>
戻り値は変数 hg に代入しているが、代入しなくてもヒストグラムは描ける。

デフォルトでは観測値がその階級の下限値より大きい、かつ上限値以下の場合に度数としてカウントする。
\[
下限値<観測値\le 上限値
\]

日本では、観測値がその階級の下限値以上、かつ上限値未満の場合に度数としてカウントするのが一般的である。
\[
下限値\le 観測値<上限値
\]
これを行うには、引数に ''right = FALSE'' を入れる（''breaks'' をベクトルで与える場合は、''include.lowest = TRUE'' が必要かもしれない）。 
<code rsplus>
> hg <- hist(data, right = FALSE)
</code>

階級数はデフォルトでは__**スタージェスの公式**__（$k$ は階級数、$n$ はデータ数）
\[
k\simeq \log_{2}n+1
\]
を使って決められるが、階級の度数次第で微調整されるようである。
スタージェスの公式は関数 ''nclass.Sturges()'' で計算できる。
<code rsplus>
> nclass.Sturges(data)
</code>


^  引数のオプション       ^  機能                                                                                                                              ^
| ''breaks = "Sturges"''  | 階級数、階級数を決める方法、境界値などを設定する。\\ 例えば、"Sturges" だとスタージェスの方法、"Scott" だとスコットの方法になる。  |
| ''freq''                | ''FALSE'' なら縦軸を確率密度にする。                                                                                               |
| ''xlab''                | 横軸のラベルを設定する。                                                                                                           |
| ''ylab''                | 縦軸のラベルを設定する。                                                                                                           |
| ''title''               | タイトルを設定する。                                                                                                               |
| ''plot''                | ''FALSE'' ならプロットしない。                                                                                                     |

''hist()'' の戻り値には以下の要素が含まれる。

^  ラベル  ^  要素         ^
| ''breaks''      | 階級の境界値  |
| ''counts''      | 度数          |
| ''density''     | 確率密度（相対度数を階級幅で割ったもの）      |
| ''mids''        | 階級値        |



===== 最大値と最小値 =====

最大値
<code rsplus>
> max(data)
</code>

最小値
<code rsplus>
> min(data)
</code>

===== 平均 =====
__**平均**__は分布の重心を表す統計量である。

\[
\bar{X}=\frac{1}{n}\sum_{i=1}^{n}X_{i}
\]

関数 ''mean()'' で計算できる。

<code rsplus>
> mean(data)
</code>

合計の関数 ''sum()'' を使うと以下のようになる。

<code rsplus>
> sum(data) / length(data)
</code>


===== 中央値 =====
__**中央値**__は分布をちょうど二等分する点の値である。
$X_{i}$ が小さい順（もしくは大きい順）に並んでいるとすると、定義は以下のようになる。
\[
\mathrm{Me}=
\left\{
\begin{array}{ll}
X_{(n+1)/2} & （n が奇数）\\
(X_{n/2}+X_{n/2+1})/2 & （n が偶数）
\end{array}
\right.
\]

例えば、$n=5$ だと
\[
X_{1},X_{2},X_{3},X_{4},X_{5}
\]
の中で中央値は $\mathrm{Me}=X_{3}$ になる。

$n=4$ だと
\[
X_{1},X_{2},X_{3},X_{4}
\]
の中で中央値は $\mathrm{Me}=(X_{2}+X_{3})/2$ になる。



関数 ''median()'' で求められる。

<code rsplus>
> median(data)
</code>

===== 最頻値 =====
__**最頻値**__は級別データで度数が最大となるところの階級値である。
ヒストグラムの最も山が高いところでもある。

''hist()'' の戻り値 hg を使って以下のようにして求められる。
<code rsplus>
> hg$mids[which.max(hg$counts)]
</code>

関数 ''which.max()'' はデータの中で最大値の要素の番号を返す。
===== 小テスト =====

[[KMS>|Moodle Server（非公式）]]で第4回の小テストを受験しなさい。

===== 四分位数 =====
__**四分位数**__ は分布を四等分したときの三つの分割点の値である。
$X_{i}$ を小さい順に並べて四等分し、その三つの分割点を小さい順に
第一四分位数 $Q_{1}$、
第二四分位数 $Q_{2}$、
第三四分位数 $Q_{3}$ という。
$Q_{1}$ は 25%パーセンタイル、$Q_{2}$ は 50%パーセンタイル、$Q_{3}$ は 75%パーセンタイルにあたる((パーセンタイルとは、データの小さい方からデータ数を数えて全体の何パーセントの位置にあるかを示す値である。))
。

データを小さい順に並べたときの分布を四分位数を使って「箱」と「ひげ」の長さで表したグラフを__**箱ひげ図**__という。
データがどの付近に集中して分布しているのかを見やすくするためのものである。

四分位数の定義には複数あるが、ここでは二つの定義を挙げる。

----

==== テューキーの五数要約の定義 ====
John Wilder Tukey が提案した定義である。

  - データの中央値を $Q_{2}$ としてデータを半分に分ける。
    * $n$ が奇数なら $Q_{2}$ を<fc #ff0000>**両半分に含める**</fc>。
    * $n$ が偶数ならちょうど半分に分けられる。
  - 小さい方の半分の中央値を $Q_{1}$ とする。
  - 大きい方の半分の中央値を $Q_{3}$ とする。

これに最小値 $Q_{0}$ と最大値 $Q_{4}$ を加えて__**五数要約**__という。

五数要約は ''fivenum()'' で計算できる（もしくは ''summary()'' や ''quantile()''）。
<code rsplus>
> fivenum(data)
</code>

このテューキーの五数要約による箱ひげ図は ''boxplot()'' で簡単に描くことができる。

<code rsplus>
> boxplot(data)
</code>
デフォルトでは縦方向の箱ひげ図だが、引数に ''horizontal = TRUE'' を入れると横方向の箱ひげ図になる。

データのプロットと並べると以下のようになる。
データのプロットは見やすいように ''jitter()'' を使って y 方向にランダムに少しずらしてある。
<code rsplus>
> layout(1 : 2)
> plot(x = data, y = jitter(numeric(length(data))), ylim = c(-0.2, 0.2), ylab = "", yaxt = "n", xlab = "オゾン濃度[ppb]", main = "データ")
> boxplot(data, horizontal = TRUE, xlab = "オゾン濃度[ppb]", main = "テューキーの方法")
> layout(1)
</code>

{{ socialip_2021:boxplot-tukey.png?nolink |箱ひげ図（テューキーの方法）}}

IQR（Inter Quantile Range） は四分位範囲（箱の長さ）である。
この中にはデータの 50% が含まれている。
\[
\mathrm{IQR}=Q_{3}-Q_{1}
\]

ひげの下端は $Q_{1}$ から $1.5\times\mathrm{IQR}$ 以内の範囲にあるデータの最小値、
ひげの上端は $Q_{3}$ から $1.5\times\mathrm{IQR}$ 以内の範囲にあるデータの最大値である。
ひげの下端より小さいデータや上端より大きいデータは__**外れ値**__として○で描かれる。

外れ値がなければ、ひげの下端は $Q_{0}$、上端は $Q_{4}$ に等しい。

----
==== 中学校指導要領の定義 ====

  - データの中央値を $Q_{2}$ としてデータを半分に分ける。
    * $n$ が奇数なら $Q_{2}$ を<fc #ff0000>**両半分に含めない**</fc>。
    * $n$ が偶数ならちょうど半分に分けられる。
  - 小さい方の半分の中央値を $Q_{1}$ とする。
  - 大きい方の半分の中央値を $Q_{3}$ とする。


R ではこの方法で四分位数を求める関数は組み込まれていないが、以下のようにして計算できる。
<code rsplus>
> q <- quantile(data, type = ifelse(length(data) %% 2 == 0, 5, 6))
</code>

中学校指導要領で定義されている箱ひげ図は、
箱の上端を $Q_{3}$、中の横線を $Q_{2}$、下端を $Q_{1}$ とし、
ひげの上端を最大値、下端を最小値にする。
この定義だと外れ値はない。

データのプロットと箱ひげ図を並べると以下のようになる。
''boxplot()'' だとテューキーの方法になってしまうので、代わりに $Q_{1},Q_{2},Q_{3}$ などを指定して''bxp()'' で箱ひげ図を描いている。
<code rsplus>
> layout(1 : 2)
> plot(x = data, y = jitter(numeric(length(data))), ylim = c(-0.2, 0.2), ylab = "", yaxt = "n", xlab = "オゾン濃度[ppb]", main = "データ")
> bp.japan <- list(stats = matrix(q, 5, 1), n = length(data))
> bxp(bp.japan, horizontal = TRUE, xlab = "オゾン濃度[ppb]", main = "中学校指導要領の方法")
> layout(1)
</code>


{{ socialip_2021:boxplot-japan.png?nolink |箱ひげ図（中学校指導要領）}}

----
==== 例 ====

''boxplot()'' は複数の種類のデータをまとめて箱ひげ図として描画できる。

次のプログラム boxplot.R は第2回に e-Stat からダウンロードした山形県の市町村データ（ファイル名は yamagata.csv にしている）において、2015年度の第1次、第2次、第3次就業者比率をまとめて箱ひげ図で表示するプログラムである。

<file rsplus boxplot.R[enable_line_numbers="true"]>
library(readr)

# データの読み込み
df <- read_csv("yamagata.csv", locale = locale(encoding = "SHIFT-JIS"))

df$`/項目` <- NULL
df$調査年 <- as.numeric(sapply(df$調査年, gsub, pattern="年度", replacement = ""))
df <- df[order(df$調査年),]

# 2015年度の就業者比率
data <- df[df$調査年 == 2015,
           c("#F01201_第1次産業就業者比率【％】",
             "#F01202_第2次産業就業者比率【％】",
             "#F01203_第3次産業就業者比率【％】")]

# ラベルの書き換え
colnames(data) <- c("第1次産業", "第2次産業", "第3次産業")

# 箱ひげ図
boxplot(data, horizontal = TRUE, main = "山形県市町村の就業者比率", xlab = "就業者比率【％】")
</file>

===== 分散 =====
__**分散**__は分布が平均からどれくらい広がっているかを表す統計量である。

----
==== 分散 ====
\[
\sigma^{2}=\frac{1}{n}\sum_{i=1}^{n}(X_{i}-\bar{X})^{2}=s^{2}\left(1-\frac{1}{n}\right)
\]

R には分散を直接計算する関数は組み込まれていないので以下で計算する。
<code rsplus>
> var(data) * (1 - 1 / length(data)) 
</code>

----
==== 不偏分散 ====
標本から計算し、母分散の推定に使われる。

\[
s^{2}=\frac{1}{n-1}\sum_{i=1}^{n}(X_{i}-\bar{X})^{2}
\]

<code rsplus>
> var(data)
</code>
===== 標準偏差 =====
__**標準偏差**__は分散の平方根をとったものである。

----
==== 標準偏差 ====
\[
\sigma=\sqrt{\sigma^{2}}=s\sqrt{1-\frac{1}{n}}
\]

<code rsplus>
> sd(data) * sqrt(1 - 1 / length(data))
</code>
----
==== 不偏標準偏差 ====
\[
s=\sqrt{s^{2}}
\]

<code rsplus>
> sd(data)
</code>
もしくは
<code rsplus>
> sqrt(var(data))
</code>

===== 課題  =====
[[KMS>|Moodle Server（非公式）]] で第4回の課題を行いなさい。

締め切り：2021年7月15日（木）20時