1 分布の推定.

  • QQプロットは、（片方の\(q\)分位数 , もう片方の\(q\)分位数）をプロットしたもの。
  • PPプロットは、それぞれの分布関数 \(F_1,F_2\)に対し、\((F_1(x),~F_2(x))\)をプロットしたもの。

注意
\(0<q<1\) は媒介変数としての役割を果たすので、曲線の媒介変数表示をしていることになります。

  • ２つの分布が一致しているかを確認したいときは、Q-QプロットまたはP-Pプロットが y=x に近いかで判断
  • ２つの分布が同じ種類の分布なのかを確認したいときは、Q-Qプロットが直線かどうかで判断

具体例で学ぶ数学様「Q-QプロットとP-Pプロットの意味と違い」参照

2 基礎事項の確認.

2.1 qnorm.

qnormは、分布関数の逆関数.

n <- 10
y1 <- (1:n)/n 
plot(y1,qnorm(y1),xlim=c(0,1),ylim=c(-2,2))

data <- data.frame("y1"=y1,"qnorm(y1)"=qnorm(y1))
data
y1 qnorm.y1.
0.1 -1.2815516
0.2 -0.8416212
0.3 -0.5244005
0.4 -0.2533471
0.5 0.0000000
0.6 0.2533471
0.7 0.5244005
0.8 0.8416212
0.9 1.2815516
1.0 Inf

y1=1 のときの値が \(\infty\) となり欠測値となってしまうので修正する必要があるだろう。.

座標を入れ替えたグラフに正規分布の分布関数の理論曲線と重ねてみると、曲線上にあることが分かる。
つまり、qnormは分布関数の逆関数になっている。

plot(qnorm(y1),y1,ylim=c(0,1),xlim=c(-2,2))
curve(pnorm(x,0,1),add=T)

次に、

  • qnormのグラフ.
  • 縦軸の目盛をy1の値にしたもの.

を並置する。

par(mfrow = c(1,2))
plot(y1,qnorm(y1),xlim=c(0,1),ylim=c(-2,2))
plot(y1,qnorm(y1),xlim=c(0,1),ylim=c(-2,2),yaxt="n")
axis(2, qnorm(y1), y1)#座標目盛

右のグラフは、縦軸の目盛にy1の値を表示しているので等間隔ではない。
右図の目盛にすると、計算なしで手書きのプロットが可能になる。

2.1.1 半整数値補正.

n <- 10
y <- (1:n-0.5)/n 

par(mfrow = c(1,2))
plot(y,qnorm(y),xlim=c(0,1),ylim=c(-2,2))
plot(y,qnorm(y),xlim=c(0,1),ylim=c(-2,2),yaxt="n")
axis(2, qnorm(y), y)#縦座標目盛にyをprint

右のグラフは、縦軸の目盛にyの値を表示しているので等間隔ではない。

data <- data.frame(X=y,Y=qnorm(y))
data
X Y
0.05 -1.6448536
0.15 -1.0364334
0.25 -0.6744898
0.35 -0.3853205
0.45 -0.1256613
0.55 0.1256613
0.65 0.3853205
0.75 0.6744898
0.85 1.0364334
0.95 1.6448536 
plot(qnorm(y),y,ylim=c(0,1),xlim=c(-2,2))
curve(pnorm(x,0,1),add=T)

かなり見やすくなっている。

2.2 分位点.

確率の分位点の取り方を 群馬大学 青木繁先生「正規確率紙」. に従って変えてみる。

probs <- c(0.01, 0.1, 1, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 95, 99, 99.9, 99.99)/100
par(mfrow = c(1,2))
plot(probs,qnorm(probs),xlim=c(0,1),ylim=c(-3,3))
plot(probs,qnorm(probs),xlim=c(0,1),ylim=c(-3,3),yaxt="n")
axis(2, qnorm(probs), probs)#縦座標目盛にyをprint

