このセクションで使う R パッケージ一覧

library(tidyverse)
library(patchwork)

1. 統計学の種類

統計学は記述統計学 (descriptive statistics) と推測統計学 (inferential statistics) の二つに分類できる

統計学の種類	特徴
記述統計学:	統計量を使ってデータ傾向や性質をつかむ
推測統計学:	母集団の母数（母平均や母分散）を検定し推定する

記述統計学

集めたデータの統計量（平均、分散、標準偏差など）を計算して、そのデータを我々が理解できるような形で記述し要約する
データの分布を明らかにすることで、データ傾向や性質をつかむための統計手法

推測統計学

入手可能な記述統計量を使って未観測の事象を予測、推定するスキルが「推測統計学」
母集団 (population) から標本 (sample) を無作為に抽出し、その標本によって得られた標本平均や標本不偏分散などの統計量を使って、母集団の母数（母平均や母分散）を検定し推定する
推測統計学の基本的な考えは、母集団からランダムに抽出した標本を増やし、無限に試行を繰り返せば、全体の一部である標本から、巨大で未知の母集団を推測できる、ということ
推測統計学では、分析対象が確率分布すると考える

「推測統計学」の肝は「統計的推測」

「統計的推測」とは、標本（サンプル）を使って母集団を推定すること
「統計的推測」は「統計的推定」と「仮説検定」から構成される

→　ここでは「統計的推定」と「仮説検定」を解説する

2. 母集団と標本

2.1 「標本数」と「標本サイズ」の違い

「標本サイズ」と「標本数」が混同されがちなので、注意が必要
標本サイズ (sample size)・・・・・・・n で表す抽出した観測数のこと
標本数 (the number of samples) ・・・ n 個の観測数で抽出した標本セット数

例：　　

矢内さんは人口328万人の郡山市に行って、東日本大震災に関して二つの標本をとった
標本 A では 1,000 人を対象に世論調査を実施し、標本 B では 2,000 人に対して世論調査を行った

→ 標本数は 2（標本 A と標本 B）であり、標本サイズはそれぞれ 2,000 と 1,000 である　　　　

2.2無作為抽出法 (`Random Sampling`)

(1) 単純ランダム・サンプリング法 (simple random sampling):

偏りのない標本を選ぶ代表的な方法
単純無作為抽出された標本は、母集団の偏りのない縮図と見なせる
→ その標本を使って、母集団を推定できる

(2) 多段抽出法:

何段階かに分けて標本の抽出を行う方法
→例: 日本の有権者から1000人を標本として抽出する
1 段階: 市町村を無作為に選ぶ
2 段階: 選ばれた市町村から無作為に標本を選ぶ
母集団分布が広範囲に分布している場合に多段抽出法を使う
→多段抽出法の一例・・・層化 2 段抽出
各市町村を大都市、中小都市、郡部と「層化」する
各層の大きさに応じたウェイ卜（例えば 3 : 2 : 1 ) をつけて市町村を無作為に選ぶ

【母平均 \((μ)\) と標本平均 \(\bar{x}\) の関係】

母集団に関する母数は未知のことが多いが、ここでは R を使って人工的に母集団を作成してみる
作成する母集団は、母平均 10、母標準偏差 5（母分散 = 25）
dnorm()関数を使うと、母集団を人工的に作成でき、その分布密度を表示できる

# 母平均 10、母標準偏差 5 と指定
# -10 から 30 の範囲を指定
curve(dnorm(x, 10, 5), from = -10, to = 30)

この母集団からランダムに標本を抽出する
抽出する標本サイズは 20 (N = 20)で、その標本に \((x1)\) と名前を付ける
ここで人工的に作成した母集団からランダムサンプリングする時は rnorm() 関数を使う

# 平均10、標準偏差 5 (sd = 5) の母集団から 20 個の数字を無作為抽出して x1 と名付け、データとヒストグラムを表示する

x1 <- rnorm(20, mean = 10, sd =5) 
x1

 [1]  2.866353 11.660236  9.974924  4.715395 15.918723 11.109943  6.555090
 [8]  8.345732 10.115209  7.905942 15.160718  6.664991 11.404893  7.057985
[15]  8.555314 17.871686  5.433826 11.915541  8.843969  9.404119

表示したい桁数は round() 関数を使って指定できる例えば、整数を表示したければ、小数点以下の表示桁数 (digits) を 0 と指定すればよい

round(x1, digits = 0)　# 小数点以下表示したい桁数を指定

 [1]  3 12 10  5 16 11  7  8 10  8 15  7 11  7  9 18  5 12  9  9

hist(x1)　　　　　　　# 取り出したサンプルをヒストグラムで表す

標本 \((x)\) の標本平均 \(\bar{x}\) は次の式で求めることができる: \[\bar{x} = \frac{\sum_{i=1}^n x_i}{n}\]

・R では mean( ) を使えば標本平均を求めることができる

mean(x1)              # 取り出したサンプル x1 の平均値

[1] 9.57403

2.3 「推定量」と「推定値」の違い

\[\bar{x} = \frac{\sum_{i=1}^n x_i}{n}\]

推定量 (estimator) :
母数の推定のために利用され、標本の関数として表現される・・・ここでは上記 \(\bar{x}\) の式
推定値 (estimate) :
実際に標本の数値を代入して計算した結果の数値・・・ここでは上記 \(\bar{x}\) の値
標本サイズ 20 の標本を無作為抽出して x2 と名付け、標本平均を求めてみる
母集団の平均は 10 (mean = 10)、標準偏差は 5 (sd = 5) と指定

x2 <- rnorm(20, mean = 10, sd =5) 
x2

 [1] 17.9553272  6.8431829  9.2014304  8.5704822  4.4044775 15.3251180
 [7] 10.2299722  7.0803491 18.7658989 17.8757231 12.2767327  8.6178387
[13]  8.0483961  0.7042173  2.1161221  7.3793859  5.2901019  5.7531278
[19]  2.5132830 13.6235880

round(x2, digits = 0)      # 小数点以下表示したい桁数を指定

 [1] 18  7  9  9  4 15 10  7 19 18 12  9  8  1  2  7  5  6  3 14

