61. データの分布とヒストグラム・密度推定


データの分布

まず,データの要素をグループ化するには,順序無し因子を生成する関数 factor() を使い,関数 table() を用いることで度数分布表を作ることが出来る.

 
 data(warpbreaks)          # 1台の織機当たりの反り破損の数
 attach(warpbreaks)        # 使うデータを固定
 fc <- factor(tension)     # 要素をグループ化
 levels(fc)                # グループ化されているかを確認

[1] "L" "M" "H"

 table(fc)                 # グループ化したデータ fc に関する度数分布表

fc
 L  M  H 
18 18 18 

他にも関数 xtabs() でクロス集計が出来る.

関数 cut() と組み合わせると,以下の様なことも出来る.

 
 tb <- factor(cut(breaks, breaks = 10+10*(0:6)))
 table(tb, fc)

         fc
tb        L M H
  (10,20] 3 6 8
  (20,30] 7 7 7
  .............
  (60,70] 2 0 0

また,データの平均や分散を計算する場合は関数 tapply() を用いる.tapply() の書式と例は以下の通り.

 
 # tapply(数値ベクトル,数値ベクトルと同じ長さの文字列ベクトル,適用したい関数) 
 ( mymeans <- tapply(breaks, fc, mean) )     # breaks の平均をグループごとに計算

       L        M        H 
36.38889 26.38889 21.66667 

 ( mysd <- tapply(breaks, fc, sd) )          # グループごとに計算することも出来る

        L         M         H 
16.446487  9.121009  8.352527 

ところで,データの集合の要約は関数 summary() と fivenum() で,度数表示は stem() で知ることが出来る.stem() は,stem(データ, x) で幹の長さが既定値(=1)の x 倍になる.

 
 stem(breaks)

  The decimal point is 1 digit(s) to the right of the |
  1 | 0234555667788899
  2 | 001111445666678889999
  3 | 00156699
  4 | 1234
  5 | 124
  6 | 7
  7 | 0

 summary(breaks)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  10.00   18.25   26.00   28.15   34.00   70.00 

 fivenum(breaks)  # 最小値・下側ヒンジ・中央値・上側ヒンジ・最大値

[1] 10 18 26 35 70 

しかし,1 次元データや 2 次元データならば視覚的に見た方がデータの特徴をつかみやすい場合が多い.そこでヒストグラムと密度推定の方法を以下で述べることにする.

ヒストグラムによる密度推定(準備)

φを平均 -1,分散 1 の正規分布に従う確率変数,ψを平均 2 ,分散 1 の正規分布に従う確率変数として

となる密度から得られる乱数を生成する関数 generator() を定義する.

 
 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 <- 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を採択するための手続
 }

関数 hist() でヒストグラムを表示することも出来る.区切り幅は R が適当に選択してくれる.probability=T と指定することによって相対度数で作図するように変更することが出来る.

 
 x <- generator(1000)
 hist(x)

ヒストグラムによる密度推定

区切り幅は『適当に選択される』が『適切に選択される』わけではない.というのも,hist() のデフォルトは 『データの範囲を log2n + 1 ( n はデータの個数) 個の階級に分割して各階級に属するデータの数を棒グラフとして作図する』という Sturges (1926!!) の方法を用いているため,まず平滑化をし過ぎる嫌いがあり,さらにデータが正規分布 (正確には二項分布) から遠ざかれば遠ざかるほど当てはめが悪くなる.

そこでパッケージ MASS にある関数 truehist() (この関数では Scott (1992) が提唱した方法を用いている) や,パッケージ KernSmooth にある関数 dpih() (この関数では Wand (1995) が提唱した方法を用いている) を用いることで,より正確なヒストグラムを描くことが出来る.

 
 library(MASS)
 x <- generator(1000)                         # hist(x, breaks="Scott") でも可
 truehist(x)                                  # hist() の breaks には他に Sturges と FD を指定することが出来る

 library(KernSmooth)

KernSmooth 2.22 installed
Copyright M. P. Wand 1997

 x <- generator(1000)
 h <- dpih(x)                                 # bin には区切り幅の点を表す x 座標を指定
 bins <- seq(min(x)-0.1, max(x)+0.1+h, by=h)  # bin が等差数列でなければ区切り幅の横幅
 hist(x, breaks=bins)                         # もバラバラな長さになる

 

 

