![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
Rの確率分布システム
かって統計の本には必ず付いていた数表はもはや必要ない。R 自身がオンライン数表である。
R における確率分布関係関数は次のような統一的命名規則を持つ。つまり、分布名を xxx とすると
| dxxx | 密度関数(離散分布なら確率関数)の値 |
| pxxx | は累積確率分布関数値 P[X <= x] もしくは P[X > x] |
| qxxx | クォンタイル値(確率 q に対する P[X <=x]=q または P[X > x]=q となる x の値) |
| rxxx | (疑似)乱数 |
例として正規分布を例にあげると
| dnorm(x, mean=0, sd=1, log = FALSE) |
| pnorm(q, mean=0, sd=1, lower.tail = TRUE, log.p = FALSE) |
| qnorm(p, mean=0, sd=1, lower.tail = TRUE, log.p = FALSE) |
| rnorm(n, mean=0, sd=1) |
ここで引数の意味は
| x | 密度関数を求める位置(のベクトル) |
| q | 確率分布関数を求める位置(のベクトル) |
| p | クォンタイル値を求める確率(のベクトル) |
| n | 発生すべき乱数の数 |
(分布毎に異なる)パラメータ引数
| mean | 平均(のベクトル) | 既定値は 0 |
| sd | 標準偏差の(ベクトル) | 既定値は 1 |
その他のオプションパラメータ
| log | 論理値 | 既定値は FALSE | もし TRUE なら確率は対数値 |
| log.p | 論理値 | 既定値は FALSE | もし TRUE なら確率は対数値 |
| lower.tail | 論理値 | 既定値は TRUE | TRUE なら下側確率 P[X <= x]、もし FALSE なら上側確率 P[X > x] |
貢献パッケージには多次元確率分布を含む、更に多くの関数が提供されている。
| 分布名 | R での名前 | パラメータ | 分布名 | R での名前 | パラメータ |
| ベータ | beta | shape1, shape2, ncp | 2項 | binom | size, prob |
| コーシー | cauchy | location, scale | カイ自乗 | chisq | df, ncp |
| 指数 | exp | rate | F | f | df1, df1, ncp |
| ガンマ | gamma | shape, scale | 幾何 | geom | prob |
| 超幾何 | hyper | m, n, k | 対数正規 | lnorm | meanlog, sdlog |
| ロジスティック | logis | location, scale | 多項 | multinom | n, size, prob |
| 負の2項 | nbinom | size, prob | 正規 | norm | mean, sd |
| ポアソン | pois | lambda | t | t | df, ncp |
| 一様 | unif | min, max | ワイブル | weibull | shape, scale |
| ウィルコクソン | wilcox | m, n | スチューデント化範囲(チューキー) | tukey | nmenas, df, nranges |
(tukey 分布は ptukey と qtukey のみ)
> ?rnorm > dnorm(c(-0.7, 0.5, 1,2, 2.1)) # N(0,1) の密度関数 [1] 0.31225393 0.35206533 0.24197072 0.05399097 0.04398360 > pnorm(c(-0.7, 0.5, 1,2, 2.1)) # N(0,1) の累積確率分布 [1] 0.2419637 0.6914625 0.8413447 0.9772499 0.9821356 > qnorm((0:10)/10) # N(0,1) のクォンタイル [1] -Inf -1.2815516 -0.8416212 -0.5244005 -0.2533471 0.0000000 [7] 0.2533471 0.5244005 0.8416212 1.2815516 Inf > rnorm(10) # N(0,1) の乱数 [1] -0.68562741 1.54190545 -0.75403197 0.25997618 -0.54360004 1.11077348 [7] 0.11983464 -0.05869832 -0.21042407 0.14462494 > mean(rnorm(100, mean=1, sd=1.3)) [1] 1.026899 # 平均 1 を推定 > sd(rnorm(100, mean=1, sd=1.3)) [1] 1.253605 # 標準偏差 1.3 を推定 > pnorm(1.3, mean=1:3, sd=1:3) # (mean,sd)=(1,1),(2,2),(3,3) で計算 [1] 0.6179114 0.3631693 0.2854703 > pnorm(1.3, mean=1:3, sd=1.3) # 平均にベクトル指定 [1] 0.59125296 0.29512923 0.09548885 > pnorm(1.3, mean=1.3, sd=1:3) # 標準偏差にベクトル指定 [1] 0.5 0.5 0.5 > qnorm((0:10)/10, lower.tail=FALSE) #上側クォンタイル P[X > x]=q となる x [1] Inf 1.2815516 0.8416212 0.5244005 0.2533471 0.0000000 [7] -0.2533471 -0.5244005 -0.8416212 -1.2815516 -Inf > pnorm((0:10)/10, lower.tail=FALSE, log.p=TRUE) # log P[X <-x] [1] -0.6931472 -0.7761546 -0.8657395 -0.9621028 -1.0654340 -1.1759118 [7] -1.2937038 -1.4189678 -1.5518513 -1.6924928 -1.8410216
> x <- 1:12 > sample(x) # x のランダムな並べ変え [1] 7 3 10 12 11 1 5 8 4 2 9 6 > sample(10) # 1:10 のランダムな並べ変え [1] 4 7 1 10 6 5 9 3 2 8 > sample(x, replace = TRUE) # x からのランダムな復元抽出(つまり離散一様分布乱数) [1] 11 6 9 9 8 9 4 9 8 8 2 5 # いま風にいえばブートストラップ・サンプリング > sample(1:10, prob=1:10, replace=TRUE) # 抽出比率を指定した復元抽出(つまり離散非一様分布) [1] 9 9 7 3 7 9 7 9 6 9
> sample(c(0, 1), 100, replace = TRUE) # 公平なコイン投げ100回のベルヌイ試行
[1] 0 1 1 0 1 0 1 0 0 0 0 1 1 1 0 0 1 1 0 0 1 0 0 0 0 1 0 1 1 1 1 0 0 1 1 1 0
[38] 0 1 0 0 1 1 0 0 1 1 0 0 1 0 0 1 0 1 0 0 0 1 1 0 1 1 0 0 0 0 0 1 0 0 1 0 1
[75] 0 1 0 0 1 0 0 1 1 1 0 0 0 1 1 0 1 0 1 1 1 1 1 1 1 1
> sample(c("T", "H"), 100, replace = TRUE)
[1] "H" "T" "H" "T" "H" "H" "H" "H" "T" "T" "H" "H" "H" "H" "T" "H" "H" "H"
[19] "H" "T" "T" "T" "T" "H" "H" "H" "H" "H" "T" "H" "T" "T" "T" "T" "T" "T"
[37] "H" "T" "T" "H" "H" "T" "H" "T" "H" "T" "T" "H" "T" "H" "H" "H" "T" "T"
[55] "T" "T" "T" "H" "T" "T" "H" "H" "T" "H" "T" "H" "H" "H" "H" "H" "H" "H"
[73] "H" "H" "T" "T" "H" "H" "H" "T" "H" "H" "H" "T" "H" "H" "H" "T" "H" "T"
[91] "H" "T" "H" "T" "H" "H" "H" "T" "T" "H"
## 単純無作為抽出のパラドックス ## n=10,000, 1000,000, 10,000,000 個の同質有限集団から100個を単純無作為抽出 ## 標本標準偏差(つまり調査精度)はほとんど同一 ## 一万個の母集団の平均情報を100個の標本で調査できるなら、 ## 一千万個の(同質)母集団の平均情報も100個の標本で十分調査できる！ > y <- rnorm(10000) > y1 <- sample(y, size=100); sd(y1)/sqrt(100) [1] 0.1033224 > y <- rnorm(1000000) > y1 <- sample(y, size=100); sd(y1)/sqrt(100) [1] 0.1066535 > y <- rnorm(10000000) > y1 <- sample(y, size=100); sd(y1)/sqrt(100) [1] 0.1021400
