1 乱数と確率変数

連続型確率変数で、最も基本的な分布は正規分布（normal distribution）です。
平均０、分散１の正規分布は、
　　　\(\phi \left( z\right) =\dfrac{1}{\sqrt{2\pi }}\exp \left( -\dfrac{z^{2}}{2}\right)\sim N(0,1)\)　　（標準正規分布）

ガウスが天文学観測データの測定誤差の研究から発見したと言われています。

関数式の発見や導出は先人の研究に任せ、ここでは乱数を用いて正規分布に関連した確率分布の性質を探究します。

1.1 乱数の生成

次のような手順を基本とします。

  • 乱数を発生させる。
  • ヒストグラムを作成する。
  • 平均・分散を乱数から計算する。
  • 研究者の発見した関数式のグラフを重ねてみる。

ヒストグラムの縦軸は確率密度（Density）としておくとグラフと比較しやすくなります。

1.1.1 正規分布

rnorm(n, mean, sd)：nは乱数の個数、sdは標準偏差

x <- rnorm(1000,0,1)
hist(x,prob=T,breaks = "Scott",xlim=c(-3,3))
curve(dnorm(x,0,1),add=T)

mean(x)
## [1] 0.02703465
var(x)
## [1] 1.018515

平均０を中心に、左右標準偏差3つ分の範囲に99%のデータが収まります。

1.1.2 カイ二乗分布

標準正規分布の二乗はどんな分布をするでしょうか。\(Z\) が -3 から 3 の値だから、\(Z^2\) は 0 から 9 くらいの値で 0が一番多いだろうというのは常識的な判断です。

x <- rnorm(1000,0,1)
hist(x^2,prob=T,breaks = "sturges",xlim=c(0,20))

標準正規分布の平方和はどのような分布になるでしょうか。
５組の正規乱数を発生させ、各要素の二乗５個の和で表される乱数sをつくります。
sの確率密度のヒストグラムを描き、理論上のカイ二乗分布のグラフを重ねます。

s = 0
for (i in 1:5) {
  x <- rnorm(1000)
  s <- s + x^2
}
hist(s,prob=T,breaks = "sturges",xlim=c(0,20))
curve(dchisq(x,5),add=T)

mean(s)
## [1] 4.955147
var(s)
## [1] 9.668334

正規乱数平方５個の和が従う分布を、自由度５のカイ二乗分布といいます。
\(Z_1,Z_2,\cdots,Z_5\sim N(0,1):i.i.d\) のとき、
　　　\(\chi^2=Z_1^2+Z_2^2+\cdots+Z_5^2\sim \chi^2(5)\) 　と書きます。
期待値と分散は、
　　　\(E(\chi^2)=5,~~V(\chi^2)=10\)　と推定されます。

正規乱数の平方k個の和は、自由度kのカイ二乗分布に従います。

2 カイ二乗独立検定

四分表の検定についてシミュレーションを行う。 四分表は２要因２水準のデータをtable形式で集計したものである。

2.1 例4-5(p125)

白血病患者140名をランダムに２群に分け A, B ２種類の治療したときの生存・死亡数を比較したものである。

　 Y=0 (A群) Y=1 (B群)
X=0 (生存) 48 24
X=1 (死亡) 32 36

この四分表に対するカイ二乗値は表の数値を、
\(a=48,~b=24,~c=32,~d=36,~総計~N=140\) として、
　　　　\(\chi^2=\dfrac{(ad-bc)^2N}{(a+b)(c+d)(a+c)(b+d)}=5.49\)　　で計算される。

「四分表のカイ二乗値は自由度１のカイ二乗分布に従う」　ことが分かっている。

無限回の調査をして、無限個の \(\chi^2\) を得たとして、それが本当にカイ二乗分布 \(\chi^2(1)\) に従うだろうか。
カイ二乗分布 \(\chi^2(1)\) に従うのなら、カイ二乗値の平均 mean(X)=1 となるはずである。
シミュレーションで確かめよう。

\((a,b,c,d)\) は多項分布乱数を用いると出力できる。
上と同じ比率を指定した上で1000回出力する。

2.2 多項分布rmultinom.