ただし,ヒストグラムは左端の端点の位置によって形が変わってしまうという欠点がある (データ数が 50 未満の場合) .

 
 library(KernSmooth)
 x <- generator(100)
 h <- dpih(x)
 bins <- seq(min(x)-2*h, max(x)+h, by=h)
 hist(x, breaks=bins)                             # 左下の図
 bins <- seq(min(x)-3*h/2, max(x)+3*h/2, by=h)
 hist(x, breaks=bins)                             # 右下の図

 

 

 

カーネルによる推定

前述のデータについて,関数 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)

データに対応する x 軸上の点に縦線で印を付ける場合は以下のようにする.

 
 plot(density(x))
 rug(x)
 rug(jitter(x, amount = .1), side = 3, col = "light blue")

関数 density() の引数 bw には以下を指定することが出来る.

引数 bw

機能

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)の方法でバンド幅を選択する.

bw 以外に以下の引数が指定できる.

引数

機能

adjust

"bw" で指定されたバンド幅は,実際には adjust*bw となる.これは 『デフォルトの半分』 のようなバンド幅を指定する時に有用となる.

kernel ,window

推定に用いるウインドウを指定する文字列で,"gaussian"(デフォルト),"rectangular" ,"triangular" ,"epanechnikov" ,"biweight" ,"cosine" ,"optcosine" が指定できる.個人的には汎用性のある "epanechnikov" がお勧め.これらの文字列入力が面倒ならば window="g" などの様に先頭の一字に略しても良い.
(注)"cosine"は S の命令であるが,これは "optcosine" よりも平滑化に優れる. 

width

S 言語との互換性の為に存在する引数で,云うなれば "bw" である.width が与えられていて bw が与えられていなければ bw に width がセットされる.

give.Rkern

これを TRUE にすれば密度推定が行われず,代わりに選ばれたカーネルの canonical bandwidth が返される.

n

密度関数の値を求める等間隔点の数.

from ,to

密度を from から to までの区間を n-1 等分した点で求める.

cut

推定された密度を,両端において 0 に下がらせるためのもの.左と右の値はデータの両極端を越えて「切られた」バンド幅となる.

na.rm

これを TRUE にすれば欠損値が x から取り除かれる.FALSEならば,欠損値が見つかればエラーを起こす.

前で定義した generator() で f(x) に従う乱数を 1000 個生成して密度推定した (バンド幅は "SJ" により選択) .真の密度を赤,density() で推定した密度を黒としてプロットしている

 
 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() ならば,もっとデータに合わせた推定を行う.他にも関数 bkde() などが用意されている.

 
 library(KernSmooth)
 data <- generator(1000);  h <- dpik(data)
 est <- bkde(data, bandwidth=h)
 plot(est, type="l", xlim=c(-6,6), ylim=c(0,0.3))

二次元データの密度推定

二次元データの当てはめを行う場合は,頻度ポリゴンを描く関数 hist2d() (パッケージ:gregmisc) ,カーネル密度推定を行う関数 bkde2D()  (パッケージ:KernSmooth) ,関数 kde2d()  (パッケージ:MASS) などが用意されている.

 
  library(gregmisc)

Attaching package 'gregmisc':
        The following object(s) are masked from package:base :
        lowess 

  x <- rnorm(2000, sd=4)           # データ:無相関な2変量正規乱数
  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)

  library(MASS)                    # 関数 width.SJ() を使うので MASS を呼び出し
  library(KernSmooth)              # 関数 bkde2D() を使うので KernSmooth を呼び出し
  x <- rnorm(2000, sd=4)           # データの準備
  y <- rnorm(2000, sd=1)           # データの準備
  f1 <- bkde2D(cbind(x, y), bandwidth=c(width.SJ(x), width.SJ(y)))
  persp(f1$fhat, phi = 30, theta = 20, d = 5)

 

 

lowess() による平滑化

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

他にも Nadaraya-Watson 推定による核関数を用いた回帰平滑化を行う関数 ksmooth() や,3次の平滑化スプライン関数を当てはめる関数 smooth.spline(),Friedman の SuperSmoother を用いて平滑化を行う関数 supsmu() や,移動中央値平滑化を行う関数 runmed() ,smooth() ,smoothEnds() などがある.