スタートページ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分布について個別説明します。

一様分布 unif

〓〓〓 一様分布・乱数発生
基本形: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


二項分布 binom

〓〓〓  二項分布・確率密度関数
基本形: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分布 t

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検定」で扱います。

t分布・有意水準→t値 qt

基本形: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

t分布・t値→p値

基本形: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)