p1=72/140;p2=80/140
d <- rmultinom(n=1000, size = 140, prob = c(p1*p2,p1*(1-p2),p2*(1-p1),(1-p1)*(1-p2)))  # 乱数を発生
data <- as.data.frame(t(d))
colnames(data) <- c("a", "b", "c", "d")
data[1:2,]
a b c d
44 40 29 27
35 39 33 33

2.2.1 カイ二乗値を計算

X = c()
for (i in 1:1000) {
  a=data[i,1]
  b=data[i,2]
  c=data[i,3]
  d=data[i,4]
  kai=(a*d-b*c)^2*140/((a+b)*(c+d)*(a+c)*(b+d))
  X <-c(X,kai)
}
head(X)
## [1] 0.004770667 0.102013778 1.015030946 0.002382898 0.715600082 0.687960688
hist(X,prob=T,breaks = "sturges")
curve(dchisq(x,1),add=T)

mean(X)
## [1] 1.000238

理論上のグラフとほぼ一致しており、平均mean(X)も１に近い。
\(\chi^2=5.49\) はあり得ない数値として帰無仮説は棄却される。

2.3 別の方法.

有意差をシミュレーションするため、１列目と２列目を別々に出力する。
１列目と２列目のprobが同じであれば四分表は独立である。

X = c()
for (i in 1:1000) {
  a1 <- rmultinom(n=1, size = 80, prob = c(72,68)/140)
  a2 <- rmultinom(n=1, size = 60, prob = c(72,68)/140)
  data <-c(a1,a2)
  a=data[1]
  c=data[2]
  b=data[3]
  d=data[4]
  kai=(a*d-b*c)^2*140/((a+b)*(c+d)*(a+c)*(b+d))
  X <-c(X,kai)
}
mean(X)
## [1] 0.9966872
hist(X,prob=T)
curve(dchisq(x,1),add=T)

上のグラフから見ると kai=6 より大きければ異常であり、A,Bの生存率に差があるとみなせると思うが、kai=3 くらいなら偶然変動でありうる値かもしれない。

2.4 有意差

１行目と２行目のprobを変えてみよう。
統計量kaiが異常な値を出力するようであれば、データが独立でないことを検出できることになる。

2.4.1 生存率 A 0.59、B 0.44

X = c()
for (i in 1:1) {
  a1 <- rmultinom(n=1, size = 80, prob = c(82,58)/140)
  a2 <- rmultinom(n=1, size = 60, prob = c(62,78)/140)
  data <-c(a1,a2)
  a=data[1]
  c=data[2]
  b=data[3]
  d=data[4]
  kai=(a*d-b*c)^2*140/((a+b)*(c+d)*(a+c)*(b+d))
  X <-c(X,kai)
}
mean(X)
## [1] 0.6106454

有意となっているだろうか？
1回の出力なので偶然変動があり、何回かシミュレーションしてみるとkaiの値は変動する。

2.4.2 生存率 A 0.62、B 0.40

X = c()
for (i in 1:1) {
  a1 <- rmultinom(n=1, size = 80, prob = c(90,50)/140)
  a2 <- rmultinom(n=1, size = 60, prob = c(55,85)/140)
  data <-c(a1,a2)
  a=data[1]
  c=data[2]
  b=data[3]
  d=data[4]
  kai=(a*d-b*c)^2*140/((a+b)*(c+d)*(a+c)*(b+d))
  X <-c(X,kai)
}
mean(X)
## [1] 16.03898

出力されたカイ二乗値kaiは大きな値となっており、異常が検出できる。
P値を算出すると、

1-pchisq(mean(X), df=1, ncp=0, lower.tail = TRUE)
## [1] 6.205161e-05

確率的にあり得ないことが起こっていると判断され、帰無仮説は棄却される。

3 RStudio Cloud の導入

初めは、インストールなしで使える”RStudio Cloud” がお勧めです。 ブラウザーで利用可能です。 無料版は容量などの制限がありますが、当面はこれでできると思います。
R Studio Cloudの使い方　が参考になります。 RStudioのサイトは　https://posit.cloud　
（Googleアカウントでも利用できます）

3.1 Rマークダウン

この資料はRマークダウンで作成しています。ワードとエクセルを1つにしたようなソフトで、論文や研究発表に使われています。出力形式を変えれば、ワード形式、パワーポイント形式、html形式で出力できるため効率的です。
Rstudioを使い始めてからは利便性が飛躍的に向上しました。 それまではpagesで作った原稿をホームページビルダーに１行ずつ貼り付けてホームページを作っていましたが、Rstudioは数式まで含めてそのままhtml出力できます。

