R言語で経験密度関数と経験分布関数を求める方法及びグラフを描画する方法を解説していきます。
この記事では経験密度関数と分布関数を計算する関数とその使用例について見ていきます。
記事で扱うスクリプトは以下からダウントードできます。
経験密度関数と経験分布関数
経験分布を計算する関数を紹介します。
経験密度関数と経験分布関数のための関数をいかにまとめました。
経験密度関数
経験密度関数を計算するには次の関数densityを用います。
デフォルトでガウスカーネルの下での経験密度を計算することができ、引数 によってカーネルを変更することもできます。
density(x, …)
# S3 method for default
density(x, bw = "nrd0", adjust = 1, kernel = c("gaussian", "epanechnikov", "rectangular", "triangular", "biweight", "cosine", "optcosine"), weights = NULL, window = kernel, width, give.Rkern = FALSE, n = 512, from, to, cut = 3, na.rm = FALSE, …)
| 引数 | 説明 |
| x | 経験密度を計算するデータ。 |
| bw | カーネル平滑化のパラメーター。 |
| adjust | カーネル平滑化のパラメーターを補正する値。パラメータは正確には、adjust × bwで計算される。 |
| kernel, window | character型。平滑化に用いるカーネルを決定する。"gaussian", "rectangular", "triangular", "epanechnikov", "biweight", "cosine", "optcosine"を選択できる。 |
| weight | 非負のnumeric型のベクトル。観測値の重みを表す。 |
| width | widthが与えられていてbwが与えられていない場合、bwがwidthの値で与えられる。 |
| give.Rkern | logical型。TRUEの場合、密度は推定されない。代わりに選択したkernelのcanonical bandwidthが返される。 |
| n | 経験密度を推定する等間隔点の数。 |
| from,to | 経験密度を計算する範囲。デフォルトだとcut × bw。 |
| cut | from,toの範囲に関係するパラメーター。経験密度の裾が0になるようにする。 |
| na.rm | logical型。TRUEのとき欠測値はxから除外される。 |
| ... | 他のメソッドに渡す引数。 |
経験分布関数
経験分布関数を計算するには関数ecdfを用います。
ecdf(x)
# S3 method for ecdf
plot(x, …, ylab="Fn(x)", verticals = FALSE, col.01line = "gray70", pch = 19)
# S3 method for ecdf
print(x, digits= getOption("digits") - 2, …)
# S3 method for ecdf
summary(object, …) # S3 method for ecdf quantile(x, …)
関数ecdfや関連する関数の引数は以下の通り。
| 引数 | 説明 |
| x, object | xはnumeric型の観測値から成るベクトル。objectはecdfから継承されるオブジェクト。 |
| ... | plotに渡す引数。 |
| ylab | y軸のラベル。 |
| verticals | logical型。経験分布関数の縦線を描くかどうか。 |
| col.01line | y=0とy=1の横線の色を表すnumericまたはcharacter型。 |
| pch | 点の形を指定する。 |
| digits | 有効数字を表す値。 |
実行例
実際にRで経験密度関数と経験分布関数を計算していきます。
計算した経験分布のグラフを描画するためにパッケージggplot2を用います。事前にインストールしておいてください。
1 2 | install.packages("ggplot2") library(ggplot2) |
経験分布の計算とプロット
次のquakes（地震）のデータについて経験分布を求めたいと思います。
1 2 3 4 5 6 7 8 9 | > dataset <- quakes > head(dataset) lat long depth mag stations 1 -20.42 181.62 562 4.8 41 2 -20.62 181.03 650 4.2 15 3 -26.00 184.10 42 5.4 43 4 -17.97 181.66 626 4.1 19 5 -20.42 181.96 649 4.0 11 6 -19.68 184.31 195 4.0 12 |
quakesのmagについてヒストグラムを描くと次のようになります。
magの経験分布を計算していきたいと思います。
経験密度関数を計算するには関数density、経験分布関数を計算するには関数ecdfを用います。
次のように関数の引数に経験分布を求めたいデータを渡すことで、そのデータの経験分布を得ることができます。
1 2 | density_quakes <- density(dataset$mag) ecdf_quakes <- ecdf(dataset$mag) |
density_quakesは次のようにxやyなどの経験密度関数の引数とそれに対応する経験密度などが格納されています。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 | > density_quakes Call: density.default(x = dataset$mag) Data: dataset$mag (1000 obs.); Bandwidth 'bw' = 0.09105 x y Min. :3.727 Min. :0.0000497 1st Qu.:4.463 1st Qu.:0.0196331 Median :5.200 Median :0.1868643 Mean :5.200 Mean :0.3390539 3rd Qu.:5.937 3rd Qu.:0.6270866 Max. :6.673 Max. :1.0294696 |
経験密度を得るには次のようにdensity_quakesの後ろにxとyをつけます。
1 | > epdf_x_quakes <- density |
_quakes$x #経験密度関数の引数 > epdf_y_quakes <- density_quakes$y #xに対応する経験密度 > > head(epdf_x_quakes) [1] 3.726836 3.732601 3.738367 3.744133 3.749899 3.755665 > > head(epdf_y_quake) [1] 0.002341667 0.002829994 0.003400496 0.004063127 0.004858902 0.005791372
ecdf_quakesは関数として与えられており、次のようにecdf_quakes(x)の形で経験分布関数を計算することができます。
1 2 3 4 5 | > x <- c(4.5, 4.6, 4.7, 4.8, 4.9) > ecdf_y_quakes <- ecdf_quakes(x) > > ecdf_y_quakes [1] 0.484 0.585 0.683 0.748 0.802 |
続いて、経験分布のプロット方法について見ていきます。
plotを用いることで簡単に経験密度関数と経験分布関数を描画することができます。
1 | plot(density_quakes) #経験密度関数 |
magの経験密度関数
階段形式の経験分布関数を描きたいときは、引数verticalsをTRUE、do.pをFALSEにします。
1 2 | plot(ecdf_quakes) #経験分布関数 plot(ecdf_quakes, verticals = TRUE,do.p = FALSE) #階段形式の経験分布関数 |
magの経験分布関数
magの経験分布関数（階段）
計算した経験密度をggplotでプロットするには、次のように関数geom_lineやgeom_ribbonを使います。
1 2 3 4 5 6 | dataset_denPlot <- data.frame(x = density_quakes$x, y = density_quakes$y) ggplot(dataset_denPlot, aes(x = x, y = y)) + geom_line() + geom_ribbon(aes(ymin = 0, ymax = y), alpha = 0.5) + labs(y = "epdf") #経験密度関数 |
magの経験密度関数（ggplot2）
また、ggplot2で用意されているgeom_densityでも経験密度をプロットすることができます。
実際のデータの最小値から最大値までの密度が描かれます。
1 2 3 | ggplot(dataset, aes(x = mag)) + geom_density(fill = "black", alpha = 0.5) + labs(y = "epdf") #geom_densityの場合 |
magの経験密度関数（geom_density）
経験分布関数をggplotでプロットするには、関数geom_stepを使います。
eom_stepは最小のxと対応するyの座標からプロットし始めるため、 xの最小値から0.1差し引いたところからプロットするようにデータを調整しています。
1 2 3 4 5 6 | dataset_ecdfPlot <- rbind(data.frame(x = dataset$mag, y = ecdf_quakes(dataset$mag)), data.frame(x = min(dataset$mag) - 0.1 , y = 0)) ggplot(dataset_ecdfPlot, aes(x = x, y = y)) + geom_step() + labs(y = "ecdf") #経験分布関数 |
magの経験分布関数（ggplot2）
また、ggplot2にはstat_ecdfという経験分布関数を描画する関数があり、こちらでも同様のグラフを描くことができます。
1 2 3 | ggplot(dataset, aes(x = mag)) + stat_ecdf() + labs(y = "ecdf") #stat_ecdfの場合 |
magの経験分布関数（geom_ecdf）
多群の場合
次に多群（複数の因子をもつ）のデータの経験分布の計算とプロットの仕方について紹介します。
データセットとして次のirisを使います。
1 2 3 4 5 6 7 8 9 | > dataset <- iris > head(dataset) 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 3 4.7 3.2 1.3 0.2 setosa 4 4.6 3.1 1.5 0.2 setosa 5 5.0 3.6 1.4 0.2 setosa 6 5.4 3.9 1.7 0.4 setosa |
各SpeciesのSepal.Lengthのヒストグラムは次のようになります。
下図のヒストグラムを持つSepal.Lengthについて、経験密度関数と分布関数を計算していきたいと思います。
1 2 | ggplot(dataset, aes(x = Sepal.Length, color = Species, fill = Species)) + geom_histogram(position = "identity", alpha = 0.75, bins = 15) |
Sepal.Lengthのヒストグラム
Speciesの値ごとにSepal.Lengthの経験分布を計算するには、次のようにtapplyを用いると便利です。
factorのレベルごとに関数を適用することで、各Speciesごとの経験密度関数と経験分布関数を計算しています。
1 2 | densities_sepalLength <- tapply(dataset$Sepal.Length, dataset$Species, density) ecdfs_sepalLength <- tapply(dataset$Sepal.Length, dataset$Species, ecdf) |
計算した経験密度関数はそれぞれ次の通りです。
例としてsetosaのSepal.Lengthの経験分布の計算結果のみ載せています。他のSpeciesについても同様に計算できます。
1 2 3 4 5 6 7 8 | > epdf_x_setosa <- densities_sepalLength$setosa$x #setosaのSepal.Lengthの経験密度関数の引数 > epdf_y_setosa <- densities_sepalLength$setosa$x #xに対応するsetosaのSepal.Lengthの経験密度 > > head(epdf_x_setosa) [1] 3.931426 3.935804 3.940182 3.944560 3.948938 3.953316 > > head(epdf_y_setosa) [1] 3.931426 3.935804 3.940182 3.944560 3.948938 3.953316 |
経験分布関数は次の通りです。
1 2 3 4 5 | > x <- c(4.5, 4.6, 4.7, 4.8, 4.9) > ecdf_y_setosa <-ecdfs_sepalLength$setosa(x) > > ecdf_y_setosa [1] 0.10 0.18 0.22 0.32 0.40 |
最後に多群の場合の経験分布のプロット方法について見ていきます。
ggplotで経験密度関数と経験分布関数をプロットするには、上で紹介したように、それぞれ関数geom_density、stat_ecdfを使います。
1 2 3 4 5 6 7 | ggplot(dataset, aes(x = Sepal.Length, color = Species, fill = Species)) + geom_density(alpha = 0.5) + labs(y = "epdf") #geom_densityの場合 ggplot(dataset, aes(x = Sepal.Length, color = Species)) + stat_ecdf() + labs(y = "ecdf") #stat_ecdfの場合 |
Sepal.Lengthの経験密度関数（geom_density）
Sepal.Lengthの経験分布関数（stat_ecdf）
また、上で計算した経験密度と経験分布を基にggplotでプロットしたい場合は、次のようにデータを加工するプロセスが必要になります。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 | dataset_denPlot <- data.frame(x = NULL, y = NULL, Species = NULL) dataset_ecdfPlot <- data.frame(x = NULL, y = NULL, Species = NULL) for (specy in names(densities_sepalLength)) { dataset_density <- data.frame(x = densities_sepalLength[[specy]]$x, y = densities_sepalLength[[specy]]$y, Species = specy) dataset_ecdf <- data.frame(x = dataset$Sepal.Length[dataset$Species == specy], y = ecdfs_sepalLength[[specy]](dataset$Sepal.Length[dataset$Species == specy]), Species = specy) dataset_denPlot <- rbind(dataset_denPlot, dataset_density) dataset_ecdfPlot <- rbind(dataset_ecdfPlot, dataset_ecdf) } dataset_ecdfPlot <- rbind(dataset_ecdfPlot, data.frame(x = max(dataset$Sepal.Length) + 0.1, y = 1, Species = levels(dataset$Species)), data.frame(x = min(dataset$Sepal.Length) - 0.1, y = 0, Species = levels(dataset$Species))) ggplot(dataset_denPlot, aes(x = x, y = y, fill = Species)) + geom_line(aes(color = Species)) + geom_ribbon(aes(ymin = 0, ymax = y), alpha = 0.5) + labs(y = "epdf") #経験密度関数 ggplot(dataset_ecdfPlot, aes(x = x, y = y, color = Species)) + geom_step() + labs(y = "ecdf") #経験分布関数 |
Sepal.Lengthの経験密度関数（ggplot2）
Sepal.Lengthの経験分布関数（ggplot2）
まとめ
R言語で経験分布を計算する関数とその使用例について紹介しました。
経験密度関数を計算するには関数density、経験分布関数を計算するには関数ecdfを用います。
