Rで統計学を学ぶ(3)

• 母集団: 対象のデータ全体　　　　　　　　　　　　　　　
• 標本: 母集団から一部のデータを取り出した(抽出した)もの　　　　　　　　　　　　　　　　
• 標本抽出: 母集団から標本をとりだすこと
• 母数: 母集団の性質をあらわす統計的指数（比率、平均、分散、相関係数など）

推測統計

手元にある標本から母集団の性質(母数)を推測することを推測統計といいます。 推測統計は推定と検定に分けられます。

「推定」は、母数の値に対し具体的な値を求める点推定、つまり一つの値で推定の結果を表すものと、区間推定、つまりある幅を持った区間で結果を表すものに分けられます。

また「検定」は母集団について述べた異なる立場の二つの主張(仮説)のうちどれを採用するかを決定するものです。例えば、日本の中学生のテスト得点が5年前と現在とで変化したという仮説と、変化していないという仮説がある場合、どちらを採択するかをきめるのが検定です。

例えば、10人の17歳男性の身長データがあるとします。これは日本人の17歳男性全体という『母集団』から抽出された『標本』とみなせます。そしてこのとき標本のデータの個数が10ですから、n = 10と書きます。ここでnはサンプルサイズを表す記号です。

身長 <- c(165.2, 175.9, 161.7, 174.2, 172.1, 163.3, 170.6, 168.4, 171.3)

この標本から、「17歳の日本人男性全体の身長」(母集団)の平均(母数の一つ)を点推定してみましょう。 それには、次のようにします。

mean(身長)

• 標本統計量：標本データから計算される平均や分散などの値
• 推定値：標本データを用いて計算した母数の推定量

• 実際の母数の値にどの位近い推定値が得られるか？
• 推定の結果はどのくらい信用できるか？

• 標本誤差：　母数の値と推定値とのズレ--- データを沢山集めれば、生じうる誤差は小さくできる
• 標本分布: どのような推定値が得られる可能性があるか、についての手がかり

1. 標本抽出の方法--- 単純無作為抽出と呼ばれる方法が前提
   • 単純無作為抽出: 母集団の中のどのデータにも平等に選ばれる可能性を持っているような標本抽出の方法
   • 無作為標本: 単純無作為抽出によって得られた標本
2. 単純無作為抽出によって得られるデータの性質としての確率変数
   
   • 確率変数: 実際に結果が得られるまでどのような値が得られるかが決まっていない変数
     例： サイコロを振った時の出目
     同じ手続でデータを得ても結果に再現性がない⇒ 単純無作為抽出によって得られるデータは確率変数といえる
3. 確率変数がどのような値をとるのかを示す確率分布
   
   • 確率変数の実現値：抽出の結果、確率変数がとる値
   • 確率分布: 確率変数がどのような値をどのような確率でとるかということを表した分布 [サイコロの出目の確率分布]
   理論的には、サイコロの出目の確率分布はつぎのようになります:

   サイコロの出る目	1	2	3	4	5	6
   確率	1/6	1/6	1/6	1/6	1/6	1/6

   
   それに対し、次のシミュレーションにより、サイコロの出る目の度数分布は次のようになりました：

    set.seed(1)
    dice6times <- ceiling(runif(6,0,6))
    table(dice6times)
    size=60000
    dice60000times  <- ceiling(runif(size,0,6))
    table(dice60000times)
    table(dice60000times)/(size/6)
   

   (1) 6回サイコロを振った場合:

   サイコロの目	1	2	3	4	5	6
   度数	0	2	1	1	0	2
   相対度数	0	2/6	1/6	1/6	0	2/6

   
   (2) 60,000回サイコロを振った場合:

   サイコロの目	1	2	3	4	5	6
   度数	10,097	9,918	9,936	9,960	10,050	10,039
   相対度数	1.0097	0.9918	0.9936	0.9960	1.0050	1.0039

   このように、非常に多くの実現値が得られれば、度数分布は確率分布とほぼ同じになります
4. 確率分布を用いた母集団の表現としての母集団分布
   • 母集団分布：ある変数の母集団における分布
     無作為抽出の場合、標本データの確率分布は母集団分布と同じになります
5. 代表的な母集団分布である正規分布
6. Rを使って正規分布の母集団から標本を抽出する方法

