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

R Note/統計/確率分布 の変更点

Top / R Note / 統計 / 確率分布

*確率分布(確率密度関数) [#l9080f24]

-Rのサンプルデータを使ってさまざまな確率分布を再現してみる(([[Wikipedia 確率分布:http://ja.wikipedia.org/wiki/%E7%A2%BA%E7%8E%87%E5%88%86%E5%B8%83]]を参考にした))
-参考
--[[R と確率分布(同志社大 金先生):http://mjin.doshisha.ac.jp/R/09.pdf]]
--[[R による統計解析の基礎(確率密度関数,分布関数,分位点関数についてR での表現の一覧):http://minato.sip21c.org/statlib/stat.pdf#page=25]](中澤先生)
--[[Rにおける確率分布(RjpWiki):http://www.okada.jp.org/RWiki/?R%A4%CB%A4%AA%A4%B1%A4%EB%B3%CE%CE%A8%CA%AC%C9%DB]]

#contents

**連続変数(continuous variable) [#j5bf8607]
***正規分布(Normal distribution)、ガウス分布(Gaussian distribution)(([[統計学:Rを用いた入門書:http://www3.imperial.ac.uk/naturalsciences/research/statisticsusingr]]、Michael J.Crawley 著、野間口謙太郎・菊池泰樹 訳、井立出版、2008 より引用([[データセット:http://www3.imperial.ac.uk/pls/portallive/docs/1/22719696.ZIP]] を c:\\temp に入れておく))) [#z52702ee]
-[0,10]区間でランダムに生成した10000個の値のヒストグラム

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

#ref(Rstudy1i.png,, 50%);
#ref(Rstudy1i.png,,70%);

-[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))

#ref(Rstudy1j.png,, 50%);
#ref(Rstudy1j.png,,70%);

--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)

#ref(Rstudy1k.png,, 50%)
#ref(Rstudy1k.png,,70%);

-標準正規分布(平均0,分散1)
> nd<-seq(-3,3,0.01)
> y<-dnorm(nd)
> plot(nd,y,type="l")

#ref(Rstudy1l.png,, 50%)
#ref(Rstudy1l.png,,70%);

--以下のようにしても作れる。
> curve(dnorm(x, 0, 1), -3, 3)

***混合正規分布 [#b2444351]
-参考:[[苦労する遊び人の覚え書き/ガウス混合分布(GMM):http://qh73xe.jimdo.com/%E3%82%AF%E3%83%A9%E3%82%B9%E3%82%BF%E3%83%AA%E3%83%B3%E3%82%B0%E3%83%91%E3%83%83%E3%82%B1%E3%83%BC%E3%82%B8/%E3%82%AC%E3%82%A6%E3%82%B9%E6%B7%B7%E5%90%88%E5%88%86%E5%B8%83/]]

***指数分布 [#yf70e89e]
***t分布 [#e4aabe66]
***カイ2乗分布 [#p596b125]
***ガンマ分布 [#bed7e603]
***ベータ分布 [#l8a50109]
***F分布 [#o9ea1711]
-分子の自由度3、分母の自由度16のF分布((Rによるやさしい統計学, オーム社, 2008. より引用))
> curve(df(x, 3, 16), 0, 5)

***コーシー分布 [#y6f41dc7]
***アーラン分布 [#o01b2832]
***三角分布 [#db52fad7]
***ラプラス分布 [#k4fdff28]
***レイリー分布 [#le32514e]
***ロジスティック分布 [#faa32112]
***パレート分布 [#ra7a9f8b]
***ワイブル分布 [#fd5d9655]

**連続変数の分布の種類を特定する [#n37e1e32]
-例として、[[Fisher のアヤメのデータ:http://speechresearch.fiw-web.net/113.html#teaa87d8]]を使います。Rではデフォルトで使用可能です。
> data(iris)
> iris

--
     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関数の使い方は[[こちら:http://speechresearch.fiw-web.net/121.html#l688e65c]]を参照)。
>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)

#ref(R_002.png,,80%);

-t検定やANOVAを使うためには、「標本分布が正規分布をしている」前提を満たす必要があります。そこで、このアヤメのデータが、''どの程度正規分布に近いのか''を検定します。まずは、標本分布と同じ平均・標準偏差の正規分布曲線を生成して、ヒストグラムに重ねてみます((統計学:Rを用いた入門書,Michael J.Crawley 著,野間口謙太郎・菊池泰樹 訳,井立出版,2008.を参考にしました。))。
>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)

#ref(R_003.png,,80%);

--標本分布と正規分布は、あまり大きくずれてはいないようです。
-続けて、'''[[コルモゴロフ-スミルノフ検定(KS検定):http://speechresearch.fiw-web.net/restricted/index.php?R%20Note%2F%E8%87%AA%E7%BF%92#l2352aae]]''' を使って、「標本分布が正規分布(もしくは、指定した他の分布)に従うのか、そうでないのか」を検定します。
--最終的にANOVAを使うなら、正規分布へ適合するかどうかだけ検定して、(正規分布へ適合していなかったら)ヒストグラムを見ながら分布の変形処理を行います。
--[[一般化線形モデル:http://speechresearch.fiw-web.net/124.html#q03d38ba]]を使うなら、何の分布へ適合するか厳密に検定する必要があります。

-setosaについて関数ks.testを実行します。((参考:[[Rでt検定 2 正規性の検定(帯広畜産大 増田先生):http://www.obihiro.ac.jp/~masuday/resources/stat/r_t-test02.html]]))
> 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 パッケージで提供される分布一覧:http://www.okada.jp.org/RWiki/?R%A4%CB%A4%AA%A4%B1%A4%EB%B3%CE%CE%A8%CA%AC%C9%DB#xb99d531]] にある分布は全て指定できます。ここでは、既知の分布に「pnorm(p(累積確率分布関数値)+norm(正規分布)((参考:[[RWiki Rにおける確率分布:http://www.okada.jp.org/RWiki/index.php?R%A4%CB%A4%AA%A4%B1%A4%EB%B3%CE%CE%A8%CA%AC%C9%DB]])))」を指定しています。pnormが要求するパラメーターmeanおよびsdは、標本分布のmeanとsdを指定しています。
--帰無仮説は「標本分布と、既知の分布に差がない」です。したがって、p値が大きくなれば、指定した既知の分布と標本分布が等しい、ということになります。''ここでは、p-value > 0.05 なので、標本分布は正規分布しているといえます。'' → &color(red){''本当?有意差がない=同じ、ではないのでは?''};
--「タイがあるため~」の警告は、まったく同じ数値のデータが複数存在している時に出るようです((参考:[[Yusuke-Memo タイがあるため・・・:http://yusuke-memo.blogspot.com/2009/10/blog-post.html]]))。

-KS検定以外の方法として、'''[[カイ二乗検定によって標本分布と正規分布のヒストグラムの類似度を検定する方法:http://speechresearch.fiw-web.net/restricted/index.php?cmd=read&page=R%20Note%2F%E7%B5%B1%E8%A8%88%2F%E6%A4%9C%E5%AE%9A&word=%E3%82%AB%E3%82%A4%E4%BA%8C%E4%B9%97#zf68e56b]]'''や、'''[[正規確率の線形性で判定する方法:http://speechresearch.fiw-web.net/restricted/index.php?%E5%88%86%E6%95%A3%E5%88%86%E6%9E%90#e9d7fd9a]]''' もあります。

***一群のヒストグラムの歪度と尖度を求める(([[統計学:Rを用いた入門書:http://www3.imperial.ac.uk/naturalsciences/research/statisticsusingr]]、Michael J.Crawley 著、野間口謙太郎・菊池泰樹 訳、井立出版、2008 を参考にしました。)) [#y4c168e7]

-[[統計解析フリーソフト R の備忘録頁 ver.3.1 59. 基本統計量の算出:http://cse.naro.affrc.go.jp/takezawa/r-tips/r/59.html]](竹澤先生)

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

-[[MATLAB Note/統計/確率分布/尖度 kurtosis と 歪度 skewness:http://speechresearch.fiw-web.net/115.html#v50f0af0]] でも議論しています。

**正規分布でない標本を正規分布に近づける [#x381a770]
-検討1. 標本数を増やす(中心極限定理)
--'''t 検定をするとき、サンプリングしたデータが正規分布に従っているか心配で仕方がないという人がいるかもしれませんが、実のところそれほど心配する必要はありません。中心極限定理によって母集団がどんな分布であっても「標本平均の分布は正規分布に従う」ことが分かっているからです。'''([[コラム『統計備忘録』第66話「m の分布 -標本平均と中心極限定理-」:http://software.ssri.co.jp/statweb2/column/column0904.html#c066]] より引用)
--'''統計法をある程度学習している者にとって,次のような反論がされるかもしれない。「間隔尺度版統計法は分析データが正規分布でなければならない。しかし順序尺度版データは正規分布とは限らない。順位化を行っても,順位化したデータが正規分布していなければ間隔尺度用統計法の使用は不適切ではないか。」結論から言えば,後述の中心極限定理により,非正規分布データに対しても,標本データ数が大きい場合には,使用可能であることが保証されている。'''([[順序尺度データにおける多様な多重比較法(広島大 林先生):http://ir.lib.hiroshima-u.ac.jp/metadb/up/74006416/Kyoikugaku_kenkyuka_3_54_197_203_hayashi.pdf]] より引用)
--'''t 検定,分散分析,多重比較法は標本データの正規分布性を仮定して開発された手法であるが,直接分析対象となるものが標本データの分布ではなく,標本平均の分布であることから,中心極限定理が作用する。そのため,少なくともデータ数が多い場合には,これらの統計法を用いても問題はない。'''([[順序尺度データにおける多様な多重比較法(広島大 林先生):http://ir.lib.hiroshima-u.ac.jp/metadb/up/74006416/Kyoikugaku_kenkyuka_3_54_197_203_hayashi.pdf]] より引用)
---もしもt検定や分散分析以外の方法、例えば&color(red){''重回帰分析や[[線形混合モデル:http://speechresearch.fiw-web.net/124.html#e5dd0c7d]]などを使う場合は、中心極限定理が作用しない''};ので、十分に注意する必要がある。

-検討2. 正規分布に近づくような変換を施す
--L字型(左(0方向)に山が寄っている)分布を正規分布に近づけるには((参考:田中敏, 山際勇一郎, ユーザーのための教育・心理統計と実験計画法 : 方法の理解から論文の書き方まで, 1992.))
---開平変換(平方根変換)(square root transformation)
---''角変換(逆正弦変換、アークサイン変換)(Angular transformation, Arcsine square root transformation, Arcsine transformation)''
---対数変換((対数をとると正規分布になるような分布のことを、[[対数正規分布:http://aoki2.si.gunma-u.ac.jp/lecture/Bunpu/log-normal.html]] という。))
---逆数変換
--J字型(右に山が寄っている)分布を正規分布に近づけるには
---Log変換  → &color(red){''逆対数変換でなくて?要確認。''};
--形が分からない分布を正規分布に近づけるには
---'''[[Box-Cox変換:http://speechresearch.fiw-web.net/restricted/index.php?cmd=read&page=R%20Note%2F%E7%B5%B1%E8%A8%88%2F%E3%83%99%E3%82%A4%E3%82%BA%E7%B5%B1%E8%A8%88&word=Box-Cox#na8930da]]'''(最適な変換方法を'''[[AIC:http://speechresearch.fiw-web.net/restricted/index.php?cmd=read&page=%E6%83%85%E5%A0%B1%E9%87%8F%E5%9F%BA%E6%BA%96&word=AIC]]'''に基づいて自動的に推定)
---形状のゆがみを統計的に解析する方法として、[[一群のヒストグラムの歪度と尖度を求める:http://speechresearch.fiw-web.net/116.html#y4c168e72]]もあります。
--'''分布の正規性と分散の等質性という2つの条件は、一般に1つの変換により同時に満足される。なぜならば、正規分布においては、分散の期待値は平均値にかかわらず一定になるからである'''((森敏昭(著), 吉田寿夫, 心理学のためのデータ解析テクニカルブック, p.38 より引用))

-ただし、検討2について、''これらの変換にかけても、データは(正確には)正規分布にならない''ことが指摘されています。
--強制選択変数(forced-choice variables)、質問の正解率(question-answer accuracy)などといった、カテゴリカルな結果変数は、角変換にかけても正規分布にならない((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.))
--比例尺度にアークサイン平方根変換を適用しても、ANOVAの結果は誤り(第一種・第二種の検定誤りのリスクが回避できない)。このようなデータでは、[[ロジットモデル(ロジスティック回帰モデル):http://software.ssri.co.jp/statweb2/column/column0711.html#c027]]を使うべきである((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.))
--ポアソン分布(L字型)を平方根変換しても、すそ野のフィットが悪くなる(分布の右端が欠けるため、左右対称の正規分布にはならない)(([[線形モデルにおけるDos&Dont's:http://homepage.mac.com/daichis/osjstat/files/2007MorYam.pdf]] p.12(立教大森本さん))) → 以下に例示します。

--はじめに、L字型のポアソン分布に従う標本集団を作ります(([[フリーソフトRをもちいた生態学データの加工と統計解析 II. GLM (+GUI)編:http://www001.upp.so-net.ne.jp/ito-hi/stat/R2.html]] より引用))。
>  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) # 頻度分布を確認

#ref(R_004.png,,50%);

--ポアソン分布への適合度検定の結果(群馬大青木先生のスクリプトを使用しています)
> 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より小さく、正規分布に適合していません。((注:ここではp値は大きいほど正規分布に適合しています(帰無仮説:正規分布である)。))

--なお、正確な正規分布を作って適合度をみると以下くらいの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)

#ref(R_006.png,,50%);

--正規分布と一致しているかどうか検定します。
> 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)

#ref(R_007.png,,50%);

--正規分布と一致しているかどうか検定します。
> 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より小さく、正規分布に適合していません。

--このように、&color(red){L字型分布を変換して、左右対称(に見える)分布に変換しても、正規性の検定を行うと正規分布には適合しません。};前述の文献(([[線形モデルにおけるDos&Dont's:http://homepage.mac.com/daichis/osjstat/files/2007MorYam.pdf]](立教大森本さん)))によれば、変換ではなく、元々の分布の形状を考慮したモデル([[一般化線形モデル:http://speechresearch.fiw-web.net/113.html#q03d38ba]])の方がより適当です。

**離散変数(categorical variable) [#p6cdf4e3]
***二項分布 [#odea2b90]
-試行回数 n = 20、各試行における成功確率 p = 0.7 の二項分布(([[Wikipedia 二項分布:http://ja.wikipedia.org/wiki/%E4%BA%8C%E9%A0%85%E5%88%86%E5%B8%83]] および [[R と確率分布(同志社大 金先生):http://mjin.doshisha.ac.jp/R/09.pdf#page=3]] を参考にしました。))
> x <- 0:40
> plot(x, dbinom(x,20,0.7),type="h")

#ref(Rstudy_dbinom1.png,,50%);
#ref(Rstudy_dbinom1.png,,70%);

***ポアソン分布 [#lb2eaa87]

***負の二項分布 [#b1e876f9]
***ベルヌーイ分布 [#lf2dc99b]
***超幾何分布 [#j8131a6a]
***多項分布 [#tf6b9bfc]
***ゼータ分布、ジップ分布 [#t33c3243]