2011-08-20

度数カウントもどき、高速フーリエ変換、畳み込み

２次元を例に
- 分布(たとえば正規分布)の重なりとして推定する
Rの関数bkde2D()の中身を追いかけよう
２次元プロットデータを作る

Npt<-300
x<-y<-rep(0,Npt)
r<-10
for(i in 1:Npt){
	t<-runif(1)*2*pi/1
	r2<-r*(sin(t*200)+10)+rnorm(1)*0.1
	x[i]<-r2*cos(t)
	y[i]<-r2*sin(t)
}
x<-cbind(x,y)

bkde2d() の中身を確認する

# 引数(実行条件)
gridsize = c(51L, 51L) # 推定値を計算する格子点の数
truncate = TRUE # 周辺部の取扱い条件
bandwidth=c(10,10) # 平滑化のときの「隣」とみなす範囲
# 処理の中身  
    n <- nrow(x) # 点の数
    M <- gridsize # 格子点数
    h <- bandwidth # 「隣」の幅
    tau <- 3.4 # どうやって決められたか(まだわからない)定数
    if (length(h) == 1L) 
        h <- c(h, h)

推定範囲の設定
- 引数指定をしなければ、各軸の最小・最大値範囲から「平滑幅」の1.5倍

    #if (missing(range.x)) {
        range.x <- list(0, 0)
        for (id in (1L:2L)) range.x[[id]] <- c(min(x[, id]) - 
            1.5 * h[id], max(x[, id]) + 1.5 * h[id])
    #}

推定範囲・推定格子点の座標指定
- 長方形の２点の座標を指定し
- 格子点の２軸の座標値のベクトルを作成

    a <- c(range.x[[1L]][1L], range.x[[2L]][1L])
    b <- c(range.x[[1L]][2L], range.x[[2L]][2L])
    gpoints1 <- seq(a[1L], b[1L], length = M[1L])
    gpoints2 <- seq(a[2L], b[2L], length = M[2L])

観察データから、格子点に値を与える
- この関数の中身はこちらで
- １次元ならは、格子が座標x1,x2にあって、観察が、 $x1 \le x \le x2$ に１サンプルなされたときには、x1の格子に $\frac{x2-x}{x2-x1}$ 、x2の格子に $\frac{x-x1}{x2-x1}$ の観察がなされたように配分する
- ２次元では格子が(x1,y1),(x1,y2),(x2,y1),(x2,y2)とあり、観察が $x1 \le x \le x2$ , $y1 \le y \le y2$ となったときには、４点にそれぞれ $\frac{x2-x}{x2-x1}\times \frac{y2-y}{y2-y1}$ 、 $\frac{x2-x}{x2-x1}\times \frac{y-y1}{y2-y1}$ 、 $\frac{x-x1}{x2-x1}\times \frac{y2-y}{y2-y1}$ 、 $\frac{x-x1}{x2-x1}\times \frac{y-y1}{y2-y1}$ のように配分する

    gcounts <- linbin2D(x, gpoints1, gpoints2)

平滑化フィルタリング行列の作成
- 行列要素の算出は、行と列とで別個に行い
- その積で行列のセルの値を出している
- 行の値、列の値の計算にdnorm()を使っているが、これが「正規分布」で平滑にする、ということに対応する
- フィルタリング行列(フィルタリング用の処理情報〜カーネル)を「掛ける」
  - 元のデータの凸凹した分布に、フィルタリング用の分布を掛けることは、「関数」ｘ「関数」で、これが畳み込み
  - 畳み込むためだけならフーリエ変換しなくたってよい(こちら)
  - 畳み込みの親切な説明記事はこちら

    L <- numeric(2L)
    kapid <- list(0, 0)
    for (id in 1L:2L) {
# 「平滑幅」がいくつの格子幅を持つかを整数で算出して、平滑化のための分布関数の確率計算点を作る(「平滑化係数計算幅」)
        L[id] <- min(floor(tau * h[id] * (M[id] - 1)/(b[id] - 
            a[id])), M[id] - 1L)
# 「中心」から「平滑化係数計算幅」の番号付け
        lvecid <- 0:L[id]
# 平滑化係数計算幅
        facid <- (b[id] - a[id])/(h[id] * (M[id] - 1L))
# 標準正規分布の確率密度を中心から離散的に計算する
        z <- matrix(dnorm(lvecid * facid)/h[id])
# 軸に沿って、正負の２方向に離散的に「平滑化係数」計算点のそれを計算するので、その分を考慮している。「距離０」の点は１点でそれ以外は２点。「積分」なので、「平滑化係数計算幅」を掛けている
# それを足し合わせたものは、「分布の積分」に近似できるから、それをtotとしたうえで、
        tot <- sum(c(z, rev(z[-1L]))) * facid * h[id]
# 全体が１になるようにtotで割っている
        kapid[[id]] <- z/tot
    }
# ２軸について、そのような係数ベクトルを作って、２次元化
    kapp <- kapid[[1L]] %*% (t(kapid[[2L]]))/n

畳み込みを考えるにあたって、ある程度の広さを確保。フーリエ変換なので、無限の領域(に周期関数がのたくっている)仮定。その一部を切り取り
- 周期的なので「循環畳み込み」→こちら

    if (min(L) == 0) 
        warning("Binning grid too coarse for current (small) bandwidth: consider increasing 'gridsize'")
    P <- 2^(ceiling(log(M + L)/log(2)))
    L1 <- L[1L]
    L2 <- L[2L]
    M1 <- M[1L]
    M2 <- M[2L]
    P1 <- P[1L]
    P2 <- P[2L]

rpは平滑化係数でできたカーネル
- 正方形の４隅に「ピーク」が来て、減衰を考慮した範囲は値があるが、中心付近は全部０

    rp <- matrix(0, P1, P2)
    rp[1L:(L1 + 1), 1L:(L2 + 1)] <- kapp
    if (L1) 
        rp[(P1 - L1 + 1):P1, 1L:(L2 + 1)] <- kapp[(L1 + 1):2, 
            1L:(L2 + 1)]
    if (L2) 
        rp[, (P2 - L2 + 1):P2] <- rp[, (L2 + 1):2]

観察データから「格子点に観察されたものとみなした値」の行列として作ったgcountsを、フーリエ変換用の広さに準備する

    sp <- matrix(0, P1, P2)
    sp[1L:M1, 1L:M2] <- gcounts

畳み込むべき２つの関数(に対応する行列)ができたが、これを真面目に掛け算するよりも、もっといい方法がある
- その方法をとってもよいことを表したのが「畳み込み定理」(こちら)
  - フーリエ変換して、スペクトルにして、スペクトル同士をかけて、逆フーリエ変換で戻してやっても、地道に(真面目に)掛け算したのと同じだよ、という定理

    rp2 <- fft(rp)
    sp2 <- fft(sp)
    rp3 <- Re(fft(rp2 * sp2, inverse = TRUE)/(P1 * P2))[1L:M1, 1L:M2]
# 0以下の値が入ってしまったセルを0にする処理
    rp4 <- rp3 * matrix(as.numeric(rp3 > 0), nrow(rp3), ncol(rp3))
    list(x1 = gpoints1, x2 = gpoints2, fhat = rp4)