正規分布の表現方法：　N (μ,σ²) --- μは平均値、σ²は分散 (σは標準偏差)

curve(dnorm(x,mean=0,sd=1), from=-4, to=4)    		# この場合mean=, sd=, from=, to= は省略可能
curve(dnorm(x,mean=1,sd=1), add=TRUE)   		# add=TRUE により出力されたグラフに重ね書きできる

正規分布 N (0,1)のグラフ	正規分布 N (1,1) のグラフを重ね書き

• 標準正規分布： 平均0 分散1の正規分布 --- N（0,1）
• 離散変数：　整数など、とびとびの値をとる変数（例：サイコロの出目)
• 連続変数：　実数値をとる変数（例：正規分布に従う確率変数）
• 確率密度 : 正規分布の図の縦軸は確率密度
  　 　確率密度×確率変数の範囲＝確率
• 確率密度関数：確率密度を確率変数の値の関数として表したもの　

• μ-σからμ+σ の範囲(偏差値で言えば、40点から60点の範囲)の確率はおよそ 0.683
• μ-2σからμ+2σ の範囲(偏差値で言えば、30点から70点の範囲)の確率はおよそ 0.954
• μ-3σからμ+3σ の範囲(偏差値で言えば、20点から80点の範囲)の確率はおよそ 0.997

μ-σからμ+σ の範囲	μ-2σからμ+2σ の範囲	μ-3σからμ+3σの範囲

hnorm <- function(range, color) {
   plot(dnorm, -4, 4)
   xvals <- seq(-range, range, length=30)
   dvals <- dnorm(xvals)
   polygon(c(xvals, rev(xvals)), c(rep(0,30), rev(dvals)), col=color)
   abline(v=0, lty=2)
   abline(v=-range, lty=3)
   abline(v=range, lty=3)
   abline(h=0, lty=2)
}
# 使用例: hnorm(2, "blue")

> 標本 <- rnorm(n=5,mean=50,sd=10)
> 標本
[1]59.87265 55.41718 34.19793 52.95355 59.10080
> hist(標本)
> 大標本 <- rnorm(n=10000,mean=50,sd=10)
> hist(大標本)

標本(n=5)	大標本(n=10000)

標本分布は標本における個々のデータの実現値を表した度数分布ではなく、標本統計量の確率分布
標本分布は母集団分布と標本統計量の種類、そしてサンプルサイズが決まると理論的（数学的）に導かれる。 実際のデータから作成されるわけではない

1. 標本分布が母数の本当の値を中心として分布しているか⇒平均
2. 標本分布が横に大きく広がっていないか⇒標準偏差

• 母集団分布がどのような分布になるかが不明 　
• 母集団分布がわからないと、標本から得られる推定値の信頼性も不明
  ⇒　母集団が正規分布であると仮定して「もし母集団分布がこのような正規分布ならば、このくらいあてになる推定値が得られる」ということを検討

> 標本 <- rnorm(n=10,mean=50,sd=10)
> 標本			　# データの確認
　[1] 43.17949  39.07841  51.59838  44.00960  42.78063  53.01545
  [7] 47.99463  48.81240  46.76994  65.39066
> mean(標本)
　[1] 48.26296

標本平均　<- numeric(length=10000)　　　　	# 推定値を格納する場所を予約
for( i  in  1:10000 ){   			# 括弧の中を10000回処理する
　　標本 <- rnorm(n=10,mean=50,sd=10)			# 正規分布に基づく標本を生成
    標本平均[i] <-  mean(標本)   			# 標本平均を計算
}

hist(標本平均, freq=FALSE)			# 度数ではなく割合をy軸として表示させる
curve(dnorm(x,mean=50,sd=sqrt(10)),add=TRUE)		# N(50,10)の正規分布のグラフを重ね合わせて表示

誤差の絶対値が5以下、すなわち母平均の推定値が45以上、55以下の範囲で得られた回数を調べてみましょう:

> 誤差絶対値5以下　<- ifelse(abs(標本平均 - 50) <= 5, 1, 0)
> table(誤差絶対値5以下)
   誤差絶対値5以下
　　　0          1
  1167        8833