Rマークダウンで学習し、htmlで書き出して保存すれば記録としても使えます。
学習ノート・記録の代わりになるのではと思います。つまり、生徒も学習記録として自宅でも復習に使えるかもしれません。Rマークダウン はPDFでも出力できますが、ブラウザで見るにはhtmlのほうが便利と思います。

ネットコモンにアップすれば、校内で共有することもできるようです。 https://www.netcommons.org/manual/ユーザーマニュアル/各プラグインについて/情報共有のための機能/キャビネット/ファイルのダウンロード・アップロード

3.2 基本文法

  • 半角2で改行（改行しても出力は改行にならない）
  • #と半角スペースで大項目、## 中項目、### 小項目（番号は自動）
  • 必要に応じてパッケージをインストールする（Toolsから）
  • パッケージを使用するときは、library(name)で読み込む

R Studio Cloudの使い方　が参考になります。 RStudio Cloud　のサイトはGoogleアカウントでも利用できます。

4 はじめての統計15講

小寺平次先生の「はじめての統計15講」の演習問題の解説です。
詳細については、本冊 を参照してください。

4.1 1 度数分布

演習問題１(p7)

参考サイト
http://tips-r.blogspot.com/2014/11/r.html histの利用
https://evrythingonmac.blogspot.com/2017/04/r.html 累積折線

a<- c(160,159,164,158,156,169,153,154,159,161,
      157,157,151,159,166,145,155,161,161,150,
      158,164,165,157,158,160,162,165,159,155,
      162,155,154,157,165,157,154,163,162,158,
      169,155,160,164,158,159,149,154,156,155)
h <- hist(a,breaks=seq(144.5, 172.5, by = 4))

str(h)
## List of 6
##  $ breaks  : num [1:8] 144 148 152 156 160 ...
##  $ counts  : int [1:7] 1 3 12 18 10 4 2
##  $ density : num [1:7] 0.005 0.015 0.06 0.09 0.05 0.02 0.01
##  $ mids    : num [1:7] 146 150 154 158 162 ...
##  $ xname   : chr "a"
##  $ equidist: logi TRUE
##  - attr(*, "class")= chr "histogram"
h$mids#階級値
## [1] 146.5 150.5 154.5 158.5 162.5 166.5 170.5
plot(ecdf(a), col="blue")

h$counts
## [1]  1  3 12 18 10  4  2
f <- cumsum(h$counts)#累積
par(family="HiraKakuProN-W3")
plot(f,type="l",xaxt="n",#デフォルトのX軸メモリを消す
    xlab="身長",ylab="頻度")
axis(1,1:7,labels=h$mids)

n <- length(h$counts) # 階級の数
class_names <- NULL # 階級の名前格納用
for(i in 1:n) {
  class_names[i] <- paste(h$breaks[i], "～", h$breaks[i+1])
}
frequency_table <- data.frame(class=class_names, frequency=h$counts)
frequency_table
class frequency
144.5 ～ 148.5 1
148.5 ～ 152.5 3
152.5 ～ 156.5 12
156.5 ～ 160.5 18
160.5 ～ 164.5 10
164.5 ～ 168.5 4
168.5 ～ 172.5 2

4.2 14 適合度・独立性の検定

参考サイト

https://data-science.gr.jp/implementation/ist_r_chi_squared_test.html カイ二乗検定

演習問題14(p101)
(1)適合度検定

vx=c(57,48,64,51,39,41)
chisq.test(x=vx,p=c(1,1,1,1,1,1)/6)
## 
##  Chi-squared test for given probabilities
## 
## data:  vx
## X-squared = 9.04, df = 5, p-value = 0.1075

(2)独立検定

vx=matrix(c(17,29,14,
            20,62,38,
            11,29,20),nrow=3,byrow=T)
chisq.test(vx)
## 
##  Pearson's Chi-squared test
## 
## data:  vx
## X-squared = 4.1889, df = 4, p-value = 0.381

4.3 終わりに

仕組みが分からなくとも、Rは検定結果を出してくれる。
しかし、「そこにどんな意味があるのか」「どこまで限界があるのか」
検定結果の妥当性を判断するには、乱数によるシミュレーションが経験として必要であると思う。