『入門・社会統計学――2ステップで基礎から〔Rで〕学ぶ』

第3章　推測統計の基礎

Author

SUGINO Isamu, Build with R4.2.2

Published

October 16, 2024

0 全体の章構成

1 母集団(population)と標本(sample)

1-1 無作為抽出(random sampling)と可能な標本の数

N個からn個を選び出す「場合の数」の公式は次の通りである。

\[{}_N \mathrm{C}_n = \frac{N!}{n!(N-n)!} = \frac{N(N-1)\cdots(N-n+1)}{n!}\]

組合せの数の例示

N=10, n=2とすると次のように計算できる。

\[{}_{10} \mathrm{C}_2 = \frac{10!}{2!(10-2)!} = \frac{10!}{2!8!} = \frac{10\cdot9\cdot8!}{2!8!} = \frac{10\cdot9}{2\cdot1} = 45\]

N=40000, n=100の場合は次のようになる。

\[{}_{40000} \mathrm{C}_{100} = \frac{40000!}{100!(40000-100)!} = \frac{40000!}{100!39900!} = \frac{40000\cdot39999\cdot39998\cdots39902\cdot39901}{100!}\]

これは手計算は無理なので，Rの関数の choose(N, n) を利用する。

choose(40000, 100)

[1] 1.521272e+302

これは，「1.521272×10の302乗」と読む。e+302は10の302乗である。

MS Excelでは？

MS-Excelでは，$combin(N, n)$ という関数を利用する。

しかしこれらの関数は，Nやnがもう少し大きくなると，RやExcelでは（そのままでは）計算不能になる。

Rでの計算結果を以下に示すので理解してほしい。

choose(40000, 100)

[1] 1.521272e+302

choose(45900, 100)

[1] 1.460143e+308

choose(46000, 100)

[1] Inf

choose(40000, 102)

[1] 2.350825e+307

choose(40000, 103)

[1] Inf

“Inf”というのは「無限大」という意味であるが，本当に無限大なのではなく，「簡単には計算できないくらい大きな数」という意味である。

対数変換を利用して莫大な組み合わせの数を計算

　非常に大きな数についての組合せの数を計算するには，「対数変換」を利用する方法がある。場合の数（組合せの数）を以下のように k と置く。

\[k = \frac{N(N-1)\cdots(N-n+1)}{n!}\]

この両辺の常用対数を取り，「積の対数は対数の和に，商（割り算）の対数は対数の差に」変形できることを利用すると，次のように変形できる。

\[\begin{align*} logk &= log\frac{N(N-1)\cdots(N-n+1)}{n!} \\ &= logN(N-1)\cdots(N-n+1)-log(n!)\\ &={logN + log(N-1) + \cdots + log(N-n+1)} - {logn + log(n-1) + \cdots + log2 + log1} \end{align*}\]

対数の和や差であればなんとかRで計算できるので，log(k)を計算してから，10のlogk乗を計算して，kを求めることができる（この計算でも多少近似的かも知れないが，おおよその大きさを把握することはできるだろう）。

N <- 100000000 # 母集団サイズを任意に指定
n <- 1000 # 標本サイズを任意に指定

# N-n+1からNまでのヴェクトルの常用対数の総和
nume <- sum(log10(c((N-n+1):N)))

# 1からnまでのヴェクトルの常用対数の総和
denom <- sum(log10(c(1:n)))

# nume-denomの整数部分を取り出す
b <- floor(nume-denom)

# nume-denomの小数部分だけを10の「べき指数」として計算する
a <- 10^(nume-denom-b)

#　結果を表示する
sprintf("%.5f×10の%d乗", a, b)

[1] "2.47279×10の5432乗"

　これが何をしているかがわからない場合は，c((N-n+1):N)だけを実行してみる，次に log10(c((N-n+1):N))を実行してみる，そして sum(log10(c((N-n+1):N))) を実行してみる，nume - denom を計算した後で，floor(nume - denom) を計算してみる，のように部分的・段階的に実行していくと分かる場合が多い。
　Nやnを書き換えるだけで，このスクリプトを利用することができる。

N <- 100000000 # 母集団サイズを任意に指定
n <- 3000 # 標本サイズを任意に指定

# N-n+1からNまでのヴェクトルの常用対数の総和
nume <- sum(log10(c((N-n+1):N)))

# 1からnまでのヴェクトルの常用対数の総和
denom <- sum(log10(c(1:n)))

# nume-denomの整数部分を取り出す
b <- floor(nume-denom)

# nume-denomの小数部分だけを10の「べき指数」として計算する
a <- 10^(nume-denom-b)

#　結果を表示する
sprintf("%.5f×10の%d乗", a, b)

[1] "2.30400×10の14869乗"

1-2 標本統計量の標本抽出分布(sampling distribution)

　前項のスクリプトを使うと，4,894人から50人を選ぶ場合の数（可能な標本の数）は “7.78123×10の119乗”となる。
　4,894人分の模擬年収データをcsvファイルで提供するので，まずは自分のRのworking directoryを getwd( ) で確認し，必要があれば setwd( ) で変更する。
　download.file( )でworking directoryにcsvファイルを保存する。
このファイルは，実際の調査データに多少の加工を施して作成している（Rで自動でダウンロードしようとするとうまく行かない事もあるので，右クリックで手動でダウンロードしても良い）。
　ファイルをworking directoryに保存したら，次のスクリプトで data01 として読み込もう。ついでに関数names( )，head( )，tail( )を例示する。

# download.file("http://sgn.sakura.ne.jp/text/mock_income.csv", "mock_income.csv")
data01 <- read.csv("mock_income.csv") #csvファイルをデータフレイムに読み込む
str(data01) # data01の構造を表示

'data.frame':   4894 obs. of  2 variables:
 $ case  : int  1 2 3 4 5 6 7 8 9 10 ...
 $ income: int  799 387 605 301 505 243 771 289 359 210 ...

母集団分布

　本来は母集団分布や母平均，母比率，母標準偏差などは未知・不可知であるが，ここではシミュレイションなので（言わば神の視点から）母集団分布と母平均を覗き見してみよう。

本来は不可知の母集団分布を覗き見

hist(data01$income, freq = TRUE, ylab = "人数",
     breaks = seq(0, 2000, by = 50 ),
     main = paste("母集団（サイズ", nrow(data01), "人）における年収のヒストグラム"),
     sub = "※ 赤の線は母平均を示す",
     xlab = paste("母平均", round(mean(data01$income), 1),
                  "万円，母標準偏差", round(sd(data01$income), 1), "万円"),
     col = "#ffcccc80", border="#ffcccc")

segments(mean(data01$income), 0, mean(data01$income), 300, col = "red")

ある一つの標本における値の分布（≠標本分布）

　この収入データ(data01$income)から，無作為に50個の数字を抽出してその平均や標準偏差を計算してみよう。
一般に「標本分布」と訳されるsampling distribution（「標本抽出分布」）はこれのことではない点に注意。

