1 + 1 |      | #数値演算 |
x <- c(2, 3, 5) |      | #ベクトルの定義 |
y <- c(10, 14, -7) |      |      |
x + y |      | #ベクトルの和 |
2x |      | #ベクトルのスカラー倍 |
sum(x) |      | #ベクトル要素の和 |
1:10 |      | #連続した値 |
x <- 1:10 |      | # |
y <- c(1,3,2.8,2.5,3,3.5,3.8,4,3.5,3.6) |      | # |
length(y) |      | # y の長さ |
plot(x, y) |      | #グラフ表示 |
plot(x, y, type="l") |      | #グラフ表示 |
length(eigo) |      | #データ数 |
mean(eigo) |      | #標本平均 |
var(eigo) |      | #標本分散 |
sd(eigo) |      | #標本標準偏差 |
boxplot(eigo, main="英語得点の箱ヒゲ図") |      | #箱ひげ図 |
boxplot.stats(eigo) |      | #箱ひげ図用統計量 |
hist(eigo, breaks=seq(0, 100, by=5), xlab="英語得点", ylab="頻度", main="") | ||
title(main = "英語得点のヒストグラム") |      | #グラフタイトル |
stem(eigo, scale=2) |      | #幹葉表示 |
boxplot(eigo,kokugo, xaxt="n") axis(1, 1:2, labels=c("英語","国語")) title(main = "英語と国語得点の箱ひげ図")
> stem(eigo,scale=2) The decimal point is 1 digit(s) to the right of the |      2 | 34 2 | 3 | 23 3 | 68 4 | 23 4 | 6677 5 | 0001112223344 5 | 5566777788899 6 | 0011122334444 6 | 5555556667777888899 7 | 011223 7 | 67 8 | 01 |
![]() |
![]() |
偏差を標本標準偏差(SD)で割ったもの. | ![]() |
標準化データ(zi)の平均は 0,標準偏差は 1,(分散も1) |
変数 X    | x1 | x2 | x3 | x4 | x5 |
---|---|---|---|---|---|
確率 P    | p1 | p2 | p3 | p4 | p5 |
変数 X    | 1 | 2 | 3 | 4 | 5 | 6 |
---|---|---|---|---|---|---|
確率 P    | 1/6 | 1/6 | 1/6 | 1/6 | 1/6 | 1/6 |
x <- 1:6 |      | # 1 から 6 まで並べたベクトル |
p <- rep(1/6, 6) |      | # 1/6 を6個並べたベクトル |
m <- sum(x*p) |      | #平均 |
s2 <- sum(p*(x-m)^2) |      | #分散 |
n <- 5 | #試行回数 |
x <- 0:n | #回数 |
p <- 0.3 | #成功確率 |
hit <- dbinom(x, size=n, prob=0.3) | #二項確率 |
plot(x, hit, type="h", ylim=c(0,0.4), xlim=c(0,n), cex.lab=0.8, xlab="ヒット数", ylab="確率密度") | |
title(main="ヒット数の分布(n=5,p=0.3)") | #タイトル |
sum(x*hit) | #平均(np = 1.5) |
sum(hit*(x - n*p)^2) | #分散(np(1 - p) = 1.05) |
   5,6の個数    |    0 |    1 |    2 |    3 |    4 |    5 |    6 |
---|---|---|---|---|---|---|---|
出た回数 |    185    |    1149    |    3265    |    5475    |    6114    |    5194    |    3067    |
   5,6の個数    |    7 |    8 |    9 |    10 |    11 |    12 |    合計 |
出た回数 |    1331    |    403    |    105    |    18    |    0    |    0    | 26306    |
   モデル    |    平均    |    分散    |
---|---|---|
   データ    |    4.052    |    2.696    |
   2項分布(p=0.338)    |    4.052    |    2.685    |
   2項分布(p=0.333)    |    4    |    2.667    |
x <- 0:12 | #個数 |
dice <- c(185,1149,3265,5475,6114,5194,3067,1331,403,105,18,0,0) | #回数データ |
sum(dice) | #試行回数 |
pdice <- dice/sum(dice) | #回数の確率 |
m <- sum(x*pdice) | #平均 |
p <- m/12 | #5,6の出る確率 |
s2 <- sum(pdice*(x-m)^2) | #分散 |
v <- 12*p*(1-p) | #二項分布のもとでの分散 |
h1 <- dbinom(x, 12, 1/3) | #正しいサイコロのもとでの二項確率分布 |
h2 <- dbinom(x, 12, p) | #推定確率からの二項確率分布 |
dicedis <- rbind(pdice,h2,h1) | #行ベクトル−>行列 |
colnames(dicedis) <- as.character(0:12) | #列の名前 |
barplot(dicedis, beside=TRUE, cex.axis=0.8, cex.lab=1.0, xlab="5,6の個数", ylab="確率", legend=c("データ","p=0.338", "p=0.333")) | |
title(main="Weldon のサイコロ実験の分布") | #グラフタイトル |
   女児数    |    0 |    1 |    2 |    3 |    4 |    5 |    6 |
---|---|---|---|---|---|---|---|
度数 |    7    |    45    |    181    |    478    |    829    |    1112    |    1343    |
   女児数    |    7 |    8 |    9 |    10 |    11 |    12 |    合計 |
度数 |    1033    |    670    |    286    |    104    |    24    |    3    | 6155    |
x <- 0:8 | #グラフのx軸の範囲 |
lam <- 2 | #λの定義 |
yp <- dpois(x,lam) | #ポアソン分布の確率密度 |
y1 <- dbinom(x, 5, 0.4) | #二項分布(n = 5,p = 0.4)の確率密度 |
y2 <- dbinom(x, 10, 0.2) | #二項分布(n = 10,p = 0.2)の確率密度 |
y3 <- dbinom(x, 20, 1/10) | #二項分布(n = 20,p = 0.1)の確率密度 |
y4 <- dbinom(x, 40, 1/20) | #二項分布(n = 40,p = 0.05)の確率密度 |
plot(x, y1, type="b", ylab="確率") | #二項分布(n = 5,p = 0.4)のプロット(黒) |
points(x, y2, type="b", col="green") | #二項分布(n = 10,p = 0.2)のプロット(赤) |
points(x, y3, type="b", col="blue") | #二項分布(n = 20,p = 0.1)のプロット(青) |
points(x, y4, type="b", col="purple") | #二項分布(n = 40,p = 0.05)のプロット(紫) |
points(x, yp, type="b", col="red") | #ポアソン分布のプロット(赤) |
title(main="二項分布(np = 2)がポアソン分布に近づく様子") | |
# 凡例の記述(locator(1)は,凡例の記述場所をクリックで指定) | |
legend(locator(1), c("p=0.4", "p=0.2", "p=0.1", "p=0.05", "ポアソン"), | |
lty=1, col=c("black", "green", "blue", "purple", "red")) |
   死亡記事件数    |    0 |    1 |    2 |    3 |    4 |    5 |    6 以上 |
