スタートページ> JavaScript> 他言語> R言語
R言語は統計分野を重点にしているので、統計分布関連の関数が豊富にあります。
分布 乱数発生関数 一様分布 runif(n, min=0, max=1) 二項分布 rbinom(n, size, prob) ポアソン分布 rpois(n, lambda) 正規分布 rnorm(n, mean=0, sd=1) カイ2乗分布 rchisq(n, df, ncp=0) t分布 rf(n, df) F分布 rf(n, df1, df2) ガンマ分布 rgamma(n, shape) ベータ分布 rbeta(n, shape1, shape2) 対数正規分布 rlnorm(n, meanlog = 0, sdlog = 1) ロジスティック rlogis(n) 指数分布 rexp(n, rate = 1) 負二項分布 rnbinom(n, size, prob, mu) 多項分布 rmultinom(n, size, prob) 幾何分布 rgeom(n, prob) 超幾何分布 rhyper(nn, m, n, k) コーシー分布 rcauchy(n,location=0,scale = 1) ワイブル分布 rweibull(n, shape, scale = 1)
発生した乱数を確認したいとき、データが膨大なので、部分的に取り出して表示するのが適切です。
head(r, n=a) tail(r, n=a) rの先頭・末尾のa件を表示。nの省略時は6
hist(r, breaks=seq(xmin, xmax, dx)) ヒストグラム描画の関数。 seq はX軸目盛り(省略可)
(これらは乱数だけでなく、一般の大量データに適用できます。)
h <- hist(r) h$breaks # ヒストグラムのX軸の区分点 h$mids # 中央値(軸の代表値) h$counts # 軸の件数 h$density # 確率密度 例: r <- runif(1000) head(r) # 最初の6件を表示 # 0.7983842 0.9224371 0.9664134 0.8658408 0.7727329 0.18132631 h <- hist(r) h$breaks # 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 h$density # 0.96 0.86 0.83 1.04 1.10 1.18 0.95 1.15 0.90 1.03
それぞれの乱数発生関数の先頭の r を d, p, q にした関数があります。
関数 一様分布の例 d:確率密度関数 dunif(x, min=0, max=1,・・・) p:下側確率 punif(q, min=0, max=1,・・・) q:分位点 qunif(p, min=0, max=1,・・・) 例:「サイコロを3回投げたとき、1の目が2回だけでる確率」 dbinom(r, n, prob) d <- dbinom(2, 3, 1/6) # 0.3230 0.2907 例:「平均0,標準偏差1の正規分布で、X<2(2σ以下)の確率」 p <- pnorm(2, mean=0, sd=1) # p = 0.977 例:「0~1の一様分布で0~qが90%となるq」 q <- qunif(0.9, min=0, max=1) # q = 0.9
runif(n, min=0, max=1)、pbinom(q, size, prob) を例にします。 第1引数(x, q, p) 省略できません 乱数発生関数以外は複数指定することができます。pbinom(c(2, 3), size=10, prob=0.5) 戻り値はベクトルになります。 第2以降の引数 全体を省略できる関数もあります。個々の関数で省略時として記述します。 どちらの記述も同じ意味です。 runif(100, min=0, max=1) と runif(100, 0, 1)、 pbinom(3, size=10, prob=0.5) と pbinom(3, 10, 0.5) 名前付きのときは、順序は自由です。runif(100, max=1, min=0)
一様分布、正規分布、二項分布、t分布について個別説明します。
〓〓〓 一様分布・乱数発生 基本形:runif(n, min, max) 省略時 min=0, max=1 r <- runif(1000) # 0~1の一様分布乱数1000個発生しれrに入れる r <- floor(runif(1000, min=1, max=7)) # サイコロの目(切り捨て整数化) 〓〓〓 一様分布・確率密度関数 基本形:dunif(x, min, max) dunif(0.6) # d = 1 dunif(3, 0, 10) # d = 0.1 (=1/10) 〓〓〓 一様分布・下側確率 基本形:punif(q, min, max) x≦q の確率 punif(0.6) # 0~1 の一様分布で 0<x<0.6 である確率 p=0.6 punif(4, 0, 10) # 0~10 の一様分布で 0<x<4 である確率 p = 0.4 〓〓〓 一様分布・分位点 基本形:qunif(p, min, max) 指定した上側確率に対応するz値(境界値)。下側確率と逆の関係がある。 qunif(0.6) # 0~1 の一様分布で 0<x<q である確率が 0,6 となるqの値 q=0.6 qunif(0.4, 0, 10) # 0~10 の一様分布で 0
〓〓〓 正規分布・乱数発生 基本形:rnorm(n, mean, sd) 省略時 mean=0, sd=1 r <- rnorm(1000) 〓〓〓 正規分布・確率密度関数 基本形:dnorm(x, mean) 右図の赤点のy座標 dnorm(-1.75) # d = 0.086 dnorm(-1) # d = 0.242 dnorm(0) # d = 0.399 dnorm(1) # d = 0.242 dnorm(2) # d = 0.054 〓〓〓 正規分布・下側確率 基本形:pnorm(q, mean, sd) -∞ < x < q の累積確率(青の面積) pnorm(-2.326) # p = 0.01 pnorm(-1.645) # p = 0.05 pnorm(0) # p = 0.5 左右対称 pnorm(1) # p = 0.841 pnorm(1.645) # p = 0.95 pnorm(2.326) # p = 0.99 pnorm(3.090) # p = 0.999 3σ以上の確率は0.01%以下 〓〓〓 正規分布・累積確率 基本形:qnorm(p, mean, sd) 累積確率(青の面積)を与えてqを得る qnorm(0.05) # q = -1.645 qnorm(0.1) # q = -1.281552 qnorm(0.5) # q = 0 qnorm(0.84) # q = 1 qnorm(0.9) # q = 1.28 qnorm(0.95) # q = 1.645 qnorm(0.99) # q = 2.326 qnorm(0.999) # q = 3.090 上位0.1%に入るには、3σ以上が必要
ある母集団の身長が平均167.6cm、標準偏差7.0cm であることがわかっています。
r <- rnorm(1000, mean=167.6, sd=7.0)
hist(r)
により、右のヒストグラムが得られました。
r <- rnorm(1000, mean=167.6, sd=7.0) ● 身長>180 の割合 X→p:下側確率(身長<180)= pnorm(180, mean=167.6, sd=7.0) # p=0.962 → 3.8% ● 身長が高い10%は何センチ以上か p→q:分位点(累積確率=0.9) qnorm(0.9, mean=167.6, sd=7.0) # 176.57cm ● 60%の人はどの範囲にあるか 下側20%と上側20%(下側80%)の分位点 qnorm(0.2, mean=167.6, sd=7.0) # 161.71 qnorm(0.8, mean=167.6, sd=7.0) # 173.49
〓〓〓 二項分布・確率密度関数 基本形:dbinom(r, size, prob) ある事象が起こる確率が plob であり、n回の試行をしたところ、ちょうどr回発生する確率 サイコロを10回投げたところ、1の目が2回でる確率 dbinom(2, 10, 1/6) # 0.29071 1の目が0回、1回、2回でる確率を個別に計算 dbinom(c(0,1,2), 10, 1/6) # 0.1615 0.3230 0.2907 〓〓〓 二項分布・下側確率 基本形:pbinom(q, size, prob) 成功回数が q 回以下である確率 サイコロを10回投げたところ、1の目が0回、1回、2回でる確率の合計 pbinom(2, 10, 1/6) # 0.775 (= 0.1615 + 0.3230 + 0.2907) 〓〓〓 二項分布・累積確率 基本形:qbinom(p, size, prob) サイコロを10回投げたところ、成功回数が q 回以下である確率は、p 以上だが、 成功回数が (n−1) 回以下である確率は、p未満となる q を求める。 qbinom(0.78, 10, 1/6) # 3 上のケースの逆算 qbinom(0.75, 10, 1/6) # 2 qbinom(0.95, 10, 1/6) # 4
t検定に関しては、「基本統計量・検定・信頼区間」statistic-testing.txt を参照してください。そこではt値とp値が重要なものになっています(関数の変数名ではt値はqとなっています)。
このpを表現するのに、信頼水準(95%の信頼度をもって~)という表現と危険水準(5%の危険はあるが~)という表現があります。 qt() は信頼水準によっているので、qt(0.95, 10) ような記述になります。 危険水準に合わせるならならば、qt(1-0.95,10) とするか、-qt(0.05,10) とします。 ここでは、信頼水準とします。
t分布での分析では、一般に母集団は正規分布に従うとされることが多いので、rt が利用されることはあまりありません。確率密度関数の df も同様です。それでここでは、検定などで重要なt値を求める qt、P値を求める pt を扱います。
t分布は、検定や信頼区間の分野で重要な分布です。それに関しては、別章「t検定」で扱います。
基本形:qt(p, df) 自由度dfにおいて、有意水準pに対応するt値を求める
p:両側検定のときはp/2とする。
t値:確率pで帰属仮説を棄却するためのtの最小値
同じ自由度 df(=10) であれば、信頼水準 p を大ににする(検定条件を厳しくする)と、
帰無仮説を棄却できる q (t値)は大になります。
qt(0.9, 10) # t値 = 1.372184
qt(0.95, 10) # t値 = 1.812461 A
qt(0.975,10) # t値 = 2.228139
qt(0.99, 10) # t値 = 2.763769 B
qt(0.995,10) # t値 = 3.169273
有意水準を逆に解釈することもあります。
「95%の信頼度で」を「5%の危険度で」と考えれば、t値は符号が逆になります、
検定などでは、このほうが通常です。
qt(0.05, 10) # t値 = -1.812461 C:Aの逆
同じ信頼水準では、自由度が小だとt値は大になります。
データが少ないとデータの値が少しでも異なると、結果に大きな影響を与えるからです。
qt(0.95, 3) # t値 = 2.353363
qt(0.95, 5) # t値 = 2.015048
qt(0.95, 10) # t値 = 1.812461
qt(0.95, 20) # t値 = 1.724718
qt(0.95, 50) # t値 = 1.675905
基本形:pt(q, df, lower.tail=FALSE) q:t値 p値:検定量t0に対する下側確率P(tz0)になる 1-下側確率P 「… より大きいか」の検定にはこれを用いる 両側検定のときは、この値を2倍する qtの逆関数ですので、いくつかの例を挙げるだけにします。 pt( 1.81, 10) # p値 = 0.950 Aの逆 pt( 2.76, 10) # p値 = 0.990 Bの逆 pt(-1.81, 10) # p値 = 0.050 Cの逆 pt( 1.81, 10, lower.tail=FALSE) # p値 = 0.050 (=1-0.950) pt(-1.81, 10, lower.tail=FALSE) # p値 = 0.950 (=1-0.050)