n1 <- 50

sample02 <- sample(data01$income, size = n1)
m1 <- mean(sample02); sd1 <- sd(sample02)
hist(sample02, freq = TRUE, ylab = "人数",
     breaks = seq(0, 2000, by = 50),
     main = paste("ある標本（サイズ",  n1,"人）における年収のヒストグラム\n"),
     xlab = paste("標本平均", round(m1, 1), "（万円），標準偏差",
                  round(sd1, 1), "（万円）"),
     col = "#00640080", border = "#006400")

　無作為抽出なので，実行するたびに結果が変わる。この3行だけを何度か繰り返し実行してみると良い。

標本抽出分布（sampling distribution, 標本分布）

　以下のスクリプトを使うと，本文のヒストグラムに類似したグラフが描ける。実行するたびにグラフは少し変化する。

times <- 10000 # 抽出回数の指定
n <- 50 # サンプルサイズの指定

# times個得られるはずの標本平均の入れ物ヴェクトルを用意
sample_mean <- vector(length = times) 

for (i in 1:times) { # times回だけ以下の操作を繰り返す
   temp_sample <- sample(data01$income, n) # サイズnの標本を抽出する
   sample_mean[i] <- mean(temp_sample) # 抽出した標本の平均を計算する
}

hist(sample_mean, freq = TRUE, ylab = "標本の数",
     main = paste("年収の標本平均のヒストグラム\n（標本サイズ", n, "人を",
                  times, "回抽出）"),
     xlab = paste("年収（万円），母平均", round(mean(data01$income), 1), "万円"),
     col = "#00008040", border = "#000080")

　変数$x$の標本平均$\bar{x}$は，可能な標本全体を考えれば確率変数(random variable)であり，近似的に，以下の様な正規分布に従う。
$\mu$は変数$x$の母平均，$\sigma$は母標準偏差，$n$はサンプルサイズである。

理論的な標本抽出分布

curve(dnorm(x, mean = 439.9, sd = 312/sqrt(50)), xlim = c(250, 650),
      bty = "n", yaxt = "n", xlab = "標本平均（万円）", ylab = "",
      col = "#000080", main = "上記の標本平均の理論的な標本抽出分布")

abline(h = 0, col = "#00000060")
.x <- seq(250, 650 ,by = 1)
segments(.x, 0, .x, dnorm(.x, mean = 439.9, sd = 312/sqrt(50)),
         col = "#00008080")

標準誤差(standard error; SE)

\[\frac{\sigma}{\sqrt{n}}\] 上式・（標本平均の）標本抽出分布の標準偏差は，標本平均が母平均からどの程度ズレてしまうのか（誤差を生じるのか）の指標である為，標準誤差(standard error, s.e., SE)と呼ぶ。

標準誤差と有限母集団修正項

　有限なサイズの母集団から標本を抽出した場合は，標準誤差はより正確には以下の様な式で計算される。

\[s.e.=\sqrt{\frac{N-n}{N-1}}\frac{\sigma}{\sqrt{n}}\] $\sqrt{\frac{N-n}{N-1}}$を，有限母集団修正項(fpc; finite population correction)と呼ぶ。
式をよく見れば分かるが，$1<n<N$の時，fpcは必ず1未満になる。

\[0<\sqrt{\frac{N-n}{N-1}}<1\] つまり，これを計算に含んだ場合，標準誤差は小さくなる。
Nに対するnの割合が小さくなければ，fpcは無視しえない値となるが，nに対してNが非常に大きければ，fpcは殆ど1となるので，計算上は省略してもほぼ同じである。

plot((n <- 1:99), fpc <- sqrt(((N <- 100)- n)/(N-1)), bty="l", 
     xlab="Nに対するnの比(%)", ylab="有限母集団修正項fpc")

具体的に示せば，$N=100,000, n=1,000$の時の$fpc$は約0.995となり，$N=100,000,000, n=3,000$の時は1である。

　fpcを明示すれば，母集団の大きさも（一応）関係すると云う事を示す事が出来，「母集団が大きくても小さくても関係ないと云うのはおかしいのではないか？」と云う素朴な疑問への回答にもなる。
　しかし，敢えて誤差を（僅かに）小さくするfpcを明示しても，現実の調査の誤差は理論的な（最も単純な場合の計算式通りの）誤差よりも大きくならざるを得ないので，fpcを省略する事は，実際の誤差は多少大きくなると云う事も含んだものと見る事も出来るだろう。

仮に母標準偏差が既知であるならば…

　標本平均$\bar{x}$が正規分布に従うならば，確率変数$\bar{x}$を標準化(standardize)したものは，標準正規分布に従うと云う事である。

標準正規分布に従う確率変数（標準正規変量）であるならば，95%の確率で以下の不等式が成り立つ。

この式を単純に変形すると以下の様にもなる。
標本平均が95%の確率で‘落ちる’区間: 未知であるが固定した区間

母平均の95%信頼区間(Confidence ;Interval): 標本が異なれば信頼区間（の計算結果）も変わる。

標準正規分布の描画

　テキストの，標準正規分布の95%区間のグラフを描くスクリプトを（少し簡素化して）示しておくので，pの値を.90や.99などに書き換えて実行してみよう。
コマンドを一つずつ実行していくと，それぞれが何をするコマンドなのかが分かり易くなる。

p <- .95 # この値を .90や.99など書き換えて実行してみよ。
x <- seq(-4,4,by=0.1)
lwr <- qnorm((1-p)/2, mean=0, sd=1)
upr <- qnorm(1-(1-p)/2, mean=0, sd=1)
plot(x, dnorm(x, 0, 1), bty="n", xlab=paste(p*100, "%区間の限界値"),
     ylab="確率密度", main="標準正規分布",
     type="l", ylim=c(0, 0.45), las=1, xaxt="n")
segments(-4, 0, 4, 0)
axis(side=1, at=c(0), las=1)

x2 <- c(seq(-4, lwr, by=.01), seq(upr, 4, by=.01))
axis(side=1, at=c(round(lwr, 3), round(upr, 3)), las=2)
segments(x2, 0, x2, dnorm(x2, 0, 1), col="gray")

1-3 母平均の区間推定(interval estimation)

母標準偏差$\sigma$は通常は未知なので，それを標本から計算された不偏分散$\hat{\sigma}^2$の平方根$\hat{\sigma}$で推定する。 \[\hat{\sigma}^2=\frac{1}{n-1}\sum_{i=1}^{n}(x_{i}-\bar{x})^2,\;\;\;\hat{\sigma}=\sqrt{\frac{1}{n-1}\sum_{i=1}^{n}(x_{i}-\bar{x})^2}\] 　確率変数としての標本平均$\bar{x}$を標準化した式の$\sigma$にこの推定値（代用品）$\hat{\sigma}$を代入すると，標準正規分布ではなく，自由度$n-1$のｔ分布に従う事になる。