> pnorm(55, mean=50, sd=sqrt(10)) - pnorm(45,mean=50,sd=sqrt(10))
[1] 0.8861537

平均μ、分散σ²の正規分布 N (μ, σ²)に従う母集団からサンプルサイズnの標本を抽出したとき、標本平均の標本分布は、理論的にN (μ, σ²/n) に従うことは、次のようにして実験的に確かめることができます。

for ( k in c(50,60,70,80,90,100,150,200,250,300) )  { # サンプルサイズを指定
  標本平均　<- numeric(length=10000)         # 推定値を格納する場所を予約
  for( i  in  1:10000 ){   			# 括弧の中を10000回処理する
　　標本 <- rnorm(n=k,mean=50,sd=10)			# 正規分布に基づく標本を生成(n=k)
    標本平均[i] <-  mean(標本)   			# 標本平均を計算
  }
  cat("サンプルサイズ=",k," 平均=",mean(標本平均), 
      "分散=",var(標本平均)*(k-1)/k, "母分散/分散=",
     100/(var(標本平均)*(k-1)/k), "\n")
}

サンプルサイズ= 50  平均= 49.99211 分散= 1.949846 母分散/分散= 51.2861 
サンプルサイズ= 60  平均= 50.00122 分散= 1.633351 母分散/分散= 61.22382 
サンプルサイズ= 70  平均= 50.02179 分散= 1.416551 母分散/分散= 70.59402 
サンプルサイズ= 80  平均= 49.97832 分散= 1.235205 母分散/分散= 80.95824 
サンプルサイズ= 90  平均= 50.01785 分散= 1.114071 母分散/分散= 89.76091 
サンプルサイズ= 100  平均= 49.99396 分散= 0.9985369 母分散/分散= 100.1465 
サンプルサイズ= 150  平均= 49.9866  分散= 0.667744 母分散/分散= 149.758 
サンプルサイズ= 200  平均= 50.00681 分散= 0.499225 母分散/分散= 200.3105 
サンプルサイズ= 250  平均= 50.00481 分散= 0.3976698 母分散/分散= 251.4649 
サンプルサイズ= 300  平均= 50.00694 分散= 0.3322331 母分散/分散= 300.9935

例: N (50, 10²)の正規母集団から n=10 の標本を抽出したとき、 標本平均の標本分布はN(50, 10²/n)に従うので、標準誤差はsqrt(n)、つまり(n=10なので) およそ3.2となります。 ということは、標本平均から母集団の平均を推定しようとすると、おおよそ68%の確率で46.8から53.2の間の値が得られますが、運が悪いと5以上の誤差が生じることがあることがわかります。

この例から、(1)母集団の分散(母分散)が大きければ大きいほど標本の平均値は母平均から離れた値を取りやすい(標準誤差が大きくなりやすい)、(2)標本のサイズが小さいほど標準誤差が大きい、ということがわかります。 言い換えれば、標準誤差を小さくするにはサンプルサイズnを大きくするのが有効な方法であると言えます。

一般に、どのような標本統計量でも、サンプルサイズが大きいほど標準誤差は小さくなります

中心極限定理: 母集団が正規分布でなくとも、標本平均の標本分布は「サンプルサイズnが大きい時」はほぼ正規分布になる

課題3-2

標本分布を求めるの項目ではn=10の標本を10000個作り、その標本平均を求めた。 この課題ではn=100として10000個の標本平均を求めてみよう。そしてn=10の場合のヒストグラムとn=100の場合のヒストグラムを比較し、標準誤差が小さくなっていることを確かめてみよ。 [ヒント]

n=10のときの標本平均が変数 SampleAverage10に入っており、 n=100のときの標本平均が変数 SampleAverage100に入っているとする。 この2つのヒストグラムを重ねあわせ表示するには、x軸とy軸をともに合わせる必要がある。 それには次のように xlim と ylim の値を指定すればよい:

hist(SampleAverage10, col="yellow", xlim=c(35, 65), ylim=c(0,2500))		# n=10の方を黄色で表示
par(new=T)						# ヒストグラム同士を重ねあわせ表示する
hist(SampleAverage100, col="green", xlim=c(35, 65), ylim=c(0,2500))		# n=100の方を緑色で表示

その結果、以下の様な図が表示される: