トップ   編集 凍結解除 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索   ヘルプ   最終更新のRSS

R Note/統計/確率分布

Last-modified: 2016-10-20 (木) 20:00:56
Top / R Note / 統計 / 確率分布

確率分布(確率密度関数)

連続変数(continuous variable)

正規分布(Normal distribution)、ガウス分布(Gaussian distribution)*2

  • [0,10]区間でランダムに生成した10000個の値のヒストグラム

hist(runif(10000) * 10, main = "")

Rstudy1i.png
  • [0,10]区間で、1回につき5個の値をランダムに生成し、その平均を求め(標本平均)、その値を10000個集めたヒストグラム

means<-numeric(10000)

for (i in 1:10000){

means[i]<- mean(runif(5)*10)

}

hist(means,ylim=c(0,1600))

Rstudy1j.png
  • 0~10からランダムに生成した5個の平均なので、値は5の近くになる。
  • このヒストグラムは、どの程度正規分布に一致しているか?
    • 平均と分散を求める

mean(means)

[1] 5.003137

sd(means)

[1] 1.286634
  • 確率密度曲線の横軸の分割点のベクトルを作る(0から10までを100分割する)

xv <- seq(0, 10, 0.1)

xv

 [1]  0.0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1.0  1.1  1.2  1.3
[15]  1.4  1.5  1.6  1.7  1.8  1.9  2.0  2.1  2.2  2.3  2.4  2.5  2.6  2.7
[29]  2.8  2.9  3.0  3.1  3.2  3.3  3.4  3.5  3.6  3.7  3.8  3.9  4.0  4.1
[43]  4.2  4.3  4.4  4.5  4.6  4.7  4.8  4.9  5.0  5.1  5.2  5.3  5.4  5.5
[57]  5.6  5.7  5.8  5.9  6.0  6.1  6.2  6.3  6.4  6.5  6.6  6.7  6.8  6.9
[71]  7.0  7.1  7.2  7.3  7.4  7.5  7.6  7.7  7.8  7.9  8.0  8.1  8.2  8.3
[85]  8.4  8.5  8.6  8.7  8.8  8.9  9.0  9.1  9.2  9.3  9.4  9.5  9.6  9.7
[99]  9.8  9.9 10.0
  • ヒストグラムの面積を求める
    • データ数 10000、ヒストグラムの一つのビンの幅が 0.5 なので、10000×0.5 = 5000 がヒストグラムの面積である。
    • 求める正規分布の面積は1.0なので、各値に5000をかける必要がある。
  • Rの関数を使って正規分布の確率密度関数を生成するには
    • 上記に加えて、以下を実行

      yv<-dnorm(xv,mean=5.003137,sd=1.286634)*5000

      lines(xv,yv)

Rstudy1k.png
  • 標準正規分布(平均0,分散1)

    nd<-seq(-3,3,0.01)

    y<-dnorm(nd)

    plot(nd,y,type="l")

Rstudy1l.png
  • 以下のようにしても作れる。

    curve(dnorm(x, 0, 1), -3, 3)

混合正規分布

指数分布

t分布

カイ2乗分布

ガンマ分布

ベータ分布

F分布

  • 分子の自由度3、分母の自由度16のF分布*3

    curve(df(x, 3, 16), 0, 5)

コーシー分布

アーラン分布

三角分布

ラプラス分布

レイリー分布

ロジスティック分布

パレート分布

ワイブル分布

連続変数の分布の種類を特定する

  •     Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
    1            5.1         3.5          1.4         0.2     setosa
    2            4.9         3.0          1.4         0.2     setosa
    (略)
    51           7.0         3.2          4.7         1.4 versicolor
    52           6.4         3.2          4.5         1.5 versicolor
    (略)
    101          6.3         3.3          6.0         2.5  virginica
    102          5.8         2.7          5.1         1.9  virginica
    (略)
    150          5.9         3.0          5.1         1.8  virginica
  • ここでは、「がくの長さ(Sepal.Length)」にのみ注目します(subset関数の使い方はこちらを参照)。

    setosa.Length = subset(iris, Species=="setosa")$Sepal.Length

    versicolor.Length = subset(iris, Species=="versicolor")$Sepal.Length

    virginica.Length = subset(iris, Species=="virginica")$Sepal.Length

  • ヒストグラムを作ります。

    par(mfrow=c(3,1)) #三段組みのプロットをする

    brk = seq(4, 8, 0.3) #分割範囲を指定、4~8まで、0.3ずつ

    hist(setosa.Length, breaks=brk)

    hist(versicolor.Length, breaks=brk)

    hist(virginica.Length, breaks=brk)

R_002.png
  • t検定やANOVAを使うためには、「標本分布が正規分布をしている」前提を満たす必要があります。そこで、このアヤメのデータが、どの程度正規分布に近いのかを検定します。まずは、標本分布と同じ平均・標準偏差の正規分布曲線を生成して、ヒストグラムに重ねてみます*4

    par(mfrow=c(3,1))

    hist(setosa.Length, breaks=brk)

    setosa.Norm <- dnorm(brk, mean=mean(setosa.Length), sd=sd(setosa.Length))*(0.3*50) #正規分布を作る

    lines(brk,setosa.Norm) #重ねてプロット

    hist(versicolor.Length, breaks=brk)

    versicolor.Norm <- dnorm(brk, mean=mean(versicolor.Length), sd=sd(versicolor.Length))*(0.3*50)

    lines(brk,versicolor.Norm)

    hist(virginica.Length, breaks=brk)

    virginica.Norm <- dnorm(brk, mean=mean(virginica.Length), sd=sd(virginica.Length))*(0.3*50)

    lines(brk,virginica.Norm)

R_003.png
  • 標本分布と正規分布は、あまり大きくずれてはいないようです。
  • 続けて、コルモゴロフ-スミルノフ検定(KS検定) を使って、「標本分布が正規分布(もしくは、指定した他の分布)に従うのか、そうでないのか」を検定します。
    • 最終的にANOVAを使うなら、正規分布へ適合するかどうかだけ検定して、(正規分布へ適合していなかったら)ヒストグラムを見ながら分布の変形処理を行います。
    • 一般化線形モデルを使うなら、何の分布へ適合するか厳密に検定する必要があります。
  • setosaについて関数ks.testを実行します。*5

    ks.test(setosa.Length,"pnorm",mean=mean(setosa.Length),sd=sd(setosa.Length))

  •        One-sample Kolmogorov-Smirnov test
    data:  setosa.Length 
    D = 0.1149, p-value = 0.5245
    alternative hypothesis: two-sided 
    
    警告メッセージ: 
    In ks.test(setosa.Length, "pnorm", mean = mean(setosa.Length), sd = sd(setosa.Length)) :
       タイがあるため、正しい p 値を計算することができません
  • ks.testのオプションとして、「比較したい既知の分布」の名称とそのパラメーターを指定すれば、自動的にその分布と比較してくれます。RWiki Rにおける確率分布 base パッケージで提供される分布一覧 にある分布は全て指定できます。ここでは、既知の分布に「pnorm(p(累積確率分布関数値)+norm(正規分布)*6)」を指定しています。pnormが要求するパラメーターmeanおよびsdは、標本分布のmeanとsdを指定しています。
  • 帰無仮説は「標本分布と、既知の分布に差がない」です。したがって、p値が大きくなれば、指定した既知の分布と標本分布が等しい、ということになります。ここでは、p-value > 0.05 なので、標本分布は正規分布しているといえます。本当?有意差がない=同じ、ではないのでは?
  • 「タイがあるため~」の警告は、まったく同じ数値のデータが複数存在している時に出るようです*7

一群のヒストグラムの歪度と尖度を求める*8

  • 尖度(kurtosis)は分布の頂点の尖り方が正規分布とどの程度一致しているか
  • 歪度(skewness)が負の値をとるときは左に、正の値をとるときは右に歪んでいる

正規分布でない標本を正規分布に近づける

  • 検討1. 標本数を増やす(中心極限定理)
    • t 検定をするとき、サンプリングしたデータが正規分布に従っているか心配で仕方がないという人がいるかもしれませんが、実のところそれほど心配する必要はありません。中心極限定理によって母集団がどんな分布であっても「標本平均の分布は正規分布に従う」ことが分かっているからです。コラム『統計備忘録』第66話「m の分布 -標本平均と中心極限定理-」 より引用)
    • 統計法をある程度学習している者にとって,次のような反論がされるかもしれない。「間隔尺度版統計法は分析データが正規分布でなければならない。しかし順序尺度版データは正規分布とは限らない。順位化を行っても,順位化したデータが正規分布していなければ間隔尺度用統計法の使用は不適切ではないか。」結論から言えば,後述の中心極限定理により,非正規分布データに対しても,標本データ数が大きい場合には,使用可能であることが保証されている。順序尺度データにおける多様な多重比較法(広島大 林先生) より引用)
    • t 検定,分散分析,多重比較法は標本データの正規分布性を仮定して開発された手法であるが,直接分析対象となるものが標本データの分布ではなく,標本平均の分布であることから,中心極限定理が作用する。そのため,少なくともデータ数が多い場合には,これらの統計法を用いても問題はない。順序尺度データにおける多様な多重比較法(広島大 林先生) より引用)
      • もしもt検定や分散分析以外の方法、例えば重回帰分析や線形混合モデルなどを使う場合は、中心極限定理が作用しないので、十分に注意する必要がある。
  • 検討2. 正規分布に近づくような変換を施す
    • L字型(左(0方向)に山が寄っている)分布を正規分布に近づけるには*9
      • 開平変換(平方根変換)(square root transformation)
      • 角変換(逆正弦変換、アークサイン変換)(Angular transformation, Arcsine square root transformation, Arcsine transformation)
      • 対数変換*10
      • 逆数変換
    • J字型(右に山が寄っている)分布を正規分布に近づけるには
      • Log変換 → 逆対数変換でなくて?要確認。
    • 形が分からない分布を正規分布に近づけるには
    • 分布の正規性と分散の等質性という2つの条件は、一般に1つの変換により同時に満足される。なぜならば、正規分布においては、分散の期待値は平均値にかかわらず一定になるからである*11
  • ただし、検討2について、これらの変換にかけても、データは(正確には)正規分布にならないことが指摘されています。
    • 強制選択変数(forced-choice variables)、質問の正解率(question-answer accuracy)などといった、カテゴリカルな結果変数は、角変換にかけても正規分布にならない*12
    • 比例尺度にアークサイン平方根変換を適用しても、ANOVAの結果は誤り(第一種・第二種の検定誤りのリスクが回避できない)。このようなデータでは、ロジットモデル(ロジスティック回帰モデル)を使うべきである*13
    • ポアソン分布(L字型)を平方根変換しても、すそ野のフィットが悪くなる(分布の右端が欠けるため、左右対称の正規分布にはならない)*14 → 以下に例示します。
  • はじめに、L字型のポアソン分布に従う標本集団を作ります*15

    set.seed(1) # 乱数系列の初期化

    n <- 1000 # 処理あたりの標本数

    m <- 4 # 無処理の場合、平均4個の実生が発生するとする。

    ex <- 6 # 処理Xの効果: 6倍とする。

    ey <- 0.25 # 処理Yの効果: 0.25倍とする。

    x <- c(rep(0, n), rep(0, n), rep(1, n), rep(1, n))

    # 処理Xの有無(なし->0,あり->1)のベクトルを作成

    y <- c(rep(0, n), rep(1, n), rep(0, n), rep(1, n))

    # 処理Yの有無(なし->0,あり->1)のベクトルを作成

    z <- c(rpois(n, m), rpois(n, m*ey), rpois(n, m*ex), rpois(n, m*ex*ey))

    # 結果zを乱数を使って生成

    data <- data.frame(x, y, z)

    # データフレームdataを作成

    data.x0y0 <- data[(data$x == 0 & data$y == 0),]

    # x==0かつy==0のデータだけとりだす

    hist(data.x0y0$z) # 頻度分布を確認

R_004.png
  • ポアソン分布への適合度検定の結果(群馬大青木先生のスクリプトを使用しています)

    pois.Hist <- hist(data.x0y0$z)$counts

    source("http://aoki2.si.gunma-u.ac.jp/R/src/poissondist.R", encoding="euc-jp")

    result.pois <- poissondist(pois.Hist)

    result.pois

  •         ポアソン分布への適合度の検定
    
    data:  pois.Hist 
    X-squared = 58.8944, df = 9, p-value = 2.188e-09
    sample estimates:
           n   lambda 
    1000.000    3.034 
    • p-value が0.05より小さいので、ポアソン分布に適合しています。
    • 注:ここではp値は小さいほどポアソン分布に適合しています(帰無仮説:ポアソン分布でない)。
  • 正規分布との適合度も判定してみます。

    ks.test(data.x0y0$z,"pnorm",mean=mean(data.x0y0$z),sd=sd(data.x0y0$z))

  •         One-sample Kolmogorov-Smirnov test
    data:  data.x0y0$z 
    D = 0.1323, p-value = 1.221e-15
    alternative hypothesis: two-sided 
    
     警告メッセージ: 
    In ks.test(data.x0y0$z, "pnorm", mean = mean(data.x0y0$z), sd = sd(data.x0y0$z)) :
       タイがあるため、正しい p 値を計算することができません 
  • 「データ内に同順位がある」という警告が出ています(とりあえず無視します)。
  • p-value は0.05より小さく、正規分布に適合していません。*16
  • なお、正確な正規分布を作って適合度をみると以下くらいのp値になります。

    means<-numeric(10000)

    for (i in 1:10000){

    means[i]<- mean(runif(5)*10)

    }

    hist(means,ylim=c(0,1600))

    ks.test(means,"pnorm",mean=mean(means),sd=sd(means))

  •         One-sample Kolmogorov-Smirnov test
    data:  means 
    D = 0.0085, p-value = 0.4595
    alternative hypothesis: two-sided 
  • 開平変換(平方根変換)をしてみます。

    z2 = sqrt(data.x0y0$z)

    hist(z2)

R_006.png
  • 正規分布と一致しているかどうか検定します。

    ks.test(z2,"pnorm",mean=mean(z2),sd=sd(z2))

  •         One-sample Kolmogorov-Smirnov test
    
    data:  z2 
    D = 0.1297, p-value = 4.774e-15
    alternative hypothesis: two-sided 
    
     警告メッセージ: 
    In ks.test(z2, "pnorm", mean = mean(z2), sd = sd(z2)) :
       タイがあるため、正しい p 値を計算することがで
  • 警告が出たものの、p-value は0.05より小さく、正規分布に適合していません。
  • 注:ここではp値は大きいほど正規分布に適合しています(帰無仮説:正規分布でない)。
  • 角変換(逆正弦変換、アークサイン変換)をしてみます。

    z3 = asin(sqrt(data.x0y0$z/length(data.x0y0$z)))

    ##ヒストグラムを作る

    hist(z3)

R_007.png
  • 正規分布と一致しているかどうか検定します。

    ks.test(z3,"pnorm",mean=mean(z3),sd=sd(z3))

  •         One-sample Kolmogorov-Smirnov test
    
    data:  z3 
    D = 0.1296, p-value = 5.218e-15
    alternative hypothesis: two-sided 
    
     (警告は省略)
    • p-value は0.05より小さく、正規分布に適合していません。
  • このように、L字型分布を変換して、左右対称(に見える)分布に変換しても、正規性の検定を行うと正規分布には適合しません。前述の文献*17によれば、変換ではなく、元々の分布の形状を考慮したモデル(一般化線形モデル)の方がより適当です。

離散変数(categorical variable)

二項分布

  • 試行回数 n = 20、各試行における成功確率 p = 0.7 の二項分布*18

    x <- 0:40

    plot(x, dbinom(x,20,0.7),type="h")

Rstudy_dbinom1.png

ポアソン分布

負の二項分布

ベルヌーイ分布

超幾何分布

多項分布

ゼータ分布、ジップ分布


*1 Wikipedia 確率分布を参考にした
*2 統計学:Rを用いた入門書、Michael J.Crawley 著、野間口謙太郎・菊池泰樹 訳、井立出版、2008 より引用(データセット を c:\\temp に入れておく)
*3 Rによるやさしい統計学, オーム社, 2008. より引用
*4 統計学:Rを用いた入門書,Michael J.Crawley 著,野間口謙太郎・菊池泰樹 訳,井立出版,2008.を参考にしました。
*5 参考:Rでt検定 2 正規性の検定(帯広畜産大 増田先生)
*6 参考:RWiki Rにおける確率分布
*7 参考:Yusuke-Memo タイがあるため・・・
*8 統計学:Rを用いた入門書、Michael J.Crawley 著、野間口謙太郎・菊池泰樹 訳、井立出版、2008 を参考にしました。
*9 参考:田中敏, 山際勇一郎, ユーザーのための教育・心理統計と実験計画法 : 方法の理解から論文の書き方まで, 1992.
*10 対数をとると正規分布になるような分布のことを、対数正規分布 という。
*11 森敏昭(著), 吉田寿夫, 心理学のためのデータ解析テクニカルブック, p.38 より引用
*12 Jaeger, T.F. (2008). Categorical data analysis: Away from ANOVAs (transformation or not) and towards logit mixed models. Journal of Memory and Language, 59, 434-446.
*13 Jaeger, T.F. (2008). Categorical data analysis: Away from ANOVAs (transformation or not) and towards logit mixed models. Journal of Memory and Language, 59, 434-446.
*14 線形モデルにおけるDos&Dont's p.12(立教大森本さん)
*15 フリーソフトRをもちいた生態学データの加工と統計解析 II. GLM (+GUI)編 より引用
*16 注:ここではp値は大きいほど正規分布に適合しています(帰無仮説:正規分布である)。
*17 線形モデルにおけるDos&Dont's(立教大森本さん)
*18 Wikipedia 二項分布 および R と確率分布(同志社大 金先生) を参考にしました。

添付ファイル: fileRstudy1l.png 961件 [詳細] fileRstudy1k.png 968件 [詳細] fileRstudy1j.png 951件 [詳細] fileRstudy1i.png 1004件 [詳細] fileRstudy_dbinom1.png 953件 [詳細] fileR_003.png 1185件 [詳細] fileR_002.png 1230件 [詳細] fileR_007.png 1217件 [詳細] fileR_006.png 1170件 [詳細] fileR_004.png 1195件 [詳細]