読者です 読者をやめる 読者になる 読者になる

Rを通じて統計学を学ぶ備忘録ブログ

SPSSからRに移行したい私のような人向けのR解説ブログ兼学習用備忘録。

第15回〜第21回 確率分布と乱数

第15回〜第21回を修正してまとめて投稿します。第15回〜第21回は「確率分布と乱数」です。

二項分布(binomial distribution)

二項分布は、結果が成功か失敗のいずれか、つまり2値の値しか取らない場合に、独立n回の試行を行ったときの成功を表す離散確率分布。各試行の成功確率pは一定、失敗確率は1-p、このような試行をベルヌーイ試行と呼びます。

確率質量関数

二項分布

二項分布の形状

plot(0, type = "n", xlim = c(1, 50), ylim = c(0, 0.2), xlab = "x", ylab = "prob", main = "binominal")
for(i in 1:5) {
  lines(dbinom(x = 1:50, size = 50, prob = i*0.1), col=i)
}
legend("topright", legend=paste0("binom(50, ", seq(0.1, 0.5, 0.1), ")"), col=1:5, lty=1)

スクリーンショット 2017-01-13 23.08.17

二項分布の乱数

rbinom(n = 20, size = 5, prob = 0.5)
[1] 2 2 4 1 3 2 2 2 3 3 4 4 3 2 2 3 3 3 3 1

ポアソン分布(Poisson distribution)

ポアソン分布は、所与の時間間隔で発生する離散的な事象を数える特定の確率変数Xを持つ離散確率分布のことです。

確率質量関数

ポアソン分布

ポアソン分布の形状

plot(0, type = "n", xlim = c(1, 15), ylim = c(0, 0.4), xlab = "x", ylab = "prob", main = "poison")
for(i in 1:5) {
  lines(dpois(x = 0:20, lambda = i*1), col=i)
}
legend("topright", legend=paste0("λ( ", seq(1, 5, 1), ")"), col=1:5, lty=1)

binom

ポアソン分布の乱数

rpois(n = 20, lambda = 5)
[1]  5  5  6  3  6  4  5 10  2  4  4  2  6  7  3  2  4  7  4  3

正規分布(normal distribution)

正規分布(ガウス分布)は、平均値の付近に散らばるデータの分布を表した連続的な変数に関する確率分布です。中心極限定理により、独立な項目の和として表される確率変数は正規分布に従う。これにより、正規分布統計学や自然科学、社会科学の様々な場面で複雑な現象を簡単に表すモデルとして用いられます。

確率密度関数

正規分布

正規分布の形状

x <- -10:10
plot(0, type = "n", xlim = c(-20, 20), ylim = c(0, 0.2),xlab = "x", ylab = "prob", main = "normal")
for(i in 1:5) {
  curve(dnorm(x,  mean = 0, sd = 2*i), col=i, add = TRUE)
}
legend("topright", legend=paste0("Norm(0, ", seq(2, 10, 2), ")"), col=1:5, lty=1)

スクリーンショット 2017-01-13 23.15.54

正規分布の乱数

round(rnorm(n = 20, mean = 0, sd = 1), 2)
 [1] -0.68 -0.63 -1.17  1.27 -1.19 -1.21 -0.57 -0.93 -1.44  0.10  0.47 -0.60
[13]  1.04  0.76 -0.92 -0.79  0.06 -1.49  2.01  1.75

F分布(F distribution)

F分布は、分散分析や標本間の等分散性の検定等に利用されます。パラメタは自由度d1およびd2。自由度はともに正の整数です。

確率密度関数

スクリーンショット 2017-01-13 23.17.52

F分布の形状

plot(0, type = "n", xlim = c(0, 5), ylim = c(0, 1),xlab = "x", ylab = "prob", main = "F")
for(i in 1:5) {
  curve(df(x, df1 = i, df2 = i), col = i , add = TRUE)
}
legend("topright", legend=paste0("F(", seq(1, 5, 1), ", ", seq(1, 5, 1), ")"), col = 1:5, lty = 1)

スクリーンショット 2017-01-13 23.21.09

F分布の乱数

round(rf(n = 20, df1 = 5, df2 = 5), 2)
 [1]  2.10  1.94 51.45  0.38  2.63  0.77  0.84  0.84  1.74  0.44  1.22  4.73
[13]  3.79  4.04  1.84  0.74  3.08  2.45  0.63  1.00 

ガンマ分布(gamma distribution)

ガンマ分布は、連続確率分布の一種です。形状母数k、尺度母数θの2つのパラメータで特徴づけられます。

確率密度関数

ガンマ分布

ガンマ分布の形状

scale = 1/rateなので,scale=2でrateで指定したい場合はrate=1/2と指定。