hist(x2)                 　# 抽出したサンプルをヒストグラムで表す

mean(x2)                  # 取り出したサンプル x2 の平均値

[1] 9.128738

R は母集団から無作為に数値を抽出しているので、サンプリングにおいて「偏り」はないが「誤差」はある
ここでは母平均を 10 に設定したにもかかわらず、一回目のサンプル x1 と二回目のサンプル x2 の標本平均がぴったり 10 にはならない
しかし、サンプリングを何度も繰り返すと（つまり x1, x2, ..... xn)、母平均より大きかったり小さかったりする標本平均が得られ、それらの値を平均すると母平均に一致する
→ 標本から得られた統計量を使って母集団のパラメタを統計的に推定できる根拠
標本平均以外にも min( ), median( ), max( )などのコードを使うと最小値、中央値、最大値などの記述統計が簡単に得られる

min(x2)　　# min = 最小値

[1] 0.7042173

median(x2) # median = 中央値

[1] 8.309439

max(x2)    # max = 最大値

[1] 18.7659

5数は次のコマンドで一度に表示できる

quantile(x2, c(0, .25, .5, .75, 1))

        0%        25%        50%        75%       100% 
 0.7042173  5.6373714  8.3094391 12.6134465 18.7658989

2.4 標本不偏分散と標本分散

例えば、学生 10 人が TOEFL (iBT) を受験したとする
10人のテストスコアに x という変数名を付ける

x <- c(22, 33, 44, 55, 66, 77, 88, 99, 100)
x

[1]  22  33  44  55  66  77  88  99 100

・標本を使った通常の分析では「標本不偏分散 \(U^2\)」を使う方が一般的
→ R で分散を計算する関数 var() では標本不偏分散 \(U^2\) をデフォルトに設定している

「標本不偏分散」と「標本分散」は使い分ける必要がある

・標本不偏分散 \(U^2\)・・・「母集団に興味がある場合」に使う散（＝標本標本不偏分散）
・標本分散 \(S^2\)・・・「抽出したサンプルにだけ興味がある場合」に使う

標本不偏分散 \(U^2\)を使う場合 (R ではデフォルト):

・母集団データの一部だけが手元にあり、入手したデータの背後にある母集団データの散らばりを知りたい場合

【標本不偏分散 \(U^2\) を求める式】:var() で計算可能

・仮に \(x\) が母集団から抽出された標本だとしよう
・母分散は未知なので、標本として得られた統計量を使って未知の母分散を推定する
・母分散を推定するための統計量が標本不偏分散 \(U^2\)
・標本不偏分散 \(U^2\)は次の式で求められる

\[U^2 = \frac{\sum_{i=1}^N (x_i - \bar{x})^2}{N-1}\]

・\(\bar{x}\): 標本平均

・Rで計算してみる

sum((x - mean(x))^2) / (length(x) - 1) # var(x) でも同じ結果を得る

[1] 808.6111

var() でも x の標本不偏分散 \((S^2)\) の値が出力される

# 標本不偏分散
var(x)

[1] 808.6111

標本分散 \(S^2\)を使う場合:

・母集団から抽出された標本(サンプル)が手元にある
・興味の対象が母集団ではなく、標本（サンプル）だけにあるとき
・そのサンプルがどの程度散らばっているかということだけを知りたい時

【標本分散 \(S^2\) を求める式】

\[S^2 = \frac{\sum_{i=1}^N (x_i - \bar{x})^2}{N}\]

・\(\bar{x}\): 標本平均

・R には標本分散 \(S^2\) を求める機能がない
→　\(S^2\) は次の式を作って計算しなければならない

# 標本分散
var(x) * (length(x) - 1) / length(x)

[1] 718.7654

標本分散 \(S^2\ \) (718.7654) より、標本不偏分散 (808.6111) \(S_x^2\) の方が常に大きい

・標本分散 \(S^2\) を求める式では分母が N なのに、標本不偏分散 \(U^2\) を求める式では分母が N - 1 だから
・詳細は「10.2 推測統計学（推定・検定）」を参照

分散の計算で２乗 (\(S^2\) や \(U^2\)) の和の平均をとる理由・１つの群における各データの数値の平均からの差（＝偏差）はぞれが平均値からどれだけ離れているかを表す指標
→　偏差はプラスとマイナス両方を取り得る
→　各データが平均値からどれだけ離れているかを表す指標として不適切
・偏差を２乗することにより、平均値からの距離の基準を絶対値に換算することで適切な指標に変換できる

3. 標本分布

3.1 統計量は分布する？

標本分布 (sampling distribution) :
ある標本サイズで標本を抽出したとき、そこで得られる統計量の分布
ランダムサンプリングにおいて、標本を選ぶ方法はたくさんある
例えば、10 人から 3 人を選ぶ方法 \((_{10}C_3)\) は 120 通りある
R では choose(　) コードを使って次のように計算できる

choose(10, 3)  # 10 人から 3 人を選ぶ組み合わせ

[1] 120

母集団が 1 億人の有権者から 10 人の標本を選ぶ方法を計算してみると、\((2.7)\) x \((10^{73})\) 通り以上あることがわかる

choose(100000000, 10)

[1] 2.755731e+73

特定の母集団について特定のサイズの標本を選ぶ方法は何通りもある
それぞれの標本について統計量の値を求めることは可能
しかし、統計量の値が全て一致するとは限らない
同じ大きさの標本を抽出する作業を繰り返せば、統計量は分布する
一定の標本サイズで全ての組み合わせで標本を抽出したとき、そこで得られる統計量の分布を標本分布と呼ぶ
標本分布・・・統計量と母数を確率でつなぐ統計的推定の肝

3.2 「母平均」と「標本平均」の関係

母集団を人工的に作る
3人の議員（Aさん、Bさん、Cさん）それぞれの当選回数

・当選回数 3, 6, 12回の母集団を作成し、母平均と母標準偏差を計算する

win <- c(3, 6, 12)  # 当選回数 3, 6, 12回の母集団を作成
mean(win)　　　　　 # 当選回数の母平均

[1] 7

分散を計算する

var1 <-  ((7-3)^2 +(7-6)^2 + (7-12)^2)/3
var1

[1] 14

