Kernel Regression

author : Kei Ikegami

RDではcut offでの推定値の差が処置による効果とされます。それゆえにRDにおいては推定の精度が求められます。

ここでいう精度は係数の精度というよりも従属変数の予測、推定の精度です。

確かに線形回帰は非常にゆるい仮定のもとでCEFのefficientな推定値を与えてくれます。しかし、その主眼は個々の値ではなくあくまでCEFの近似でした。RDにおいて最も重要なのはcut offにおける一つの値であるので、CEFの矜持をする意味はありません。むしろ一つ一つのデータを重視して説明力を高めることが求められます。

線形回帰の説明力（予測性能）を大きく制限しているのは関数形の仮定です。

ということなので、関数形を仮定しない回帰分析が求められます。今回紹介するのはその一つの手法であるKernel Regressionです。

（日本語が適当ですが許してください）

Nadaraya - Watson Estimator

カーネルについてはwikipediaを参照。

wikipedia)

カーネルの例がたくさんでてますね。

カーネル回帰においてカーネルは重みとして用いられます。すなわち各観測値の情報をどれほど重視するかを決定するものです。

コードをみましょう。

$y_{i} = {\rm sin}\ x_{i} + \epsilon_{i}$で発生させたデータから、カーネル回帰による推定値の一つであるNadaraya - Watson Estimatorを使って関数を推定するコードです。



In [4]:

    
# nadaraya - watson estimator 

using Gadfly
using Distributions
srand(1234)

# データの生成
n = 500
start = 0
stop = 15
x = Array(linspace(start, stop, n))
y = Array(Float64, n)
for i in 1:n
    y[i] = sin(x[i]) + randn()/2
end

# set band width
band = 0.5

# nw est
nw = Array(Float64, n)
for (i,t) in enumerate(x)
    local_pop = x[(t-band) .< x .< (t +band)]
    local_dep = y[(t-band) .< x .< (t +band)]
    weight = pdf(Normal(t), local_pop)
    weighted = weight' * local_dep
    est = weighted / sum(weight)
    nw[i] = est[1]
end

t_layer = layer(x = x, y = nw, Geom.line, Theme(default_color = colorant"black"))

plot(x = x, y = y,t_layer ,
    Geom.point,
    Guide.xlabel("x"), Guide.ylabel("y"), Guide.title("Nadaraya - Watson"),
    Theme(default_point_size = 1.5pt)
)









    Out[4]:

Nadaraya - Watson Estimatorは以下の手順で推定を行います。

band width (=h)を設定します
xから一点を選びます($x_{0}$とします)
$x_{0}$のh近傍に存在するデータをlocal dataとして保存します
local dataに存在するxに対応するyの値に対して以下のように構成した重みをつけて合計し、その和を重みの合計で割る
これによって$x_{0}$における推定値を得ることができる
これを得られている全てのxについて行う

重みの作り方

$x_{0}$を1次モーメントとする分散1の正規分布のPDFを考える

local dataに存在する各xの値に対応する上記のPDFの値を当該xに対応するyの重みとする

ここでカーネルの正体が明らかになります。

すなわち、カーネルは重みづけの際に使用される分布みたいなものです。

今回用いているものはガウスカーネルと呼ばれるもので、正規分布の形を思い出すとわかるように基準にしている点に近いデータを重視して平均をとることができます。

そのほかにも例えばrectangular kernelなどでは、周辺のデータに等しく重みをつけることができます。

とはいえ基本はガウスカーネルらしいですので、正直カーネルは重要な点ではありません。

上の Nadaraya - Watson Estimator においてカーネルを差し置いて最も重要なのは実はband width selectionなのです。

上の重み付けの説明からもわかると思いますが、カーネル回帰に置けるband widthはある一点の推定を行う際にどこまでのデータを使用するかを決めるものです。

ゆえに当然band widthの設定によってカーネル回帰の結果は大きく変わります。

下でその様子を見てみましょう。