plot(0, type = "n", xlim = c(0, 15), ylim = c(0, 1),xlab = "x", ylab = "prob", main = "Gamma")
shape.k <- c(1, 2, 3, 5, 9, 7.5, 0.5)
scale.Θ <- c(2, 2, 2, 1, 0.5, 1, 1)
for(i in 1:7) {
  curve(dgamma(x, shape = shape.k[i], scale = scale.Θ[i]), from=-1, to=15, n=2002, col=i, add=T)
}
legend("topright", legend=paste0("Gamma(", shape.k, " , ", scale.Θ, ")"), col=1:5, lty=1)

スクリーンショット 2017-01-13 23.24.06

ガンマ分布の乱数

round(rgamma(n = 20, shape = 1, scale = 2), 2)
 [1] 0.78 0.92 1.54 0.60 0.84 1.35 0.69 1.89 1.29 3.20 2.07 1.99 1.14 8.47
[15] 0.55 2.22 1.25 0.18 0.42 0.78

カイ二乗分布(Chi-square distribution)

カイ二乗分布は確率分布の一種です。図が示すように、どの自由度kでも、一定以上Zが大きいならば、Zが大きいほどその確率が低くなるので、「正規分布でランダムで値をとったば場合、その値を用いて高々二乗和をとった程度の数値Zがとてつもなく大きくなる確率は少ない」と解釈できます。

確率密度関数

カイ二乗分布

カイ二乗分布の形状

plot(0, type = "n", xlim = c(0, 20), ylim = c(0, 0.5),xlab = "x", ylab = "prob", main = "Chi-square")
x <- 0:20
for(i in 1:5) {
  curve(dchisq(x, df = i), col = i, add = TRUE)
}
legend("topright", legend=paste0("Chisq(", 1:5, ")"), col = 1:5, lty = 1)

スクリーンショット 2017-01-13 23.26.16

カイ二乗分布の乱数

round(rchisq(n = 20, df = 5, ncp = 0), 2)
 [1]  5.70  5.51  7.47  2.58  3.29  4.37  2.24  2.98  6.39  4.99  4.64 10.99
[13]  7.20  4.48  3.55  7.87 11.64  3.07  1.84 17.90

