データ解析
第十回「ノンパラメトリック推定法（１）」鈴木 大慈
理学部情報科学科
西八号館 W707 号室
s-taiji@is.titech.ac.jp

1 / 36休講情報

6/24 は休講

2 / 36

今日の講義内容
カーネル密度推定
スプライン回帰3 / 36

構成

1

カーネル密度推定2

ノンパラメトリック回帰：B-スプライン

4 / 36

カーネル密度推定の目的

ある分布の密度関数を推定したい．
いままではパラメトリックモデルの上での推定手法を紹介してきた．
(最尤推定，ベイズ推定)
もし分布を適切なパラメトリックモデルで記述できなかったら？
→ ノンパラメトリック推定
カーネル密度推定はその代表的な方法．

5 / 36

カーネル密度推定で何が得られる？

0.3
0.2
0.1
0.0p(x)

0.4

0.5

0.6

Kernel Density Estimation (gauss)−3−2−10123x6 / 36

カーネル密度推定量{Xi }ni=1 : データ (一次元とする)
カーネル密度推定量: あるカーネル関数 K を用いて，
1∑p̂(x) =K
nni=112(x − Xih).

h > 0 のことをバンド幅と呼ぶ．適切に選択する必要がある．
K は次の性質を満たすものとする:∫∫∫
K (x)dx = 1,xK (x)dx = 0,x 2 K (x) > 0.

7 / 36

−101x2
348 / 36

0.0

0.1

0.2

p(x)
0.3

0.40.5

カーネルの種類
次のようなカーネル関数がよく用いられる.
1

Gaussian:

2