ここから求められる95%信頼区間は以下の通りである。

$t_{.025}(n-1)<0$であり，$t_{.975}(n-1)>0$であるので注意。
但し，自由度$n-1$が100を超えれば，標準正規分布の2.5%点，97.5%点の値（±1.96）を用いてほぼ問題ない。

信頼係数の求め方

　99.9%信頼区間を求めるための信頼係数を本文では3.29と書いているが，p <- .999; qnorm(1-(1-p)/2, mean=0, sd=1) とすれば，3.2905267と自分で正確に求められる。
　自由度6のt分布の場合は，p <- .999; qt(1-(1-p)/2, df = 6)とすればrp <- .999; qt(1-(1-p)/2, df = 6)と求められる。自由度100だと，qt(1-(1-p)/2, df = 100)`から3.3904913となる。自由度が大きいと標準正規分布における値と殆ど変わらない。

t分布の描画

　標準正規分布とt分布を重ねて描いているグラフは，以下のスクリプトでほぼ再現できる（本文より少し簡略化してある）。

.df <- 6 # t分布の自由度を書き換えることができる
x <- seq(-4, 4, by=.01)
plot(x, dnorm(x, 0, 1), bty="n", xlab="95%区間の限界値", ylab="確率密度",
     main=paste("標準正規分布（青・実線）と自由度", .df, "のt分布（赤・点線）"),
     type="l", ylim=c(0, 0.45), las=1, xaxt="n", col="#6666FF")
segments(-4, 0, 4, 0)
axis(side=1, at=c(0), las=1)
par(new=T)
plot(x, dt(x, df=.df), bty="n", xlab="", ylab="", axes=F,
     main="", type="l", lty=2, ylim=c(0, 0.45), col="#FF6666")
x2 <- c(seq(-4, qt(.025, .df), by=.01), seq(qt(.975, .df), 4, by=.01))
x3 <- c(seq(-4, qnorm(.025, 0, 1), by=.01), seq(qnorm(.975, 0, 1), 4, by=.01))
segments(x2, 0, x2, dt(x2, .df), col="#FF666660")
axis(side=1, at=c(round(qt(.025, .df),2), round(qt(.975, .df),2)), las=2)
segments(x3, 0, x3, dnorm(x3, 0, 1), col="#6666FF60")
axis(side=1, at=c(round(qnorm(.025, 0, 1),2), round(qnorm(.975, 0, 1),2)), las=2)

　自由度6のt分布の95%限界値は，qt(.975, df=6)より 2.4469119 と求められる。詳細は以下のページを是非参照して欲しい。

確率分布を読み取ったりグラフに描いたりする方法

1-4 区間推定シミュレイション

母平均$\mu$の95%信頼区間の式（「きっとこの式が成り立っているだろう」と考える式。絶対ではない。） \[\bar{x}-1.96\frac{\hat{\sigma}}{\sqrt{n}}<\mu<\bar{x}+1.96\frac{\hat{\sigma}}{\sqrt{n}}\]

　N=4894の母集団からn=200の標本を20回抽出して，その都度信頼区間を求め，それを本文のようなグラフに作図するスクリプトを示す。長いので，前半と後半に分ける。

区間推定シミュレイション（前半）

times <-  20 # 抽出回数の指定
n     <- 200 # サンプルサイズの指定

N <- length(data01$income) # 母集団サイズを算出

sample_mean <- vector(length = times) # 標本平均（群）のヴェクトル
sample_sd   <- vector(length = times) # 標準偏差（群）のヴェクトル
uplimit     <- vector(length = times) # 信頼区間上限のヴェクトル
lowlimit    <- vector(length = times) # 信頼区間下限のヴェクトル

# times回だけ以下の標本抽出と標本統計量の計算の操作を繰り返す
for (i in 1:times) {
     temp_sample    <- sample(data01$income, n)
     sample_mean[i] <- mean(temp_sample)
     sample_sd[i]   <- sd(temp_sample)           # 11/Aug./2017修正
}

　ここまでの前半で，サイズnの標本をtimes回抽出し，それぞれの標本における標本平均と標本標準偏差を計算している。nとtimesは自由に設定出来る。
　以下の後半は，この標本平均と標本標準偏差から信頼区間を計算してグラフにするスクリプトである。
alphaを書き換えてこの後半だけを繰り返すと，同じ20組の標本群について，異なる信頼区間のグラフを得る事が出来る。
　ここでまず，n=200，times=20，alpha＝.95として，前半と後半を実行したグラフを示す。

区間推定シミュレイション（後半）

# ここから先だけを，alphaの値を変えて繰り返す事が出来る
alpha   <- .95 # 信頼水準の設定
t_alpha <- qt((1-alpha)/2, df = n-1, lower.tail = FALSE) # 限界値の計算

# 信頼区間の上限と下限の計算
uplimit  <- sample_mean + t_alpha * sqrt((N-n)/(N-1)) * sample_sd/sqrt(n)
lowlimit <- sample_mean - t_alpha * sqrt((N-n)/(N-1)) * sample_sd/sqrt(n)

# 色の設定
GorR <- rep("#006600", times)
GorR[uplimit < mean(data01$income) | lowlimit > mean(data01$income)] <- "#FF3300"

結果の描画

# 以下はグラフウィンドウの大きさの指定。実行時は行頭の # を取る。
# windows(width=7, height=5)

plot(sample_mean, c(1:times),
     bty = "l", las = 1, xlim = c(340, 540),
     main = sprintf("年収（万円）の区間推定シミュレイション\nN=%d, n=%d", N, n),
     sub = "真ん中の垂直の点線は母平均を表す", 
     xlab = sprintf("各標本から計算された母平均の%.0f%%信頼区間群", alpha*100),
     ylab = "標本抽出回数", 
     col = GorR)

segments(lowlimit, c(1:times), uplimit, c(1:times), col = GorR)
segments(mean(data01$income), 0, mean(data01$income), times, lty = 2)

信頼係数を下げる効果

次に，スクリプトの後半だけを，alpha=.90として実行したグラフを示す。

区間推定シミュレイション（後半）

alpha   <- .90 # 信頼水準の設定
t_alpha <- qt((1-alpha)/2, df=n-1, lower.tail=FALSE) # 限界値の計算

# 信頼区間の上限と下限の計算
uplimit  <- sample_mean + t_alpha * sqrt((N-n)/(N-1)) * sample_sd/sqrt(n)
lowlimit <- sample_mean - t_alpha * sqrt((N-n)/(N-1)) * sample_sd/sqrt(n)

# 色の設定
GorR <- rep("#006600", times)
GorR[uplimit < mean(data01$income) | lowlimit > mean(data01$income)] <- "#FF3300"

90%信頼区間の結果

# windows(width=7, height=5) # グラフウィンドウの大きさの指定

plot(sample_mean, c(1:times),
     bty = "l", las = 1, xlim = c(340, 540),
     main = sprintf("年収（万円）の区間推定シミュレイション\nN=%d, n=%d", N, n),
     sub = "真ん中の垂直の点線は母平均を表す", 
     xlab = sprintf("各標本から計算された母平均の%.0f%%信頼区間群", alpha*100),
     ylab = "標本抽出回数", 
     col = GorR)

segments(lowlimit, c(1:times), uplimit, c(1:times), col = GorR)
segments(mean(data01$income), 0, mean(data01$income), times, lty = 2)

信頼区間がそれぞれ少し狭くなっている事が見て確認出来るだろう。

標本サイズを大きくする効果

　最後に，n=800，alpha＝.95として前半・後半を実行したグラフを示す。
標本抽出からやり直すので上記の二つのグラフとの対応関係はないが，全体的に標本平均のバラつきがぐっと小さくなり，かつ信頼区間が狭くなっている事が見て取れるだろう。
これが，標本サイズを大きくする事で得られる効果である。

区間推定シミュレイション

times <-  20 # 抽出回数の指定
n     <- 800 # サンプルサイズの指定

N <- length(data01$income) # 母集団サイズを算出

sample_mean <- vector(length=times)
sample_sd   <- vector(length=times)
uplimit     <- vector(length=times)
lowlimit    <- vector(length=times)

# times回だけ以下の標本抽出と標本統計量の計算の操作を繰り返す
for (i in 1:times) {
     temp_sample    <- sample(data01$income, n)
     sample_mean[i] <- mean(temp_sample)
     sample_sd[i]   <- sd(temp_sample)         # 11/Aug./2017修正
}

alpha   <- .95 # 信頼水準の設定
t_alpha <- qt((1-alpha)/2, df = n-1, lower.tail = FALSE) # 限界値の計算

# 信頼区間の上限と下限の計算
uplimit  <- sample_mean + t_alpha * sqrt((N-n)/(N-1)) * sample_sd/sqrt(n)
lowlimit <- sample_mean - t_alpha * sqrt((N-n)/(N-1)) * sample_sd/sqrt(n)

# 色の指定
GorR <- rep("#006600", times)
GorR[uplimit < mean(data01$income) | lowlimit > mean(data01$income)] <- "#FF3300"

nを4倍にした95%信頼区間の結果

# windows(width=7, height=5) # グラフウィンドウの大きさの指定
plot(sample_mean, c(1:times),
     bty = "l", las = 1, xlim = c(340, 540),
     main = sprintf("年収（万円）の区間推定シミュレイション\nN=%d, n=%d", N, n),
     sub = "真ん中の垂直の点線は母平均を表す", 
     xlab = sprintf("各標本から計算された母平均の%.0f%%信頼区間群", alpha*100),
     ylab = "標本抽出回数", 
     col = GorR)

segments(lowlimit, c(1:times), uplimit, c(1:times), col = GorR)
segments(mean(data01$income), 0, mean(data01$income), times, lty = 2)

発展1　母比率の区間推定

母比率$\pi$の95%信頼区間の式（ここでも有限母集団修正項は省略）。$p$は標本比率。

\[p-1.96\sqrt{\frac{p(1-p)}{n}}<\pi<p+1.96\sqrt{\frac{p(1-p)}{n}}\]

　二項分布に関する関数を紹介しよう。仮に母比率を pi としたとき，サイズ n の標本において内閣支持者が w 人となる確率は，choose(n, w)*pi^w*(1-pi)^(n-w) として計算出来る。pi = .473，n=10，w=5 とした時の実行結果を以下に示す。

n  <- 10
pi <- .473
w  <- 5
prob <- choose(n, w) * pi^w * (1-pi)^(n-w)
prob

[1] 0.2425266

ここで，以下のように w を少し工夫して，5人以上となる確率をそれぞれ求め，その合計も求めてみる。

n  <- 10
pi <- .473
w   <- 5:n; w

[1]  5  6  7  8  9 10

prob <- choose(n, w) * pi^w * (1-pi)^(n-w)
prob

[1] 0.2425265677 0.1813963734 0.0930338678 0.0313128696 0.0062454090
[6] 0.0005605462

sum(prob)

[1] 0.5550756

この結果は，母比率が47.3%であるとき，10人の標本中5人が支持者である確率が24.3%，6人が支持者である確率が18.1%，10人全員が支持者である確率が0.056%，それらを全て合計した5人以上が支持者である確率が55.5%である事を示している。
　w <- 0:n とすると，有り得る全ての場合について計算される。確率の合計は必然的に1となる。

n  <- 10
pi <- .473
w  <- 0:n; w

 [1]  0  1  2  3  4  5  6  7  8  9 10

prob <- choose(n, w) * pi^w * (1-pi)^(n-w)
prob

 [1] 0.0016523656 0.0148305302 0.0598990200 0.1433636253 0.2251788252
 [6] 0.2425265677 0.1813963734 0.0930338678 0.0313128696 0.0062454090
[11] 0.0005605462

sum(prob)

[1] 1

　同じ計算は，dbinom(w, n, pi) という関数を使って行う事が出来る。

dbinom(w, n, pi)

 [1] 0.0016523656 0.0148305302 0.0598990200 0.1433636253 0.2251788252
 [6] 0.2425265677 0.1813963734 0.0930338678 0.0313128696 0.0062454090
[11] 0.0005605462

sum(dbinom(w, n, pi))

[1] 1

plot(w, prob, type="h") # 最低限の，飾り気のないプロット

plot(w, dbinom(w, n, pi), type="h") # これでも同じグラフ

　これではあまりに素っ気ないので，以下のようにオプションを追加して描いてみる。

plot(w, prob, 
     type = "h", lwd = 5, las = 1, col = "#009999ff",
     bty = "n", ylab = "確率", xlab = "支持者数w", family = "serif",
     main = paste("n=", n, ", π=", pi, "の二項分布グラフ"))

　テクスト本文ではn=25，π=.40のグラフを掲載しているが，以下にはn=50，π=.473の例を掲載しよう。オプション lwd は”line width”で以下ではlwd=2としてある。　

n  <- 50
pi <- .473
w  <- 0:n
prob <- choose(n, w) * pi^w * (1-pi)^(n-w)
plot(w, prob, 
     type = "h", lwd = 2, las = 1, col = "#009999ff",
     bty = "n", ylab = "確率", xlab = "支持者数w", family = "serif",
     main = paste("n=", n, ", π=", pi, "の二項分布グラフ"))

　本文同様，正規分布に非常に似通っている事が確認出来る。
　以下は，x軸を支持者数wではなく，標本比率w/nにしてグラフを描くスクリプトである。

n  <- 100
pi <- .4
w  <- 0:n
prob <- choose(n, w) * pi^w * (1-pi)^(n-w)

plot(w/n, prob, 
     type = "h", lwd = 2, las = 1, col = "#009999ff",
     bty = "n", ylab = "確率", family = "serif", xlab = "標本比率w/n",
     main = paste("標本サイズn=", n, ", 母比率π=", pi, "の二項分布グラフ"))

二値変数の分散

ここで少し脱線して，二値変数(dichotomous variable)の分散の求め方を示しておこう。
　母集団サイズを$N$人とする。或人が現在の内閣を支持していれば$u=1$，そうでなければ$u=0$の値をとる2値変数$u$を考える。これが$N$人全てについてある。$N$人全てについて$u$の平均を計算すると，それが母集団における内閣支持率$π$になる。支持者の数を$N_1$，非支持者の数を$N_0$とすると，$N_1 + N_0 = N$である（旧ヴァージョンの間違いを修正）。　 \[\begin{align*} σ_u^2&=\frac{1}{N}\sum_{i=1}^{n}(u_i-\bar{u})^2\\ &=\frac{1}{N}\left\{\sum_{i=1}^{N_1}(1-π)^2+\sum_{i=1}^{N_0}(0-π)^2\right\}\\ &=\frac{1}{N}\left\{N_1(1-π)^2+N_0π^2\right\}\\ &=\frac{1}{N}\left\{N_1(1-π)^2+(N-N_1)π^2\right\}\\ &=\frac{1}{N}\left\{N_1-2N_1π+N_1π^2+Nπ^2-N_1π^2\right\}\\ &=\frac{N_1-2N_1π+Nπ^2}{N}\\ &=\frac{N_1}{N}-2π\frac{N_1}{N}+π^2\\ &=π-2π^2+π^2\\ &=π-π^2\\ &=π(1-π) \end{align*}\]

二値変数の母分散が$σ_u^2=\frac{1}{N}\sum_{i=1}^{N}(u_i-\bar{u})^2=π(1-π)$である事が示された。標本分散も同様にして示すことが出来る。
　実際に二値変数を発生させて標本分散と不偏分散を比較してみる。

n <- 1248 # 標本サイズ
pi <- .45 # 母比率
u <- rbinom(n, 1, pi)
table(u) # 標本における変数値の分布

u
  0   1 
685 563

mean(u) # 標本平均

[1] 0.4511218

var(u) # Rによる分散（不偏分散）

[1] 0.2478095

mean((u-mean(u))^2) # 定義から計算した標本分散

[1] 0.2476109

n*mean((u-mean(u))^2)/(n-1) # 定義から計算した不偏分散

[1] 0.2478095

標本分散を用いても不偏分散を用いても殆ど違いがない。

　以下，本文の世論調査（内閣支持率p=.412，標本サイズn=1248）の区間推定を行う。
本文よりも少し工夫したスクリプトを挙げるので見比べてみると良い。

# 標本比率と標本サイズの指定
p <- .412; n <- 1248; alpha <- .95

# 標準正規分布を用いて(alpha*100)%信頼区間の上限と下限の計算
lwr <- p - qnorm(1-(1-alpha)/2, 0, 1) * sqrt(p*(1-p)/n)
upr <- p + qnorm(1-(1-alpha)/2, 0, 1) * sqrt(p*(1-p)/n)

lwr; upr # 計算結果の表示

[1] 0.3846927

[1] 0.4393073

有限母集団修正項(fpc)，不偏分散とt分布を用いた場合は下記となる。

# t分布を用いて(alpha*100)%信頼区間の上限と下限の計算
N <- 104000000
lwr2 <- p - qt(1-(1-alpha)/2, n-1) * sqrt((N-n)/(N-1)) * sqrt(p*(1-p)/(n-1))
upr2 <- p + qt(1-(1-alpha)/2, n-1) * sqrt((N-n)/(N-1)) * sqrt(p*(1-p)/(n-1))

lwr2; upr2

[1] 0.3846554

[1] 0.4393446

小標本では近似の性能は悪い

　本文で，n=10，p=.20のような場合は標準正規分布やt分布を用いて計算した信頼区間が不適切である事を述べた。F分布を用いた計算も含めて以下で示す。

# 標本比率と標本サイズの指定
x <- 2; n <- 10; p <- x / n
alpha <- .95

標準正規分布での近似

# 標準正規分布を用いて(alpha*100)%信頼区間の上限と下限の計算
lwr <- p - qnorm(1-(1-alpha)/2, 0, 1) * sqrt(p*(1-p)/n)
upr <- p + qnorm(1-(1-alpha)/2, 0, 1) * sqrt(p*(1-p)/n)
lwr; upr # 計算結果の表示

[1] -0.04791801

[1] 0.447918

t分布での近似

# t分布を用いて(alpha*100)%信頼区間の上限と下限の計算
N <- 104000000
lwr2 <- p - qt(1-(1-alpha)/2, n-1) * sqrt((N-n)/(N-1)) * sqrt(p*(1-p)/(n-1))
upr2 <- p + qt(1-(1-alpha)/2, n-1) * sqrt((N-n)/(N-1)) * sqrt(p*(1-p)/(n-1))
lwr2; upr2 # 計算結果の表示

[1] -0.1016209

[1] 0.5016209

F分布での近似

# F分布で求める信頼区間下限
2*x/(2*x+2*(n-x+1)*qf(.975, 2*(n-x+1), 2*x))

[1] 0.02521073

# F分布で求める信頼区間上限
2*(x+1)*qf(.975, 2*(x+1), 2*(n-x))/(2*(n-x)+2*(x+1)*qf(.975, 2*(x+1), 2*(n-x)))

[1] 0.5560955

二項分布での正確な計算

本来は二項分布を用いて求めるべきものであるので，以下の様に求めてみる。

― pi <- seq(0, 1, by = 0.000001) として，母比率（の候補）pi を0から1までの間で，0.0001%単位で変化させたヴェクトルを作る。
- 標本において，二値変数uが1となる人数をxとする時，関数 pbinom(x, size = n, prob = pi) は，母比率がpiである時にn人の標本においてu=1となる人数がx人以下となる確率を意味する。 - 例えば，pbinom(4, size = 10, prob = .473) の結果は 0.4449244 となる。

95%信頼区間とは，「母比率がpi = pi0である時に，標本比率が左裾2.5%にも右裾2.5%にも含まれない」と云う範囲に由来する。
母比率が定まればLとUも定まるが，母比率は未知である。
調査実施後はp (= x/n)が定まる。 ― 本来はp（またはx）が可変的でpiが定数と考えてpが95%で入る区間を導くが，それをpを固定してpiを可変的と考えて，その不等式を充たすpiの範囲を求めるのが「区間推定」の手法である。
母比率を仮にpi0とした時，その条件の二項分布で右裾の2.5%範囲と左裾の2.5%範囲に標本比率の実現値が含まれなければ，pi0は条件を満たす。
左裾の2.5%範囲に含まれないと云う事は，xまでの累積確率が少なくとも2.5%以上あると云う事。そうでなければxが不等式を充たさない。またx-1までの累積確率は非常に小さい可能性もあるので条件に出来ない。
右裾の2.5%範囲に含まれないと云う事は，x - 1までの累積確率が97.5%よりは小さいと云う事。x - 1までで97.5%に達していれば，xは明らかに右裾の中にある。ただ，右裾の累積確率が非常に小さい可能性もあるので，xまでの累積確率がどの程度に収まるかは不定である。
この条件を表現すると，pbinom(x, n, pi0) >= .025 & pbinom(x - 1, n, pi0) <= .975 となる。
より一般化して，
pbinom(x, n, pi) >= (1-alpha)/2 & pbinom(x - 1, n, pi) <= 1-(1-alpha)/2 とする。
この条件を充たすpiの集合 pi[pbinom(x, n, pi) >= (1-alpha)/2 & pbinom(x - 1, n, pi) <= 1-(1-alpha)/2]が(alpha*100)%信頼区間となるので，その下限minと上限maxを求めればよい。

pi <- seq(0, 1, by = 0.000001) # 母比率を0から1まで，.0001%単位で変化させる
lower <- min(pi[pbinom(x, n, pi) >= (1-alpha)/2 & pbinom(x - 1, n, pi) <= 1-(1-alpha)/2]) # 下限値
upper <- max(pi[pbinom(x, n, pi) >= (1-alpha)/2 & pbinom(x - 1, n, pi) <= 1-(1-alpha)/2]) # 上限値
lower; upper

[1] 0.025211

[1] 0.556095

F分布を使用して求めた区間とほぼ一致している事が確認出来る。もっと精度を粗くして良いなら，by=0.000001の部分をby = 0.00001やby = 0.0001にすると，計算時間がぐっと短縮される。

　標本1,248人中514人が支持者であった場合の信頼区間を，正規分布近似，t分布近似，F分布使用，二項分布スクリプトのそれぞれで求めたのが下記である。

# 標本比率と標本サイズ，信頼係数などの指定
n <- 1248; x <- 512; p <- x/n; alpha <- .95

# 標準正規分布を用いて(alpha*100)%信頼区間の上限と下限の計算
lwr <- p - qnorm(1-(1-alpha)/2, 0, 1)*sqrt(p*(1-p)/n)
upr <- p + qnorm(1-(1-alpha)/2, 0, 1)*sqrt(p*(1-p)/n)
lwr; upr

[1] 0.3829666

[1] 0.4375462

# t分布を用いて(alpha*100)%信頼区間の上限と下限の計算
N <- 104000000
lwr2 <- p - qt(1-(1-alpha)/2, n-1)*sqrt((N-n)/(N-1))*sqrt(p*(1-p)/(n-1))
upr2 <- p + qt(1-(1-alpha)/2, n-1)*sqrt((N-n)/(N-1))*sqrt(p*(1-p)/(n-1))
lwr2; upr2

[1] 0.3829293

[1] 0.4375835

# F分布で求める信頼区間下限と下限
2*x/(2*x+2*(n-x+1)*qf(.975, 2*(n-x+1), 2*x)); 2*(x+1)*qf(.975, 2*(x+1), 2*(n-x))/(2*(n-x)+2*(x+1)*qf(.975, 2*(x+1), 2*(n-x)))

[1] 0.3828048

[1] 0.4381318

# 二項分布スクリプト
pi <- seq(0, 1, by=0.00001)
lower <- min(pi[pbinom(x, n, pi) >= (1-alpha)/2 & pbinom(max(0, x-1), n, pi) <= 1-(1-alpha)/2])
upper <- max(pi[pbinom(x, n, pi) >= (1-alpha)/2 & pbinom(max(0, x-1), n, pi) <= 1-(1-alpha)/2])
lower; upper

[1] 0.38281

[1] 0.43813

いずれを用いても実用上はほぼ問題無い程度の違いと言えるだろう。

発展2　無回答誤差(non-response error)

　本書ではユニット無回答(unit non-response)や項目無回答(item non-response)については殆ど扱っていないので，ここでだけほんの少し無回答の影響について検討してみよう。

回収率（回答率；response rate）と無回答バイアス(non-response bias)

　ごく簡単な例として，架空の単純なモデルを考えてみよう。対象となる社会が高収入群と低収入群１：１であるとする。簡便化の為に，高収入群の平均年収を500（万円）に，高収入群の回収率（ユニット回答率）を50%に固定する。収入格差をdM（万円），回収率の差をdr（ポイント）とし，この二つの値だけを変化させて比較する。詳しくコメントを付けたスクリプトは下記である。

Mh <- 500; rh <- .50 # 高収入群
dM <- 200; dr <- .20 # 二群の差

# 母集団に於ける平均年収は，当然中間値になる。
M <- (Mh + (Mh -dM))/2; M # 母平均

[1] 400

# 二群の母平均差と全体母平均の比, 及び高所得群平均の比
dM/M; dM/Mh

[1] 0.5

[1] 0.4

# 標本平均。それぞれの設計標本サイズをn（同比率なので共通）とすると，
# (Mh*n*rh + (Mh-dM)*n*(rh-dr))/(n*rh + n*(rh-dr))
# (Mh*rh + (Mh-dM)*(rh-dr))/(rh + (rh-dr))
m <- (Mh*rh + (Mh-dM)*(rh-dr))/(rh + (rh-dr)); m

[1] 425

# 回収率差，標本平均と真の母平均の差，及びその差の母平均に対する比
dr; m - M; (m - M)/M

[1] 0.2

[1] 25

[1] 0.0625

この数値例では，二群の平均収入差は全体平均の50%，高収入群の平均の40%に相当し，回収率の差は20%ポイント，母平均400万円に対して標本平均は425万円となり，バイアスは25万円，バイアスの母平均に対する比は6.25%となる。　収入格差と回収率差を変化させるとどうなるかを以下に示すので自分で読み解いてみよう。

### 繰り返し比較用
dM <- 200; dr <- .20 # 収入格差も回答率差も中程度
M <- (Mh + (Mh -dM))/2
m <- (Mh*rh + (Mh-dM)*(rh-dr))/(rh + (rh-dr))
M; dM/M; dM/Mh

[1] 400

[1] 0.5

[1] 0.4

dr; m; m - M; (m - M)/M

[1] 0.2

[1] 425

[1] 25

[1] 0.0625

dM <- 100; dr <- .40 # 収入格差小，回答率差大
M <- (Mh + (Mh -dM))/2
m <- (Mh*rh + (Mh-dM)*(rh-dr))/(rh + (rh-dr))
M; dM/M; dM/Mh

[1] 450

[1] 0.2222222

[1] 0.2

dr; m; m - M; (m - M)/M

[1] 0.4

[1] 483.3333

[1] 33.33333

[1] 0.07407407

dM <- 400; dr <- .10 # 収入格差大，回答率差小
M <- (Mh + (Mh -dM))/2
m <- (Mh*rh + (Mh-dM)*(rh-dr))/(rh + (rh-dr))
M; dM/M; dM/Mh

[1] 300

[1] 1.333333

[1] 0.8

dr; m; m - M; (m - M)/M

[1] 0.1

[1] 322.2222

[1] 22.22222

[1] 0.07407407

項目無回答と欠損値補定（補綴；imputation）

　やや分かりにくいが，教育年数と年収の模擬データを発生させ，そこから人為的に欠損値を発生させ，どの様な結果が得られるかを例示してみる。　まず教育年数eduを適当に尤もらしい確率で乱数発生させ，次にそのeduと強く順相関する年収データincome0を生成する。更に，年収が低い程高い値を取りやすい0～1の範囲の変数missを生成させて，missの上位10%のケースについて，年収だけをNAにした変数incomeを生成する。最後に，それらの要約統計量と相関係数行列を示している。

n <- 200 # ケース数は適当に指定
# 模擬教育年数
edu <- sample(c(9, 12, 14, 16, 18), n, replace=TRUE, c(.05, .40, .15, .35, .05))
# 教育年数に相関する模擬年収データ（試行錯誤で工夫したので気にしなくて良い）
income0 <- round(scale(edu)+abs(rnorm(n, 0, 3)), 1)*50 + 300 

# ここで，低収入群からのみ，収入の項目欠損をランダムに発生させる
miss0 <- (scale(income0 + rnorm(n, 0, 50)))
miss <- (max(miss0) - miss0)/(max(miss0)-min(miss0))
income <- income0

income[miss > quantile(miss, probs=c(.9))] <- NA # 全体の10%について年収欠損の生成
# 関連する4つの変数をデータフレイムに
data01 <- data.frame(edu, "income_full"=income0, "income_NA"=income, "missing"=miss)
summary(data01)

      edu         income_full      income_NA        missing      
 Min.   : 9.00   Min.   :200.0   Min.   :215.0   Min.   :0.0000  
 1st Qu.:12.00   1st Qu.:342.5   1st Qu.:358.8   1st Qu.:0.5206  
 Median :14.00   Median :400.0   Median :410.0   Median :0.6286  
 Mean   :13.79   Mean   :411.1   Mean   :426.6   Mean   :0.6109  
 3rd Qu.:16.00   3rd Qu.:470.0   3rd Qu.:481.2   3rd Qu.:0.7211  
 Max.   :18.00   Max.   :765.0   Max.   :765.0   Max.   :1.0000  
                                 NA's   :20

cor(data01, use="pairwise") # 相関係数行列

                   edu income_full  income_NA    missing
edu          1.0000000   0.4920285  0.3855618 -0.4250250
income_full  0.4920285   1.0000000  1.0000000 -0.9049359
income_NA    0.3855618   1.0000000  1.0000000 -0.8882862
missing     -0.4250250  -0.9049359 -0.8882862  1.0000000

　当然，欠損を含んだ年収incomeは分布が全体に高い方にずれ，教育年数eduとの積率相関は完備データ(income_full)より低下している。本来は教育年数と年収はもっと強く相関しているのに，低収入層が不釣り合いに分析から除外される為に，相関が弱まって見えてしまっていると云う（仮想的）状況である。これは「完備データ分析(complete case analysis)」の問題を示している。

最も簡単な欠損値補定の例

　ここで，最も簡単な欠損値補定(imputation)の例を示して，より複雑な補定が必要である事だけ示唆しておこう。　まずは有効ケースの平均でNAを置き換える方法の結果を示す。スクリプトと結果から解読してみよう。

# 有効ケースの平均値によるimputation
income_M <- income # income_NAをcopy
M <- mean(income, na.rm=T) # 有効ケースの平均値
income_M[is.na(income)] <- M # NAならMを代入
data_M <- data.frame(edu, "income_full"=income0, "income_NA"=income, "income_M"=income_M)
summary(data_M); cor(data_M, use="pairwise")

      edu         income_full      income_NA        income_M    
 Min.   : 9.00   Min.   :200.0   Min.   :215.0   Min.   :215.0  
 1st Qu.:12.00   1st Qu.:342.5   1st Qu.:358.8   1st Qu.:363.8  
 Median :14.00   Median :400.0   Median :410.0   Median :420.0  
 Mean   :13.79   Mean   :411.1   Mean   :426.6   Mean   :426.6  
 3rd Qu.:16.00   3rd Qu.:470.0   3rd Qu.:481.2   3rd Qu.:470.0  
 Max.   :18.00   Max.   :765.0   Max.   :765.0   Max.   :765.0  
                                 NA's   :20

                  edu income_full income_NA  income_M
edu         1.0000000   0.4920285 0.3855618 0.3505670
income_full 0.4920285   1.0000000 1.0000000 0.8870631
income_NA   0.3855618   1.0000000 1.0000000 1.0000000
income_M    0.3505670   0.8870631 1.0000000 1.0000000

var(income0, na.rm=T); var(income, na.rm=T); var(income_M, na.rm=T)

        [,1]
[1,] 10845.2

         [,1]
[1,] 9487.383

         [,1]
[1,] 8533.877

よく考えれば当然だが，分布はより平均に集中してバラつきが小さくなり，教育年数eduとの相関も更に小さくなっている。主に低い年収であった筈のケースが全体平均年収で補定されているのであるから当然そうなる。欠測がランダムではない上の例の限りでは，平均値で補定する利点は無い。
　次に，収入の欠測を，収入を教育年数に回帰させた結果を用いた予測値で置換する方法を示す。偶然誤差（攪乱項）を含まない回帰補定である。回帰分析の結果を利用する為に，summary( )関数，names( )関数などで確認をしている部分も敢えて表示させている。それぞれの意味をしっかり解読しよう。
　まずは回帰分析の結果を確認する部分。セミコロン”;“で複数のコマンドを一行に書いて一度に結果を表示させているので区切りが分からないかも知れない。そう云う場合はセミコロンを使わずコマンドを一つ一つ実行して確認する。

# 誤差を含まない回帰補定
income_R <- income
summary(reg <- lm(income ~ edu)); names(reg); reg$coefficients


Call:
lm(formula = income ~ edu)

Residuals:
    Min      1Q  Median      3Q     Max 
-127.80  -64.47  -22.59   52.61  306.50 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  194.041     42.258   4.592 8.28e-06 ***
edu           16.529      2.965   5.575 9.05e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 90.12 on 178 degrees of freedom
  (20 observations deleted due to missingness)
Multiple R-squared:  0.1487,    Adjusted R-squared:  0.1439 
F-statistic: 31.08 on 1 and 178 DF,  p-value: 9.045e-08

 [1] "coefficients"  "residuals"     "effects"       "rank"         
 [5] "fitted.values" "assign"        "qr"            "df.residual"  
 [9] "na.action"     "xlevels"       "call"          "terms"        
[13] "model"

(Intercept)         edu 
  194.04083    16.52888

income_R[is.na(income)]<- reg$coefficients[2]*edu[is.na(income)] + reg$coefficients[1]

　次がここでの本来の目的の回帰補定の結果と比較である。平均値補定程はバラつきは過小にならないが，それでも本来の分布からすると集中度が高く，平均値も未だ高い。教育年数との相関も，平均値補定程は弱くならないが，本来の相関からするとまだ弱い。欠損値生成（miss0が登場する行以降）をやり直すと結果は変化するが，何度か繰り返してもこの傾向は変わらない。自分で確認してみると良い。

data_R <- data.frame(edu, "income_full"=income0, "income_M"=income_M, "income_R"=income_R)
summary(data_R); cor(data_R, use="pairwise")

      edu         income_full       income_M        income_R    
 Min.   : 9.00   Min.   :200.0   Min.   :215.0   Min.   :215.0  
 1st Qu.:12.00   1st Qu.:342.5   1st Qu.:363.8   1st Qu.:355.0  
 Median :14.00   Median :400.0   Median :420.0   Median :400.0  
 Mean   :13.79   Mean   :411.1   Mean   :426.6   Mean   :421.9  
 3rd Qu.:16.00   3rd Qu.:470.0   3rd Qu.:470.0   3rd Qu.:470.0  
 Max.   :18.00   Max.   :765.0   Max.   :765.0   Max.   :765.0

                  edu income_full  income_M  income_R
edu         1.0000000   0.4920285 0.3505670 0.4175936
income_full 0.4920285   1.0000000 0.8870631 0.9491640
income_M    0.3505670   0.8870631 1.0000000 0.9847747
income_R    0.4175936   0.9491640 0.9847747 1.0000000

var(income0, na.rm=T); var(income, na.rm=T); var(income_M, na.rm=T); var(income_R, na.rm=T)

        [,1]
[1,] 10845.2

         [,1]
[1,] 9487.383

         [,1]
[1,] 8533.877

         [,1]
[1,] 8799.796

　このほか，回帰補定に誤差（攪乱項）を交える方法や，複数の値で補定して異なるデータセット・分析結果を作り出して，その結果を統合する多重代入法(multiple imputation)などの方法がある。注目度が高いのは最後の多重代入法であるが，高度な内容になるのでここでは割愛する。

第3章の【練習問題】の解答

1) の答え

pt(2.00, df=200) - pt(-1.00, df=200)

[1] 0.817314

2) の答え

qt(.025, df=1000); qt(.975, df=1000)

[1] -1.962339

[1] 1.962339

3) の答え

m <- 60; v <- 120; n <- 285
m - qt(.95, df=n-1)*sqrt(v/n)

[1] 58.92919

m + qt(.95, df=n-1)*sqrt(v/n)

[1] 61.07081

4) の答え

　標準正規分布近似，F分布計算，二項分布スクリプトの実施例を示す。

n <- 800; p <- .45; alpha <- .95
# 標準正規分布を用いて(alpha*100)%信頼区間の上限と下限の計算
lwr <- p - qnorm(1-(1-alpha)/2, 0, 1)*sqrt(p*(1-p)/n)
upr <- p + qnorm(1-(1-alpha)/2, 0, 1)*sqrt(p*(1-p)/n)
lwr; upr

[1] 0.4155261

[1] 0.4844739

# t分布を用いて(alpha*100)%信頼区間の上限と下限の計算
N <- 106200000
lwr2 <- p - qt(1-(1-alpha)/2, n-1)*sqrt((N-n)/(N-1))*sqrt(p*(1-p)/(n-1))
upr2 <- p + qt(1-(1-alpha)/2, n-1)*sqrt((N-n)/(N-1))*sqrt(p*(1-p)/(n-1))
lwr2; upr2

[1] 0.4154523

[1] 0.4845477

# F分布で求める信頼区間下限と下限
2*x/(2*x+2*(n-x+1)*qf(.975, 2*(n-x+1), 2*x)); 2*(x+1)*qf(.975, 2*(x+1), 2*(n-x))/(2*(n-x)+2*(x+1)*qf(.975, 2*(x+1), 2*(n-x)))

[1] 0.6056455

[1] 0.6733184

# 二項分布スクリプト
x <- n*p
pi <- seq(0, 1, by=0.00001)
lower <- min(pi[pbinom(x, n, pi) >= (1-alpha)/2 & pbinom(max(0, x-1), n, pi) <= 1-(1-alpha)/2])
upper <- max(pi[pbinom(x, n, pi) >= (1-alpha)/2 & pbinom(max(0, x-1), n, pi) <= 1-(1-alpha)/2])
lower; upper

[1] 0.41515

[1] 0.48522

いずれを用いるかによって，41.6%～48.4% もしくは 41.5%～48.5%と微妙にずれるが，実際に必要な／期待される精度からすれば問題とする差ではない。

5) の答え

問題文に書かれたスクリプトを実際に実行してみるのが課題である。以下に少し変えて再掲する。

population <- rnorm(n=10000, mean=50, sd=10)

mean(population); sd(population); quantile(population)

[1] 50.05046

[1] 9.973484

       0%       25%       50%       75%      100% 
 6.655032 43.388484 49.954925 56.884866 88.692345

hist(population)

boxplot(population)

n <- 100
s1 <- sample(population, n)

mean(s1) # 標本平均を求める
mean(s1)-qnorm(.975, 0, 1)*sd(s1)/sqrt(n) # 母平均の95%信頼区間の下限
mean(s1)+qnorm(.975, 0, 1)*sd(s1)/sqrt(n) # 母平均の95%信頼区間の上限

[1] 48.61579
[1] 46.81086
[1] 50.42072

問題文にある様に，サンプリングの部分からだけ繰り返してみる。

s1 <- sample(population, n)

mean(s1) # 標本平均を求める
mean(s1)-qnorm(.975, 0, 1)*sd(s1)/sqrt(n) # 母平均の95%信頼区間の下限
mean(s1)+qnorm(.975, 0, 1)*sd(s1)/sqrt(n) # 母平均の95%信頼区間の上限

[1] 48.64357
[1] 46.65244
[1] 50.63471

s1 <- sample(population, n)

mean(s1) # 標本平均を求める
mean(s1)-qnorm(.975, 0, 1)*sd(s1)/sqrt(n) # 母平均の95%信頼区間の下限
mean(s1)+qnorm(.975, 0, 1)*sd(s1)/sqrt(n) # 母平均の95%信頼区間の上限

[1] 49.43602
[1] 47.51688
[1] 51.35516

s1 <- sample(population, n)

mean(s1) # 標本平均を求める
mean(s1)-qnorm(.975, 0, 1)*sd(s1)/sqrt(n) # 母平均の95%信頼区間の下限
mean(s1)+qnorm(.975, 0, 1)*sd(s1)/sqrt(n) # 母平均の95%信頼区間の上限

[1] 50.49131
[1] 48.5508
[1] 52.43182