ヒストグラムと密度の推定

(グラフィックス参考実例集に戻る。Rのグラフィックスパラメータを参照する。)




真の密度とデータの準備

密度 f(x) = 0.6φ(x)+0.4ψ(x)

φを平均-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)
 center

f(x) に従う乱数を生成する関数

次に, 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])) 
}  

ヒストグラムの作成

単純なヒストグラム:Sturges (1926) の方法

関数 hist() でヒストグラムを表示することも出来る.区切り幅は R が適当に選択してくれる.

> data <- generator(1000)
> hist(data)  
 center

『適当に選択』するが,適切な選択とは云えない.上のグラフを見て真の密度を思い浮かべるのは難しいだろう.Sturgesの方法(1926年の方法!!)を用いている為,データが正規分布(正確には二項分布)から遠ざかれば遠ざかるほど当てはめが悪くなる.そこでまず,パッケージ MASS にある関数 truehist() を用いる.この関数では Scott (1992) が提唱した方法を用いている.

> library(MASS)
> x <- generator(1000)
> truehist(x)
 center

追加:ユーザが自分で階級数を適切に(!)指定してやれば 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)
falsehist.png

プラグイン法によるヒストグラム:Wandの方法

次に,パッケージ 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)
 center

ところで,データに対応する 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")
 center

Average Shifted Histogram

ヒストグラムには区切り幅選択の他にもう一つの重大な欠陥がある.それはヒストグラムの端点 (データを描く始点) によってヒストグラムの形状が変わってしまうことである.この現象はデータ数が少ない場合に,より顕著に現れる.

> library(KernSmooth)
> x <- generator(100)
> h <- dpih(x)
> bins <- seq(min(x)-2*h, max(x)+h, by=h)
> hist(x, breaks=bins)
 center

上の例において,端点を (区切り幅)/2 だけ右にずらしてみる.

> bins <- seq(min(x)-3*h/2, max(x)+3*h/2, by=h)
> hist(x, breaks=bins)
 center

一つの改善案として,いくつかヒストグラムを描いたものの平均をとる方法がある.以下に関数 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
 center

密度の推定

前述のデータについてヒストグラムを描いたが,デコボコすぎてデータの分布があまりよく分からない.離散分布に従っている(と仮定された)母集団からのデータならば問題は無いが,大抵は連続分布に従っている(と仮定された)母集団から取られたデータならば曲線でデータの概観を捉えたいものである.

カーネル推定

実は,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)
 center

関数 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="")
 center

プラグイン法による密度推定

パッケージ 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)
 center

二変量データのヒストグラムと密度推定

二変量データのヒストグラム:イメージ風

二変量データをイメージ的なヒストグラムで表すには以下のようにする.

> library(fields)
> look<- image.count( precip$x, nrow=32, ncol=32)
> image.plot( look)
 center

二変量データのヒストグラム:marginal

二変量データの 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)
 center

二変量データのヒストグラム:頻度ポリゴン

二変量データのヒストグラム(頻度ポリゴン)を描く場合は以下のようにする.

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

二変量データのカーネル密度推定

二変量データのカーネル密度推定を行う場合は関数 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)
 center
> 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)
 center

カーネル密度推定のわかりやすい文献(なんでも掲示板2003-10-03の投稿より)

R言語は多くの種類のカーネル密度推定法をサポートしている. この手法をわかりやすく解説した書物は以下のようなものがある(上3つは竹澤さんの記事より抜粋)

lowess() による平滑化

関数 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)  
 center

添付ファイル: filehist04.png 3517件 [詳細] filesmooth-00.png 3655件 [詳細] filehist-image.png 4178件 [詳細] filehist-polygon.png 3743件 [詳細] filekern-00.png 3736件 [詳細] filehist06.png 3542件 [詳細] filehist-02.png 3947件 [詳細] filehist-01.png 4436件 [詳細] filehist07.png 3467件 [詳細] filekern-01.png 3443件 [詳細] filekern-04.png 4277件 [詳細] filekern-03.png 3634件 [詳細] filekern-02.png 3698件 [詳細] filehist-marginal.png 3947件 [詳細] filehist05.png 3522件 [詳細] filefalsehist.png 3116件 [詳細] filehist-00.png 3845件 [詳細] filehist-03.png 4832件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2023-03-25 (土) 11:19:17