data <- data.frame("probs"=probs,"qnorm(probs)"=qnorm(probs))
data
probs qnorm.probs.
0.0001 -3.7190165
0.0010 -3.0902323
0.0100 -2.3263479
0.0500 -1.6448536
0.1000 -1.2815516
0.2000 -0.8416212
0.3000 -0.5244005
0.4000 -0.2533471
0.5000 0.0000000
0.6000 0.2533471
0.7000 0.5244005
0.8000 0.8416212
0.9000 1.2815516
0.9500 1.6448536
0.9900 2.3263479
0.9990 3.0902323
0.9999 3.7190165 
plot(qnorm(probs),probs,ylim=c(0,1),xlim=c(-3,3))
curve(pnorm(x,0,1),add=T)

かなり見やすくなっている。

3 経験分布関数.

群馬大学 青木繁先生「正規確率紙」 をみると、データから分布関数を描く方法が示されている。

  • データを昇順にソートする.
  • \(k\)番目のデータの確率を \(k/n\) とする（nはデータの総数）.

3.1 正規乱数.

#x <- runif(500, min=0, max=1)
x <- rnorm(200,0,1)
n <- length(x)               # データの個数
x <- sort(x)         # 昇順に並べ替える
#x
y <- (1:n)/n　　 # 累積相対度数
#y
plot(x, y)
curve(pnorm(x, mean=0, sd=1),add=T,col="red")

データの数が200個になると、ほぼ理論曲線上に乗ってくる。

3.2 一様乱数.

x <- runif(200, min=0, max=1)
n <- length(x)               # データの個数
x <- sort(x)         # 昇順に並べ替える
#x
y <- (1:n)/n　　 # 累積相対度数
#y
plot(x, y)
curve(punif(x,min=0, max=1),add=T,col="red")

一様分布の分布関数は、直線 \(y=x\) である。

4 正規確率紙.

分布の比較には、

  • 分布関数を比較することが考えられる.

4.1 qnorm.

4.2 経験確率分布.

\(N(0,1)\) に従う正規乱数500個から、分布関数を構成する。

probs <- c(0.01, 0.1, 1, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 95, 99, 99.9, 99.99)/100

x <- rnorm(500,0,1)
n <- length(x)               # データの個数
x <- sort(x)         # 昇順に並べ替える
y <- (1:n-0.5)/n　　 # 累積相対度数
plot(x, y)
curve(pnorm(x, mean=0, sd=1),add=T,col="red")

第２座標に経験確率分布に対応する理論値を置く。
データが理論分布通りなら、第１座標、第２座標は一致し、プロットは \(y=x\) 上に並ぶはずである。

plot(x, qnorm(y))
abline(v=0,h=0,col="lightblue")

4.2.1 平均・標準偏差を変える.

x <- rnorm(500,5,1)
n <- length(x)               # データの個数
x <- sort(x)         # 昇順に並べ替える
y <- (1:n-0.5)/n　　 # 累積相対度数
plot(x, y)

plot(x, qnorm(y))
abline(v=0)

4.3 理論値.

x <- rnorm(50,0,1)
x <- sort(x)         # 昇順に並べ替える
n <- length(x) 

y <- pnorm(x)#理論値
y1 <- (1:n-0.5)/n　　# 経験値

par(mfrow = c(1,2))
plot(x,qnorm(y),yaxt="n")
axis(2, qnorm(y1), y1)#座標目盛

plot(x,qnorm(y1),yaxt="n")
axis(2, qnorm(y1), y1)#座標目盛

4.4 整理.

x <- rnorm(50,0,1)
y <- pnorm(x)#理論値
x1 <- sort(x)         # 昇順に並べ替える
y1 <- (1:n-0.5)/n　　# 累積相対度数
par(mfrow = c(1,2))
plot(x,y)
plot(x1,y1)

