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

// 舟尾
COLOR(red){SIZE(25){ヒストグラムと密度の推定}}

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

#contents
~
-混合分布からのサンプリングコードは鮮やかですね。data3 <- (runif(n) <= 0.6) の部分が最初何をしているのか理解しにくかったので括弧を入れさせてもらいました。論理値ベクトルを生成しているのですね。 --  &new{2004-01-29 (木) 18:00:21};
-修正ありがとうございます。ついでにコメントも入れました。 --  &new{2004-01-30 (金) 20:13:39};
-多(2)変量密度推定に関するパッケージや関数についての情報はありますか？ --  &new{2004-02-23 (月) 00:08:02};
-library(MASS) にある関数 kde2d() を参照してみてください． -- &new{2004-02-23 (月) 20:48:24};
-kde2d と bkde2D について例を挙げてみました． --  &new{2004-03-06 (土) 05:29:44};
-あらら…，２変量のヒストグラムについて書くのを忘れていました・・・ので，あわてて UP しました． --  &new{2004-03-25 (木) 01:11:11};
-関数 generator() を最適化して頂いたので，早速掲載しました．ありがとうございます．&new{2004-07-01 (木) 01:06:38};
-[[Rコードの最適化:http://www.okada.jp.org/RWiki/?R%A5%B3%A1%BC%A5%C9%A4%CE%BA%C7%C5%AC%B2%BD%CE%E3%A1%A7%BA%AE%B9%E7%C0%B5%B5%AC%CD%F0%BF%F4%A4%CE%C8%AF%C0%B8%A5%B3%A1%BC%A5%C9]] より関数 rmixnorm4() の定義を転載しました．この関数定義を御教授下さった方，ありがとうございます．&new{2004-07-04 (日) 00:00:00};

#comment
~

*真の密度とデータの準備 [#yb7a0eb5]

**密度 f(x) = 0.6φ(x)+0.4ψ(x) [#s5aac18f]
φを平均-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)

#ref(ヒストグラムと密度の推定/hist-00.png, center)

**f(x) に従う乱数を生成する関数 [#xc7aa680]
次に， 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)

【[[参考:http://www.okada.jp.org/RWiki/?R%A5%B3%A1%BC%A5%C9%A4%CE%BA%C7%C5%AC%B2%BD%CE%E3%A1%A7%BA%AE%B9%E7%C0%B5%B5%AC%CD%F0%BF%F4%A4%CE%C8%AF%C0%B8%A5%B3%A1%BC%A5%C9]]】φを 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])) 
 }

*ヒストグラムの作成 [#x4565714]

**単純なヒストグラム：Sturges (1926) の方法 [#n3c9cc0e]
関数 hist() でヒストグラムを表示することも出来る．区切り幅は R が適当に選択してくれる．

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

#ref(ヒストグラムと密度の推定/hist-01.png, center)

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

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

#ref(ヒストグラムと密度の推定/hist-03.png, 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)

#ref(falsehist.png)
----

**プラグイン法によるヒストグラム：Wandの方法 [#x16f8b96]
次に，パッケージ ''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)

#ref(ヒストグラムと密度の推定/hist-02.png, 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")

#ref(ヒストグラムと密度の推定/hist04.png, center)

**Average Shifted Histogram [#u14adf44]

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

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

#ref(ヒストグラムと密度の推定/hist05.png, center)

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

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

#ref(ヒストグラムと密度の推定/hist06.png, 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

#ref(ヒストグラムと密度の推定/hist07.png, center)

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

**カーネル推定 [#r04a6f52]
実は，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)

#ref(ヒストグラムと密度の推定/kern-00.png, 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 
で確認することをお勧めする．

-x ： 観測値のベクトル．この観測値の従う確率分布の密度を推定する．
-bw ： 用いられる平滑化バンド幅を指定する．指定された(もしくは既定の) bw の値は以下の引数で与える adjust で adjust 倍される．また，文字列でも指定することが出来，"nrd0", "nrd", "ucv", "bcv", "SJ-ste", "SJ-dpi" が指定できる．詳しくは help(bw.nrd) にて．
--nrd0 ： 平滑化カーネル（ガウスカーネル）の標準偏差になるようにスケール化される．その既定値は標本の大きさ（サンプルサイズ）の 1/5 乗の 1.34 倍で割った標準偏差と四分偏差の小さい方を 0.9 倍したもの（= Silverman (1986) の 『経験則』）でバンド幅を選択する．ただし四分偏差が0の場合は除かれ，このときは bw > 0 が保証される．（Silverman, B. W. (1986) Density Estimation. London: Chapman and Hall.）
--nrd ： nrd0 を Scott (1992) の方法により一般的したものを用いてバンド幅を選択．（Scott, D. W. (1992) Multivariate Density Estimation: Theory, Practice, and Visualization. Wiley.）
--ucv ： バイアス無しのクロスバリデーション規準によりバンド幅を選択．
--bcv ： バイアス付きのクロスバリデーション規準によりバンド幅を選択．
--SJ-ste ： パイロット評価を使用して（方程式を解くことにより）帯域幅を選択する Sheather & Jones (1991)の方法でバンド幅を選択する．（Sheather, S. J. and Jones, M. C. (1991) A reliable data-based bandwidth selection method for kernel density estimation. Journal of the Royal Statistical Society series B, 53, 683-690.）
--SJ-dpi ： パイロット評価を使用して（プラグイン法により）帯域幅を選択する Sheather & Jones (1991)の方法でバンド幅を選択する．
-adjust ： "bw" で指定されたバンド幅は，実際には adjust*bw となる．これは 『デフォルトの半分』 のようなバンド幅を指定する時に有用となる．
-kernel ，window ： 推定に用いるウインドウを指定する文字列で，"gaussian"（デフォルト），"rectangular" ，"triangular" ，"epanechnikov" ，"biweight" ，"cosine" ，"optcosine" が指定できる．これらの文字列入力が面倒ならば window="g" などの様に先頭の一字に略しても良い．
（注）"cosine"は S の命令であるが，これは "optcosine" よりも平滑化に優れる．
-width ： S 言語との互換性の為に存在する引数で，云うなれば "bw" である．width が与えられていて bw が与えられていなければ bw に width がセットされる．それ以外の場合は bw を width のカーネルに依存する倍数にセットする．
-give.Rkern ： これを TRUE にすれば密度推定が行われず，代わりに選ばれたカーネルの canonical bandwidth が返される． プロットをしない場合，もしくは give.Rkern が TRUE の場合の返り値の説明は以下の通り． 
--x ： 密度関数が推定される n 個の点の座標値
--y ： 推定された密度関数値
--bw ： 用いられるバンド幅
--N ： 欠損値を除いた後の標本の大きさ（サンプルサイズ）
--call ： 結果を生み出した呼出し
--data.name ： 引数 x の deparsed name
--has.na ： 互換性のための論理値（常に FALSE ）
-n ： 密度関数の値を求める等間隔点の数．
-from ，to ： 密度を from から to までの区間を n-1 等分した点で求める．
-cut ： 推定された密度を，両端において 0 に下がらせるためのもの．左と右の値はデータの両極端を越えて「切られた」バンド幅となる．
-na.rm ： これを TRUE にすれば欠損値が ｘ から取り除かれる．FALSEならば，欠損値が見つかればエラーを起こす．

前で定義した 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="")

#ref(ヒストグラムと密度の推定/kern-01.png, center)

**プラグイン法による密度推定 [#q0caa129]
パッケージ ''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)

#ref(ヒストグラムと密度の推定/kern-02.png, center)

*二変量データのヒストグラムと密度推定 [#ta85f364]

**二変量データのヒストグラム：イメージ風 [#hc6b12fd]
二変量データをイメージ的なヒストグラムで表すには以下のようにする．

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

#ref(ヒストグラムと密度の推定/hist-image.png, center)

**二変量データのヒストグラム：marginal [#w52316be]
二変量データの 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)

#ref(ヒストグラムと密度の推定/hist-marginal.png, center)

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

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

#ref(ヒストグラムと密度の推定/hist-polygon.png, center)

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

#ref(ヒストグラムと密度の推定/kern-04.png, 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)

#ref(ヒストグラムと密度の推定/kern-03.png, center)

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

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

-Haerdle, W. (1991) : Smoothing Techniques with Implementation in S. （New York: Springer-Verlag.）
-A. Bowman and A. Azzalini(1997): Applied smoothing techniques: The kernel approach with S-Plus illustrations. (Oxford Statistical Science Ser (OSSS), No. 18), Oxford University Press ISBN 0-19-852396-3
-W.N. Venables and B.D. Ripley(2002):Modern applied statistics with S-PLUS, 4th ed. (Statics and Computing), Springer-Verlag.ISBN 0-387-95457-0
-B. W. Silverman(1986): Density estimation for statistics and data analysis. Chapman&Hall, ISBN 0-412-24620-1
-ジェフリー S. シモノフ 著，竹澤邦夫・大森宏 訳(1999) 「平滑化とノンパラメトリック回帰への招待」(農林統計協会，ISBN 4-541-02398-9) 
-竹澤邦夫 著(2003) みんなのためのノンパラメトリック回帰(第２版)上・下 (吉岡書店 ISBN 4-8427-0313-X (上)，4-8427-0314-8 (下))

*lowess() による平滑化 [#g4e374c0]
関数 lowess() によってデータの平滑化を行うことも出来る．
lowess() は平滑結果の座標を与える成分 x と y のリストを返す．
平滑化は関数 lines() を用いて元の散布図に追加することが出来る．
書式は以下の通り．

lowess(x, y, f=2/3, iter=3, delta=.01*diff(range(x)))

引数の説明は以下の通り． 
-x ，y ： 散布図中のプロットの座標を与えるベクトル．単一のプロット構造を指定しても良い．
-f ： 平滑幅．これは各位置での平滑に影響を及ぼすプロットの点の割合を与える．値が大きい程より滑らかになる．
-iter ： 実行されるべき頑健化繰り返しの数．小さな iter の値は lowess の実行を速くする．
-delta ： 互いに距離 delta 以内に位置する x の値は lowess の出力で単一の値に置き換えられる．

以下に使用例を示す．

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

#ref(ヒストグラムと密度の推定/smooth-00.png, center)

タイムスタンプを変更しない