はじめに
東京大学/株式会社Nospare リサーチャーの栗栖です．
この記事では前回の局所線形回帰についての記事から引き続き，ノンパラメトリックな回帰分析の統計手法の一つである局所多項式回帰について紹介します．
モデル
ノンパラメトリックな回帰モデル
Y_i = m(\boldsymbol{X}_i) + U_i,\ E[U_i|\boldsymbol{X}_i] = 0
から独立同分布 (independent and identically distributed, i.i.d.) なデータ $(Y_i, \boldsymbol{X}_i)$, $i = 1,\dots,n$ が得られたとします．ただし，$Y_i \in \mathbb{R}$, $\boldsymbol{X}_i=(X_{1,i},\dots,X_{d,i})' \in \mathbb{R}^d$とします．
以下ではデータから回帰関数 $m(x) = E[Y_i|\boldsymbol{X}_i = x]$, $x = (x_1,\dots,x_d) \in \mathbb{R}^d$ を推定するという問題を考えましょう．推定にあたり，上記のモデルに対して以下を仮定します．
(A1:LP)
[1] $X_i$は密度関数 $f$ をもち，$f$ は $x$ のある近傍 $N_x$ で連続．
[2] $m$ は $N_x$ で $p+1$ 回連続偏微分可能．
[3] $\sigma^2(x) = E[U_i^2|\boldsymbol{X}_i = x]$ は $x$ で連続．
また $x \in \mathbb{R}^d$ で $p+1$ 回連続偏微分可能な関数 $g: \mathbb{R}^d \to \mathbb{R}$ に対して，
\partial_{j_1 \dots j_L} g(x) = {\partial \over \partial x_{j_1} \dots \partial x_{j_L}}g(x),\ 1\leq j_1 \leq \dots \leq j_L \leq d,\ L=1,\dots, p+1
とします．
局所多項式回帰
回帰関数 $m$ の推定にあたり，以下の最適化問題を考えましょう．
\begin{align*}
\hat{\beta} &= (\hat{\beta}_0, \hat{\beta}_{11},\dots, \hat{\beta}_{dd}, \hat{\beta}_{111},\dots, \hat{\beta}_{ddd},\dots, \hat{\beta}_{1\dots 1}, \dots, \hat{\beta}_{d \dots d})\\
&= \underset{\beta \in \mathbb{R}^{D}}{\text{argmin}}\!\sum_{i=1}^n \!\!\left(\! Y_i \!-\! \beta_0 \!-\! \sum_{L=1}^{p}\sum_{1 \leq j_1 \leq \dots \leq j_L \leq d}\!\!\!\!\!\!\!\!\!\!\! \beta_{j_1 \dots j_L}\! \prod_{\ell=1}^{L}(X_{j_\ell, i} - x_{j_\ell})\!\!\right)^2 \!\!\! K \!\! \left(\!{\boldsymbol{X}_i - x \over h} \!\right),\\
K\left({\boldsymbol{X}_i - x \over h}\right) &= K\left({X_{1,i} - x_1 \over h},\dots, {X_{d,i} - x_d \over h}\right).
\end{align*}
ここで，
D = 1 + \#\{(j_1,\dots, j_L): 1\leq j_1 \leq \dots \leq j_L, 1\leq L\leq p\},
$h = h_{n} \to 0$, $n \to \infty$ はバンド幅，$K:\mathbb{R}^d \to \mathbb{R}$ はカーネル関数と呼ばれ，以下を満たすとします．
(A2:LP)
[1] $\int K(u)du=1$
[2] $K(u)$ は有界かつ $[-1,1]^d$ に台をもつ ($K(u) = 0$, $\|u\| >1$)
[3] $\tilde{u} = (u_{j_1\dots j_L})_{1\leq j_1 \leq \dots \leq j_L\leq d, 1 \leq L \leq p}$, $\check{u} = (1, \tilde{u})'$ として，$S=\int \check{u}\check{u}'K(u)du$ は正則(逆行列をもつ)
上記の最適化問題は $\boldsymbol{X}_i = x$ の近傍で「局所」的に $Y_i$ を「多項式関数」に「回帰」していると解釈できるので，「局所多項式回帰」と呼ばれます．また得られる解は「局所多項式推定量」(local polynomial estimator, LP 推定量) と呼ばれることもあります．特にこの方法では，$\boldsymbol{X}_i = x$ の近傍で
\begin{align*}
Y_i &\approx m(\boldsymbol{X}_i) \approx m(x) + \sum_{L=1}^{p}\sum_{1 \leq j_1 \leq \dots \leq j_L\leq d}{\partial_{j_1 \dots j_L} m(x) \over s_{j_1\dots j_L}!}\prod_{\ell=1}^{L}(X_{j_\ell, i} - x_{j_\ell})
\end{align*}
と近似できることを利用して $m(x)$ (とその偏微分 $\partial_{j_1\dots j_L} m(x)$) を推定しています．ここで，
\begin{align*}
s_{j_1\dots j_L}! &= \prod_{k=1}^{d}s_{j_1\dots j_Lk}!,\\
s_{j_1\dots j_Lk} &= \#\{j_\ell: j_\ell = k, 1\leq \ell \leq L\}.
\end{align*}
です．上記の最適化問題を解くと，$m(x)$ の LP推定量が以下で与えられることが分かります．
\begin{align*}
\hat{\beta} &= \left(\hat{m}(x),\partial_1 \hat{m}(x) ,\dots, \partial_d \hat{m}(x), {\partial_{11}\hat{m}(x) \over s_{11}!},\dots,{\partial_{dd}\hat{m}(x) \over s_{dd}!} \right.,\\
&\left. \quad {\quad \partial_{111}\hat{m}(x) \over s_{111}!},\dots,{\partial_{ddd}\hat{m}(x) \over s_{ddd}!},\dots,{\partial_{1\dots 1}\hat{m}(x) \over s_{1\dots1}!},\dots,{\partial_{d \dots d}\hat{m}(x) \over s_{d \dots d}!}\right) \\
&= \left(\sum_{i=1}^{n}K\left({\boldsymbol{X}_i - x \over h}\right)\tilde{\boldsymbol{X}}_i\tilde{\boldsymbol{X}}'_i \right)^{-1}\sum_{i=1}^{n}K\left({\boldsymbol{X}_i - x \over h}\right)\tilde{\boldsymbol{X}}_iY_i
\end{align*}
ここで $\tilde{\boldsymbol{X}}_i = (1, \check{\boldsymbol{X}}_i)'$, $\check{\boldsymbol{X}}'_i = ((\boldsymbol{X}_i - x)'_1,\dots, (\boldsymbol{X}_i - x)'_p)'$,
(\boldsymbol{X}_i - x)_L = \left(\prod_{\ell=1}^{L} (X_{j_\ell,i} - x_{j_\ell})\right)_{1 \leq j_1 \leq \dots \leq j_L \leq d}
です．このとき，$\boldsymbol{X}_i = x$ の近傍で
\begin{align*}
Y_i &\approx m(\boldsymbol{X}_i) \approx \tilde{\boldsymbol{X}}'_iM(x),\\
M(x) &= \left(m(x), \check{M}(x)\right)',\ \check{M}(x) = \left({\partial_{j_1 \dots j_L} m(x) \over s_{j_1\dots j_L}!}\right)'_{1 \leq j_1 \leq \dots \leq j_L\leq d, 1 \leq L \leq p}
\end{align*}
と書けることに注意しておきます．従って，$\hat{\beta}$ は $M(x)$ の推定量とみることができます．
LP推定量の漸近正規性
最後にLP推定量の漸近的性質について紹介しておきます．まず(A1:LP), (A2:LP) に加えて以下を仮定します．
(A3:LP)
[1] $n \to \infty$ として，$h \to 0$, $nh^{d+2p} \to \infty$, $nh^{d+2(p+1)} \to c \in (0,\infty]$．
[2] また $x \in \mathbb{R}^d$ に対して $f(x)>0$，
[3] ある $\delta>0$ と $x$ の近傍 $N_x$ に対して $\sup_{z \in N_x}E[|U_i|^{2+\delta}|\boldsymbol{X}_i = z] \leq U(x)<\infty$ ．
$H = \text{diag}(1,h,\dots,h, h^2,\dots, h^2, \dots, h^p, \dots, h^p) \in \mathbb{R}^D$ とします．
仮定(A1:LP)-(A3:LP)の下で，$n \to \infty$ として以下が成り立ちます．
\begin{align*}
\sqrt{nh^d}\left(H\left(\hat{\beta} - M(x)\right)- S^{-1}B_p M_{p,n}(x)\right) &\stackrel{d}{\to}N
\left( \boldsymbol{0}, {\sigma^2(x) \over f(x)}S^{-1} V S^{-1}
\right).
\end{align*}
ここで，$V = \int \check{u}\check{u}'K^2(u)du$,
\begin{align*}
B_p &= \int \check{u}(u)'_{p+1}K(u)du,\ (u)_{p+1} = \left(\prod_{\ell=1}^{p+1}u_{j_\ell}\right)'_{1 \leq j_1 \leq \dots \leq j_{p+1} \leq d},\\
M_{p,n}(x) &= \left({\partial_{j_1 \dots j_{p+1}} m(x) \over s_{j_1\dots j_{p+1}}!}h^{p+1}\right)'_{1 \leq j_1 \leq \dots \leq j_{p+1}\leq d}.
\end{align*}
LP推定量の性質
- 多項式の保存：$e_1 = (1,0,\dots,0)' \in \mathbb{R}^{D}$, $\breve{\boldsymbol{X}}_i = ((\boldsymbol{X}_i)'_1,\dots, (\boldsymbol{X}_i)'_p)'$, $\breve{x} = ((x)_1,\dots,(x)_p)'$ として
W(\boldsymbol{X}_i) = e'_1\left(\sum_{i=1}^{n}K\left({\boldsymbol{X}_i - x \over h}\right)\tilde{\boldsymbol{X}}_i\tilde{\boldsymbol{X}}'_i \right)^{-1}K\left({\boldsymbol{X}_i - x \over h}\right)\tilde{\boldsymbol{X}}_i
とおくと，
\hat{\beta}_0 = \sum_{i=1}^{n}W(\boldsymbol{X}_i)Y_i
と書けることから，LL推定量と同様に $W(\boldsymbol{X}_i)$ は $m(x)$ を推定する際の $Y_i$ の $\boldsymbol{X}_i$ 貢献度と解釈することができます．またLP推定量は多項式関係を復元することが分かります．実際，任意の $(u_0, \tilde{u})' \in \mathbb{R}^D$ に対して $Y_i = u_0 + \tilde{u}'\breve{\boldsymbol{X}}_i$ とすると，
\begin{align*}
\sum_{i=1}^{n}W(\boldsymbol{X}_i)(u_0 + \tilde{u}'\breve{\boldsymbol{X}}_i) &= \sum_{i=1}^{n}W(\boldsymbol{X}_i)(u_0 + \tilde{u}'\breve{x} + \check{u}'\check{\boldsymbol{X}}_i)\\
&= \sum_{i=1}^{n}W(\boldsymbol{X}_i)\tilde{\boldsymbol{X}}_i \left(
\begin{array}{c}
u_0 + \tilde{u}'\breve{x} \\
\tilde{u}
\end{array}
\right) \\
&= u_0 + \tilde{u}'\breve{x}
\end{align*}
が成り立ちます．これは LC, LL推定量にはない性質です．
まとめ
この記事ではノンパラメトリック回帰分析の１つとして局所多項式回帰を紹介しました．またこの記事ではデータが独立同分布の場合を考えましたが，次の記事では空間データ分析への応用として，データが空間相関をもつ場合への局所多項式推定量の拡張について紹介する予定です．
株式会社Nospareには統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータ分析につきましては株式会社Nospareまでお問い合わせください．