「データが正規分布に従っているか作図によって検証する方法」を整理する。

  • xを昇順ソートしたx1から経験確率分布pを得る.
    逆対応：\(p\rightarrow x1\).
  • 同じpに対し、正規分布累積関数の逆関数qnormによりx2を対応させる. 逆対応：\(p\rightarrow x2\).
  • x1が正規分布に従っていれば x1=x2 となり、\((x1,x2)\)は直線 y=x 上に乗る.
x <- rnorm(10,0,1)
n <- length(x) 
x1 <- sort(x)         # 昇順に並べ替える
p <- (1:n-0.5)/n　　# 経験累積確率分布
x2 <- qnorm(p)#理論値
plot(x1,x2)
abline(a=0,b=1,col="lightblue")

#par(mfrow = c(1,2))
#plot(x,qnorm(y),yaxt="n")
#axis(2, qnorm(y1), y1)#座標目盛

#plot(x,qnorm(y1),yaxt="n")
#axis(2, qnorm(y1), y1)#座標目盛

4.5 手書きプロットの工夫.

qnorm(p)を手書きでプロットするには、

  • 目盛をqnorm(probs)で構成し.
  • 目盛に従ってpを記入すれば、qnorm(p)の高さになる.

4.5.1 \(y=x^2\) のプロット.

x <- c(1:6)
par(mfrow = c(1,2))
plot(x,x^2)
plot(x,x^2,yaxt="n")
axis(2, x^2, x)#座標目盛

目盛を \(y=x^2\) で描いておけば、
　　xの値をプロットが \(y=x^2\) の高さになっている。

5 演習15(2).

度数分布表から、正規確率紙を用いて、

  • 正規性の判定.
  • 母平均、母標準偏差の推定.

を行う。

  • データから経験累積分布関数の値を計算し.
  • 経験累積分布関数の分位点から理論分布関数によりxを逆算する.

簡単のため、mean=53, sd=13　を用いて理論分布関数を描くと、ほぼ一致していることがわかる（左図）。

  • xの逆算で qnorm(y1,50,10) を用いると、(x1,x2) 直線上に位置する.
  • 回帰直線を描くと、50%点に x=53 が対応し、これが平均の予測値となる（右図）.
y <- c(3,12,18,35,15,5)
x1 <- c(29,39,49,59,69,79)
y1 <- cumsum(c(3,12,18,35,15,5))/90#probsに相当
par(mfrow = c(1,2))
plot(x1,y1)
curve(pnorm(x, mean=53, sd=13),add=T,col="red")
#予測点
x2 <- qnorm(y1,50,10)
data <- data.frame(x1=x1,x2=x2)
model <- lm(x2 ~ x1, data = data)

plot(x1,x2,yaxt="n")
abline(model,col="lightblue")
abline(h=50,col="red",lty = 3)
abline(v=53,col="red",lty = 3)
axis(2, x2, round(x2,0))#座標目盛

qnorm(y1,53,13) を用いれば、回帰直線は y=x になる。 平均・標準偏差を変えても直線性は保たれる。
（一次変換で移るからと考えられる）.

右図を大きく描いてみよう.

probs <- c(0.01, 0.1, 1, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 95, 99, 99.9, 99.99)/100
x3 <- qnorm(probs,50,10)
x3
##  [1] 12.80984 19.09768 26.73652 33.55146 37.18448 41.58379 44.75599 47.46653
##  [9] 50.00000 52.53347 55.24401 58.41621 62.81552 66.44854 73.26348 80.90232
## [17] 87.19016
plot(x1,x2,yaxt="n")
abline(model,col="lightblue")
#abline(a=1,b=1,col="lightblue")
axis(2, x3, round(x3,0))#座標目盛
abline(h=50,col="red",lty = 3)
abline(v=53,col="red",lty = 3)

注意. 度数分布表でなく元のデータがある場合は、青木先生の関数を用いるのが良いと思う。群馬大学 青木繁先生「正規確率紙」 参照.