標準偏差は分散の平方根なので、標準偏差は

sqrt(var1)

[1] 3.741657

平均当選回数は 7 回、当選回数の標準偏差は 3.7
従って、母平均 \((μ = 7)\)、母標準偏差 \((σ = 3.7)\)
この母集団からサイズ = 2 の標本を全て選ぶ
A, B, C の 3 人から 2 人が選ばれるのは次の 3 通り：「A と B」、「A と C」、「B と C」
例えば「選ばれる人物」が「A と B」の当選回数に関する標本平均と標本標準偏差を R で計算してみると、

ab <- c(3, 6)  # A と B が選ばれたサンプルを ab と名付ける
mean(ab)　　　 # サンプル ab の標本平均

[1] 4.5

分散を計算する

var2 <-  ((4.5-3)^2 + (4.5-6)^2)/2
var2

[1] 2.25

標準偏差は分散の平方根なので、標準偏差は

sqrt(var2)

[1] 1.5

上の表の「標本平均」を見ると、3 つの標本の標本平均 (4.5, 7.5, 9) が全て異なっている
しかも、事前に設定した母集団の母平均 \((μ = 7)\)とも異なる
標本から母数について推定できること・・・標本が母集団の偏りのない縮図であること
「偏りのない」＝ひとつひとつの統計量の値が母数と異なっても、得られた値の平均値が母数に一致する
3 つの標本平均の平均値 \((μ)\) を求めてみると

\[{μ} = \frac{4.5+7.5+9}{3} = 7\]

抽出可能なすべての組み合わせで標本を抽出し、各標本から得られた平均値の平均値を求めると、母平均に一致する
→ 不編性 (unbiasedness)　　
→ 不編推定量 (unbiased estimator): 不偏性のある推定量のこと
→ 「標本平均 \(\bar{x}\) は母平均 \((μ)\) の不偏推定量であるといえる」
→ 母平均 \((μ)\) は標本平均 \(\bar{x}\) を使って推定する

3.3 「母標準偏差」と「標本標準偏差」の関係

しかし、3 つの標本標準偏差の平均値を求めてみると、母標準偏差 (3.7) に一致しない

\[\frac{1.5+4.5+3}{3} = 3 < {3.7}.\]

標本標準偏差の平均は、常に母標準偏差より小さい
標本標準偏差（統計量）は母標準偏差（母数）を過小評価するという体系的な bias（偏り）をもつ
→ 「標本標準偏差 \((s)\)は母標準偏差 \((σ)\) の不偏推定量とはいえない」
→ 標本不偏標準偏差 \((u_x)\)を使って母標準偏差 \((σ)\) を推定する
→ 母標準偏差 \((σ)\) は標本不偏標準偏差 \((u_x)\) を計算できる
sd( ) 関数を使って推定する

sd(win) # win の不偏標準偏差を求める

[1] 4.582576

母分散 \((\sigma^2)\) を求める式・・・（母標準偏差 \((σ)\) は \((σ^2)\) の平方根）

\[\sigma^2 = \frac{\sum_{i=1}^N (x_i - \mu)^2}{N}\]

標本不偏分散 \((u_x^2)\) を求める式・・・（不偏標準偏差 \((u_x)\) は \((u_x^2)\) の平方根）

\[u_x^2 = \frac{\sum_{i=1}^n (x_i - \bar{x})^2}{N-1}\]

R を使って計算する時に使うコード:

母分散 \((\sigma^2)\) ・・・var(x) - (length(x) - 1) / length(x)

母標準偏差 \((\sigma)\) ・・・sqrt(var(x) - (length(x) - 1) / length(x))

標本不偏分散 \((u_x^2)\) ・・・var(x)

不偏標準偏差 \((u_x)\) ・・・sd(x)

3.4 大きな母集団を想定したシミュレーション

人工的に日本人成人男性の身長に関する正規分布する母集団（n = 10万人、母平均 = 170cm、母標準偏差 = 6）を作り、height と名付ける

height <- rnorm(100000, 170, 6)
hist(height)

10 万人から 10 人選ぶ方法は、\((2.7)\) x \((10^{43})\) 通り以上あることがわかる

choose(100000, 10)

[1] 2.754492e+43

\((2.7)\) x \((10^{43})\) 通りの組み合わせで標本平均を求めるのは不可能
代替策 → 10 万人の母集団から無作為に 10 人を選び sample1 と名付け、標本平均を求め、ヒストグラムで示す

sample1 <- sample(height, 10, replace = FALSE) # 標本サイズ 10 で無作為抽出する
hist(sample1)　　　　　　　　　　　　　　　　　# ヒストグラムで表す

mean(sample1)　　　　　　　　　　　　　　　　  # sample1 の平均値を表示する

[1] 167.6417

母集団から 10 人を無作為に選ぶという作業を 500 回行い、その分布を調べる
標本サイズを 10 と 90 の二種類 (n = 10, 90)、標本数は 2 つ作成し、アニメーションで表示する

アニメーションを使ったシミュレーション (n = 10, 90)

人工的に日本人成人男性の身長に関する母集団（n = 10万人、母平均 = 170cm、母標準偏差 = 6）を作り、height と名付ける

height <- rnorm(100000, 170, 6)  
## Parameters: n = 標本サイズ
##              T = 標本数
height_experiment <- function(n = 10, T = 500, seed = 2016-01-24) {
    means <- rep(NA, T)
    s_mat <- matrix(NA, ncol = n, nrow = T)
    set.seed(seed)
    height <- rnorm(100000, 170, 6)
    for (i in 1:T) {
        s <- sample(height, n, replace = TRUE)
        s_mat[i,] <- s
        means[i] <- mean(s)
        par(family = "HiraKakuPro-W3", cex.lab = 1.4, cex.axis = 1.2, cex.main = 1.6, mex = 1.4)
        hist(means[1:i], axes = FALSE, freq = FALSE, col = "gray",
             xlim = c(min(height), max(height)),
             xlab = "sample mean", ylab = "probability density",
             main = paste0("sample size = ", n, ", sample id = ", i))
        axis(1, min(height):max(height))
        abline(v = 170, lwd = 2, col = "red")
    }
    axis(2)
    return(s_mat)
}