---|---|---|---|---|---|---|---|
日数 |    484    |    391    |    164    |    45    |    11    |    1    | 0    |
死亡記事件数データの平均は0.8239,分散は0.8294,であった.このデータがポアソン分布に従っている
と考える.ポアソン分布の平均は λ なので,λ = 0.8239 のポアソン分布にあてはめてみたところ,非常によく
一致していた.
また,ポアソン分布は,平均と分散が等しいという特徴がある.データの平均と分散の値が
近いことから,データはポアソン分布によく適合していることを示している.
![]() |
|
x <- 0:6 | #グラフのx軸の範囲 |
y <- c(484, 391, 164, 45, 11, 1, 0) | #死亡記事件数データ |
s <- sum(y) | #データ総数 |
m <- sum(x*y/s) | #データ分布の平均 |
v <- sum((x-m)^2*y/s) | #データ分布の分散 |
yp <- dpois(x, m) | #平均 m のポアソン分布確率密度 |
plot(x, y/s, type="h", ylab="確率") | #データの棒グラフ表示 |
points(x, yp, type="b", col="red") | #ポアソン分布の重ねがき(赤) |
title(main="死亡記事件数へのポアソン分布のあてはめ") | |
legend(3.5, 0.4, c("データ", "ポアソン分布"), lty=1, col=c("black","red")) |
7 日間で売れる商品の個数は,平均 2×7 = 14 のポアソン分布に従うと考えられる.ポアソン分布 の95%(分位)点を求めればよい.これより,必要な在庫は20個とわかる.
m <- 14 | # ポアソン分布のパラメータ(平均) |
x <- 0:30 | # グラフのx軸の範囲 |
yp <- dpois(x, m) | # ポアソン分布確率密度 |
cyp <- ppois(x,m) | # ポアソン分布累積確率 |
stok <- qpois(0.95, m) | # 95%(分位)点 |
stok | # 答えの表示 |
op <- par(mfrow = c(1, 2)) | # 横に2つのグラフを並べる |
plot(x, yp, type="h", ylab="確率密度") | # 確率密度グラフ |
points(x[2:21],yp[2:21], type="h", col="red") | # x = 20 まで赤色表示 |
plot(x,cyp,type="l", ylab="累積確率") | # 累積確率グラフ |
arrows(stok,0.95,stok,0, length=0.1, col="red") | # 赤矢印 |
segments(0,0.95, stok,0.95, col="red") | # |
par(op) | # グラフ表示もとに戻す |
title(main="ポアソン分布(λ = 14)の95%点") | # グラフタイトル |
平均: E[X ] = (a + b)/2,分散: Var[X ] = (b - a)2/12
x <- runif(10000) | #(0,1)一様乱数1000個列 |
hist(x, main="(0, 1) 一様乱数 10000個") | #ヒストグラム表示 |
mean(x) | # 一様乱数の平均(理論値は 1/2) |
var(x) | # 一様乱数の分散 (理論値は 1/12) |
n <- 10000 | #一様乱数の個数 |
x <- runif(n, -1, 1) | #(-1, 1) の範囲の一様乱数 n 個生成 |
y <- runif(n, -1, 1) | # |
r <- x^2 + y^2 | #原点からの距離の2乗 |
plot(x,y, type="n", xlim=c(-1,1), ylim=c(-1,1)) | #グラフの表示範囲の指定 |
abline(h=0) | # x 軸の表示 |
abline(v=0) | # y 軸の表示 |
segments(-1, 1, 1, 1) | #(-1, 1)から(1, 1)までの直線 |
segments(1, 1, 1, -1) | # |
segments(-1, -1, 1, -1) | # |
segments(-1, -1, -1, 1) | # |
pin <- (1:n)[r<1] | #乱数のうち単位円内に入る乱数の番号 |
points(x[-pin], y[-pin], pch=".", col="green") | #単位円の外の乱数を緑点で表示 |
points(x[pin], y[pin], pch=".", col="red") | #単位円内の乱数を赤い点で表示 |
s <- 0:360 | # 0 度から 360 度 |
theta <- s*pi/180 | #度をラジアンに変換 |
xp = sin(theta) | #単位円の x 座標 |
yp = cos(theta) | #単位円の y 座標 |
points(xp,yp, type="l") | #単位円を表示 |
title(main="(-1,1)一様乱数による点列と単位円") | |
length(pin) | #単位円内に入った乱数の個数 |
4*length(pin)/n | #πの近似値 |
![]() |
![]() |
# 平均の異なる正規分布 | |
curve(dnorm(x, 40, 4), 30, 70, ylim=c(0,0.2), xlab="",ylab="確率密度") | # 平均:40,標準偏差:4 |
curve(dnorm(x, 50, 4), add=TRUE, col="red") | # 平均:50,標準偏差:4 |
curve(dnorm(x, 60, 4), add=TRUE, col="blue") | # 平均:60,標準偏差:4 |
title(main="正規分布(異なる平均(μ))\n平均=40, 50, 60,σ = 4") | # タイトル |
legend(52,0.19,c("μ = 40","μ = 50","μ = 60"), lty=1, col=c("black","red","blue")) | # 凡例 |
# 標準偏差(分散)の異なる正規分布 | |
curve(dnorm(x, 50, sd=2), 30, 70, ylim=c(0,0.2), xlab="",ylab="確率密度") | # 平均:50,標準偏差:2 |
curve(dnorm(x, mean=50, sd=4), add=TRUE, col="red") | # 平均:50,標準偏差:4 |
curve(dnorm(x, mean=50, sd=6), add=TRUE, col="green") | # 平均:50,標準偏差:6 |
title(main="正規分布(異なる標準偏差(σ))\n平均=50,σ = 2, 4, 6") | # タイトル |
legend(55,0.19,c("σ = 2","σ = 4","σ = 6"), lty=1, col=c("black","red","green")) | # 凡例 |
![]() |
![]() |
pmorm(-1) | # = 0.16(赤矢印) |
pmorm(1) - pnorm(-1) | # = 0.683 |
pmorm(2) - pnorm(-2) | # = 0.954 |
qnorm(0.975) | # = 1.96(青矢印),両側 5 %点 |
   身長    |    57 |    58 |    59 |    60 |    61 |    62 |    63 |    64 |    65 |    66 |
---|---|---|---|---|---|---|---|---|---|---|
人数 |    2    |    4    |    14    |    41    |    83    |    169    |    394    |    669    |    990    |    1223    |
   67 |    68 |    69 |    70 |    71 |    72 |    73 |    74 |    75 |    76 |    77 |
   1329    |    1230    |    1063    |    646    |    392    |    202    |    79    |    32    |    16    |    5    |    2    |
