2013-03-03

どのカイ分布・カイ二乗分布を選ぶか

決断 Decision Theory カイ分布カイ二乗分布超幾何関数ベッセル関数距離 R

Rのdensity()関数の中身をedit(density.default)として表示させると、ガウシアンでのカーネル密度推定では、bwというパラメタを「バンド幅」としていることがわかる
さらにこのbwは個々の観測値についてそれを平均とする正規分布を作って、全観測値について足し合わせるときの、その単位正規分布の標準偏差として使われるものである
このbwは結局のところどうやっているか、というとbw.nrd0(x)という処理で得ている。このxは密度関数を推定しようとしているその数値列である
さらにこのbw.nrd0()関数は…というと

> bw.nrd0
function (x) 
{
    if (length(x) < 2L) 
        stop("need at least 2 data points")
    hi <- sd(x)
    if (!(lo <- min(hi, IQR(x)/1.34))) 
        (lo <- hi) || (lo <- abs(x[1L])) || (lo <- 1)
    0.9 * lo * length(x)^(-0.2)
}
<bytecode: 0x000000000aba7178>
<environment: namespace:stats>

となっており、また、help(bw.nrd0)でもわかるとおり

bw.nrd0 implements a rule-of-thumb for choosing the bandwidth of a Gaussian kernel density estimator. It defaults to 0.9 times the minimum of the standard deviation and the interquartile range divided by 1.34 times the sample size to the negative one-fifth power (= Silverman's ‘rule of thumb’, Silverman (1986, page 48, eqn (3.31)) unless the quartiles coincide when a positive result will be guaranteed.

とあってSilverman's role of thumbによって(元の分布に正規分布を想定した(？）上でのリスク関数の最小化をもたらす幅を「簡易的に求める式」で求めている。ちなみに、このリスク関数は、推定分布と真の分布のずれの二乗を全体にわたって積分したものをを最小にしようとしたもので、Wikiの記事(こちら)でのmean integrated squared error(MISE)であり、それを漸近的に求めるAMISE(Asymptotic MISE)の「簡易計算」になっている
一応、Rのdensity()関数がそうやっていることを、べたに確認しよう

d <- c(runif(10),runif(20)*3,runif(5)*20,rexp(15))

dd <- density(d)

bw <- bw.nrd0(d)

x <- seq(from=min(dd$x),to=max(dd$x),length=500)

y <- rep(0,length(x))
for(i in 1:length(d)){
	y <- y + dnorm(x,d[i],bw)/length(d)
}

xlim <- range(x)
ylim <- c(0,max(y))
plot(density(d),xlim=xlim,ylim=ylim)
par(new=TRUE)
plot(x,y,type="l",xlim=xlim,ylim=ylim,col=2)

じゃあ、「幅」は「漸近的な解」を「簡易的に求め」ているんだから、カイ分布・カイ二乗分布のときもそれを使わせてもらうことにしちゃってもそんなに悪くないのでは…(分布が正規分布でないときにも正規分布を仮定しているようだし…)
さらに、観測値ごとの小さい分布を作るときに、その小分布の標準偏差をそのbwにしているのなら、カイ分布・カイ二乗分布でもそうしてやるというのは、一つのやり方として通りそうだ
- ちなみに、カイ分布の分散は $df+m^2-\mu^2$ (mは形を決めるパラメタ、 $\mu$ は平均。カイ二乗分布は $2(df+2M)$ (Mは形を決めるパラメタ)
さて、やってみよう。bwはバンド幅(であって、それを分布の標準偏差にしようとしている)
カイ分布の２つのパラメタ(dfは自由度、mは非心性のパラメタ)
- このほかに、(平均と分散)を変数で表しておこう
  - ここで $bw=sd(x)$ :xは観測値のセット
- 関係は
  - $df+m^2-\mu^2=v=bw^2$
  - $\sqrt{\frac{\pi}{2}} L_{\frac{1}{2}}^{(\frac{df}{2}-1)}(-\frac{m^2}{2})=\mu$
- $\mu,bw$ が与えられ、これを満足する $df,m$ の２変数について、解く、と
カイ二乗分布の２つのパラメタ(dfは自由度、Mは非心性のパラメタ)
- 同様に、(平均と分散)を変数で表しておこう
  - ここで $bw=sd(x)$ :xは観測値のセット
- 関係は
  - $2(df+2M)=bw^2$
  - $df+m = \mu$
- が与えられ、これを満足するの２変数について、解く、と
  - こちらは平均が大きいときに分散を小さくできなそうなので、このままではだめだろう…(一定以上に大きかったら、負の領域の存在は無視して正規分布にする…とかで対処する…)