医療職からデータサイエンティストへ

統計学、機械学習に関する記事をまとめています。

R の 確率密度関数 ( rnorm, pnorm, qnorm, dnorm ) とは何なのか

統計学を学んでいてRでシミュレーションをしたくなったある日のこと、確率密度関数の使い方がよく分からず困った....
今回は、確率密度関数の使い方を正規分布を例にまとめてみます!!

正規分布の確率密度

Rの確率密度関数を調べてみると、

  • rnorm
  • dnorm
  • pnorm
  • qnorm

があるようです。
最初の1文字が違うだけで、何が違うんだ!と私は思ったので、一つずつ検証してみました。

まずはdnromから

まずはdnormから試してみます。 helpをみると第一引数にquantiles、第二引数にmean、第三引数にsdとあります。第一引数がよく分かりませんが、とりあえず平均0、分散1の標準正規分布で試してみます。

> dnorm(3,mean = 0,sd = 1)
[1] 0.004431848
> dnorm(0,mean = 0,sd = 1)
[1] 0.3989423

出ました!どうやらdnormのdはdensityのdで、正規分布の確率密度の値を返してくれるようです。 ちなみに確率密度関数はこんなのです!

  f(x)=\dfrac{1}{\sqrt{2\pi}\sigma}\exp\left\{-\dfrac{(x-\mu)^{2} }{ 2\sigma^{2}}\right\}

ここのuに0、σに1が入っています。

試しにマイナス無限大からプラス無限大まで積分してみると、、、

> integrate(dnorm,-Inf,Inf)
1 with absolute error < 0.000094

無事に1になりましたー

グラフでみてみると f:id:h-wadsworth02:20181218215751p:plain:w400

確かにX軸が3のところは0に近く、0のところは0.4付近にあります!
これがdnormです。しかし、グラフを書いたりする以外にあまり使い道はなさそうです.....

rnorm

次はrnomです。dnormと違って第一引数はsizeなので、個数のようです。
とりあえず実行してみます。

> rnorm(10,0,1)
 [1]  0.4626736 -0.1871801  0.4449611  0.8936388
 [5] -0.6656489  0.7289715 -2.0133530 -0.7423946
 [9] -0.6219160  0.4033237

出ました!これは標準正規分布に従う10個の乱数のようですね。 適当に1000個ぐらい発生させて、平均と分散を計算してみます。

> rm <-rnorm(1000,0,1)
> mean(rm)
[1] 0.06051489
> var(rm)
[1] 1.016525

平均と分散も理論上の値とだいたいあってます。

qnormとpnorm

続いてqnormとpnormです。この二つは同時に学んだ方が理解が深まります。

> qnorm(0.95,0,1)
[1] 1.644854
> pnorm(1.6448,0,1)
[1] 0.9499945

qnormは第一引数に確率を受けて、確率密度を返し、
pnormは第一引数に確率密度を受けて、確率を返すようです。
グラフでみてみましょう! f:id:h-wadsworth02:20181218225746p:plain:w400
qnormは累積確率が95%になる確率密度の値(1.6448)を返して、pnormは確率密度がマイナス無限大から1.6448までの累積確率(94.9%)を返します。

最後にマイナス無限大から1.6448までの範囲で積分してみます。

> integrate(dnorm , -Inf, 1.6448)
0.9499945 with absolute error < 0.000000086

予想通りほぼ95%になりました!

まとめ

  • dnorm : xの値を受けて、確率密度関数の値を返す。あまり使わなさそう。
  • rnorm : 第一引数は個数を受けて、正規分布に従う乱数を返す。
  • qnorm : 第一引数に確率を受けて、確率密度を返す。
  • pnorm : 第一引数に確率密度を受けて、確率を返す。

使いながら慣れるのが一番ですね!
最後までお読み頂き、ありがとうございました!

参考にさせて頂いたサイト

http://www2.hak.hokkyodai.ac.jp/fukuda/lecture/SocialLinguistics/Rshagen/05distributionR.html
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/60.html
http://cse.naro.affrc.go.jp/minaka/R/R-normal.html