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では，\(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), "人）における年収のヒストグラム"),
     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")

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")

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

sample_mean <- vector(length=times) # 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")

curve(dnorm(x, mean=439.9, sd=312/sqrt(50)), from=250, to=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")

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

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\)のｔ分布に従う事になる。
\[\frac{\color{red}{\bar{x}}-\color{blue}{\mu}}{\frac{\color{red}{\hat{\sigma}}}{\sqrt{\color{green}{n}}}} \sim t_{(df=n-1)}\] ここから求められる95%信頼区間は以下の通りである。

\[t_{.025}(n-1)<\frac{\color{red}{\bar{x}}-\color{blue}{\mu}}{\frac{\color{red}{\hat{\sigma}}}{\sqrt{\color{green}{n}}}} < t_{.975}(n-1)\] \[\color{red}{\bar{x}}-t_{.975}(n-1)\frac{\color{red}{\hat{\sigma}}}{\sqrt{\color{green}{n}}} < \color{blue}{\mu} < \color{red}{\bar{x}} - t_{.025}(n-1)\frac{\color{red}{\hat{\sigma}}}{\sqrt{\color{green}{n}}}\] \(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と求められる。
　標準正規分布と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 と求められる。詳細は以下のページを是非参照して欲しい。

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

標準正規分布と，自由度によって変わるt分布

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),
     sub="真ん中の垂直の点線は母平均を表す", ylab="標本抽出回数", col=GorR,
     xlab=sprintf("各標本から計算された母平均の%.0f%%信頼区間群", alpha*100),
     main=sprintf("年収（万円）の区間推定シミュレイション\nN=%d, n=%d", N, n))
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"

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

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

二値変数の分散