x <- 57:77 | # 身長(x)の範囲 |
y <- c(2, 4, 14, 41, 83, 169, 394, 669, 990, 1223, 1329, | # 身長ごとのデータ |
1230, 1063, 646, 392, 202, 79, 32, 16, 5, 2) | # |
s <- sum(y) | # データ総数 |
m <- sum(x*y/s) | # データの平均 |
v <- sum(y/s*(x-m)^2) | # データの分散 |
plot(x, y/s, type="h", xlab="身長(インチ)", ylab="頻度") | # データの棒グラフ表示 |
curve(dnorm(x, m, sqrt(v)), 57, 77, add=T, col="red") | # 正規密度のグラフ表示 |
title(main="身長データに対する正規分布のあてはめ") | # タイトル |
m <- 67.02; s <- 2.57 | # 平均と標準偏差の指定 |
1 - pnorm(70, mean=m, sd=s) | # 1. 70 までの累積確率を 1 から引く |
qnorm(0.9, mean=m, sd=s) | # 2. 累積確率が 0.9 となる身長 |
pnorm(70, mean=m, sd=s) - pnorm(65, mean=m, sd=s) | # 3. (70 までの累積確率) − (60 までの累積確率) |
n <- 5 | # 打数 |
x <- 0:n | # xの範囲 |
p <- 0.3 | # 打率 |
hit <- dbinom(x, size=n, prob=0.3) | # 二項確率 |
y <- pbinom(x, size=n, prob=0.3) | # 二項累積確率 |
m <- n*p | # 平均 |
sd <- sqrt(n*p*(1-p)) | # 標準偏差 |
op <- par(mfrow = c(1, 2)) | # |
plot(x, hit, type="h", ylim=c(0,0.4), xlim=c(0,7), xlab="ヒット数", ylab="確率密度") | |
curve(dnorm(x, mean=m, sd=sd), add=TRUE, col="red") | # 確率密度 |
plot(x, y, type="s", ylim=c(0,1), xlim=c(0,7), xlab="ヒット数", ylab="累積確率") | |
curve(pnorm(x, mean=m, sd=sd), add=TRUE, col="red") | # 累積確率 |
par(op) | # |
title(main ="二項分布:n = 5 打数,打率 p = 0.3; 正規分布:N(1.5, 1.05) ") |
課題:打数 n を大きくして,二項分布が正規分布に近づく様子を確かめよ.
![]() |
![]() |
x <- 57:77 | # 身長の範囲 |
y <- c(2, 4, 14, 41, 83, 169, 394, 669, 990, 1223, 1329, 1230, 1063, 646, 392, 202, 79, 32, 16, 5, 2) | |
hei <- NULL | # 身長ごとのカウント |
for(i in 1:length(x)) | # |
hei <- c(hei, rep(x[i], y[i])) | # カウントデータを個人データに変換 |
qqnorm(hei, xlab="正規分布分位点", ylab="データ分位点", main="") | # 正規 Q - Q プロット表示 |
qqline(hei, col="red") | # 正規分布の直線表示 |
title(main="英国人成人男子身長の正規 Q-Q プロット") | # グラフタイトル |
N <- 10000 | # 乱数列の長さ |
hit <- rbinom(N, size=100, prob=0.3) | # 100打席でのヒット数の乱数列生成 |
qqnorm(hit, xlab="正規分布分位点", ylab="データ分位点", main="") | # 正規 Q - Q プロット表示 |
qqline(hit, col="red") | # 正規分布の直線表示 |
title(main="100打席ヒット数分布の正規 Q-Q プロット") | # グラフタイトル |
![]() |
![]() |
![]() |
N <- 10000 | # 乱数列の長さ |
n <- 2 | # 標本平均のサイズ |
u <- matrix(data=runif(n*N), ncol=n) | # N×n の一様乱数行列 |
um <- apply(u, 1, mean) | # 行ごとの平均 |
op <- par(mfrow = c(1, 2)) | # 標本平均のヒストグラム |
hist(um, breaks=seq(0,1,by=0.02), freq=FALSE, ylim=c(0, 2.5), xlab="", ylab="頻度", main="") | |
m <- mean(um) | # 標本平均列の平均 |
s <- sd(um) | # 標本平均列の標準偏差 |
curve(dnorm(x, m, s), 0, 1, add=TRUE, col="red") | # 正規分布の重ねがき |
qqnorm(um, xlab="正規分布分位点", ylab="データ分位点", main="") | # 正規 Q - Q プロット |
qqline(um, col="red") | # 正規分布の直線表示 |
par(op) | # |
title(main="一様乱数 2 個の標本平均分布に正規分布 N(0.5, 0.042) を重ね書き") |
課題:サンプルサイズ n を大きくして,標本平均の分布が正規分布に近づく様子を確かめよ.
curve(dchisq(x, 1), 0, 20) | # 自由度 1 の χ2 分布のグラフの表示 |
abline(v=0, h=0) | # x 軸と y 軸の表示 |
curve(dchisq(x, 4), add=T, col=2) | # 自由度 4 の χ2 分布のグラフを色 2(赤)で追加 |
curve(dchisq(x, 10), add=T, col=3) | # 自由度 10 の χ2 分布のグラフを色 3(緑)で追加 |
legend(10, 0.7, c("n = 1", "n = 4", "n = 10"), lty=1, col=c(1, 2, 3)) | |
title(main="χ2 分布の自由度 n による形状の違い") | # タイトル |
N <- 10000 | # 乱数列の長さ |
n <- 2 | # 自由度 |
u <- matrix(rnorm(n*N), ncol=n) | # N×n の標準正規乱数乱数行列 |
u2 <- u^2 | # 行列の要素の2乗 |
un <- apply(u2, 1, sum) | # 行ごとの和 |
umx <- ceiling(max(un)) | # 最大値を超える整数 |
op <- par(mfrow = c(1, 2)) | # |
hist(un, breaks=seq(0,umx,by=0.5), freq=FALSE, xlim=c(0,15), xlab="標準正規乱数(n = 2)の2乗和", ylab="頻度", main="") | |
curve(dchisq(x, 2), 0, 15, add=T, col=2) | # 自由度 2 の χ2 分布の重ね合わせ |
plot(ecdf(un), do.points=F, verticals=T, xlim=c(0,12), ylab="累積確率", main="") | |
curve(pchisq(x,2), 0, 15, add=T, col=2) | # 自由度 2 の χ2 累積分布関数の重ね合わせ |
par(op) | # |
title(main="標準正規乱数の2乗和に χ2 分布の重ね合わせ") | # タイトル |
課題:2乗和する数 n を大きくした場合も,標準正規乱数の2乗和の分布が 自由度 n の χ2 分布に従うことを確かめよ.
op <- par(mfrow = c(1, 2)) | # |
curve(df(x, 1, 10), 0, 5, ylim=c(0,1.5), ylab="確率密度", xlab="n = 10") | # m = 1,n = 10 の F 分布 |
abline(v=0, h=0) | # x 軸,y 軸 |
curve(df(x, 2, 10), 0, 5, add=T, col=2) | # m = 2,n = 10 の F 分布 |
curve(df(x, 4, 10), 0, 5, add=T, col=3) | # m = 4,n = 10 の F 分布 |
curve(df(x, 8, 10), 0, 5, add=T, col=4) | # m = 8,n = 10 の F 分布 |
legend(2.5, 1.4, c("m = 1", "m = 2", "m = 4", "m = 8"), lty=1, col=1:4) | # 凡例 |
curve(df(x, 4, 50), 0, 5, col=2, ylab="確率密度", xlab="m = 4") | # m = 4,n = 50 の F 分布 |
abline(v=0, h=0) | # x 軸,y 軸 |
curve(df(x, 4, 10), 0, 5, add=T) | # m = 4,n = 10 の F 分布 |
legend(2.5, 0.7, c("n =10", "n = 50"), lty=1, col=c("black","red")) | # 凡例 |
par(op) | # |
title(main="F 分布の分子と分母の自由度の違いによる形状") | # タイトル |
N <- 10000 | # 乱数列の長さ |
m <- 4 | # 分子自由度 |
n <- 10 | # 分母自由度 |
um0 <- matrix(rnorm(m*N), ncol=m) | # N×m の標準正規乱数乱数行列 |
um2 <- um0^2 | # 行列の要素の2乗 |
um <- apply(um2, 1, sum) | # 行ごとの和 |
um <- um/m | # 自由度で割る |
un0 <- matrix(rnorm(n*N), ncol=n) | # |
un2 <- un0^2 | # |
un <- apply(un2, 1, sum) | # |
un <- un/n | # |
fv <- um/un | # χ2 分布乱数の比 |
fmx <- ceiling(max(fv)) | # fv の最大値を超える整数 |
hist(fv, breaks=seq(0,fmx,by=0.2), freq=FALSE, xlim=c(0,6), main="") | |
curve(df(x, 4, 10), 0, 6, add=T, col=2) | # 自由度 4,10 の F 分布の重ね合わせ |
title(main="独立な χ2 分布乱数の比に自由度 4,10 の F 分布のあてはめ", cex.main=0.9) |
課題:分子,分母の自由度を変えて,χ2 分布乱数の比が F 分布に従うことを確かめよ.
curve(dnorm(x), -6, 6) | # 標準正規分布表示 |
abline(v=0, h=0) | # x 軸と y 軸の表示 |
title(main="t 分布の自由度による形状の違い") | # タイトル |
課題: t 分布の確率密度関数は,自由度を n として dt(x, n) で与えられる.n の値を変えることにより, 上の図のようなグラフを描け.
y <- rt(1000, df=10) | # 自由度 10 の t 分布乱数1000個生成 |
qqnorm(y, xlab="正規分布分位点", ylab="t 分布分位点", main="") | # 正規 Q - Q プロット |
qqline(y, col=2) | # 正規分布の四分位範囲直線表示 |
title(main="t(10) の正規 Q - Q プロット") | # |
課題: 自由度を n の値を変えて,標準正規分布とのずれの様子を正規 Q - Q プロット で確かめよ.
N <- 10000 | # シミュレーション回数 |
n <- 10 | # 1回のサンプルサイズ |
un0 <- matrix(rnorm(n*N), ncol=n) | # N×n の標準正規乱数行列 |
un2 <- un0^2 | # 標準正規乱数行列の要素の2乗 |
un <- apply(un2, 1, sum) | # 要素の2乗の各行の和(自由度 n の χ2 分布乱数 N 個) |
unr <- sqrt(un/n) | # 自由度 n の χ2 分布乱数を n で割った平方根 |
z <- rnorm(N) | # 標準正規乱数 N 個 |
tv <- z/unr | # t 値 N 個 |
tmx <- ceiling(max(abs(tv))) | # t 値 の絶対値の最大 |
hist(tv, breaks=seq(-tmx,tmx,by=0.2), freq=FALSE, xlim=c(-5,5), main="") | |
curve(dt(x, 10), add=T, col=2) | # 自由度 10 の t 分布の重ねがき |
title(main="標準正規乱数から生成した t 値に t 分布の重ね合わせ") |
N <- 10000 | # 標本抽出の回数 |
n <- 20 | # 1回のサンプルサイズ |
m <- 50; s <- 10 | # 平均,標準偏差 |
r <- rnorm(N, m, s) | # 母集団からの標本 |
mean(r) | # 標本平均 |
var(r) | # 標本分散 |
sd(r) | # 標本標準偏差 |
rn <- matrix(rnorm(N*n, m, s), ncol=n) | # 大きさ n の標本抽出 N 回 |
rm <- apply(rn, 1, mean) | # 標本平均 |
mean(rm) | # 標本平均の平均 |
var(rm) | # 標本平均の分散 |
sd(rm) | # 標本平均の標準偏差 |
op <- par(mfrow = c(1, 2)) | # |
hist(r, breaks=seq(0,100, by=1), xlim=c(15,85),ylim=c(0,2500),xlab="",ylab="頻度",main="") | |
title(main="正規分布 N(50,100) のヒストグラム", cex.main=0.9) | |
hist(rm, breaks=seq(30,70, by=1), xlim=c(15,85),ylim=c(0,2500), xlab="",ylab="頻度",main="") | |
title(main="大きさ 20 の標本平均のヒストグラム", cex.main=0.9) | |
par(op) | # |
N <- 10000 | # 標本抽出の回数 |
m1 <- 172.5; s1 <- 6 | # 集団 A の平均と標準偏差 |
m2 <- 168; s2 <- 4.5 | # 集団 B の平均と標準偏差 |
a <- rnorm(N, m1, s1) | # 集団 A からの標本 |
b <- rnorm(N, m2, s2) | # 集団 B からの標本 |
op <- par(mfrow = c(1, 2)) | # |
hist(b, breaks=seq(140,200, by=2), ylim=c(0, 2000), col="gray",ylab="頻度", xlab="身長(cm)", main="") | |
title(main="2つの集団の重ね合わせ") | # |
par(new=T) | # グラフの重ね合わせ |
hist(a, breaks=seq(140,200, by=2), ylim=c(0,2000),density=0.1, ylab="", xlab="",main="", col="red") | |
d <- a - b | # 2つの母集団からの標本の差 |
hist(d, breaks=seq(-30,40,by=2), ylab="頻度", xlab="身長(cm)",main="") | # |
title(main="2つの標本の差の分布") | # |
par(op) | # |
length(d[d>0])/length(d) | # シミュレーションによる確率 |
dm <- m1-m2; dv <- s1^2 + s2^2 | # 差の分布の平均と分散 |
1 - pnorm(0, mean=dm, sd=sqrt(dv)) | # 差が0以上の確率 |
課題:上記の確率を求めよ.
v <- 2; n <- 10 | # 母分散と標本の大きさ |
d <- qnorm(0.975)*sqrt(v/n) | # 95%信頼区間の幅 |
x <- c(-2.5, -2.5, 2.5, 2.5) | # グラフの x の範囲 |
y <- c(0, 100, 100, 0) | # グラフの y の範囲 |
plot(x, y, type="n", xlab="", ylab="") | # グラフ領域確保 |
segments(0, 0, 0, 100, col="red") | # 母平均のライン |
for(i in 0:100){ | # 100回の繰り返し |
m <- mean(rnorm(n, mean=0, sd=sqrt(v))) | # N(0, v)からの大きさ n の標本平均 |
segments(m-d, i, m+d, i) | # 信頼区間の表示 |
points(m, i, pch=4, col="red", cex=0.8) | # 標本平均の赤×表示 |
if(m-d>0 || m+d<0) text(-2.5, i, "*") | # 信頼区間が母平均を含まない(失敗)した場合 |
} | # |
title(main="N( 0, 2 )からの大きさ 10 の標本から得られた | |
平均 μ の 95 %信頼区間を 100 回作成", cex.main=1.0) | # |
m <- 20; v <- 4 | # 正規母集団の平均と分散 |
n <- 20 | # 標本の大きさ |
x <- c(0, 0, 20, 20) | # グラフの x の範囲 |
y <- c(0, 100, 100, 0) | # グラフの y の範囲 |
plot(x, y, type="n", xlab="", ylab="") | # グラフ領域確保 |
segments(v, 0, v, 100, col="red") | # 母分散のライン |
for(i in 0:100){ | # 100回の繰り返し |
r <- rnorm(n, m, sqrt(v)) | # 母集団からの無作為標本 |
s <- (n-1)*var(r) | # 標本の偏差平方和 |
s1 <- s/qchisq(0.975, df=n-1) | # 95 %信頼区間の下限 |
s2 <- s/qchisq(0.025, df=n-1) | # 95 %信頼区間の上限 |
segments(s1, i, s2, i) | # 信頼区間の表示 |
points(var(r), i, pch=4, col="red", cex=0.8) | # 標本分散の赤×表示 |
if(s1 > v || s2 < v) text(0, i, "*") | # 信頼区間が母分散を含まない(失敗)した場合 |
} | # |
title(main="N( 20, 4 )からの大きさ20の標本から得られた | |
分散の95%信頼区間を100回作成", cex.main=1.0) | # |
|
![]() |
統計学で扱う仮説とは,母集団に対する断定や推測.たとえば,
統計的仮説検定で用いられる仮説は,まず,帰無仮説という形式で与えられる.
帰無仮説は棄却されることに意味がある仮説である.
帰無仮説と反対の仮説を対立仮説という.
上の3番目の例でみると,
母集団 A と母集団 B は異なる処理(薬の投与など)をしているので,実験の目的
は,母集団 A と母集団 B の平均は異なる(処理効果がある)ことを言いたい
(対立仮説が正しいことを望む)のだが,まずは
「等しい(処理効果無し)」と仮定してみようという考え方.
統計的仮説検定では,たとえば2つの母集団平均が等しいという帰無仮説を考えると,
この帰無仮説のもとで,検定統計量(標本平均の差に基づく t 値など)以上
(もしくは未満)の値が得られる確率を求める.
くだけた言い方をすれば,帰無仮説が正しいとしたときに,標本のようなデータが得られる確率
を求める.
これが十分小さい(ほとんどありえない)ときは,平均が等しいと仮定したことが誤りであったと判断して
帰無仮説を棄却し,2つの母集団平均には差があると結論づける.
この確率がそれほど小さく
ない場合は,このような統計量が得られることもありえると考え,帰無仮説を採択し,平均が等しいと考え
てもよいとする.
棄却か採択かの判断の基準となる確率を有意水準といい,
5 % や 1 % がよく用いられる.
実験状況によっては,薬投与などの処理を行った集団(処理群)平均 μA が,薬を投与しない 集団(対照群)の平均 μB より小さくなることはないことが事前に わかっているような場合が ある.このようなとき,
母集団平均に対する両側検定は,母集団平均に対する信頼区間と大きな関係がある. いま,帰無仮説(H0)と対立仮説(H1)が,
検定は,仮説を棄却するか採択するかのいずれかであるが,
統計量は分布をもつので,この判定には間違いが起こることがある.
以下のように,この過誤には
2 種類がある.
  | 仮説の棄却 | 仮説の採択 |
---|---|---|
仮説が真のとき | 第1種の過誤 | 正解 |
仮説が偽のとき | 正解 | 第2種の過誤 |
第1種の過誤が有意水準である.また,第2種の過誤の確率を β としたとき, 仮説が偽のとき正しく仮説を棄却する確率,1 - β,を検出力という. よい検定は,第1種の過誤を固定したもとで検出力の高い検定方式である.
なお,この検定は,対のある標本に適用できる.対のある標本とは,n 組のペアー標本,
文献
engei <- read.csv("engei.csv") | # csv データ読み込み |
engei | # データの表示 |
x <- engei[,3] | # 園芸療法前の MSE 得点を x に格納(列数で指定) |
x <- engei$療法前 | # 園芸療法前の MSE 得点を x に格納(変数名で指定) |
y <- engei$療法後 | # 園芸療法後の MSE 得点を x に格納(変数名で指定) |
t.test(y, x, paired=TRUE) | # 1母集団 t 検定,paired = TRUE で対標本を指定 |
d <- y - x | # 療法後 − 療法前で療法の効果をみる |
d | # 療法の効果の表示 |
t.test(d) | # 1母集団 t 検定,先ほどと同じ検定 |
n <- length(d) | # 標本の大きさ(サンプルサイズ) |
mean(d) | # 標本平均 |
sd(d) | # 標本標準偏差 |
dv <- n - 1 | # 標本の自由度 |
t <- sqrt(n)*mean(d)/sd(d) | # 効果がないとの帰無仮説のもとでの t 値 |
t | # 検定統計量 t 値の表示 |
2*(1 - pt(t, df=dv)) | # 両側検定の p 値 |
t0 <- qt(0.975, df=dv) | # 両側 5 %検定の閾値 |
dw <- t0*sd(d)/sqrt(n) | # 95%信頼区間の幅 |
mean(d)-dw | # 95%信頼区間の下限 |
mean(d)+dw | # 95%信頼区間の上限 |
t.test(d, alternative="greater") | # 片側検定 |
このデータでは,クライアント全体では園芸療法の効果は認められなかった.また,MSE の満点が30点なので, 元々心理機能に問題がなく満点近い得点のクライアントでは,効果が認められないのも当然といえる. その上,園芸療法処方の前後で半年 程度のタイムラグがあるので,何もしなくても認知症の症状が進行する場合もあり,効果が検出しずらい こともある程度予測できる.
課題: クライアントのうち若年層( 85 歳以下)を取り出し,園芸療法の効果を検定せよ.
N <- 10000 | # シミュレーション回数 |
n <- 10 | # サンプルサイズ |
n1 <- rnorm(N*n, mean=10, sd=2) | # N(10, 4) から大きさ n の標本を N 回シミュレーション |
n1.mat <- matrix(data=n1, ncol=n) | # データ行列 |
n1.mean <- apply(n1.mat, 1, mean) | # 各サンプルの平均 |
n1.var <- apply(n1.mat, 1, var) | # 各サンプルの分散 |
n1.var3 <- n1.var/n | # 各標本平均の分散 |
n1.td <- (n1.mean - 10)/sqrt(n1.var3) | # 各サンプルの t 値 |
mean(n1.td) | # t 値の平均 |
sd(n1.td) | # t 値の標準偏差 |
mt <- ceiling(max(abs(n1.td))) | # t 値の絶対値の最大値 |
xq1 <- qt(0.025, df = (n-1)) | # t(9) の 2.5% 点 |
xq2 <- qt(0.975, df = (n-1)) | # t(9) の 97.5% 点 |
length( n1.td[abs(n1.td) > xq2]) | # 5% 検定で有意となった個数 |
hist(n1.td, breaks=seq(-mt, mt, by=0.2), probability=TRUE, xlab="t 値", ylab="密度", main="") | |
abline(v=xq1, col="red") | # 採択域の下限 |
abline(v=xq2, col="red") | # 採択域の上限 |
#curve(dt(x, df=(n-1)), -6, 6, add=T, col="red") | # t(9) の表示 |
title(main="t 値のヒストグラムと有意となった個数") | # タイトル |
m <- 11 | # 母集団平均(帰無仮説が偽の場合) |
n1 <- rnorm(N*n, mean=m, sd=2) | # N(m, 4) から大きさ n の標本を N 回シミュレーション |
課題:母集団平均を m = 12 としたときの検出力を求めよ.
2つの母集団 A,B があり,それぞれが平均を μA,μB, 分散を σA2,σB2 の正規分布に従って いるが,その値は未知であるとする.いま,両集団の分散の値が等しく, σA2=σB2=σ2,と仮定 できるとしよう.このとき,
母集団 A から大きさ nA,母集団 B から大きさ nB の標本を抽出した. 母集団 A からの標本の標本平均が x-A, 標本分散が sA2 であり,母集団 B の 標本平均が x-B, 標本分散が sB2 であった.母集団 A,B が共通の 分散 σ2 をもつとすると,その推定値 s2 は 以下のように推定できる.
母集団 A からの標本の偏差平方和: | SA=(nA−1)sA2 |
母集団 B からの標本の偏差平方和: | SB=(nB−1)sB2 |
母集団 A,B 全体での偏差平方和: | S = SA + SB =(nA−1)sA2+ (nB−1)sB2 |
母集団 A,B 共通の標本分散: |
![]() |
帰無仮説(H0: μA = μB)のもとでは,μA−μB=0,なので, 標本平均の差は,
母集団A,Bで分散の同等性が疑われるときは,ウェルチ(Welch)の検定を用いる.
engei <- read.csv("engei.csv") | # csv データ読み込み |
x <- engei$療法前 | # 園芸療法前の MSE 得点を x に格納(変数名で指定) |
y <- engei$療法後 | # 園芸療法後の MSE 得点を y に格納(変数名で指定) |
boxplot(x, y, names= c("対照区","処理区"), ylab="MSE 得点", cex.axis=0.8) | |
title(main="園芸療法処理区と対照区の MSE 得点分布の箱ひげ図", cex.main=1.0) | |
t.test(y, x, var.equal=T) | # 2標本 t 検定 |
mx <- mean(x) | # 対照区平均 |
my <- mean(y) | # 処理区平均 |
nx <- length(x) | # 対照区サンプルサイズ |
ny <- length(y) | # 処理区サンプルサイズ |
dfv <- nx+ny-2 | # データの自由度 |
vx <- var(x) | # 対照区標本分散 |
vy <- var(y) | # 処理区標本分散 |
v <- ((nx-1)*vx + (ny-1)*vy)/dfv | # こみにした共通の分散 |
d <- my - mx | # 処理区平均 − 対照区平均 |
vd <- v/nx + v/ny | # 処理平均の差の分散 |
tv <- d/sqrt(vd) | # t 値 |
2*(1 - pt(abs(tv), df=dfv)) | # p 値 |
t0 <- qt(0.975, df=dfv) | # 自由度 dfv の t 分布の97.5%点 |
dw <- t0*sqrt(vd) | # 95 %信頼区間の幅 |
d - dw | # 平均の差の 95%信頼区間の下限 |
d + dw | # 平均の差の 95%信頼区間の上限 |
t.test(x, y) | # ウェルチ(Welch)の検定 |
課題: クライアントのうち若年層( 85 歳以下)を取り出し,2標本 t 検定を行え.
![]() |
![]() |
![]() |
![]() |
![]() |
1   | 2   | 3   | 4   | R+    | R-    | T    |
---|---|---|---|---|---|---|
+ | + | + | + | 10 | 0 | 0 |
− | + | + | + | 9 | 1 | 1 |
+ | − | + | + | 8 | 2 | 2 |
+ | + | − | + | 7 | 3 | 3 |
+ | + | + | − | 6 | 4 | 4 |
− | − | + | + | 7 | 3 | 3 |
− | + | − | + | 6 | 4 | 4 |
− | + | + | − | 5 | 5 | 5 |
+ | − | − | + | 5 | 5 | 5 |
+ | − | + | − | 4 | 6 | 4 |
+ | + | − | − | 3 | 7 | 3 |
− | − | − | + | 4 | 6 | 4 |
− | − | + | − | 3 | 7 | 3 |
− | + | − | − | 2 | 8 | 2 |
+ | − | − | − | 1 | 9 | 1 |
− | − | − | − | 0 | 10 | 0 |
![]() |
   確率:γ    |    0.99    |    0.95    |    0.90    |    0.85    |    0.80    |
---|---|---|---|---|---|
   臨界値:kγ    |    1.63    |    1.36    |    1.22    |    1.14    |    1.07    |
有意水準 5 %の検定は,自由度 n - 1 の χ2 分布の 2.5%点と 97.5%点をそれぞれ χ2(n - 1)0.025, χ2(n - 1)0.975 とすると,
母集団 A から大きさ nA,母集団 B から大きさ nB の標本を抽出した. 母集団 A からの標本の標本平均が x-A, 標本分散が sA2 であり,母集団 B の 標本平均が x-B, 標本分散が sB2 であるとする.すると, 標本分散に関係した量はそれぞれ
ところで,帰無仮説が正しいとする と,σA2 = σB2 とおけるので, 母集団の分散比は,γ0 = σA2/σB2 = 1, となる.このとき,標本分散の分散比の統計量 γ が,
engei <- read.csv("engei.csv") | # csv データ読み込み |
x <- engei$療法前 | # 園芸療法前の MSE 得点を x に格納(変数名で指定) |
y <- engei$療法後 | # 園芸療法後の MSE 得点を y に格納(変数名で指定) |
var(x) | # 対照区の標本分散 |
var(y) | # 処理区の標本分散 |
var.test(x, y) | # 分散の等質性の検定 |
課題: 分散の等質性検定で表示される数値を R の基本統計関数を用いて出せ.
N <- 10000 | # シミュレーション回数 |
n <- 20 | # サンプルサイズ |
s1 <- 1; s2 <- 2 | # 母集団分散(分散比 = 2) |
n1 <- rnorm(N*n, mean=0, sd=sqrt(s1)) | # 正規乱数 |
n2 <- rnorm(N*n, mean=0, sd=sqrt(s2)) | # 正規乱数 |
n1.mat <- matrix(data=n1, ncol=n) | # N×n データ行列(母集団1) |
n1.var <- apply(n1.mat, 1, var) | # 各行の分散 |
n2.mat <- matrix(data=n2, ncol=n) | # N×n データ行列(母集団2) |
n2.var <- apply(n2.mat, 1, var) | # 各行の分散 |
ratio <- n2.var/n1.var | # 分散比 |
f0 <- qf(0.025, df1=(n-1), df2=(n-1)) | # F 分布2.5%点 |
f1 <- qf(0.975, df1=(n-1), df2=(n-1)) | # F 分布97.5%点 |
m0 <- length(ratio[ratio < f0]) | # 2.5%以下の個数 |
m1 <- length(ratio[ratio > f1]) | # 97.5%以下の個数 |
m0+m1 | # 帰無仮説を棄却した個数 |
mf <- ceiling(max(ratio)) | # 分散比の最大値 |
hist(ratio, breaks=seq(0, mf, by=0.2), probability=TRUE, xlab="分散比", ylab="密度", main="") | |
abline(v=f0, col="red") | # F 分布2.5%点の表示 |
abline(v=f1, col="red") | # F 分布97.5%点の表示 |
title(main="分散比のヒストグラムと有意となった個数") | # |
このように,中心極限定理を利用して,標準正規近似を行って検定を行うやり方を大標本(large sample)理論 といい,コンピュータが発達する以前はもっぱら大標本理論に基づいた検定を行っていたが,現在では正確な確率 が短時間で計算できるので,大標本理論の重要性は大きく低下したといえる.
ここで,さらに近似を加えて,z02 の項を消去すると, p の2次方程式の根は,
ところで,正規近似による信頼区間の構成では,場合により信頼区間が負になったり 1 を超えることがあるが, このときは,0 と 1 で切り詰める.
r <- 7; n <- 10 | # 成功回数と試行回数 |
prop.test(r, n, p=0.5, correct=F) | # H0:p = 0.5 の正規近似検定(補正なし) |
prop.test(r, n, p=0.5) | # H0:p = 0.5 の正規近似検定(連続性の補正) |
binom.test(r, n, p=0.5) | # 二項確率に基づく正確な検定 |
binom.test(r, n, p=0.4) | # p = 0.4 の検定 |
課題: 試行回数と成功回数をともに増やして(n = 30),同様な検定を行え.
課題:ビデオリサーチ社によると,NHK 大河ドラマの関東地区世帯視聴率は26.2%であった. 真の世帯視聴率の 95 %信頼区間を求めよ.
成 功 | 失 敗 | |
---|---|---|
観測度数 | X | n - X |
期待度数 | np0 | n(1 - p0) |
となる.ここで,ピアソン(Pearson)のχ2 値,
r <- 7; n <- 10 | # 成功回数と試行回数 |
x <- c(r, n-r) | # 成功回数と失敗回数のベクトル |
p0 <- c(0.5, 0.5) | # 成功確率と失敗確率の帰無仮説 |
chisq.test(x, p=p0) | # ピアソン χ2 適合度検定 |
セル1 | セル2 | セル3 | セル4 | セル5 |   計   | |
---|---|---|---|---|---|---|
観測度数 | n1 | n2 | n3 | n4 | n5 | n |
想定確率分布 | p1 | p2 | p3 | p4 | p5 | 1 |
成功回数 |   0 回   |   1 回   |   2 回   |   3 回   |   4 回   |   5 回   |   6 回   |   計   |
---|---|---|---|---|---|---|---|---|
期待観測度数 | 8.8 | 26.3 | 32.9 | 22.0 | 8.2 | 1.7 | 0.1 | 100 |
N <- 10000 | # シミュレーション回数 |
m <- 6 | # 試行1回でのベルヌイ試行の数 |
n <- 100 | # 二項確率の試行数 |
p0 <- dbinom(x=0:6, size=6, p=1/3) | # 二項確率 B(6, 1/3) |
w <- rbinom(n, m, p=1/3) | # n = 100 回の試行での成功回数(観測値)の系列) |
tw <- table(factor(w,levels=0:m)) | # 成功回数の表(データ) |
r <- chisq.test(tw, p=p0) | # データと二項確率との適合度検定 |
p1 <- c(p0[1:4], p0[5]+p0[6]+p0[7]) | # 少ない確率をまとめた二項確率 |
y <- NULL | # χ2 値のベクトル |
for(i in 1:N){ | # |
w <- rbinom(n, m, p=1/3) | # n = 100 回の試行での成功回数(観測値)の系列) |
fw <- table(factor(w,levels=0:m)) | # 成功回数の表(データ) |
tw <- c(fw[1:4], fw[5]+fw[6]+fw[7]) | # 成功回数 5,6,7 をまとめる |
r <- chisq.test(tw, p=p1) | # 適合度検定の χ2 値 |
y <- c(y, r$statistic) | # |
} | # |
my <- ceiling(max(y)) | # y の最大値を超える最小の整数 |
mean(y) | # 平均(真値 = 4) |
var(y) | # 分散(真値 = 8) |
hist(y, breaks=seq(0,my,by=0.5),freq=FALSE,xlim=c(0,20),xlab="χ2 値",ylab="頻度",main="") | |
x <- seq(0,20, by=0.1) | # |
curve(dchisq(x, 4), 0, 20, add=T, col=2) | # 自由度 4 χ2 分布 |
title(main="二項分布の適合度と自由度 4 カイ自乗分布") | # |
q0 <- qchisq(0.95,df=4) | # χ2(4) の 95%点 |
segments(q0,0, q0,0.1, col="red") | # 棄却域 |
pv <- length(y[y>q0])/N | # |
s <- paste("有意水準 = ", pv) | # 帰無仮説の棄却率(真値 = 0.05) |
text(15,0.06, s) | # |
セル1 | セル2 | セル3 | セル4 | セル5 |   計   | |
---|---|---|---|---|---|---|
観測度数 | n1 | n2 | n3 | n4 | n5 | n |
推定確率分布 | p^1 | p^2 | p^3 | p^4 | p^5 | 1 |
成功回数 |   0 回   |   1 回   |   2 回   |   3 回   |   4 回   |   5 回   |   6 回   |   計   |
---|---|---|---|---|---|---|---|---|
観測度数 | 12 | 18 | 25 | 32 | 13 | 0 | 0 | 100 |
N <- 10000 | # シミュレーション回数 |
m <- 6 | # 試行1回でのベルヌイ試行の数 |
n <- 100 | # 二項確率の試行数 |
y <- NULL | # χ2 値のベクトル |
for(i in 1:N){ | # |
w <- rbinom(n, m, p=1/3) | # n = 100 回の試行での成功回数(観測値)の系列) |
pd <- mean(w)/m | # データのもとでの推定成功確率 |
p0 <- dbinom(x=0:6, size=6, p=pd) | # 推定成功確率のもとでの二項確率分布 |
p1 <- c(p0[1:4], p0[5]+p0[6]+p0[7]) | # 少ない確率をまとめた二項確率 |
fw <- table(factor(w,levels=0:m)) | # 成功回数の表(データ) |
tw <- c(fw[1:4], fw[5]+fw[6]+fw[7]) | # 成功回数 5,6,7 をまとめる |
r <- chisq.test(tw, p=p1) | # 適合度検定の χ2 値 |
y <- c(y, r$statistic) | # |
} | # |
my <- ceiling(max(y)) | # y の最大値を超える最小の整数 |
mean(y) | # 平均(真値 = 4) |
var(y) | # 分散(真値 = 8) |
hist(y, breaks=seq(0,my,by=0.5),freq=FALSE,xlim=c(0,20),xlab="χ2 値",ylab="頻度",main="") | |
x <- seq(0,20, by=0.1) | # |
curve(dchisq(x, 3), 0, 20, add=T, col=2) | # 自由度 3 χ2 分布 |
title(main="二項分布の適合度と自由度 3 カイ自乗分布") | # |
q0 <- qchisq(0.95,df=3) | # χ2(3) の 95%点 |
segments(q0,0, q0,0.1, col="red") | # 棄却域 |
pv <- length(y[y>q0])/N | # |
s <- paste("有意水準 = ", pv) | # 帰無仮説の棄却率(真値 = 0.05) |
text(15,0.06, s) | # |
   5,6の個数    |    0 |    1 |    2 |    3 |    4 |    5 |
---|---|---|---|---|---|---|
出た回数 |    185    |    1149    |    3265    |    5475    |    6114    |    5194    |
   5,6の個数    |    6 |    7 |    8 |    9 | 10 以上 |    合計 |
出た回数 |    3067    |    1331    |    403    |    105    |    18    | 26306    |
x <- 0:10 | # 5,6の個数 |
y <- c(185,1149,3265,5475,6114,5194,3067,1331,403,105,18) | # 出た回数 |
m <- sum(x*y)/sum(y) | # 出た回数の平均 |
pd <- m/12 | # 5,6 の出る確率の推定値 |
p <- dbinom(x=0:12, size=12, p=1/3) | # 正しいサイコロ(p = 1/3)のときの確率分布 |
p0 <- c(p[1:10], p[11]+p[12]+p[13]) | # 出る回数 10 回以上をまとめた確率分布 |
rbind(y, p0*sum(y)) | # 観測度数とモデルのもとでの期待度数 |
chisq.test(y, p=p0) | # データが p = 1/3 の二項分布に従っていることの検定 |
q <- dbinom(x=0:12, size=12, p=pd) | # データから推定された確率に基づく二項分布 |
q0 <- c(q[1:10], q[11]+q[12]+q[13]) | # 出る回数 10 回以上にまとめた確率分布 |
rbind(y, q0*sum(y)) | # 観測度数とモデルのもとでの期待度数 |
r <- chisq.test(y, p=q0) | # データが p = pd の二項分布に従っていることの検定 |
pchisq(r$statistic, df=r$parameter-1, lower.tail=F) | # 有意確率(自由度1落とす) |
   女児数    |    0 |    1 |    2 |    3 |    4 |    5 |    6 |
---|---|---|---|---|---|---|---|
度数 |    7    |    45    |    181    |    478    |    829    |    1112    |    1343    |
   女児数    |    7 |    8 |    9 |    10 |    11 |    12 |    合計 |
度数 |    1033    |    670    |    286    |    104    |    24    |    3    | 6155    |
課題:「Weldon のサイコロ実験」と同様に「12人の兄弟中の女児数のデータ」に対して 二項分布に対する適合度検定を行え.
   得点    |   40 未満   |   41 〜 50   |   51 〜 55   |   56 〜 60   |   61 〜 65   |   66 〜 70   |   70 以上   |
---|---|---|---|---|---|---|---|
人数 |    6     |    9     |    12     |    13     |    17     |    14     |    9     |
期待度数 |    4.33     |    13.95     |    11.91     |    13.60     |    12.93     |    10.24     |    13.06     |
しかしながら,英語得点データのように正規分布などの連続分布をデータにあてはめた場合の あてはまりの良さを調べる場合,χ2 適合度検定は,離散化に恣意性が入るのであまり 薦められない.連続型データの場合は,次に述べる Kolmogorov-Smirnov 検定を用いるのが普通である.
eigo <- c( 36,70,56,68,76,60,50,63,62,42,64,60,50,68,71,67, | # 英語得点データ |
50,65,67,57,72,64,61,66,46,80,46,51,59,32,55,65, 65,52,57,64,23,57,53,54,38,71,57,69,77,61,51,64, | |
63,43,65,61,51,69,72,68,53,66,68,58,73,65,62,67, 47,81,47,52,59,33,56,66,67,52,58,65,24,58,54,55) | |
d <- 5 | # ヒストグラムの階級幅 |
op <- par(mfrow = c(1, 2)) | # グラフを横に2つ並べて表示 |
hist(eigo, breaks=seq(0, 100, by=d), xlab="英語得点", ylab="頻度", main="") | |
n <- length(eigo) | # データ数 |
m <- mean(eigo) | # 平均 |
s <- sd(eigo) | # 標準偏差 |
x <- 0:100 | # |
curve(n*d*dnorm(x, m, s), 0, 100, add=TRUE, col="red") | # 推定正規分布重ねて表示 |
title(main="英語得点のヒストグラム") | # |
qqnorm(eigo, xlab="正規分布分位点", ylab="英語得点分位点", main="") | |
qqline(eigo, col=2) | # Q-Q プロット |
title(main="正規 Q - Q プロット") | # |
par(op) | # |
a <- c(0,40,50,55,60,65,70,100) | # 階級の区切り点の定義 |
b <- hist(eigo, breaks=a) | # 各セルに入る人数の計算 |
y <- b$counts | # セルの観測度数 |
l <- length(a) | # |
p <- pnorm(a[2:(l-1)], mean=m, sd=s) | # 階級の区切りまでの累積確率 |
ps <- c(0, p) | # |
pe <- c(p, 1) | # |
q <- pe-ps | # セルの確率分布 |
rbind(y, sum(y)*q) | # セルの観測度数と期待度数 |
r <- chisq.test(y, p=q) | # χ2 適合度検定 |
pchisq(r$statistic, df=r$parameter-2, lower.tail=F) | # 有意確率(自由度 2 落とす) |
plot(ecdf(eigo), xlab="英語得点", do.points=F, ylab="累積確率", main="") | |
curve(pnorm(x, m, s), 0, 100, add=TRUE, col="red") | # 正規分布の累積分布 |
title(main="英語得点の標本累積分布") | # |
ks.test(eigo, "pnorm", m, s) | # Kolmogorov-Smirnov 1標本検定 |
課題: 英語得点データでは 25 点以下の 2 名が他の集団から離れているようにみえる. この 2 名を異常値(out-lier)としてデータから除き,正規分布との適合度検定を行い, あてはまりがどう変化したかを考察せよ.