スチューデントのt分布(Student's t distribution)

t分布は、連続確率分布の一つです。正規分布する母集団の平均と分散が未知でサンプルサイズが小さい場合に平均を推定する問題に利用される。また、 2つの平均値の差の統計的有意性を検討するt検定で利用される。t分布は、一般化双曲型分布の特別なケースでもあります。

確率密度関数

スチューデントのt分布

スチューデントのt分布の形状

plot(0, type = "n", xlim = c(-5, 5), ylim = c(0, 0.4),xlab = "x", ylab = "prob", main = "Student's t")
x <- -5:5
for(i in 1:5) {
  curve(dt(x, df = i), col = i, add = TRUE)
}
legend("topright", legend=paste0("t(", c(1:5), ")"), col=1:6, lty=1)

スクリーンショット 2017-01-13 23.28.21

スチューデントのt分布の乱数

round(rt(n = 20, df = 5), 2)
 [1]  0.90 -0.13 -0.99 -2.15  0.31  0.53 -0.87 -0.63  0.37  0.13 -1.88  0.59
[13] -0.26  2.73 -0.29 -0.29  1.12 -1.01 -0.88  1.61

ベータ分布 (beta distribution)

ベータ分布は、連続確率分布のひとつです。ベータ分布のパラメーターはα(>0)およびβ(>0)であり、このパラメタを変えることで、様々な形状の確率密度を表現できる。そのような性質から、ベイズ統計学においては事前分布のモデル分布として多用される。

確率密度関数

ベータ分布

ベータ分布の形状

plot(0, type = "n", xlim = c(0, 1), ylim = c(0, 3),xlab = "x", ylab = "prob", main = "Beta")
shape.α <- c(0.5, 5, 1, 2, 2)
shape.β <- c(0.5, 1, 3, 2, 5)
for(i in 1:5) {
  curve(dbeta(x, shape1 = shape.α[i], shape2 = shape.β[i]), col = i, add = TRUE)
}
legend("topright", legend=paste0("Beta(", shape.α, " , ", shape.β, ")"), col=1:5, lty=1)

スクリーンショット 2017-01-13 23.31.43

ベータ分布の乱数

round(rbeta(n = 20, shape1 = 0.5, shape2 = 0.5), 2)
 [1] 0.01 0.02 0.29 0.51 0.80 0.24 0.67 0.28 0.86 0.46 0.77 0.68 0.10 0.81
[15] 1.00 0.57 0.03 0.98 0.89 0.38

一様分布(Uniform distribution)

一様分布は、離散型あるいは連続型の確率分布である。 サイコロを振ったときの、それぞれの目の出る確率など、すべての事象の起こる確率が等しい現象のモデルです。

確率密度関数

一様分布

一様分布の形状

plot(0, type = "n", xlim = c(-6, 6), ylim = c(0, 0.6),xlab = "x", ylab = "prob", main = "Uniform")
for(i in 1:5) {
  curve(dunif(x, min = -i, max = i), col = i, add = TRUE)
}
legend("topright", legend=paste0("Unif(", -1:-5, " , ", 1:5, ")"), col=1:6, lty=1)

スクリーンショット 2017-01-13 23.33.56

一様分布の乱数

round(runif(n = 20, min = 0, max = 1), 2)
 [1] 0.91 0.54 0.32 0.09 0.99 0.91 0.22 0.66 0.68 0.50 0.66 0.51 0.09 0.76
[15] 0.96 0.84 0.47 0.17 0.44 0.47

対数正規分布(log-normal distribution)

対数正規分布は、連続確率分布の一種です。この分布に従う確率変数の対数をとったとき、対応する分布が正規分布に従うものとして定義される。そのため中心極限定理の乗法的な類似が成り立ち、独立同分布に従う確率変数の積は漸近的に対数正規分布に従うとされています。

確率密度関数

対数正規分布

対数正規分布の形状

plot(0, type = "n", xlim = c(0, 4), ylim = c(0, 2),xlab = "x", ylab = "prob", main = "Logarithmic")
mu <- c(0)
sigma <- c(10, 3/2, 1, 1/2, 1/4, 1/8)
for(i in 1:6) {
  curve(dlnorm(x, mu, sigma[i]), col = i, add = TRUE)
}
grid()
legend("topright", legend=paste0("Lnorm(0, ", sigma, ")"), col = 1:6, lty = 1)

スクリーンショット 2017-01-13 23.36.01

対数正規分布の乱数

round(rlnorm(n = 20, meanlog = 0, sdlog = 1), 2)
 [1] 0.14 0.10 1.72 0.80 0.67 0.86 0.42 3.07 1.61 0.74 3.70 0.22 5.72 0.17
[15] 0.67 0.63 0.51 2.61 0.47 0.35

指数分布(exponential distribution)

指数分布は、ある離散的な事象に対して、ポアソン分布が単位時間当たりの生起確率を示すのに対し、指数分布は生起期間の確率を示します。

確率密度関数

スクリーンショット 2017-01-13 23.39.33

指数分布の形状

plot(0, type = "n", xlim = c(0, 4), ylim = c(0, 2.0), xlab = "x", ylab = "prob", main = "Exponential")
rate <- c(0.5, 1.0, 1.5, 2.0, 2.5)
for(i in 1:5) {
  curve(dexp(x, rate = rate[i] ), col = i, add = TRUE)
}

スクリーンショット 2017-01-13 23.45.34

指数分布の乱数

round(rexp(n = 20, rate = 1), 2)
 [1] 0.81 2.14 0.75 0.11 1.04 0.23 2.79 0.17 0.45 0.08 4.78 0.87 0.12 0.19
[15] 0.84 1.24 0.56 0.97 1.39 0.02

負の二項分布(negative binomial distribution)

負の二項分布は、二つのパラメータを持ちます。成功回数を表す定数rと、おのおのの試行で成功する確率pです。rは、正の整数で、pは、0 から1までの実数です。

確率質量関数

スクリーンショット 2017-01-13 23.47.47

負の二項分布の形状

plot(0, type = "n", xlim = c(1, 500), ylim = c(0, 0.05), xlab = "x", ylab = "prob", main = "nbinominal")
for(i in 1:5) {
  lines(dnbinom(x = 1:500, size = 50, prob = i*0.1), col=i)
}
legend("topright", legend=paste0("nbinom(50, ", seq(0.1, 0.5, 0.1), ")"), col=1:5, lty=1)

スクリーンショット 2017-01-14 0.09.21

負の二項分布の乱数

rnbinom(n = 20, size = 5, prob = 0.5)
 [1] 0 3 4 4 7 2 8 4 2 7 1 7 2 3 3 8 8 5 5 5

多項分布(multinomial distribution)

二項分布を考える際、ベルヌーイ試行によって得られる結果は2通りと考えましたが、これがm通りだとするならば、結果はm個とななります。 この試行をn回繰り返す試行の結果を表すのが多項分布です。

確率質量関数

スクリーンショット 2017-01-13 23.51.21

多項分布の質量関数と乱数

rmultinom(n = 10, size = 5, prob = c(0.1, 0.5, 0.9))
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]    1    0    0    0    0    0    0    0    0     0
[2,]    1    2    3    2    1    1    2    2    2     1
[3,]    3    3    2    3    4    4    3    3    3     4

p1<- 0.2; p2<- 0.3; p3<- 0.5;    # 3 possible outcomes
dmultinom(c(5,5,5), prob=c(p1,p2,p3))  # prob. of 5 each
[1] 0.01838917
dmultinom(c(0,0,9), prob=c(p1,p2,p3))  # prob. of 9 3's
[1] 0.001953125
dmultinom(c(0,2,7), prob=c(p1,p2,p3))  # prob. of 2 2's and 7 3's.
[1] 0.0253125
広告を非表示にする