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. メタブレーン) という本に詳しく載っています.
ちょっと時間がないので,これは別の機会に.