Rectangular:( 2)
1xK (x) = √ exp −.
22π{K (x) =3Triangular:

4

Epanechnikov:

120(|x| ≤ 1),(otherwise).

{|x| (|x| ≤ 1),
K (x) =0(otherwise).
{K (x) =324 (1 − x )0(|x| ≤ 1),(otherwise).

9 / 36

0.8

1.0

カーネルの種類

kernel functions

0.0

0.2

0.4K(x)

0.6

Rectangular
Triangular
Gaussian
Epanechnikov−3−2−10123x
10 / 36

カーネル密度推定量の様子

0.3
0.2
0.1
0.0p(x)

0.4

0.5

0.6

Kernel Density Estimation−3−2−10123x
11 / 36

カーネル密度推定量の様子

0.3
0.2
0.1
0.0p(x)

0.4

0.5

0.6

Kernel Density Estimation−3−2−10123x
11 / 36

カーネル密度推定量の様子

0.3
0.2
0.1
0.0p(x)

0.4

0.5

0.6

Kernel Density Estimation−3−2−10123x11 / 36

カーネル関数の種類と推定結果
n = 100

−3−2

−1

01x

rectangular

2

3

0.6
0.5
0.0

0.1

0.2

0.3

p(x)

0.4

0.5
0.0

0.1

0.2

0.3

p(x)

0.4

0.5
0.4
0.3
0.2
0.1
0.0p(x)

Kernel Density Estimation (gauss)

0.6

Kernel Density Estimation (triangular)

0.6

Kernel Density Estimation (rectangular)−3−2−101xtriangular23−3−2−10123xgauss

12 / 36

カーネル関数の種類と推定結果
n = 10000−3−2−101xrectangular

230.6
0.5
0.0

0.1

0.2

0.3p(x)

0.4

0.5
0.0

0.1

0.2

0.3p(x)

0.4

0.5
0.4
0.3
0.2
0.1
0.0p(x)

Kernel Density Estimation (gauss)

0.6

Kernel Density Estimation (triangular)

0.6

Kernel Density Estimation (rectangular)−3−2−101x

triangular23−3−2−10123xgauss

12 / 36

カーネル関数の種類と推定結果
n = 1000000.3
0.2
0.1
0.0p(x)0.4

0.50.6

Kernel Density Estimation (gauss)−3−2−10123xgauss
12 / 36

0.50.6

カーネル密度推定量とバンド幅

band width

0.3
0.2
0.1
0.0p(x)

0.4

true
0.024
0.24
2.4

−3−2−10123x

バンド幅は適切に選ぶ必要がある．

13 / 36

最小二乗クロスバリデーション
最小二乗クロスバリデーション: LSCV, Least Squares Cross Validation.
真の密度との二乗距離 (p̂h をバンド幅 h のカーネル密度推定量とする)∫(p̂h (x) − p(x))2 dx∫∫∫= p̂h (x)2 dx − 2 p̂h (x)p(x)dx + p 2 (x)dx.
|{z}=:J(h)

J(h) を最小化すれば良い．しかし p(x) による積分がわからない→サンプルで代用．
ただし，手元にあるサンプルと p̂h は相関がある のでクロスバリデーションする．
p̂h,(−i) : i 番目のサンプル Xi を抜いて推定した密度関数.
∫n2∑bJ(h) = p̂h (x)2 dx −p̂h,(−i) (Xi ).
ni=1bJ(h) を最小にする h を採用すればよい．
14 / 36ヒューリスティクス
Silverman の方法
12

サンプル標準偏差: ŝ =√ ∑n1n2i=1 (Xi − X̄ ) .

サンプル分位点 q(0.25): サンプルの 0.25 分位点．
{}
q̂(0.75) − q̂(0.25)σ̂ = min ŝ,1.34として，次のようにする:

1.06σ̂.
n1/5
1.06 は人によって別の定数に置き換えたりする (1.06 を Scott のルール，0.9 をSilverman のルールと言う)．
ĥ =特に理論的根拠はないが，二乗誤差を漸近展開すると，[]1/5ĥ = n(∫ f C′′K(x)dx),∫K (x)2 dx(ただし CK = (∫ x 2 K (x)dx)2 ) が漸近的に最小二乗誤差を与えることが示せて，上∫のヒューリスティクスは未知の値 f ′′ (x)dx に当たりを付ける経験則とみなせる．
15 / 36

0.60.5band width

0.3
0.0

0.1

0.2

p(x)

0.4

trueCV
Silverman (1.06)−3−2−10123x

ヒューリスティクスも良い推定結果を出している．
16 / 36

多変量カーネル密度推定i
カーネル密度推定は多変量へ拡張できる． h1 K ( x−Xh ) の代わりに，d∏1j=1(j)のようにする．なお，Xihj(K(j)

x − Xihj)

はサンプル Xi の j 番目の座標．

library(MASS)
d=kde2d(x,y,c(bandwidth.nrd(x),bandwidth.nrd(y)),n=80)
image(d,xlab="latitude",ylab="longitude")17 / 36

20

40

area

60

80

100

世田谷区中古マンション価格データ

2e+07

4e+07

6e+07

8e+07

price

18 / 36

100

120

世田谷区中古マンション価格データ

1e−10
1.5e−1080

2.5e−1010404e−−10

20

3.5e−

4e10

10

3e−
2e−100area

60

3.5e−105e−110e+00

2e+07

4e+07

6e+07

8e+07

1e+08

price

18 / 36

世田谷区中古マンション価格データ

nsity

estimated deareaepric

18 / 36

カーネル密度推定を行える R 関数
使い方はスクリプトを参照

1 次元のカーネル密度推定
density(x,bw,kernel):
バンド幅: bw = ”nrd0”がデフォルト (シルバーマンの方法で 0.9 を採用),
bw=”’nrd’ で 1.06. bw=”ucv” で (バイアス修正した) クロスバリデーション,bw=”bcv” でクロスバリデーション．
カーネル関数: kernel=”gaussian” がデフォルト．他にも”epanechnikov”,”rectangular”, ”triangular”, ”biweight”,”cosine”, ”optcosine”が指定可能．

bkde: ‘KernSmooth’ パッケージに入っている．
二次元以上

kde2d: ‘MASS’ パッケージに入っている．2 次元用．
bkde2d: ‘KernSmooth’ パッケージに入っている．２次元用．
kde: ‘ks’ パッケージに入っている．3 次元以上の密度推定も行える．
MASS パケージの bandwidth.nwd は bw.nwd の 4 倍の値を返す．

19 / 36

構成

1

カーネル密度推定2

ノンパラメトリック回帰：B-スプライン

20 / 36

ノンパラメトリック回帰

4y18

6810

cubic spline fitting

method2splinelinear51015Index

21 / 36

スプライン回帰

ノンパラメトリック回帰の基本的構造:

f (x) =q∑αj Bj (x).

j=1

Bj (x) はある非線形な基底関数．ここでは，「B-スプライン基底」を考える．

22 / 36B-スプライン基底
B-スプライン基底関数は局所的な多項式関数である．
例えば，３次 B-スプラインの場合，Bk は次のような局所３次多項式である:{
ak + bk x + Ck x 2 + dk x 3(tk ≤ x ≤ tk+4 ),Bk (x) =0(otherwise).
ここで，tk は節点 (節点) と言う．
節点 (節点): t1 ≤ · · · ≤ tq+1３次スプラインの場合，tk の候補としてサンプル点 Xi を用いて次のようにした
りする:

t1 ≤ t2 ≤ t3 ≤ t4 = X1 < · · · < Xn = tn+3 ≤ tn+4 ≤ tn+5 ≤ tn+6 .
((X1 , . . . , Xn ) の外側に 6 個余分に節点を適当に設定)

23 / 36B-スプライン基底の具体的な形
係数 ak , bk , ck , dk の決め方は本講義の範疇を超えるが，j 次 B-スプライン基底は{1 (tk ≤ x ≤ tk+1 )(1)Bk (x) =,0 (otherwise)(j)

Bk (x) =x − tktk − x(j−1)(j−1)
B(x) +B(x),tk+j − tk k
tk+j+1 − tk+1 k+1なる漸化式で決まる．j = 3 の時に 3 次 B-スプライン基底を得る．

要は非線形な関数を局所的に 3 次多項式で近似しましょうということ.
3 次 B-スプラインを 3 次-スプライン，キュービック-スプラインと呼んだりする．

24 / 36

回帰係数の決定と平滑化
以後，3 次スプライン基底考え，各データ点 Xi は節点になっているものとする．
節点が Xi であるような 3 次スプライン基底を Bi (x) (i = 1, . . . , n) と書き直し，Bn+1 (x) = x, Bn+2 (x) = 1 とする．
この合計 n + 2 個の基底を用い，次のような関数を構成する:

f (x; α) =n∑αi Bi (x) + αn+1 x + αn+2 .

i=1α ∈ Rn+2 をデータから推定したい．3 次スプラインでは次のようにして推定する．
)2∫ Xn ( 2n∑d f (x; α)dx .
minn+2(Yi − f (Xi ; α))2 + λdx 2α∈RX1
i=1|{z}正則化項※ 正則化項を加えることで，推定した関数が滑らかになるように調整 →平滑化．
これがなければデータに過適合してしまう．
※ λ > 0 は適切に選ぶ必要がある (クロスバリデーション)．
※ リッジ回帰と対応．
25 / 36二次関数への展開
二乗損失は次のように展開される:
nnn+2∑∑∑(Yi − f (Xi ; α))2 =(Yi −αj Bj (Xi ))2i=1=i=1n+2 ∑n+2
∑αj αj ′(Bj (Xi )Bj ′ (Xi )) −j=1 j ′ =1i=1⊤⊤⊤

j=1n∑n+2∑αjn∑

Yi Bj (Xi ) + Y ⊤ Y

i=1

j=1
⊤

=: α B Bα − Y Bα + Y Y .
正則化項は次のように展開される:

2)2∫ Xn ( 2∫ Xn ∑n+22
d f (Xi ; α)
dB(x)
j

 dx
dx =αjdx 2dx 2
X1X1
j=1=n+2 ∑n+2∑∫ Xnαj αj ′j=1 j ′ =1=:n+2n+2 ∑∑X1d2 Bj (x) d2 Bj ′ (x)dxdx 2dx 2αj αj ′ Gj,j ′ = α⊤ G α.

j=1 j ′ =126 / 36

二次式の最小化∫ XnBi,j = Bj (Xi ), Gj,j ′ =X1d2 Bj (x) d2 Bj ′ (x)dx,dx 2dx 2として，

minn+2α∈R⇔⇒)2∫ Xn ( 2n∑d f (x; α)dx(Yi − f (Xi ; α))2 + λ
dx 2
X1i=1
⊤min α B⊤ Bα − Y ⊤ Bα + λα⊤ G αα∈Rn+2α̂ = (B⊤ B + λG )−1 B⊤ Y ,で回帰係数が求まる．

27 / 36平滑化パラメータの決定: CV, GCV
さて，λ はどう選ばばよいか？→ 例のごとくクロスバリデーション (CV).

∑n
CV =i=1 (Yi − f(−i) (Xi ))ˆ2nfˆ(−i) (Xi ) は (Xi , Yi ) を抜いて推定した 3 次スプライン関数．
∑n

GCV =2ˆ
i=1 (Yi − f (Xi )),n(1 − Tr[A(λ)]/n)2

ただし，A(λ) := B(B⊤ B + λG )−1 B⊤ . GCV は CV の計算を楽にするために提案されたが，統計的に良い性質があることが知られている．
R (smooth.spline) のデフォルトは GCV.

28 / 36

スプラインを実行
y18 <- c(1:3, 5, 4, 7:3, 2*(2:5), rep(10, 4))artdata <- data.frame(x=1:18,y=y18)
s01 <- smooth.spline(artdata)で自動的に λ も GCV で選択．

24y186810

cubic spline fitting

510Index1529 / 36

正則化パラメータの影響

4y18

6810

cubic spline fitting

lambda2CV
smallLarge

51015Index

正則化パラメータが小さいと過適合，大きすぎると線形回帰になる．
ちょうど良いパラメータは GCV で選ぶ．
s02 <- smooth.spline(artdata, spar = 0.02)
s03 <- smooth.spline(artdata, spar = 1)
30 / 36

スプラインを実行: gam

同様のことが mgcv パッケージに入っている gam という関数でも可能．


gam.s01 <- gam(y~s(x),data=artdata)s(x) と書くことでスプラインを使うことを指定．細かい次数や節点の設定も可能．
これで勝手に GCV で正則化パラメータを選んでフィッティングしてくれる．

31 / 36

gam の GCV 値

gam を使って，GCV 値を計算してみる．


sp<-seq(from=0.0001,to=0.005,length=100)
GCV_art<-numeric()for(i in 1:length(sp)){g.m<-gam(y~s(x),sp=sp[i],data=artdata)
GCV_art[i]<-g.m$gcv.ubre}plot(sp,GCV_art,type="l",lwd=2,xlab="lambda",ylab="GCV")

32 / 36

0.70

0.75

GCV

0.80

0.85

gam の GCV 値

0.0000.001

0.002

0.003

0.004

0.005

lambda33 / 36

smooth.spline と gam の比較

4y18

68

10smooth.spline and gam

method2smooth.spline(GCV)
overfitgam(GCV)

51015Index

ほぼ同じ結果34 / 36

まとめ２つのノンパラメトリックな推定手法を紹介した．
1

密度推定→カーネル密度推定

2回帰→ B-スプライン回帰

どちらも，バンド幅や正則化パラメータといったパラメータを CV などを用いてうまく調整する必要があった．

35 / 36

講義情報ページ
http://www.is.titech.ac.jp/~s-taiji/lecture/dataanalysis/dataanalysis.html

36 / 36