ヒストグラムと密度の推定
(グラフィックス参考実例集に戻る。Rのグラフィックスパラメータを参照する。)
φを平均-1,分散1の正規分布に従う確率変数,ψを平均2,分散1の正規分布に従う確率変数として
f(x) = 0.6φ(x)+0.4ψ(x)
となる密度関数 truedensity() を定義する.
truedensity <- function (x) { 0.6/sqrt(2*pi)*exp(-(x+1)^2/2) + 0.4/sqrt(2*pi)*exp(-(x-2)^2/2) }
> curve(truedensity, xlim=c(-6,6), ylim=c(0,0.3), col=2)
次に, f(x) に従う乱数を生成する関数を定義する.
generator_tmp <- function(n) { data1 <- rnorm(n)-1 # φ(x) に従う乱数 data2 <- rnorm(n)+2 # ψ(x) に従う乱数 data3 <- (runif(n) <= 0.6) # 確率0.6でdata1を採択し, data1*data3+data2*(1-data3) # 確率0.4でdata2を採択するための手続 }
上と同じ動作をする関数で,最適化された関数を以下に挙げる(速度は約2倍,メモリ使用量は約3/4).
generator <- function(n) { (rnorm(n)-1)+3*(runif(n) > 0.6) } > data <- generator(1000)
【参考】φを N1(μ1, σ1^2) に従う確率変数,ψを N2(μ2, σ2^2) に従う確率変数として
f(x) = p*φ(x) + (1-p)*ψ(x)
に従う乱数を生成する関数 rmixnorm4() を定義すると以下のようになる.
rmixnorm4 <- function(n, p1, N1, N2) { n1 <- sum(runif(n) < p1) c(rnorm(n1, mean=N1[1], sd=N1[2]), rnorm(n-n1, mean=N2[1], sd=N2[2])) }
関数 hist() でヒストグラムを表示することも出来る.区切り幅は R が適当に選択してくれる.
> data <- generator(1000) > hist(data)
『適当に選択』するが,適切な選択とは云えない.上のグラフを見て真の密度を思い浮かべるのは難しいだろう.Sturgesの方法(1926年の方法!!)を用いている為,データが正規分布(正確には二項分布)から遠ざかれば遠ざかるほど当てはめが悪くなる.そこでまず,パッケージ MASS にある関数 truehist() を用いる.この関数では Scott (1992) が提唱した方法を用いている.
> library(MASS) > x <- generator(1000) > truehist(x)
追加:ユーザが自分で階級数を適切に(!)指定してやれば truehist を使うまでもない。(何が true か?)それに,
The default for breaks is "Sturges": see nclass.Sturges. Other names for which algorithms are supplied are "Scott" and "FD" / "Freedman-Diaconis" (with corresponding functions nclass.scott and nclass.FD). Case is ignored and partial matching is used. Alternatively, a function can be supplied which will compute the intended number of breaks as a function of x.
ということになっているのだが。
> set.seed(777) > x <- generator(1000) > layout(matrix(1:2, 2)) > library(MASS) > truehist(x) > hist(x, nclass=20, freq=FALSE) > layout(1)
次に,パッケージ KernSmooth にある関数 dpih() を用いる.この関数では Wand (1995) が提唱した方法を用いている.
> library(KernSmooth) KernSmooth 2.22 installed Copyright M. P. Wand 1997 > x <- generator(1000) > h <- dpih(x) > bins <- seq(min(x)-0.1, max(x)+0.1+h, by=h) > hist(x, breaks=bins)
ところで,データに対応する x 軸上の点に縦線で印を付ける場合は以下のようにする.
> data(faithful) > attach(faithful) > plot(density(eruptions, bw=0.15)) > rug(eruptions) > rug(jitter(eruptions, amount = .01), side = 3, col = "light blue") > detach("faithful")
ヒストグラムには区切り幅選択の他にもう一つの重大な欠陥がある.それはヒストグラムの端点 (データを描く始点) によってヒストグラムの形状が変わってしまうことである.この現象はデータ数が少ない場合に,より顕著に現れる.
> library(KernSmooth) > x <- generator(100) > h <- dpih(x) > bins <- seq(min(x)-2*h, max(x)+h, by=h) > hist(x, breaks=bins)
上の例において,端点を (区切り幅)/2 だけ右にずらしてみる.
> bins <- seq(min(x)-3*h/2, max(x)+3*h/2, by=h) > hist(x, breaks=bins)
一つの改善案として,いくつかヒストグラムを描いたものの平均をとる方法がある.以下に関数 myASH(データ,平均するヒストグラムの数) を定義する.
# my Average Shifted Histogram myASH <- function(data,breaks) { # Scott (1992) ( ≒ Haedle (1991)) library(KernSmooth) # パッケージの呼び出し h <- dpih(data) # 区切り幅の推定 delta <- h/breaks # ヒストグラムをズラす幅 binnum <- length(seq(min(data)-h-delta/2, max(data)+h+delta/2, by=delta)) # ヒストグラムの棒の数 pts <- rep(0, binnum) # 生データを選別して入れるための変数 binnum <- c() # メモリを開放(されているのかは疑問) for (i in 1:breaks) { # breaks 回くり返し tmpbins <- seq(min(data)-i*delta, max(data)+(breaks-i-1)*delta+h, by=h) # データをそれぞれの区切り幅に選別 tmppts <- hist(data, breaks=tmpbins, prob=T, plot=F)$counts for (j in 0:(length(tmppts)-1)) { # 選別は合計 breaks 回行われ, pts[breaks*j+i] <- tmppts[j+1] # それぞれの結果は breaks だけ } # ズラされて配列に格納される } bins <- seq(min(data)-3*h/2-delta, max(data)+3*h/2+delta, by=h/breaks) pts <- append(rep(0,breaks), pts) # bin : plot の x 座標 pts <- append(pts, rep(0,breaks)) # pts : 両端を 0 で埋める meanpts <- rep(0, length(pts)-breaks+1) for (i in 1:(length(pts)-breaks+1)) { for (j in 1:breaks) { meanpts[i] <- meanpts[i]+ pts[i+j-1] } # meanpts : plot の y 座標 meanpts[i] <- meanpts[i]/breaks # 周辺 breaks 個の平均が plot } # の y 座標となる plot(bins, meanpts, type="l") # 結果を plot する }
> myASH(x,5) # 5 つの histogram の平均を plot
前述のデータについてヒストグラムを描いたが,デコボコすぎてデータの分布があまりよく分からない.離散分布に従っている(と仮定された)母集団からのデータならば問題は無いが,大抵は連続分布に従っている(と仮定された)母集団から取られたデータならば曲線でデータの概観を捉えたいものである.
実は,myASH() の引数 breaks を無限大に飛ばすとカーネル推定という方法を行うことになる.その方法を用いた推定を行う関数 density() を用いることで密度推定することが出来る.この関数は与えられたカーネル関数とバンド幅を用いて密度関数推定を行う.ここでは,上の方で描いた真の密度に上書きすることにする.
> data <- generator(1000) > plot(density(data), xlim=c(-6,6), ylim=c(0,0.3)) > lines(density(data), xlim=c(-6,6), ylim=c(0,0.3)) # ここでは plot と同じ働き >par(new=T) curve(truedensity, xlim=c(-6,6), ylim=c(0,0.3), col=2)
関数 density() の書式は以下の通りである.
density(x, bw = "nrd0", adjust = 1, kernel = c("gaussian", "epanechnikov", "rectangular", "triangular", "biweight", "cosine", "optcosine"), window = kernel, width, give.Rkern = FALSE, n = 512, from, to, cut = 3, na.rm = FALSE)
引数の説明は以下の通り.ただし R のバージョンによって異なる場合があるので help で確認することをお勧めする.
前で定義した generator() で f(x) に従う乱数を 1000 個生成して密度推定した (バンド幅は "SJ" により選択).真の密度を赤,density() で推定した密度を黒としてプロットしている.個人的にはSJで推定するのが一番良い気がするが,滑らかになり 過ぎる嫌いがある.滑らかになり過ぎないようにするにはbcvの方が良いかもしれない.
> data <- generator(1000) > plot(density(data,bw="SJ"), xlab="x", ylab="y", xlim=c(-6,6)) > truedensity <- function (x) { + 0.6/sqrt(2*pi)*exp(-(x+1)^2/2)+ 0.4/sqrt(2*pi)*exp(-(x-2)^2/2) + } > curve(truedensity, xlim=c(-6,6), ylim=c(0,0.3), col=2) > par(new=T) > plot(density(data,bw="SJ"), xlim=c(-6,6), ylim=c(0,0.3), xlab="", ylab="")
パッケージ KernSmooth にある関数 dpik() を用いることでもデータに合わせた推定を行うことが出来る.
> library(KernSmooth) > h <- dpik(data) > est <- bkde(data, bandwidth=h) > plot(est, type="l", xlim=c(-6,6), ylim=c(0,0.3)) > par(new=T) > curve(truedensity, xlim=c(-6,6), ylim=c(0,0.3), col=2)
二変量データをイメージ的なヒストグラムで表すには以下のようにする.
> library(fields) > look<- image.count( precip$x, nrow=32, ncol=32) > image.plot( look)
二変量データの marginal をヒストグラムで表すには以下のようにする.
> library(ade4) > data(rpjdl) > coa1 <- dudi.coa(rpjdl$fau, scannf = FALSE, nf = 4) > s.hist(coa1$li) > s.hist(coa1$li, cgrid = 2, cbr = 3, adj = 0.5, clab = 0)
二変量データのヒストグラム(頻度ポリゴン)を描く場合は以下のようにする.
> library(gregmisc) > # データ:無相関な2変量正規乱数 > x <- rnorm(2000, sd=4) > y <- rnorm(2000, sd=1) > # 遠近法プロット (persp) のためのデータを hist2d() を使用して作成 > h2d <- hist2d(x,y,show=FALSE, same.scale=TRUE, nbins=c(20,30)) > persp( h2d$x, h2d$y, h2d$counts, + ticktype="detailed", theta=60, phi=30, + expand=0.5, shade=0.5, col="cyan", ltheta=-30)
二変量データのカーネル密度推定を行う場合は関数 kde2d() や 二次元データの当てはめ関数 bkde2D() などを用いればよい.
> library(MASS) > data(geyser) > attach(geyser) > f1 <- kde2d(duration, waiting, n = 50, lims = c(0.5, 6, 40, 100)) > f2 <- kde2d(duration, waiting, n = 50, lims = c(0.5, 6, 40, 100), + h = c(width.SJ(duration), width.SJ(waiting)) ) > persp(f2, phi = 30, theta = 20, d = 5)
> library(KernSmooth) KernSmooth 2.22 installed Copyright M. P. Wand 1997 > data(geyser, package="MASS") > x <- cbind(geyser$duration, geyser$waiting) > est <- bkde2D(x, bandwidth=c(0.7,7)) > persp(est$fhat)
R言語は多くの種類のカーネル密度推定法をサポートしている. この手法をわかりやすく解説した書物は以下のようなものがある(上3つは竹澤さんの記事より抜粋)
関数 lowess() によってデータの平滑化を行うことも出来る. lowess() は平滑結果の座標を与える成分 x と y のリストを返す. 平滑化は関数 lines() を用いて元の散布図に追加することが出来る. 書式は以下の通り.
lowess(x, y, f=2/3, iter=3, delta=.01*diff(range(x)))
引数の説明は以下の通り.
以下に使用例を示す.
> data(cars) > plot(cars, main = "lowess(cars)") > lines(lowess(cars), col = 2) > lines(lowess(cars, f = 0.2), col = 3) > legend(5, 120, c(paste("f = ", c("2/3", ".2"))), lty = 1, col = 2:3)