Rで密度推定

ヒストグラムカーネル密度推定についてすごく簡単にまとめました.Kashiwa.R#4で発表する内容を文章にしたものです.


使用するデータセット

まずここでは,faithful データセットの,eruption を例に用います.
イエローストーン国立公園にあるOld Faithful間欠泉の噴出時間 (分) だそうです.

# データセットの確認
head(faithful)
# 噴出時間データのみを抜き出す
fe <- faithful$eruption


まず度数分布を見てみると以下のようになります.

stem(fe)
  The decimal point is 1 digit(s) to the left of the |

  16 | 070355555588
  18 | 000022233333335577777777888822335777888
  20 | 00002223378800035778
  22 | 0002335578023578
  24 | 00228
  26 | 23
  28 | 080
  30 | 7
  32 | 2337
  34 | 250077
  36 | 0000823577
  38 | 2333335582225577
  40 | 0000003357788888002233555577778
  42 | 03335555778800233333555577778
  44 | 02222335557780000000023333357778888
  46 | 0000233357700000023578
  48 | 00000022335800333
  50 | 0370

>


ヒストグラム

これをヒストグラム (相対度数) で表すのですが,区切り幅や端点の選択によって,形が大きく変化してしまいます.

# Sturges (1926) の方法にしたがった区切り幅 (デフォルト)
hist(fe, freq = FALSE)

# Scott (1992) の方法にしたがった区切り幅
# KernSmoothパッケージのインストールが必要です
library(KernSmooth)
fe.h <- dpih(fe)
fe.bins <- seq(min(fe), max(fe) + fe.h, by = fe.h)
hist(fe, freq = FALSE, breaks = fe.bins)

# 区切り幅の1/2だけ端点をずらしてみる
fe.bins.2 <- seq(min(fe) - fe.h / 2, max(fe) + fe.h * 3/2, by = fe.h)
hist(fe, freq = FALSE, breaks = fe.bins.2)


困りました…
なにが問題かというと,階級内に含まれるデータの分布を考慮していない点です.
たとえば,c(0.2, 0.6, 0.8) という3つの標本があったとき,示した図はどれも,データを表すヒストグラムになり得てしまうわけです.

data <- c(0.2, 0.7, 0.9)
oldpar <- par(mfrow = c(1, 3))
hist(data, breaks = seq(0, 1, 1/2))
points(data, rep(-0, 3), pch = 16)

hist(data, breaks = seq(0, 1, 1/3))
points(data, rep(-0, 3), pch = 16)

hist(data, breaks = seq(0, 1, 1/6))
points(data, rep(-0, 3), pch = 16)


カーネル密度推定

カーネル密度推定では,この問題を回避して,境界の設定に関係なく分布を推定できるようにするための方法です,ひとつひとつの標本に,カーネル関数から決定された分布を与えて,それを積み上げて全体の分布とするようなイメージです.


カーネル密度推定 > 定義 - Wikipedia
カーネル密度推定 > 直感的説明 - Wikipedia


Rでは,density()関数を使えば,簡単にカーネル密度推定ができます.

fe.d <- density(fe)
hist(fe, freq = FALSE, breaks = fe.bins)
lines(fe.d)


ただし,個々の標本に与える分布の形 (カーネル関数) とその裾野 (バンド幅) をどう設定するかによって,推定の結果が変わってきます.デフォルトでは,カーネル関数に正規分布確率密度関数,バンド幅に bw.nrd0 関数で計算された値が使用されますが,たとえば以下のように,バンド幅を変えると,滑らかさ (平滑化の程度) が変化します.

fe.d <- density(fe)
fe.d.1 <- density(fe, bw = 0.1)
fe.d.2 <- density(fe, bw = 1)
hist(fe, freq = FALSE, breaks = fe.bins)
lines(fe.d)
lines(fe.d.1, col = "red")
lines(fe.d.2, col = "blue")

と,このようなわけで,絶対的な解はないので気をつけましょうね,というところで終わってしまいます.はい.


ヒストグラムの平滑化

まずヒストグラムを作成して,それをノンパラメトリック回帰を利用して平滑化する,という方法もあるそうです.

Rによるノンパラメトリック回帰の入門講義』 (竹澤邦夫. 2009. メタブレーン) という本に詳しく載っています.

ちょっと時間がないので,これは別の機会に.