シミュレーションの結果を表示する際には、次のコマンドを入力する

source("<path-of-the-directory>/height_sim.R")
sim_1 <- height_experiment(n = 10) # サイズ 10 の標本を 500 個抽出
sim_2 <- height_experiment(n = 90)　# サイズ 90 の標本を 500 個抽出

母集団からサイズ 10 の標本を 500 個抽出した時の標本平均 \(\bar{x}\) の分布

母集団からサイズ 90 の標本を 500 個抽出した時の標本平均 \(\bar{x}\) の分布

どちらの標本においても、ヒストグラムの中心が母平均 \(({μ} = 170)\) に一致 → どんなに母集団が大きくても、標本平均は母平均の普遍推定量である
母集団のばらつきより、標本平均のばらつきが小さいその理由 → 標本平均では極端な値を取りにくいから

標本平均は母集団の分布の両端付近の値をとりにくく、母集団よりも狭い範囲に分布する
分布のばらつき・・・標準偏差で測る
標本平均のばらつき・・・次の式で表す

\[SD (\bar{x})= \frac{\sigma}{\sqrt{n}}\]

式の分母に \(n\) が含まれているので、標本サイズ \(n\) が大きくなるほど、標本平均はますます極端な値を取りにくくなる
母集団の身長の母標準偏差 \(({\sigma = 6})\)、抽出した 500 個の標本サイズは 10 だから、

\[SD (\bar{x})= \frac{\sigma}{\sqrt{n}}\]

\[= \frac{6}{\sqrt{10}}\] \[= 1.9\]

母集団の身長の母標準偏差 \(({\sigma = 6})\)、抽出した 500 個の標本サイズを 90 に増やすと、

\[SD (\bar{x})= \frac{\sigma}{\sqrt{n}}\]

\[= \frac{6}{\sqrt{90}}\] \[= 0.63\]

標本サイズは 10 のときの標準偏差 1.90 のおよそ3 分の 1
一致性 (consistency): 標本サイズを大きくするほど、ひとつの標本から計算される標本平均が母平均近くの値をとる確率が大きくなる
→ 一致推定量 (consistent estimator)

4. 母平均の推定と信頼区間

4.1 \(z\) 値と \(t\) 値