ここで少し脱線して，二値変数(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 
704 544 

mean(u) # 標本平均

[1] 0.4358974

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

[1] 0.246088

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

[1] 0.2458909

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

[1] 0.246088

標本分散を用いても不偏分散を用いても殆ど違いがない。
　以下，本文の世論調査（内閣支持率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分布を用いた計算も含めて以下で示す。

# 標本比率と標本サイズの指定
p <- .20; n <- 10; 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分布を用いて(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

x <- n*p
# 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%単位で変化させたヴェクトルを作る。
x=n*p（標本において，二値変数uが1となる人数）である時，関数 pbinom(x, size=n, prob=pi) は，母比率がpiである時にn人の標本においてu=1となる人数がx人以下となる確率を意味する。
例えば，pbinom(4, size=10, prob=.473) の結果は 0.4449244 となる。
　95%信頼区間というのは，「母比率がpiである時に，標本比率がx/n以下（つまり該当者がx人以下）になる確率が2.5%以上であり，かつ標本比率がx/n以上になる確率も2.5%以上になる」範囲に由来する。
　該当者がx人以下になる確率が2.5%以上という条件は pbinom(x, n, pi) >= .025 と書ける。該当者がx人以上になる確率が2.5%以上という条件は，該当者がx人未満（x-1人以下）になる確率が97.5%以下になる事と等しいので，
pbinom(x-1, n, pi) <= .975 と書ける。ただし，x-1がマイナスになるとエラーが生じるので，
pbinom(max(0, x-1), n, pi) <= .975 としなければならない。
この条件を合わせると，pbinom(x, n, pi) >= .025 & pbinom(max(0, x-1), n, pi) <= .975 となる。
より一般化して，
pbinom(x, n, pi) >= (1-alpha)/2 & pbinom(max(0, x-1), n, pi) <= 1-(1-alpha)/2 とする。
この条件を充たすpiの集合 pi[pbinom(x, n, pi) >= (1-alpha)/2 & pbinom(max(0, x-1), n, pi) <= 1-(1-alpha)/2]が(alpha*100)%信頼区間となるので，その下限minと上限maxを求めればよい。

p <- .20; n <- 10; alpha <- .95
pi <- seq(0, 1, by=0.000001) # 母比率を0から1まで，.0001%単位で変化させる
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.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.000001)
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.382805

[1] 0.438131

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

回収率（回答率；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.   :205.0   Min.   :255   Min.   :0.0000  
 1st Qu.:12.00   1st Qu.:340.0   1st Qu.:355   1st Qu.:0.4640  
 Median :12.00   Median :390.0   Median :400   Median :0.5821  
 Mean   :13.46   Mean   :404.2   Mean   :416   Mean   :0.5741  
 3rd Qu.:16.00   3rd Qu.:460.0   3rd Qu.:470   3rd Qu.:0.7119  
 Max.   :18.00   Max.   :825.0   Max.   :825   Max.   :1.0000  
                                 NA's   :20                    

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

                   edu income_full  income_NA    missing
edu          1.0000000   0.4913260  0.4282038 -0.4177008
income_full  0.4913260   1.0000000  1.0000000 -0.8758243
income_NA    0.4282038   1.0000000  1.0000000 -0.8612508
missing     -0.4177008  -0.8758243 -0.8612508  1.0000000

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

# 有効ケースの平均値による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.   :205.0   Min.   :255   Min.   :255  
 1st Qu.:12.00   1st Qu.:340.0   1st Qu.:355   1st Qu.:360  
 Median :12.00   Median :390.0   Median :400   Median :415  
 Mean   :13.46   Mean   :404.2   Mean   :416   Mean   :416  
 3rd Qu.:16.00   3rd Qu.:460.0   3rd Qu.:470   3rd Qu.:460  
 Max.   :18.00   Max.   :825.0   Max.   :825   Max.   :825  
                                 NA's   :20                 

                  edu income_full income_NA  income_M
edu         1.0000000   0.4913260 0.4282038 0.3950529
income_full 0.4913260   1.0000000 1.0000000 0.9102541
income_NA   0.4282038   1.0000000 1.0000000 1.0000000
income_M    0.3950529   0.9102541 1.0000000 1.0000000

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

         [,1]
[1,] 8588.889

         [,1]
[1,] 7911.564

         [,1]
[1,] 7116.432

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

Call:
lm(formula = income ~ edu)

Residuals:
    Min      1Q  Median      3Q     Max 
-118.24  -64.86  -13.26   45.14  370.14 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  188.368     36.505   5.160 6.54e-07 ***
edu           16.656      2.635   6.322 2.02e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 80.61 on 178 degrees of freedom
  (20 observations deleted due to missingness)
Multiple R-squared:  0.1834,    Adjusted R-squared:  0.1788 
F-statistic: 39.97 on 1 and 178 DF,  p-value: 2.019e-09

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

(Intercept)         edu 
  188.36823    16.65598 

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

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.   :205.0   Min.   :255   Min.   :255.0  
 1st Qu.:12.00   1st Qu.:340.0   1st Qu.:360   1st Qu.:355.0  
 Median :12.00   Median :390.0   Median :415   Median :390.0  
 Mean   :13.46   Mean   :404.2   Mean   :416   Mean   :412.6  
 3rd Qu.:16.00   3rd Qu.:460.0   3rd Qu.:460   3rd Qu.:460.0  
 Max.   :18.00   Max.   :825.0   Max.   :825   Max.   :825.0  

                  edu income_full  income_M  income_R
edu         1.0000000   0.4913260 0.3950529 0.4568697
income_full 0.4913260   1.0000000 0.9102541 0.9561850
income_M    0.3950529   0.9102541 1.0000000 0.9843434
income_R    0.4568697   0.9561850 0.9843434 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,] 8588.889

         [,1]
[1,] 7911.564

         [,1]
[1,] 7116.432

         [,1]
[1,] 7344.615

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] 49.94621

[1] 10.16731

      0%      25%      50%      75%     100% 
13.10416 43.04637 49.93063 56.82562 89.15872 

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] 50.60007
[1] 48.61015
[1] 52.58998

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

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.05203
[1] 48.13919
[1] 51.96487

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.6906
[1] 48.83094
[1] 52.55026

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.26789
[1] 48.30576
[1] 52.23002