#はじめに
東京工業大学/株式会社Nospare の栗栖です．この記事では関数データに対するノンパラメトリック回帰分析について紹介します．この記事の内容は主に Ferraty and Vieu(2006), Ferraty and Nagy (2022) の内容をもとにしています．

この記事の内容は関数データ分析に関する記事
「ノンパラメトリック関数回帰(1)」
の続編です．

前回の記事ではノンパラメトリックな回帰関数の推定量として局所定数推定量 (local constant estimator, LC推定量) の性質について解説しましたが，この記事ではその拡張である局所線形推定量 (local linear estimator, LL推定量) について紹介し，その性質について紹介します．

#モデル

前回の記事と同様に以下のノンパラメトリック回帰モデルを考えます．

\begin{align*}
Y_i &= m(X_i) + e_i,
\end{align*}

$i=1,\dots, n$. ここで，$\{(Y_i,X_i)\}_{i=1}^{n}$ は観測データで $\{e_i\}$ は $E[e_{i}|X_{i}]=0$ を満たす観測誤差，$Y_i$ は実数値，$X_i$ は内積 $\langle \cdot, \cdot \rangle$ をもつ(可分)ヒルベルト空間 $H$ に値をとるとします．また $n$ 個のデータは独立同分布であるとします．

推定量の構成

以下ではデータ$\{(Y_i,X_i)\}_{i=1}^{n}$ を用いて回帰関数 $m$ をノンパラメトリックに推定する問題を考えてみましょう．

局所定数推定量 (LC推定量)

ここでは特に以下の推定量(Nadaraya-Watson 推定量, NW推定量)を用いて $m$ を推定することを考えます．

\begin{align*}
\widehat{m}(x) &= {\sum_{i=1}^{n}Y_iK(\|X_i-x\|/h) \over \sum_{i=1}^{n}K(\|X_i-x\|/h)},\ x \in H.
\end{align*}

ただし，$K:[0,1] \to [0,\infty)$ は適当な条件を満たすカーネル関数，$h=h_{T} \to 0$ ($T \to \infty$) はバンド幅です．
上記のNW推定量は，$X_{t}$ が実数値 ($\mathbb{R}$) や実数値ベクトル ($\mathbb{R}^{p}$) の場合に，推定量の定義に登場するノルム $\|\cdot\| = \sqrt{\langle \cdot, \cdot \rangle}$ を絶対値やユークリッドノルムに置き換えれば通常のNW推定量と一致します．この意味で $\widehat{m}$ は通常のNW推定量を関数データの枠組みに自然に拡張したものとみることができます．さらに $\widehat{m}(x)$ は以下のように定義することもできます．

\begin{align*}
\widehat{m}(x) &= \text{argmin}_{\theta}\sum_{i=1}^{n}\left(Y_i - \theta\right)^2 K\left(\|X_i-x\|/h\right)
\end{align*}

この定義から，$\widehat{m}(x)$ は局所定数推定量 (local constant estimator, LC推定量) とも呼ばれます．NW推定量の性質については前回の記事「ノンパラメトリック関数回帰(1)」で紹介しています．

局所線形推定量 (LL推定量)

Ferraty and Vieu (2006) では LC推定量の漸近的性質を議論していましたが，Ferraty and Nagy (2022) はそこでの議論を拡張して局所線形推定量 (local linear estimator, LL推定量) の性質を調べています．
LC推定量の性質を調べる際，回帰関数 $m$ に対して以下の条件を仮定しました．任意の $x' \in B(x,1)=\{y \in H: \|x-y\|\leq 1\}$ に対して，$\beta>0$, $C>0$ が存在して

