【Python】kNNを使った分散関数の推定【統計学入門】

こんにちは。たくまろです。
今回は最近勉強した分散関数の推定について関連論文と手法について思います。
ここでは、direct methodとresidual-based mthodの2つについて、紹介します。

論文の内容をまとめつつPythonのコードを貼り付けました。
分散関数の推定について、気になった方は是非ご覧ください！

たくまろのプロフィール

大学・大学院の数学専攻で統計学の勉強をしていました。
現在はデータサイエンスとして働いています。

Twitter >>> @takumaroblog

問題設定

$X$: $d$次元特徴ベクトル
$Y$: 1次元の目的変数
$W$: ノイズ s.t. $E[W|X] = 0, E[W^2]<∞$

$$Y = f(X) + W$$

任意の$x∈R^d$について、

$$f(x) = E[Y|X=x]$$

$$σ^2(x) = Var(Y|X=x) = E[(Y-f(x))^2|X=x]$$

を満たすような$σ^2(x)$を見つけることが今回のモチベーションです。

$σ^2(x)$は大きく分けると2つの方法に分類されます。

方法1: The direct method

$$σ^2(x) = E[Y^2|X=x] – (E[Y|X=x])^2$$

の形で書くことができる。この性質を利用して次のステップで分散関数を予測します。

目的変数$Y$に関する回帰関数$f(X)$を求める
目的変数$Y^2$に関する回帰関数$g(X)$を求める
分散関数$σ^2(x) = g(X) – (f(X))^2$で計算する

ダイレクトに$E[Y^2|X=x]$と$E[Y|X=x]$の推定値を求めて、その差分から$σ^2(x)$を推定しています。

単純で分かりやすい反面、推定値$σ^2(x)$が非負にならないような問題があります。

Härdle and Tsybakov(1977)では多項式カーネルを用いた回帰でアプローチを行なっています。

方法2: The residual-based method

文字通り、残差に対する回帰で分散関数の予測を考えます。次のステップで分散関数を予測します。

目的変数Yに関する回帰関数$f(X)$を求める
残差$r=(Y-f(X))^2$を計算
目的変数$r$に関する回帰関数$σ^2(x)$を求める

この方法を用いることで上記のdirect methodで問題となっていた非負性を解消することが可能です。

パラメトリックな設定では、Ruppert et al.(1997)やFan and Yao(1998)ではKernel関数を用いた重み付けをして、局所的な回帰を考えています。また、Xu and Phillips(2011)よって重みの付けし直されたKernel関数を提案しています。また、Lawrence et al.(2007)などでは、目的変数の差分に注目した方法を考えています。

ノンパラメトリックな設定では、Wang et al.(2008)やKulik and Wichelhaus(2011)ではKernel密度推定的を使ったノンパラメトリックな回帰を考え、与えられた分散関数の漸近正規性について言及しています。またDenis et al.(2020)ではknnを用いた手法を基に、理論的に保証できないような分散関数をrejectする仕組みを考えています。

以下knnを使ったresidual methodの（適当に書いた）Pythonコードです。

from sklearn.neighbors import KNeighborsRegressor
from sklearn.model_selection import KFold

#Fold数
NFold = 5
    
def var_estimation(X_train,y_train):
  #N-fold
  def stratified(X_train, y_train):
    trn_idx_list = []
    val_idx_list = []
    folds = KFold(
        n_splits=NFold,
        shuffle=False,
    )
    for trn_idx, val_idx in folds.split(X_train.values, y_train.values):
        trn_idx_list.append(trn_idx)
        val_idx_list.append(val_idx)
    return trn_idx_list, val_idx_list
    
  trn_idx,val_idx = stratified(X_train, y_train)
  sigma_trains = []
  
  for id in range(NFold):
      
      #平均の学習用
      X_train1 = X_train.iloc[trn_idx[id]]
      y_train1 = y_train[trn_idx[id]]
      #残差の計算用
      X_train2 = X_train.iloc[val_idx[id]]
      y_train2 = y_train[val_idx[id]]
  
      #kNNによる平均関数のモデルを作成
      model1 = KNeighborsRegressor()
      model1.fit(X_train1, y_train1)
      pred_train = model1.predict(X_train2)
      sigma_train = (pred_train - y_train2)**2 
      sigma_trains.append(sigma_train)
  
  #予測した残差sigma_trainsを結合してsigma_pred作成
  sigma_pred = sigma_trains[0]
  for k in range(1,5):
      sigma_pred = sigma_pred.append(sigma_trains[k])
      
  #kNNによる分散関数のモデル作成
  model2 = KNeighborsRegressor()
  model2.fit(X_train, sigma_pred)
  return model2

残差rの計算の注意点として、N-foldを用いて回帰関数f(・)と実際に代入する(X,Y)が独立になるように設定している。下記図を参照のこと。