点推定 (point estimation): 例）「母平均は 45 点」
区間推定 (interval estimation): 例）「母平均は 43 点以上、47 点以下」
信頼区間 (confidence interval: CI)・・・推定に伴う不確実性を表すために利用される
標本平均の標準偏差 \(({SD (\bar{x}})\)）がわかる
→ 標本平均の不確実性を捉えることができる
→ 推定の不確実性を明らかにできる
標本平均の標準偏差 \(({SD (\bar{x}})\)）を求める式は、

\[SD (\bar{x})= \frac{\sigma}{\sqrt{n}}\]

母標準偏差 \({\sigma}\) と標本サイズ \(n\) の値がわかれば、標本平均の標準偏差 \(({SD (\bar{x}})\)）がわかる
標本サイズ \(n\) はわかるが、母標準偏差 \({\sigma}\) は通常わからない
→ 母標準偏差 \({\sigma}\) の代わりに不偏標準偏差 \(u_x\) を使って、標本平均の標準偏差 \(({SD (\bar{x}})\)）を推定
→ このページの「標本分布」の項目を参照

\[SE = \frac{u_x}{\sqrt{n}}\]

この \(SE\) が標準誤差 (standard Error: Std. Errとも書く)
母標準偏差 \({\sigma}\) の値がわからない時には、標本平均の標準偏差 \({SD (\bar{x}})\)の推定値として標準誤差 \(SE\) を使う
\(t\) 値は、標本平均のバラツキを表す標準誤差 \(SE\) を使って、次の式で求めることができる

\[t = \frac{\bar{x} - μ}{SE} = \frac{\bar{x} - μ}{u_x / \sqrt{n}} \]
・ \(\bar{x}\) : 標本平均
・ \(μ\) : 母平均
・ \(n\) : 標本サイズ
・ \(u_x\) : 不偏標準偏差（= 標本標準偏差）
・ \(SE\) : 標準誤差 (standard Error)

参考

母集団の標準偏差 \({\sigma}\) がわかっている場合
\(t\) 分布ではなく標準正規分布 \(z\) 分布)を利用できる
\(z\) 値は、次の式で求めることができる

\[Z = \frac{\bar{x} - μ_0}{σ / \sqrt{n}} \]

\(z\) 値は、\({\sigma}\) を使って\(\bar{x}\) を変形（＝標準化）した値
しかし現実のデータ分析で母集団の標準偏差 \({\sigma}\) が既知であることはほとんどない 　　
実際には \(t\) 値を使って分析する
自由度 (df: degree of freedom) の大きさごとに 4 種類の \(t\) 分布図を描いてみる

x <- seq(-4, 4, length=100)
hx <- dnorm(x)

df <- c(1, 3, 8, 30)
colors <- c("blue", "red", "green", "gold", "black")
labels <- c("df=1", "df=3", "df=8", "df=30", "N.Distribution")

plot(x, hx, type="l", lty=2, xlab="t value",
  ylab="Density", main="t Distribution and degree of freedom")

for (i in 1:4){
  lines(x, dt(x,df[i]), lwd=2, col=colors[i])
}

legend("topright", inset=.05, title="Degree of Freedom",
  labels, lwd=2, lty=c(1, 1, 1, 1, 2), col=colors)

\(t\) 分布の特徴は次のとおり

0 を中心とした左右対称な分布
自由度 (degree of freedom: df) の数値によって分布の形が変わる
自由度が小さい（サンプルサイズが小さいと）\(t\) 値が大きく変動する
→ \(t\) 分布表に自由度が小さい時の\(t\) 値を詳しく載せている理由
自由度が大きくなるにつれて（サンプルサイズが大きくなるにつれて）分布が中央に集まる
\(t\) 値を求める式の分子は \(({\bar{x} - μ})\) なので、\(({\bar{x} = μ})\) の時（つまり、標本平均が母平均と同じ時） \(t\) 値は 0 になる
標本平均の推定では、自由度が \(({n - 1})\) の \(t\) 分布表を使う

例えば、標本サイズ n = 100 なら、自由度 99 の \(t\) 分布を使う
「標本平均 \(\bar{x}\) を変形した \(t\) が自由度 99 の \(t\) 分布に従う」ことを利用して統計的推定を行う
\(t\) 分布図の黒い実線（自由度 = 99）が示していること
→ 標本を何度も取り直した時、自由度 99 の分布に従う \(t\) の分布の状態
→ 確率密度曲線 (probability density curve)

自由度99の \(t\) 分布の確率密度曲線

確率密度曲線と横軸（確率密度 = 0）に囲まれる部分の「面積」が「確率」を表す

自由度 \(({n - 1 = 99})\) のとき、\(t\) の95％は \(({-1.98})\) 以上、\(({1.98})\) 以下の範囲に収まる（この範囲を\(({[-1.98, 1.98]})\) と表す）
\(t\) 分布の形は自由度によって変わる
→ 自由度が変われば、\(({-1.98})\) や \(({1.98})\) 以外の値を使う
一般的に、\(t\) が自由度 \(({n - 1})\) の自由度 \(t\) 分布に従うとき
→ \(t\) の 95% は区間を\(({[-t_{n-1,0.025}, t_{n-1,0.975}]})\) に収まる
→ 曲線左右両端の空白部分の面積はそれぞれ \(({0.025 (0.025*2=0.05: 有意水準（α）= 0.05) })\)
→ \(t\) の 90% は区間を\(({[-t_{n-1,0.05}, t_{n-1,0.975}]})\) に収まる
→ 曲線左右両端の空白部分の面積はそれぞれ \(({0.05 (0.05*2=0.10: 有意水準（α）= 0.10) })\)
有意水準 (\(α\): level of significance)
どれくらい低い確率であれば「現実には起こりえない」と言えるのかという基準
\(α = 0.05 (5 パーセント)\)を使うのが一般的
qt( ) 関数を使うと、有意水準を指定して棄却域を求めることができる
自由度が 99 で、有意水準 5%、両側検定の時の棄却域は、

qt(0.025, 99)  # t分布で下側確率が0.025、自由度99のt値

[1] -1.984217

qt(0.975, 99)  # t分布で上側確率が0.975、自由度99のt値

[1] 1.984217

4.2 Excercise

Q1: 自由度が 99 で、有意水準 10%、両側検定の時の棄却域を R を使って求めなさい

Q2: 自由度が 99 で、有意水準 1%、両側検定の時の棄却域を R を使って求めなさい

5. 信頼区間

5.1 信頼区間の計算方法

日本の有権者の民主党に対する感情温度を調査
単純無作為抽出によって n = 100人の標本を入手
不偏標準偏差 \(u_x\) = 21.36
民主党に対する感情温度: 平均感情温度 \(\bar{x}\) = 55.13 度
有権者全体における母平均感情温度は 50 度だといえるのか？

標本平均 \(\bar{x}\)、標準誤差 \((SE = \frac{u_x}{\sqrt{n}})\) の時、母平均 \((μ)\) の 95％信頼区間は次の式で求められる

\(({[\bar{x}-t_{n-1,0.025}・SE, \bar{x} + t_{n-1,0.975} ・SE]})\)

平均感情温度 \(\bar{x}\) = 55.13度
\(({t_{99,0.025}})\) = 1.98

\[SE = \frac{u_x}{\sqrt{n}}= \frac{21.36}{\sqrt{100}} = 2.136\]

以上を代入すると、母平均 \((μ)\) の 95％信頼区間は

\(({[\bar{x}-t_{n-1,0.025}・SE < μ < \bar{x} + t_{n-1,0.975}・SE]})\)

55.13 - 1.98・2.136 ≦ \((μ)\) ≦ 55.13 + 1.98 ・136

50.9 ≦ \((μ)\) ≦ 59.4

結論・有権者全体の民主党に対する感情温度の 95% 信頼区間は 50.9度から 59.4度
→　「50度」は感情温度の 95% 信頼区間に含まれない
→　母集団での感情温度の平均は 50度だとはいえない

信頼区間のまとめ

複数の標本を抽出すると、標本平均 \(\bar{x}\) を変形した \(t\) が自由度 \(({n - 1})\) の \(t\) 分布に従う

\[t = \frac{\bar{x} - μ}{SE} = \frac{\bar{x} - μ}{u_x / \sqrt{n}} \]

\(t\) 分布の特性により、すべての標本の 95%について次の式が成り立つ

\[-t_{n-1,0.025} ≦ \frac{\bar{x} - μ}{SE} ≦ t_{n-1,0.975}\]

\(μ\) を中心に変形すると、

\(({[\bar{x}-t_{n-1,0.025}・SE < μ < \bar{x} + t_{n-1,0.975}・SE]})\)

得られた標本数全体の 95% については \(({[\bar{x}-t_{n-1,0.025}・SE]})\) と \(({[\bar{x}+t_{n-1,0.025}・SE]})\) の区間に、真の母平均 \((μ)\) の値が含まれるということ
こうして求められた区間を「95％信頼区間」と呼ぶ

5.2 信頼区間の解釈

「95％信頼区間」の意味：
「母数の値がこの区間に含まれている確率が95%」は間違い
それぞれの標本から得られた95%信頼区間に母数が含まれている確率は 0% か 100%のどちらか

→ 上から1 番目のサンプルの95%信頼区間に母数\((μ = 45)\) が含まれる確率・・・100%

→ 上から11番目のサンプルの95%信頼区間に母数\((μ = 45)\)が含まれる確率・・・ 0%
標本サイズ100で、標本数20個をとり、20個の「95％信頼区間」が得られたとする
そのうちの95%（つまり19個）は母数をとらえるが、5%（つまり1個）は母数を捉え損ねる

R を使って感情温度の問題を確かめてみる
日本の有権者の民主党に対する感情温度を調査
単純無作為抽出によって n = 100人の標本を入手
不偏標準偏差 \(u_x\) = 21.36
民主党に対する感情温度: 平均感情温度 \(\bar{x}\) = 55.13 度
有権者全体における母平均感情温度は 50 度だといえるのか？
・母集団における平均感情温度を 55 度、標準偏差を 20 と想定し、100人の標本をランダムに抽出し、記述統計を表示

set.seed(2017)
emotion <- rnorm(100, mean=55, sd=20)
summary(emotion)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   3.14   39.56   55.20   55.13   69.00  117.75

hist(emotion)

sd(emotion)

[1] 21.36289

R を使って母平均 = 50 という条件で t 検定する

t.test(emotion, mu=50)


    One Sample t-test

data:  emotion
t = 2.4016, df = 99, p-value = 0.01819
alternative hypothesis: true mean is not equal to 50
95 percent confidence interval:
 50.89171 59.36943
sample estimates:
mean of x 
 55.13057

下から4 〜 5 行目の結果「95 percent confidence interval: 50.89171 − 59.36943」に注目
→これが「感情温度の 95% 信頼区間」

結論・「有権者全体の民主党に対する感情温度の 95% 信頼区間は 50.9度から 59.4度」
→「50度」は感情温度の 95% 信頼区間に含まれない
→ 母集団での感情温度の平均は 50度だとはいえない

p-値を使った検定の解説は「11. t 検定（1グループ）」を参照

6. 確率変数・確率分布・確率密度関数

6.1 基本事項の確認

コイン投げ試行のように、ある確率である事象が起こるか起こらないかが定まっている試行をベルヌーイ試行という　　
コインの表がでる確率を小文字の \((p)\) で表すと、ベルヌーイ試行は次の式で表せる

\[P(X = Head) = p\] \[P(X = Tail) = 1 - p\]

\((P (X))\) は確率を表し 0 から 1 の範囲の実数をとる
コインの表裏を表す変数\((X)\) を確率変数と呼ぶ
コインの表が出るか裏が出るか、コインを投げ終わらないと確定しない → このような変数を「確率変数」という
確率変数のとる各値（この場合は “Head”” と “Tail”）に対応して、それらが起こる確率が与えられる時、その対応を確率分布と呼ぶ
確率変数 \((X)\) とは「確率分布に従って特定の値になる変数」
確率変数 \((X)\) は、確率 \((p)\) で表（“Head”）になり、確率 \((1 - p)\) で裏（“Tail”）になる
ベルヌーイ試行を繰り返す
コインの表(“Head”)が出る回数の分布を二項分布という
R は二項分布の確率密度関数 dbinom( ) を備えている
10回コインを投げるとき、表が 5 回でる二項分布の確率密度（＝確率）は、

dbinom(5, 10, 0.5)  # dbinom(x, size, probability)

[1] 0.2460938

10回コインを投げたとき、表が出る回数が 0 から 10 それぞれの確率を求めると、

dbinom(0:10, 10, 0.5)

 [1] 0.0009765625 0.0097656250 0.0439453125 0.1171875000 0.2050781250
 [6] 0.2460937500 0.2050781250 0.1171875000 0.0439453125 0.0097656250
[11] 0.0009765625

[1] のすぐ右隣に示されている数字がコインが 0 回出る確率、その右隣の数字が 1 回出る確率と続き、最後の[11] のすぐ右隣に示されている数字が 10 回出る確率
予想通り、5 が出る確率が最も大きく (0.246)、4 と 6 がほぼ同じ程度 (0.205)、3 と 7 がほぼ同じ程度 (0.117)
5 から左右に離れるにつれて確率が小さくなっていることがわかる
コインの表が出る確率 0.5 でベルヌーイ試行を 10 回繰り返すと、表が 5 回出る確率が最も高いことがわかる
10回コインを投げたとき、表が出る「回数」ではなく、表が出る回数の「確率」をプロットすると、

# コインが出る回数を 0 から 10 まで指定、試行は 10 回、表が出る確率は 0.5と指定
coin10 <- dbinom(0:10, 10, 0.5) 
plot(0:10, coin10, type = "h", lwd = 5) # "h"は図に垂直な線を指定、lwd では線の太さを指定

縦棒の高さが確率密度（＝確率）
（連続量の場合には確率密度＝確率ではなく、区間の面積＝確率）
ここに示されているデータの出現確率を確率分布という
その確率分布を計算する式を確率密度関数 (probability density function) といい、次の式で表す

\[f(X) = _NC_XP^X{(1 - P)^{N-X}}\]

この式に試行回数 \((N = 10)\)、成功回数 \((X = 5)\) を代入してみる
10回サイコロを投げて 5 回表が出る組み合わせ \((_NC_X)\) は choose(　) コードを使って計算できる

choose(10, 5)  # 10回中 5 回表が出る組み合わせ

[1] 252

\((P^X{(1 - P)^{N-X}})\)の値は、

(0.5)^5*(1-0.5)^5

[1] 0.0009765625

従って、\((f(X))\)の値は、

(252)*(0.0009765625)

[1] 0.2460938

この値は上記ヒストグラムの縦棒の高さ（確率密度＝確率）
この確率密度（＝確率）は R では次のように求めることができる

dbinom(5, 10, 0.5)  # dbinom(x, size, probability)

[1] 0.2460938

6.2 コイントスのシミュレーション

前節のシミュレーションでは、コインを 10 回投げて表が出た回数（0 から 10 までの整数）を数えた（＝離散値）
ここでは連続量を使ったシミュレーションを行う　　
連続量のデータの確率は、ある特定の値ではなく、その値を含むある「区間」を取り、その面積で表す
一例として、最小値が 140、最大値が 185 で一様分布している連続変数（早稲田大学の学生の身長）から 10 人の身長を無作為抽出してみる
ここで使うコードはrunif(n, a, b)で \((a)\) が最小値、\((b)\) が最大値、そして \((n)\) が抽出数である

runif(10, 140, 185)

 [1] 162.8133 170.4660 181.2077 177.7272 161.2677 171.0087 176.0036 156.1114
 [9] 147.5133 150.7378

runifの r は random（無作為）、unif は uniform（一様分布）を表している
これは無作為抽出なので、実行する度に、異なる数値が抽出される

runif(10, 140, 185)

 [1] 141.0996 142.8775 171.5061 140.8381 167.1321 170.8180 177.8217 150.8456
 [9] 151.2468 175.7331

学術論文を仕上げる際、無作為抽出の結果として同一の結果を得たいのために、無作為抽出する際に set.seed( ) を使って seed を設定できる

set.seed(2016-01-09)
runif(10, 140, 185)

 [1] 179.9503 179.6794 149.5198 162.7430 148.9673 175.3110 171.7291 168.4825
 [9] 154.0057 140.8939

同じ値を seed 設定すると、実行する度に、同一の数値が抽出される

set.seed(2016-01-09)
runif(10, 140, 185)

 [1] 179.9503 179.6794 149.5198 162.7430 148.9673 175.3110 171.7291 168.4825
 [9] 154.0057 140.8939

無作為抽出で得られた10個の数値をヒストグラムにすると

set.seed(2016-01-09)
hist(runif(10, 140, 185))

無作為抽出する数値の数を100個に増やして、それらの数値をヒストグラムにすると

set.seed(2016-01-09)
hist(runif(100, 140, 185))

無作為抽出する数値の数を1000個に増やして、それらの数値をヒストグラムにすると

set.seed(2016-01-09)
hist(runif(1000, 140, 185))

母集団は最小値が 140、最大値が 185 で一様分布している連続変数だから、抽出する標本サイズを増やす程、母集団の形と似てくることがわかる
一様分布だけではなく、R では次のコマンドを使うことで様々な分布から無作為抽出することが可能
よく使われる分布は次のとおりである

`rnorm()`	: 正規分布
`rt()`	: \(t\) 分布
`rchisq()`	: \((\chi^2)\)（カイ二乗）分布

それぞれの分布から無作為抽出する際、それぞれ特有の引数 (arguments) があるため注意が必要
例えば、正規分布では rnorm() のカッコの中に平均値 (mean) と標準偏差 (sd) を、\(t\) 分布と \((\chi^2)\)分布ではカッコの中に自由度 (df) を指定しなければならない
rnorm() を使って無作為抽出してみる
正規分布する母集団 N \((\mu)\), \((\sigma^2)\) から n 個の数値を無作為抽出する時：
→ rnorm(n, mean = mu, sd = sigma) に必要な情報を入力
母平均 (mu) 160、母分散 16 (sigma = 4) で正規分布する母集団から n = 100 個の数値を無作為抽出する場合、次のように入力すればよい

rnorm(100, 160, 4)

  [1] 159.3848 159.4990 164.7559 155.6251 153.1917 158.3015 163.9680 153.7831
  [9] 158.1598 160.8625 161.2840 161.1273 164.9580 164.3071 157.1850 152.9570
 [17] 157.9721 165.1655 161.9859 157.9363 159.3190 153.6560 160.8890 160.0436
 [25] 152.7490 156.1923 166.4202 159.1252 161.6708 158.8871 158.8819 163.8937
 [33] 154.8815 163.1115 159.2016 152.5330 157.0690 152.9178 158.2415 162.3454
 [41] 160.8212 160.1693 162.7420 157.0018 162.6093 156.4365 163.4589 157.0134
 [49] 161.0999 163.0928 154.6236 167.2826 165.1174 152.7759 156.5711 156.8450
 [57] 160.3522 162.3860 159.9502 162.2307 159.0250 163.0924 164.7416 161.5864
 [65] 157.7696 167.6563 159.5414 163.9407 164.6144 156.0849 156.9475 164.7797
 [73] 156.1790 157.2057 157.2734 159.0295 157.1015 157.6920 159.8148 160.7410
 [81] 164.8714 157.5850 160.9043 163.3799 157.6321 156.9417 154.8651 154.0923
 [89] 161.5911 163.2804 165.0042 159.4423 163.6247 161.3810 159.2093 160.0146
 [97] 153.1990 164.3284 150.2262 158.1045

# ヒストグラムを描くと
hist(rnorm(100, 160, 4))

母集団が正規分布しているため、標本サイズが大きくなるほど、無作為抽出される標本も正規分布の形に似てくる

# 5数は
quantile(rnorm(100, 160, 4), c(0, .25, .5, .75, 1))

      0%      25%      50%      75%     100% 
150.0572 157.2044 159.8816 163.0106 168.9840

# boxplotを描くと
boxplot(rnorm(100, 160, 4))

7. 期待値

7.1 期待値 = 平均値

期待値 \((E(X_i))\) とは「ある事柄が起こる確率 \((P(X))\)」と「その結果の数値 \((X_i)\)」を全て掛け合わせ、合計した値のこと　　
期待値は「平均値」とも呼ばれる
数式で表すと、

\[E(X) = \frac{\sum_{i=1}^N p(X_i)X_i}{N}\]

事例）
サイコロを振ったときの期待値は「それぞれの目がでる確率 (1/6)」と「その結果の数値 (1, 2, 3, 4, 5, 6)」を掛け合わせ、合計する

\[E(X) = \frac{\sum_{i=1}^N p(X_i)X_i}{N}\]

\[= \frac{(1/6)*1 + (1/6)*2 + (1/6)*3 + (1/6)*4 + (1/6)*5 + (1/6)*6}{N}\]

\[= 3.5\]

R を使って、サイコロを 10 回振り、その期待値を求めてみると、

x <- 1:6
y <- sample(x, 10, replace = TRUE)
mean(y)

[1] 2.7

理論的にはサイコロを振る際の期待値は 3.5 だが、実際にはぴったり 3.5 が得られるわけではない
もう 10 回振ってみる

y <- sample(x, 10, replace = TRUE)
mean(y)

[1] 4.3

さらにサイコロを振る回数を 10000 回に増やして、期待値を求めてみると

y <- sample(x, 10000, replace = TRUE)
mean(y)

[1] 3.5204

サイコロを振る回数を限りなく増やすと、期待値が 3.5 に限りなく近づく（＝大数の法則）

7.2 期待値のシミュレーション

確率分布・・・データの出現状況を確率に置き換えること
ex) ①サイコロを振る、②コインを投げる、etc
ある確率である事象が起こる時、その事柄を「事象」(event)という
R には一定の確率で分布する母集団から無作為に標本を抽出する機能が備えられている 【二項分布のシミュレーション】
R を使ってコインを 10 回投げるシミュレーションをしてみる
まず、コインの表を “H”、裏を “T” と定義して架空のコイン (coin) を作る
コインを 10 回投げて、表と裏が出た回数（整数）をそれぞれ記録する
このような整数を離散値と呼ぶ

coin <- c("H", "T") 
flip <- sample(coin, 10, replace = TRUE) # 10回コインを投げる
table(flip)　# 結果をテーブルに表す

flip
H T 
2 8

コインの表 (“H”) がでる確率は理論的には 50 % のはずだが、実際に10 回コインを投げると、ぴったり 5 回表が出るとは限らない
コインを 10 回投げる試みを 20 回やってみて、その結果をヒストグラムで表してみる

x <- numeric(20) # ベクトル x を作り、そこに 20 個の 0 を入れる
for(i in 1:20) { # コインを 10 回投げる度、変数iを 1 から 20 まで増やす
  flip <- sample(coin, 10, rep = TRUE) # コインを 10 回投げるコマンド
  x[i] <- sum(flip == "H") # 右側の結果"H"を左の x に代入するコマンド
}
table(x) # 結果 x を表にするコマンド

x
4 5 6 7 8 
2 5 9 3 1

hist(x, breaks = 0:10)　# ヒストグラムの区間を 0 から 10 に限定

シミュレーションをやる度に得られる結果は異なるが、10 回のコイン投げを 20 回試みると、必ずしも 5 回が最頻値であるとは限らず、いびつな分布になることが多い
10 回のコイン投げを 10000 回やってみる

x <- numeric(10000) # ベクトル x を作り、そこに 10000 個の 0 を入れる
for(i in 1:10000) { # コインを 10 回投げる度、変数iを 1 から 10000 まで増やす
  flip <- sample(coin, 10, rep = TRUE) # コインを 10 回投げるコマンド
  x[i] <- sum(flip == "H") # 右側の結果"H"を左の x に代入するコマンド
}
table(x) # 結果 x を表にするコマンド

x
   0    1    2    3    4    5    6    7    8    9   10 
  13  110  427 1213 2021 2429 2072 1139  454  114    8

hist(x, breaks = 0:10)　# ヒストグラムの区間を 0 から 10 に限定 
abline(v = mean(x), col = "red") # 平均値に赤い縦線を引く

mean(x)

[1] 4.9975

結果をヒストグラムで表してみると、平均値の 5 に確実に近づいていることが明らか

8. Exercise

Q1: R を使って、母平均 50、母標準偏差 10 の母集団から sample size 5 の sample を 3 つ (number of samples = 3) 無作為抽出し、3つの sample それぞれのヒストグラム、平均、標準偏差を示しなさい

Q2: Rを使って、母平均 50、母標準偏差 10 の母集団からsample size500 の sample を 3 つ (number of samples = 3`) 無作為抽出し、3つの sample それぞれのヒストグラム、平均、標準偏差を示しなさい

Q3: 上の二種類のシミュレーションからどのようなことが言えるか

参考文献

宋財泫 (Jaehyun Song)・矢内勇生 (Yuki Yanai)「私たちのR: ベストプラクティスの探究」

土井翔平（北海道大学公共政策大学院）「Rで計量政治学入門」

矢内勇生（高知工科大学）統計学2

浅野正彦, 矢内勇生.『Rによる計量政治学』オーム社、2018年

浅野正彦, 中村公亮.『初めてのRStudio』オーム社、2018年

Winston Chang, R Graphics Cookbook, O’Reilly Media, 2012.

Kieran Healy, DATA VISUALIZATION, Princeton, 2019

Kosuke Imai, Quantitative Social Science: An Introduction, Princeton University Press, 2017

10. 推測統計学（基礎）

Masahiko Asano

2024-04-14

1. 統計学の種類

記述統計学

推測統計学

「推測統計学」の肝は「統計的推測」

2. 母集団と標本

2.1 「標本数」と「標本サイズ」の違い

2.2無作為抽出法 (`Random Sampling`)

2.3 「推定量」と「推定値」の違い

2.4 標本不偏分散と標本分散

「標本不偏分散」と「標本分散」は使い分ける必要がある

標本不偏分散 \(U^2\)を使う場合 (R ではデフォルト):

標本分散 \(S^2\)を使う場合:

3. 標本分布

3.1 統計量は分布する？

3.2 「母平均」と「標本平均」の関係

3.3 「母標準偏差」と「標本標準偏差」の関係

3.4 大きな母集団を想定したシミュレーション

4. 母平均の推定と信頼区間

4.1 \(z\) 値と \(t\) 値

4.2 Excercise

5. 信頼区間

5.1 信頼区間の計算方法

5.2 信頼区間の解釈

6. 確率変数・確率分布・確率密度関数

6.1 基本事項の確認

6.2 コイントスのシミュレーション

7. 期待値

7.1 期待値 = 平均値

7.2 期待値のシミュレーション

8. Exercise

10. 推測統計学（基礎）

Masahiko Asano

2024-04-14

1. 統計学の種類

記述統計学

推測統計学

「推測統計学」の肝は「統計的推測」

2. 母集団と標本

2.1 「標本数」と「標本サイズ」の違い

2.2無作為抽出法 (Random Sampling)

2.3 「推定量」と「推定値」の違い

2.4 標本不偏分散と標本分散

「標本不偏分散」と「標本分散」は使い分ける必要がある

標本不偏分散 \(U^2\)を使う場合 (R ではデフォルト):

標本分散 \(S^2\)を使う場合:

3. 標本分布

3.1 統計量は分布する？

3.2 「母平均」と「標本平均」の関係

3.3 「母標準偏差」と「標本標準偏差」の関係

3.4 大きな母集団を想定したシミュレーション

4. 母平均の推定と信頼区間

4.1 \(z\) 値と \(t\) 値

4.2 Excercise

5. 信頼区間

5.1 信頼区間の計算方法

5.2 信頼区間の解釈

6. 確率変数・確率分布・確率密度関数

6.1 基本事項の確認

6.2 コイントスのシミュレーション

7. 期待値

7.1 期待値 = 平均値

7.2 期待値のシミュレーション

8. Exercise

2.2無作為抽出法 (`Random Sampling`)