|m(x) - m(x')| \leq C\|x- x'\|^{\beta}

を満たす．

回帰関数の滑らかさ

ここでは回帰関数 $m$ が $x\in H$ の近傍で以下の2次の Taylor 展開をもつと仮定します：
$0 \in H$ のある近傍 $N_{0}$ が存在して任意の $u \in N_{0}$ に対して $\zeta = x + tu$，$t \in (0,1)$ が存在して以下の展開が成り立つ．

\begin{align*}
m(x+u) = m(x) + \langle m'_{x}, u \rangle + {1 \over 2}\langle m''_{\zeta} u, u\rangle．  
\end{align*}

ここで $m'_{x} \in H$，$m''_{\zeta}$ は $H$ から $H$ への線形ヒルベルト-シュミット作用素，$\nu \mapsto m''_{\nu}$ は $x \in H$ の近傍 $N_{x}$ でリプシッツ連続とします．
上記の仮定は通常の場合 ($H = \mathbb{R}^{p}$) における Taylor 展開を一般化したものになっています．

関数シングルインデックスモデル　

Ferraty and Nagy (2022) では Berkeley growth data を用いて，
$Y_{i} =$18歳での身長，$X_{i}=$ 10歳までの発育速度 (growth velocity)
の場合における関数シングルインデックスモデル

Y_i = g(\langle X_i, \beta \rangle) + e_i

分析結果 (ノンパラメトリックなリンク関数 $g$ と $\beta$ の推定) が紹介されています．
特にここでは $g: \mathbb{R} \to \mathbb{R}$ が2階連続微分可能で，$g'(x)>0$, $\|\beta\|=1$ としておきます．$g$ の Taylor 展開を考えると

\begin{align*}
g(\langle x+u, \beta \rangle) &= g(\langle \beta, x\rangle) + g'(\langle x,\beta \rangle)\langle \beta, u \rangle + {1 \over 2}g''(\langle x+tu, \beta \rangle)\langle u,\beta\rangle^2\\
&= g(\langle \beta, x\rangle) + \langle g'(\langle x,\beta \rangle)\beta, u \rangle + {1 \over 2}\langle g''(\langle x+tu, \beta \rangle)\langle \beta, u \rangle \beta ,u\rangle 
\end{align*}

と書けるので，この場合，$m'_{x} = g'(\langle x,\beta \rangle)\beta$，$m''_{\nu}(\cdot) = g''(\langle \nu, \beta \rangle)\langle \beta, \cdot \rangle \beta$ となります．特に

$\beta$ については

\beta = {m'_{x} \over \|m'_{x}\|} = {g'(\langle x,\beta \rangle)\beta \over g'(\langle x,\beta \rangle)}

という関係が成り立つので， $m'_{x}$ をLL推定量で推定すれば $\beta$ を $\widehat{\beta} = \widehat{m}'_{x}/\|\widehat{m}'_{x}\|$ で推定することができます．

$g(x)$ については $Z_i = \langle X_i, \widehat{\beta} \rangle$，$i=1,\dots,n$ を用いて通常の $H = \mathbb{R}$ の場合のノンパラメトリック回帰 (LC or LL推定量) で推定可能です．

最適化問題

$m$ の Taylor 展開における $m(x)$，$m'_{x}$ を推定するため，以下の最適化問題を考えましょう．

\begin{align*}
(\widehat{a}, \widehat{\beta}) &= \text{argmin}_{a, \beta}\sum_{i=1}^{n}\left(Y_i - a - \langle \beta, X_i-x \rangle \right)^2 K(\|X_i-x\|/h).
\end{align*}

ここで，$\widehat{a}$, $\widehat{\beta}$ はそれぞれ $m(x)$, $m'_{x}$ の推定量です．
上記の最適化問題は $X_{i} \approx x$ を与えたときの $Y_{i}$ の条件付き期待値を局所的に線形近似することに対応しています．

\begin{align*}
E[Y_i|X_i] &\approx a + \langle \beta, X_i - x \rangle.  
\end{align*}

このような定式化から，$\widehat{a}$, $\widehat{\beta}$ は局所線形推定量と呼ばれます．
実際に $m(x)$, $m'_{x}$ の推定量を計算する際には内積 $\langle \beta, X_i - x \rangle$ を $H$ の正規直行基底 $\{\phi_j\}_{j=1}^{\infty}$ で有限次元近似したバージョン

\beta \approx \sum_{j=1}^{J}\langle \phi_j, \beta \rangle \phi_j

を用いて，以下の最適化問題を解くことを考えます．

\begin{align*}
(\widehat{a}, \widehat{b}_1,\dots, \widehat{b}_J) &= \text{argmin}_{a, b}\sum_{i=1}^{n}\left(Y_i - a - \sum_{j=1}^{J}b_j\langle \phi_j, X_i-x \rangle \right)^2 K(\|X_i-x\|/h).
\end{align*}

ここで $\widehat{a} = \widehat{m}(x)$，$\widehat{b}_j = \widehat{\langle \phi_j, \beta \rangle}$，$j=1,\dots,J$ です．

推定量の計算

まずいくつか記号を導入しておきます．
$Y = (Y_{1},\dots,Y_{n})'$, $K = \text{diag}(K(\|X_1-x\|/h),\dots, K(\|X_n-x\|/h))$,

\Phi_{X} = \left(
\begin{matrix}
1 & \langle \phi_1, X_1-x \rangle & \cdots & \langle \phi_J, X_1-x \rangle \\
\vdots & \vdots & \ddots & \vdots \\
1 & \langle \phi_1, X_n-x \rangle & \cdots & \langle \phi_J, X_n-x \rangle
\end{matrix}
\right).

$V = (\widehat{a}, \widehat{b}_1,\dots, \widehat{b}_J)'$ とすると，上記の最適化問題は以下のように書き換えられます．

V = \text{argmin}_{v}(Y - \Phi_{X} v)'K(Y - \Phi_{X}v).

右辺の関数を $v$ について微分して $=0$ とすると，

\begin{align*}
&{\partial \over \partial v}(Y - \Phi_{X} v)'K(Y - \Phi_{X}v) = 0\\
&\Leftrightarrow - \Phi'_{X}KY +  \Phi'_{X}K\Phi_{X}v = 0. 
\end{align*}

従って $V = (\Phi'_{X}K\Phi_{X})^{-1}\Phi'_{X}KY$ が得られます．$\Phi = (\phi_{1},\dots,\phi_{J})'$, $e = (1,0,\dots,0)' \in \mathbb{R}^{J+1}$，$[0|I]$を以下で定義される $J \times (J+1)$ 行列とします．

[0|I] = \left(
\begin{matrix}
0 & 1 & 0 &\cdots & 0 \\
0 & 0 & 1 & \cdots & 0 \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
0 & 0 & 0 & \cdots & 1
\end{matrix}
\right).

このとき，$V$ の定義から，

\begin{align*}
\widehat{m}(x) &=  e'(\Phi'_{X}K\Phi_{X})^{-1}\Phi'_{X}KY\\
&= \widehat{a},\\
\widehat{m}'_{x} &= \Phi'[0|I](\Phi'_{X}K\Phi_{X})^{-1}\Phi'_{X}KY\\
&= \sum_{j=1}^{J}\widehat{b}_j \phi_j = \sum_{j=1}^{J}\widehat{\langle \phi_j, \beta \rangle}\phi_j
\end{align*}

と推定量が計算されます．

推定量の性質

以下では $n \to \infty$ の場合における $\widehat{m}(x)$ の漸近的性質を紹介します．$x$ を中心とする $X_i$ の small ball probability を $\pi_{x}(h)=P(X_t \in B(x,h))$，
$S_{J}$ を $\phi_1, \dots, \phi_{J}$ が張る $H$ の部分空間，
$S_{J}^{\top}$ を $S_{J}$ の直行補空間，
$P_{S_{J}^{\top}}: H \to S_{J}^{\top}$ を $H$ から $S_{J}^{\top}$ 上への射影とします．
Ferraty and Nagy (2022) では適当な条件のもとで $\widehat{m}(x)$ の収束レートが導出されています (Theorem 1)．

\widehat{m}(x) - m(x) = O(h^2) + O\left(\|P_{S_{J}^{\top}}m'_{x}\|h\right) + O_{p}\left(\sqrt{1 \over n\pi_{x}(h)}\right).

記号が複雑になるため省略しますが，さらに Ferraty and Nagy (2022) では $\widehat{m}'_{x}$ の収束レートも導出されています (Theorem 2)．

#まとめ
この記事では関数データに対するノンパラメトリック回帰分析 (局所線形推定量) について解説しました．「ノンパラメトリック関数回帰(3)」では非定常関数時系列データに対するノンパラメトリック回帰 (局所定数推定量) について紹介する予定です．

株式会社Nospareには統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータの分析につきましては株式会社Nospare までお問い合わせください．

参考文献
[1] Ferraty, F. and Nagy, S. (2022). Scalar-on-function local linear regression and beyond. Biometrika 109, 439-455.
[2] Ferraty, F. and Vieu, P. (2006). Nonparametric Functional Data Analysis: Theory and Methods. Springer.