In [5]:

    
# set band widths
bands = Float64[0.25, 0.5, 1.0, 5.0,10.0]
colors = ["black","blue","red", "green","purple"]

p = plot(x = x, y = y,
    Geom.point,
    Guide.xlabel("x"), Guide.ylabel("y"), Guide.title("Nadaraya - Watson"),
    Theme(default_point_size = 1.5pt))

for (s,band) in enumerate(bands)
    nw = Array(Float64, n)
    for (i,t) in enumerate(x)
        local_pop = x[(t-band) .< x .< (t +band)]
        local_dep = y[(t-band) .< x .< (t +band)]
        weight = pdf(Normal(t), local_pop)
        weighted = weight' * local_dep
        est = weighted / sum(weight)
        nw[i] = est[1]
    end
    color = colors[s]
    t_layer = layer(x = x, y = nw,Geom.line, Theme(default_color = colorant"black"))
    push!(p, t_layer)
end

p









    Out[5]:

これではどのbandを使えばいいのかわかりません。

しかしRDと同様に機械的にband widthを設定することができます。

やり方は4つぐらいあるようですが、今回は平均2乗誤差を最小化する(Least squares cross validation)を基準として採用します。

（すいませんこれ少し間違ってるかもしれないです）



In [6]:

    
# band width selection

using Gadfly
using Distributions
srand(1234)

# データの生成
n = 500
start = 0
stop = 15
x = Array(linspace(start, stop, n))
y = Array(Float64, n)
for i in 1:n
    y[i] = sin(x[i]) + randn()/2
end

bands = Float64[0.25, 0.5, 1.0, 5.0,10.0]
square = Array(Float64, length(bands))
nws = Dict()

for (s,band) in enumerate(bands)
    nw = Array(Float64, n)
    for (i,t) in enumerate(x)
        local_pop = x[(t-band) .< x .< (t +band)]
        local_dep = y[(t-band) .< x .< (t +band)]
        weight = pdf(Normal(t), local_pop)
        weighted = weight' * local_dep
        est = weighted / sum(weight)
        nw[i] = est[1]
    end
    nws[s] = nw
    square[s] = ((nw - y)' * (nw - y))[1]
end

band = bands[indmin(square)]

println("selected band width is $band")
println("when band width is set to $band, Nadaraya - Watson Estimator is shown as follows")

t_layer = layer(x = x, y = nws[indmin(square)], Geom.line, Theme(default_color = colorant"black"))

p = plot(x = x, y = y, t_layer,
    Geom.point,
    Guide.xlabel("x"), Guide.ylabel("y"), Guide.title("Nadaraya - Watson"),
    Theme(default_point_size = 1.5pt))









    



selected band width is 0.25
when band width is set to 0.25, Nadaraya - Watson Estimator is shown as follows






    Out[6]:

今回はband として0.25が選択されました。

Nadaraya - Watson estimatorはlocal constant estimatorとも呼ばれます。

その他の手法

カーネル回帰ではほかにもlocal linear estimationやlocal polynomial estimationなどの手法も用いられます。

まずはlocal linear regression (LL) です

LLはRDでもやったlocal linear regressionです。ただし、回帰においては以下の手順で推定を行います。

band width (=h)を設定します
xから一点を選びます($x_{0}$とします)
$x_{0}$のh近傍に存在するデータをlocal dataとして保存します
local dataのみで、残差の二乗にカーネルで重みづけをして、その和を最小にするような線形方程式の係数を推定します
推定した係数をもとに$x_{0}$での推定値を算出し、それをlocal linear estimateとします
以上の作業を全てのxについて行います。

イメージとしては、もとの関数を$x = x_{0}$でテイラー展開した際の係数を推定していると考えるといいと思います。

local poynomial regressionは、テイラー展開を2次近似以上までやって、その係数を推定するものです。

こっちは実装するのが大変そうなので言葉での説明で勘弁してください。時間があったら実装します。