突然終わるかもしれないブログ

確率や統計の内容について記事を書く予定です.

ブートストラップ法

ブートストラップ法についてまとめたいと思います.

ブートストラップ法とは

モチベーション

シンプルな例を通して,ブートストラップ法のモチベーションを説明します.

今,平均 $\theta$ 分散1の正規分布 $\mathcal{N}(\theta,1)$ から以下のような5つの標本を得たとします:
$x_1= -1.34$
$x_2= -0.81$,
$x_3= 0.73$,
$x_4= -0.66$,
$x_5= -0.36$.

このとき$\theta$の最尤推定量は標本平均ですので,

$\hat{\theta}_{\rm MLE} = \frac{x_1+\cdots +x_5}{5}=-0.488$.

再生性から標本平均は分散$1/\sqrt{5}$の正規分布に従いますので,
95%信頼区間は
\begin{align*}
-0.488-1.96\times \frac{1}{\sqrt{5}} \leq \hat{\theta}_{\rm MLE} \leq -0.488+1.96\times \frac{1}{\sqrt{5}}
\end{align*}

このように,モデルを特定していれば推定量$\hat{\theta}$の統計的性質(信頼区間等)が計算できます.

ここからが本題ですが,以下のような問題を考察します.

1. 正規性を仮定できない場合にはどうするか?

``標本が $f(x-\theta)$ ($f$:ある確率密度関数)というモデルから得られた"
という情報しかない(というモデルを仮定した)場合に,
上述のように$\theta$の推定量(算術平均)の信頼区間を作れるでしょうか?

2. 複雑な推定量の場合にはどうするか?

例えば,positive part James--Stein 推定量のような場合に信頼区間を構成できるでしょうか?*1
(参考: James–Stein estimator - Wikipedia, the free encyclopedia)


ブートストラップ法はこのような問題に計算機ベースで答えを出します.

大事なことですが,ブートストラップの応用は広く,
上述の信頼区間の問題は応用例の一つです.
そのほかにもバイアス補正や検定,予測といった問題に応用できます.

ブートストラップ法

ブートストラップ法とは,一言で述べれば
``復元抽出に基づく計算機ベースの統計的推測"
と言えるのではないでしょうか.

具体的に上述の1の問題に答える形で,
ブートストラップ法について説明したいと思います.

まず次の1と2の手順を$B$回繰り返します.

  1. $\{x_1,\ldots,x_5\}$から復元抽出で$\{x_1^{*},\ldots,x_5^{*}\}$を作る.(ブートストラップ標本と呼ばれる)
  2. $\{x_1^{*},\ldots,x_5^{*}\}$を基に最尤推定量(算術平均)$\hat{\theta}_{\rm MLE}^*$を作る.

こうして$\hat{\theta}_{\rm MLE}^*$を$B$個作ります.

信頼区間を求めたければ,$B$個の$\hat{\theta}_{\rm MLE}^*$の95%区間を計算します.
つまり$B$個をソートして,上下2.5%のとこにある値を見れば良いです

以上です.非常にシンプルです.
実際に最初の具体例でブートストラップ法によって構成された
$B$個の$\hat{\theta}_{\rm MLE}^*$のヒストグラムを見てみます.

f:id:mkprob:20150524155725p:plain

ブートストラップ法の良さ

良いところは以下の二点ではないかと思います.

  1. シンプル
  2. 計算機ベースでノンパラメトリックな方法

2つ目について.
ノンパラ・セミパラメトリックモデルの場合,得てして理論的解析が複雑です.
そうしたモデルでも,計算機ベースでシンプルに統計的推測ができるので
ブートストラップは重宝されるのではないでしょうか.

ブートストラップは何故``うまく"いくのか?

ブートストラップが何をしているのかを最後に説明したいと思います.

Efron & Tibshirani の本の絵による説明が一番わかりやすいと思います.

f:id:mkprob:20150524165404p:plain

真のモデル $P$ が分かっていれば,そこから標本を抽出し推定量を構成することで,
推定量の従う分布がわかります.

実際は真のモデル $P$ がわからないので,経験分布 $\hat{P}$ を $P$ だと思って,
同じことをしているのがブートストラップ法ということになります.

経験分布とは観測値$x_1,\ldots,x_n$のそれぞれの値を確率$1/n$でとる確率分布の事です.
従って経験分布からの標本を得たければ,復元抽出ということになります.

また $P$ を $\hat{P}$ に置き換えてることからもわかるように,
観測数$n$が大きい方がブートストラップ法の精度も上がります.

参考文献

以下の本がわかりやすいです.
今回はブートストラップの基本的なところをまとめたので,
もう少し突っ込んだ話題も,そのうちまとめたいです.


計算統計学の方法―ブートストラップ・EMアルゴリズム・MCMC (シリーズ予測と発見の科学 5)

計算統計学の方法―ブートストラップ・EMアルゴリズム・MCMC (シリーズ予測と発見の科学 5)

An Introduction to the Bootstrap (Chapman & Hall/CRC Monographs on Statistics & Applied Probability)

An Introduction to the Bootstrap (Chapman & Hall/CRC Monographs on Statistics & Applied Probability)

*1:できたらごめんなさい.

Ridge回帰とLasso

Ridge回帰とLassoの簡単な方法の紹介です(メモ).

Rのパッケージglmnetを使って簡単にできます*1

Ridge回帰とLasso

線形回帰モデルとは
\[
y=X\beta+\varepsilon,\quad y\in\mathbf{R}^n,\quad X\in\mathbf{R}^{n\times p},\quad \beta\in\mathbf{R}^p,\quad \varepsilon\in\mathbf{R}^n,
\]
というモデルです.$y$は目的変数,$X$は説明変数(計画行列)です.誤差項$\varepsilon$には平均0の正規分布を仮定するのが基本的です.

さてデータ$(X,y)$からパラメータ$\beta$を推定することが線形回帰モデルの目標ですが,Ridge回帰は
\[
\hat{\beta}_{\mathrm{ridge}}:=\mathrm{argmin}_{\beta} ~\{(y-X\beta)^2 +\lambda |\beta|^2\},
\]
を推定量として得る方法です.ここで$|\cdot |$はユークリッドノルム,$\lambda\geq 0$は正則化項(自分で決めるもの)です.$\lambda=0$のときは最小二乗推定量となります.$\lambda>0$のとき$\beta$の大きさに罰則をつけることで過学習を防ぎます.特に説明変数の間に相関が強い場合(共線性)有効な手法です.共線性があるときに二乗ノルムの罰則を与えればよさそうという簡単な理由は,最小二乗推定量(OLS)の分散が$X^{\top}X$の行列式に比例するからです:
\[
\mathrm{E}_{\beta}|\hat{\beta}_{\mathrm{OLS}}|^2\propto \det(X^{\top}X)+|\beta|^2
\]
Hoerl and Kennard (1970)によって提案されたよく知られた手法です.

Lassoは
\[
\mathrm{argmin}_{\beta} ~\{(y-X\beta)^2 +\lambda |\beta|\},
\]
を推定量として得る方法です.Ridgeと違い$L^1$の罰則をつけることにより疎(スパース)な解を得ることができます.つまりLassoで得られた推定量$\hat{\beta}$の成分の多くが真に0をとりやすい,という性質があります.Tibshirani(1996)により提案され,非常によく研究されている手法の一つです.


どちらの手法も最小二乗法に比べて良い推定や予測ができるということがあります.しかし実際にridgeとlassoをやるには正則項$\lambda$を決める必要があります.正則化項の基本的な決め方の方法に,交差確認による決め方があります.

まず$\lambda$の候補の値(複数個,たくさん)を用意します.学習データを$n$個に分割し,そのうちの$n-1$個のパートと候補の$\lambda$の値を基に$\beta$を推定します($\hat{\beta}$).残しておいた分割に対して,学習した$\hat{\beta}$を用いた予測値をだし,実際の値と予測値との誤差を計算します.これを$n$回繰り返し,一番誤差が小さかった$\lambda$を選ぶ方法です.手間がかかります.


実際あるデータでやってみた結果です:

f:id:mkprob:20141224000310p:plain

横軸は候補の$\lambda$の値(のlog)を表し,縦軸はその$\lambda$の下での予測誤差をプロットしたものです.
赤線は予測誤差の平均,それと上下のバウンドがついています.
縦の値が一番小さい$\lambda$を選択します.この図だと2.225(logだと0.799)となっています.

実データでの例

RのパッケージSMPracticalsのairpollution dataを使った実験をします.
このデータは$n=60$,$p=15$のデータです.変数の数を増やすために二次の項(変数$i$と変数$j$の積)を加えて$p=135$としたデータでRidgeとLassoを比較します.100回ランダムに学習データとテストデータを半分ずつに分けて予測誤差を比較しました.$\lambda$を決める際の分割数はデフォルトの10としています.

#Lasso and Ridge
set.seed(1) #random partition

#input data
airpollution <- as.matrix(read.table("air_pollution.txt"))
airp.X <- airpollution[,-16]
airp.Y <- airpollution[,16]

#high dimensional
#二次の項(X_i*X_j)の追加
airp.X.low <- airp.X
airp.X <- matrix(0,ncol=135,nrow=nrow(airp.X))
airp.X[,1:15] <- airp.X.low
airp.X[,16:30] <- airp.X.low^2
for(i in 1:105){
	airp.X[,30+i] <- apply(airp.X[,combn(15,2)[,i]],1,prod)
}

#grid: 正則化の探す範囲(グリッドサーチ)
grid <- 10^seq(8,-3,length=1000)
#MSE
num.crossval <- 100
mse.mat <- matrix(0,nrow=num.crossval,ncol=2)

for(i in 1:num.crossval){
	train <- sample(nrow(airp.X),30) #ランダムに学習データ(30個)とテストデータ(30個)に分割
	test <- -train

	#ridge
	cv.out.ridge <- cv.glmnet(airp.X[train,],airp.Y[train],alpha=0,lambda=grid) #交差確認
	bestlam.ridge <- cv.out.ridge$lambda.min #最適な正則化項の値
	ridge.mod <- glmnet(airp.X[train,],airp.Y[train],alpha=0,lambda=bestlam.ridge)
	ridge.pred <- predict(ridge.mod,s=bestlam.ridge,newx=airp.X[test,]) #予測
	mse.mat[i,1] <- sqrt(mean((ridge.pred-airp.Y[test])^2))

	#lasso
	cv.out.lasso <- cv.glmnet(airp.X[train,],airp.Y[train],alpha=1,lambda=grid) #交差確認
	bestlam.lasso <- cv.out.lasso$lambda.min #最適な正則化項の値
	lasso.mod <- glmnet(airp.X[train,],airp.Y[train],alpha=1,lambda=bestlam.lasso)
	lasso.pred <- predict(lasso.mod,s=bestlam.lasso,newx=airp.X[test,]) #予測
	mse.mat[i,2] <- sqrt(mean((lasso.pred-airp.Y[test])^2))
}

#plot
boxplot(mse.mat,names=c("Ridge","Lasso"))


f:id:mkprob:20141224000800p:plain

Lassoの方が少し予測精度が悪いです.

*1:James, G., Witten, D., Hastie, T., Tibshirani, R. (2013). An Introduction to Statistical Learning: with Applications in R. Springer.

An Introduction to Statistical Learning: with Applications in R: 103 (Springer Texts in Statistics)

An Introduction to Statistical Learning: with Applications in R: 103 (Springer Texts in Statistics)

混合正規モデルに対するジェフリーズ事前分布

ある論文を読んでいたらジェフリーズ事前分布に基づく事後分布がimproperな例が載っていたので.


設定はシンプルです.

以下のような一次元パラメトリックモデル$\{p(x|\theta)|\theta\in\mathbf{R}\}$を考えます.

\[
p(x|\theta)=\frac{1}{2\sqrt{2\pi}} \exp\bigg(-\frac{x^2}{2}\bigg)+
\frac{1}{2\sqrt{2\pi}} \exp\bigg(-\frac{(x-\theta)^2}{2}\bigg).
\]

このモデルのジェフリーズ事前分布を求めます.
そのためにフィッシャー情報量を計算します.スコアは

\begin{eqnarray*}
\frac{\partial}{\partial\theta} \log p(x|\theta)=(x-\theta)-\frac{x-\theta}{1+\exp(\theta x-\theta^2/2)}.
\end{eqnarray*}

従ってフィッシャー情報量$I(\theta)$は
\begin{eqnarray*}
I(\theta)&=&\int_{-\infty}^{\infty}\bigg\{(x-\theta)-\frac{x-\theta}{1+\exp(\theta x-\theta^2/2)}\bigg\}^2p(x|\theta) {\rm d}x\\
&=&-\frac{\theta^2}{2}+\int_{-\infty}^{\infty}\frac{(x-\theta)^2\exp(-x^2/2)}{2\sqrt{2\pi}(1+\exp(\theta x-\theta^2/2))} {\rm d}x.
\end{eqnarray*}

最後の積分を計算できるのかよくわかりませんでしたので,数値的に計算させます.
ジェフリーズ事前分布$p(\theta)\propto \sqrt{I(\theta)}$は下のようになります.

f:id:mkprob:20140905211643j:plain

ソースは以下の通りです.ちょっとggplot2を使ってみました.

#Jeffreys prior for mixture normal model

library(ggplot2)

#Fisher情報量の計算の最後の積分の被積分関数
Integrand_Jeffreys <- function(x,theta){
  (x-theta)*(x-theta)*exp(-x^2/2)/(1+exp(theta*x-theta^2/2))*
    (1/(2*sqrt(2*pi)))
}

Jeffreys <- function(theta){
  FisherInfo <- sapply(theta,function(eta){integrate(function(x){Integrand_Jeffreys(x,eta)},
                          lower = -Inf,upper = Inf)$value - eta^2/2})
  sqrt(FisherInfo)
}

#plot
ggplot(data.frame(x=c(-2, 2)), aes(x)) + stat_function(fun=Jeffreys) +
  xlab(expression(theta)) + ylab("Probability Density")+
  ggtitle("Jeffreys Prior")+
  theme(axis.title.x=element_text(size=20),
        axis.title.y=element_text(size=20),
        axis.text.x=element_text(size=20),
        axis.text.y=element_text(size=20),
        plot.title=element_text(size=20))

このとき事後分布はimproperです.というのも

\[
p(x|\theta)>\frac{1}{2\sqrt{2\pi}}\exp\bigg(-\frac{x^2}{2}\bigg)>0, \quad
I(\theta) \geq 0.5,
\]
なので
\[
\int p(x|\theta)\sqrt{I(\theta)}{\rm d}\theta =\infty,\quad \forall{x},
\]
となるからです.

原点でとがっているように見えますが,積分微分の交換を認めれば,
$I(\theta)>0$は原点で微分可能です.

実際図を原点周りで拡大すると以下のようになります.

f:id:mkprob:20140905212137j:plain

ジェフリーズ事前分布の形について

どうしてこのような形(0で最小値をとり,逆ベル型)になるのか考えてみます.
(説明が正しいか保障はありませんが)

次のような通信を考えます.

送信者をA,受信者をBとします.Aはメッセージ$\theta$をある通信路を通してBに送ります.
この通信路は粗悪で,確率0.5で$\theta$に標準正規分布のノイズを足した値$x$をBに送り,
また確率0.5で全く関係ない標準正規分布からの乱数$x$をBに送り付けます.

フィッシャー情報量$I(\theta)$はBが$x$というメッセージを受け取った時に得られる
真のメッセージ$\theta$に関する情報量であると考えられます.

そう考えると,$\theta$が$0$に近い場合は,観測$x$は全くノイズであっても$\theta$にノイズが加わったものでも,$0$に近い値をとります.

従ってこの場合,$x$を観測したとき真のメッセージ$\theta$にノイズが乗ったものを観測したのか,それとも全くの雑音なのか判別がつきづらいです.よって観測$x$のもつ真のメッセージに関する情報量は小さいと考えられます.

しかし,$|\theta|$がある程度大きい場合には(たとえば$\theta=10$),$x$を観測しそれが大きい値の場合には$x$は$\theta$に近い値だなということになります.逆に$x$が$0$に近ければ雑音だなということがわかりやすくなります.

従って$|\theta|$が$0$に近いときと$|\theta|$が大きい時では観測$x$から得られる情報量に差があり,また$|\theta|$が大きくなればなるほど$\theta=0$のときに得られる情報量との差が大きくなる(判別しやすくなる)ことがわかります.

混合正規モデルではなく,平均が未知,分散が既知の正規分布の場合を考えますと,フィッシャー情報量が一定です.これは$\theta$がどんな値だろうと,観測$x$の周りに$\theta$がいるのは変わりないからだとも解釈できます(一般にlocation modelのlocation parameter).


$|\theta|=7$付近で$I(\theta)>0.7$となるのですが,$|\theta|=7$くらいから,もし$|x|<3$ならば観測$x$は雑音から,$|x|>3$ならば$x$は真のメッセージの付近,と思えるからでしょうか(自信がありませんが).



最後に参考文献を.この論文を読んでいたわけではないですが,上の話は下の論文に載っています.

Wasserman, L. (2000). Asymptotic inference for mixture models using data-dependent priors. JRSS-B, vol. 62, pp. 159-180.

この論文は客観ベイズの話(probability matching prior)です.客観ベイズも面白そうです.

Contrastive Divergence 法

Contrastive Divergence法 (Hinton, 2002)について少し勉強したので,そのまとめです.

Contrastive Divergence 法とは

(確率的な)最適化方法です.正確には正規化定数が分からない(求めるのが困難)確率分布のためのパラメータの最尤法です.特にBoltzmann分布(マルコフ確率場)における最尤法を指します.

何がうれしいか

正規化定数が分からない確率分布に対しても最尤推定量に近い推定量が得られることが利点です.
Boltzmann分布の例を挙げます.Boltzmann分布とは
\[
p(\boldsymbol{x};\boldsymbol{\theta},\boldsymbol{W})=\frac{1}{Z(\boldsymbol{\theta},\boldsymbol{W})}\exp(-\boldsymbol{\theta}^{\top}\boldsymbol{x}-\boldsymbol{x}^{\top}\boldsymbol{W}\boldsymbol{x})
\]
という確率分布です.ここで$\boldsymbol{x}=(x_1,\ldots,x_d)\in\{0,1\}^d$, $\boldsymbol{\theta}\in\mathbf{R}^d$, $\boldsymbol{W}\in\mathbf{R}^{d\times d}$であり,
\[
Z(\boldsymbol{\theta},\boldsymbol{W})=\sum_{\boldsymbol{x}\in\{0,1\}^d}\exp(-\boldsymbol{\theta}^{\top}\boldsymbol{x}-\boldsymbol{x}^{\top}\boldsymbol{W}\boldsymbol{x})
\]
は正規化定数です.また
\[
E(\boldsymbol{x};\boldsymbol{\theta},\boldsymbol{W})=\boldsymbol{\theta}^{\top}\boldsymbol{x}+\boldsymbol{x}^{\top}\boldsymbol{W}\boldsymbol{x}
\]
はエネルギー関数と呼ばれます.注意すべきは$d$が比較的多きときに$Z(\boldsymbol{\theta},\boldsymbol{W})$の計算は非常に大変であるということです.例えば$d=10$のときには$
Z(\boldsymbol{\theta},\boldsymbol{W})$を計算するのに$\exp(-\boldsymbol{\theta}^{\top}\boldsymbol{x}-\boldsymbol{x}^{\top}\boldsymbol{W}\boldsymbol{x})$を1024回足さなければなりません.$d=20$なら約100万回足さなければなりません.

Contrastive Divergence 法はどういう方法なのか

この確率分布からの標本$\boldsymbol{x}_1,\ldots,\boldsymbol{x}_n$が得られたときに,最尤推定量$\hat{\boldsymbol{\theta}}_{\rm mle}, \hat{\boldsymbol{W}}_{\rm mle}$を数値的に求めるには対数尤度 $\log p(\boldsymbol{x};\boldsymbol{\theta},\boldsymbol{W})$の勾配が必要です.例えば$\boldsymbol{\theta}=(\theta_1,\ldots,\theta_d)^{\top}$の勾配を計算すると
\[
\frac{\partial}{\partial \theta_i}\log p(\boldsymbol{x};\boldsymbol{\theta},\boldsymbol{W})=-x_i-\frac{\partial}{\partial \theta_i}\log Z(\boldsymbol{\theta},\boldsymbol{W})
\]
となり$Z(\boldsymbol{\theta},\boldsymbol{W})$の計算が必要となります.最尤推定量の計算には勾配法が基本的ですが,Boltzamnn分布の対数尤度の勾配を計算するのは困難であることが上の式からわかります.


そこでContrastive Divergence 法(CD法)が考え出されました(Hinton, 2002).Contrastive Divergence法の考え方はシンプルです.
表記の簡便のため,パラメータ$\boldsymbol{\xi}=(\boldsymbol{\theta},\boldsymbol{W})$とまとめておきます.$\xi_i$方向の勾配は
\[
\frac{\partial}{\partial \xi_i}\log p(\boldsymbol{x};\boldsymbol{\theta},\boldsymbol{W})=-\frac{\partial}{\partial \xi_i}E(\boldsymbol{x};\boldsymbol{\theta},\boldsymbol{W})-\frac{\partial}{\partial \xi_i}\log Z(\boldsymbol{\theta},\boldsymbol{W})
\]
移項して
\[
\frac{\partial}{\partial \xi_i}\log Z(\boldsymbol{\theta},\boldsymbol{W})=-\frac{\partial}{\partial \xi_i}\log p(\boldsymbol{x};\boldsymbol{\theta},\boldsymbol{W})-\frac{\partial}{\partial \xi_i}E(\boldsymbol{x};\boldsymbol{\theta},\boldsymbol{W})
\]
左辺は$\boldsymbol{x}$に依らないです.両辺を$\boldsymbol{x}$ ($\sim p(\boldsymbol{x};\boldsymbol{\theta},\boldsymbol{W})$)で期待値をとって
\[
\frac{\partial}{\partial \xi_i}\log Z(\boldsymbol{\theta},\boldsymbol{W})=-\mathrm{E}_{\tilde{\boldsymbol{x}}\sim p(\tilde{\boldsymbol{x}};\boldsymbol{\theta},\boldsymbol{W})}\bigg[\frac{\partial}{\partial \xi_i}E(\tilde{\boldsymbol{x}};\boldsymbol{\theta},\boldsymbol{W})\bigg]
\]
この関係式を代入して,勾配は
\[
\frac{\partial}{\partial \xi_i}\log p(\boldsymbol{x};\boldsymbol{\theta},\boldsymbol{W})=-\frac{\partial}{\partial \xi_i}E(\boldsymbol{x};\boldsymbol{\theta},\boldsymbol{W})+\mathrm{E}_{\tilde{\boldsymbol{x}}\sim p(\tilde{\boldsymbol{x}};\boldsymbol{\theta},\boldsymbol{W})}\bigg[\frac{\partial}{\partial \xi_i}E(\tilde{\boldsymbol{x}};\boldsymbol{\theta},\boldsymbol{W})\bigg]
\]
となります.このままでは$\mathrm{E}_{\tilde{\boldsymbol{x}}\sim p(\tilde{\boldsymbol{x}};\boldsymbol{\theta},\boldsymbol{W})}$の計算の際に$2^d$回の和が必要なので,計算が困難なことには変わりありません.
そこでCD法では期待値の計算をモンテカルロ近似します.つまり$p(\tilde{\boldsymbol{x}};\boldsymbol{\theta},\boldsymbol{W})$に従うサンプルを$k$個$(\tilde{\boldsymbol{x}}_1,\ldots,\tilde{\boldsymbol{x}}_k)$サンプリングし,
\[
\mathrm{E}_{\tilde{\boldsymbol{x}}\sim p(\tilde{\boldsymbol{x}};\boldsymbol{\theta},\boldsymbol{W})}\bigg[\frac{\partial}{\partial \xi_i}E(\tilde{\boldsymbol{x}};\boldsymbol{\theta},\boldsymbol{W})\bigg]\approx
\frac{1}{k}\sum_{m=1}^k \frac{\partial}{\partial \xi_i}E(\tilde{\boldsymbol{x}}_m;\boldsymbol{\theta},\boldsymbol{W})
\]
とします.これにより勾配の計算に必要な計算コストををおさえます.
サンプリング方法はマルコフ連鎖モンテカルロ法(MCMC),特にギブスサンプリングです.*1
CD法はMCMCが収束するように十分遷移させてサンプルするのではなく,少ない回数だけ遷移させます.
初期値は経験分布です.毎回初期値を経験分布にリセットするのはもったいない気がします.Persistent Cotrastive Divergenceという方法もあるみたいで,これはパラメータを更新するたびにMCMCの初期値を前のステップでMCMCで得た分布を使うものです.


まとめるとCD法は以下のような手順で最尤推定量を求める方法です.

  1. 初期値 $\boldsymbol{\xi}_0=(\boldsymbol{\theta}_0,\boldsymbol{W}_0)$ を決める
  2. 勾配 $\frac{\partial}{\partial \boldsymbol{\xi}}\log p(\boldsymbol{x};\boldsymbol{\theta}_k,\boldsymbol{W}_k)$ の近似値をモンテカルロ法で計算

(初期値は経験分布,MCMCの遷移回数は少ない)

  1. $\boldsymbol{\xi}_{k+1} \leftarrow \boldsymbol{\xi}_{k} + \gamma_k \times \frac{\partial}{\partial \boldsymbol{\xi}}\log p(\boldsymbol{x};\boldsymbol{\theta}_k,\boldsymbol{W}_k)$ ($\gamma_k$: ステップ幅)
  2. 停止条件を満たしてなければステップ2へ

少し考えると,この方法では勾配の計算の際,少ないMCMCの遷移回数で近似することにしています.
単純に考えると,MCMCは十分多く遷移させると所望の確率分布に収束する方法なので少数の遷移回数でよいのかは疑問です.
しかし,実は近似のための遷移回数が少なくてもCD法はうまく動くということが経験的に知られています.

自分としては近似のための遷移回数が少なくてもCD法はうまく動くというのが不思議だったので,少し勉強&検証してみました.

数値実験

まずは数値実験をして検証してみました.

設定としては$d=3$の小さい次元の場合で検証しました.サンプルサイズ$10,100,200$ごとに100回
$\boldsymbol{\theta}=(\theta_1,\theta_2,\theta_3)$を平均3,分散1の正規分布から,
$\boldsymbol{W}$の上三角かつ非対角成分$(W_{12},W_{13},W_{23})$を平均-3,分散1の正規分布から,
発生させ,それを真の値とするBoltzmann分布からサンプルサイズ分だけ標本を得ます($\boldsymbol{W}$の他のパラメータの値は全て0というモデルを考えています).
各回ごとに,最尤推定量,およびCD1,CD3を計算します.
ここで,CD1,CD3は共にCD法で,それぞれモンテカルロ近似の際に用いた遷移の数が1または3であることを意味します.
CD1,CD3で得られた推定量と最尤推定量との距離を平均二乗誤差$\|\hat{\boldsymbol{\xi}}_{\rm mle}-\hat{\boldsymbol{\xi}}_{\rm cd}\|_2$で計算しました
(各サンプルサイズごとに100回CD1とCD3の平均二乗誤差が計算されます).

結果は以下の通りです.コードは最後に回します.


f:id:mkprob:20140720031109j:plain

f:id:mkprob:20140720031115j:plain

この結果からCD3の方がCD1よりも最尤推定量に近い推定量を得られることが分かります.
CD3の方が最尤推定量からのばらつきが小さいです.
またサンプルサイズが大きくなるにつれてCD1もCD3も最尤推定量との距離が近くなっています.
個人的にはCD3でも結構最尤推定量に近いのが意外でした.


なぜ少ないサンプルでも良いのか?

やはり経験的に知られているように,今回の数値実験でも
少ない遷移回数による勾配の近似計算でも最尤推定量に近い推定量が得られました.

調べている内に次の文献にたどり着きました:

S. Akaho and K. Takabatake. (2008).
Information geometry of contrastive divergence.
International Conference on Information theory and statistical learning.

情報幾何では確率分布を多様体の一点としてとらえます.
沢山の遷移させなければMCMCで得られる確率分布(1点)はモデルの確率分布(1点)には収束しないですが,
経験分布からMCMCを数ステップ回して得られる確率分布への"方向"が
経験分布からモデル(多様体)への射影の方向に近いためうまくいくということになります.


CD法ではないですが

K. Takabatake. (2004).
Information geometry of Gibbs sampler.
Proceedings of the WSEAS International Conferences.

も面白いです.Gibbs samplerが貪欲的に経験分布からの距離を最小化しているというのは知りませんでした.

コード

#contrastive divergence
#maximum likelihood estimation
#Boltzmann distribution

#######################
#integer to binary
dec2bin <- function(num, digit=0){
  if(num <= 0 && digit <= 0){
    return(NULL)
  }else{
    return(append(Recall(num%/%2,digit-1), num%%2))
  }
}


########################
#Boltzmann distribution
energy_boltz <- function(state, theta, W){
	#state: binary vector (length=dim(theta))
	#Energy function: Energy(state)=state%*%theta+state%*%W.mat%*%state
	#prob of state: p(state)=exp(-Energy(state))/Z
	#Z: normalization const
	#W: W[1]=W_{12}, W[2]=W_{23},...
	
	vardim <- length(theta)
	W.mat <- matrix(0, ncol=vardim, vardim)
	W.mat[upper.tri(W.mat)] <- W
	exp(-state%*%theta-state%*%W.mat%*%state)
}

dboltz <- function(theta, W){
	#distribution function of Boltzmann 
	#p(x|theta,W) = exp(-Energy(state))/Z

	vardim <- length(theta)
	states <- 2^vardim
	probs <- sapply(0:(states-1),dec2bin,digit=vardim)
	probs <- rbind(probs, apply(probs,2,energy_boltz,theta=theta,W=W))
	probs[(vardim+1),] <- probs[(vardim+1),]/sum(probs[vardim+1,]) #normalization
	t(probs)
}

rboltz <- function(n, theta, W){
	#random n samples from Boltzmann

	probs <- dboltz(theta, W)
	vardim <- length(theta)
	states <- 2^vardim
	samples <- matrix(0, ncol=vardim, nrow=n)
	sample10 <- sample(x=0:(states-1),size=n,prob=probs[,vardim+1],replace=T) #sampling
	samples <- sapply(sample10,dec2bin,digit=vardim) #integer to binary
	t(samples)
}


########################
#Maximum Likelihood Estimation
#gradient ascent
likegrad <- function(samples,theta,W,gradtheta.emp,gradW.emp){
	thetadim <- length(theta)
	Wdim <- length(W)
	offdiag <- t(combn(thetadim,2)) #off-diagonal index of W
	gradtheta <- rep(0, thetadim)
	gradW <- rep(0, Wdim)
	probs <- dboltz(theta,W)
	gradtheta <- -gradtheta.emp+apply(probs[,1:thetadim]*probs[,thetadim+1],2,sum) #gradient of theta
	gradW.model <- sapply(1:Wdim, function(i) sum(apply(probs[,c(offdiag[i,],(thetadim+1))],1,prod))) #mean of gradient of W 
	gradW <- -gradW.emp + gradW.model
	return(c(gradtheta, gradW))
}

findmle <- function(samples,theta.ini,W.ini,lrate=1,maxiter=100){
	pars <- c(theta.ini, W.ini)
	thetadim <- length(theta.ini)
	offdiag <- t(combn(thetadim,2))
	Wdim <- length(W.ini)
	theta <- theta.ini
	W <- W.ini
	gradtheta.emp <- apply(samples,2,mean) #empirical mean of gradient of theta
	gradW.emp <- sapply(1:Wdim, function(i) mean(apply(samples[,offdiag[i,]],1,prod))) #empirical mean of gradient of W

	for(i in 1:maxiter){
		pars <- pars + likegrad(samples,theta,W,gradtheta.emp,gradW.emp)*lrate #gradient ascent update
		theta <- head(pars,thetadim)
		W <- head(pars,Wdim)
	}
	rtrn.ls <- list(theta, W)
	names(rtrn.ls) <- c("theta.mle", "W.mle")
	return(rtrn.ls)
}


########################
#Contrastive Divergence Method
gibbsupdate <- function(state,theta,W,mcmciter){
	dimstate <- length(state)
	state.now <- state

        #Gibbs sampler update
        #update only one component of state
	for(i in 1:mcmciter){
		state.new <- state.now
		updateindex <- sample(dimstate,1)
		state.new[updateindex] <- 1-state.now[updateindex]
		accept <- energy_boltz(state.new,theta,W)
		reject <- energy_boltz(state.now,theta,W)
		acceptflag <- sample(x=1:0,size=1,prob=c(accept,reject))
		state.now <- if(acceptflag) state.new else state.now
	}
	state.now
}

gibbslikegrad <- function(samples,theta,W,gradtheta.emp,gradW.emp,mcmciter){
	thetadim <- length(theta)
	Wdim <- length(W)
	offdiag <- t(combn(thetadim,2))
	gradtheta <- rep(0, thetadim)
	gradW <- rep(0, Wdim)
	mcmcsamples <- t(apply(samples,1,function(state)gibbsupdate(state,theta,W,mcmciter))) #samples from Gibbs sampler
	gradtheta <- -gradtheta.emp+apply(mcmcsamples,2,mean)
	gradW.model <- sapply(1:Wdim, function(i) mean(apply(mcmcsamples[,offdiag[i,]],1,prod)))
	gradW <- -gradW.emp + gradW.model
	return(c(gradtheta, gradW))
}

findCD <- function(samples,theta.ini,W.ini,lrate=1,maxiter=100,mcmciter=1){
	pars <- c(theta.ini, W.ini)
	thetadim <- length(theta.ini)
	offdiag <- t(combn(thetadim,2))
	Wdim <- length(W.ini)
	theta <- theta.ini
	W <- W.ini
	gradtheta.emp <- apply(samples,2,mean)
	gradW.emp <- sapply(1:Wdim, function(i) mean(apply(samples[,offdiag[i,]],1,prod)))

	for(i in 1:maxiter){
		pars <- pars + gibbslikegrad(samples,theta,W,gradtheta.emp,gradW.emp,mcmciter)*lrate #gradient ascent
		theta <- head(pars,thetadim)
		W <- tail(pars,Wdim)
	}
	rtrn.ls <- list(theta, W)
	names(rtrn.ls) <- c("theta.cd", "W.cd")
	return(rtrn.ls)
}


mle_to_cd <- function(n, n.sample, dim){
	print(n)
	dimW <- (dim-1)*dim/2
	risk <- rep(0, 2)
	theta.true <- rnorm(dim, mean=3)
	W.true <- rnorm(dimW, mean=-3)
	par.true <- c(theta.true, W.true)
	names(par.true) <- c("theta.true", "W.true")
	samples <- rboltz(n.sample,theta.true,W.true)
	result.mle <- findmle(samples,theta.ini=rep(0,dim),W.ini=rep(0,dimW),maxiter=500)
	result.cd1 <- findCD(samples,theta.ini=rep(0,dim),W.ini=rep(0,dimW),maxiter=500)
	result.cddim <- findCD(samples,theta.ini=rep(0,dim),W.ini=rep(0,dimW),maxiter=500,mcmciter=dim)

	#Mean Squared Errors
	risk[1] <- sqrt(sum((unlist(result.mle)-unlist(result.cd1))^2)/(dim+dimW))
	risk[2] <- sqrt(sum((unlist(result.mle)-unlist(result.cddim))^2)/(dim+dimW))

	risk
}

*1:MCMCとCD法は全く別物です.MCMCは確率分布からのサンプリング方法であり,CD法はパラメータの推定方法です.

再生核ヒルベルト空間(その2)

Reproducing Kernel Hilbert Spaces in Probability and Statistics

Reproducing Kernel Hilbert Spaces in Probability and Statistics

の勉強まとめ(4章)です.

この章の主な目的は(符号付)測度$\mu$を
\[
\mu\mapsto \int K(\cdot, x)\mathrm{d}\mu(x)\in\mathcal{H}_K,
\]
によって再生核ヒルベルト空間に埋め込むことです.

特に$\mu$が確率測度のとき,
\[
\int K(\cdot, x)\mathrm{d}\mu(x),
\]
はkernel meanとよばれ,統計・機械学習で応用があります*1*2

また$x_1,\cdots, x_n$を確率分布$P$からの独立な標本するとき経験分布
\[
\frac{1}{n}\sum_{i=1}^n\delta_{x_i},
\]
を再生核ヒルベルト空間に埋め込むことで,$\mathcal{H}_K$上の確率分布を誘導します.この経験分布は$\mathcal{H}_K$上のGaussian measureに弱収束することなどもわかります.収束レートなども調べることができます(Berry--Essen).


まずは測度をヒルベルト空間に埋め込むために,ヒルベルト空間値変数の積分の存在条件を考える必要があります.弱・強可積分性が存在します.

定義(Pettis integrability)
$X:(\Omega, \mathcal{F}, P)\to (\mathcal{H}, \mathcal{B}(\mathcal{H}))$を$\mathcal{H}$-値確率変数とする.$A\in\mathcal{A}$に対し
\[
\bigg|\int_A \langle X,f\rangle_{\mathcal{H}} \mathrm{d}P\bigg| < \infty,\quad \forall{f}\in\mathcal{H},
\]
が成立するとき,$X$は$A$上,weakly integrableという.$X$が$\Omega$上でweak integrableのとき,$X$をPettis integrableという.またこの時 $\exists{x}_A\in\mathcal{H}$で
\[
\int_A \langle X,f\rangle_{\mathcal{H}} \mathrm{d}P =\langle x_A, f\rangle,
\]
を満たすものが存在する*3.この$x_A$を$X$のPettis integralといい
\[
x_A\stackrel{\rm def}{=}\oint_A X\mathrm{d}P,
\]
と表記する.このとき
\[
\int_A \langle X, f\rangle \mathrm{d}P=\bigg\langle \oint_A X\mathrm{d}P, f\bigg\rangle
\]

定義(Bochner integrability)
$X:(\Omega, \mathcal{F}, P)\to (\mathcal{H}, \mathcal{B}(\mathcal{H}))$を$\mathcal{H}$-値確率変数とする.$A\in\mathcal{A}$に対し
\[
\int_A \|X\|_{\mathcal{H}} \mathrm{d}P< \infty,
\]
が成立するとき,$X$は$A$上,strongly integrableという.$X$が$\Omega$上strongly integrableであるとき$X$をBochner integrableという.


$X:(\Omega,\mathcal{F},P)\to (E,\mathcal{A})$を確率変数とし$\mu$を$E$上の$X$の分布とすると
\[
\int_E K(x,x)\mathrm{d}\mu(x)<\infty,
\]
であれば$\|K(\cdot, x)\|^2=K(x,x)$であり,$\|K(\cdot, x)\|\in L^2(\mu)\subset L^1(\mu)$なのでBochner integrable.よってkernel meanが存在する.

kernel meanの性質として任意の$f\in\mathcal{H}_K$に対し
\[
E[f(X)]=\int_E f(x)\mathrm{d}\mu(x)=\langle f, \mu_X\rangle\quad (\mu_X\stackrel{\rm def}{=}\int_E K(\cdot, x)\mathrm{d}\mu(x)),
\]
が成立します.定義から明らかですが,積分を和だと思えば内積積分が交換できるのは直感的納得がいきます.


またこの埋め込みにより,$E$上の符号付測度の空間に,特に$E$上の確率分布に内積が定義できます*4:
\[
\forall{\mu},\forall{\nu},\quad \langle\langle\mu,\nu\rangle\rangle\stackrel{\rm def}{=}\int_{E^2} K(x,y)\mathrm{d}(\mu\otimes\nu)(x,y).
\]
内積が定義できるので,位相が入るわけですが,面白いのはその位相が測度の弱収束位相と一致する場合があるということです.
たとえば,$E$がコンパクト距離空間で$K$が連続,かつ測度同士の内積が上のように定まるとき,一致します.

*1: T. Hofmann, B. Schölkopf, A. Smola. (2008). Kernel methods in machine learning. The Annals of Statistics, 36(3):1171–1220.

*2:Kernel embedding of distributions. In Wikipedia, http://en.wikipedia.org/wiki/Kernel_embedding_of_distributions

*3:実際 $T_X:\mathcal{H}\ni f\mapsto \langle X, f\rangle \in L^1(A, P)$は閉作要素であり,閉グラフ定理から連続である.何故ならば$f_n\to f$ in $\mathcal{H}$,$T_X(f_n)\to \exists{g}\in L^1(A, P)$とすると,$\{f_n\}$の部分列で$T_X(f_n)\to g$ a.e. となるものが取れる.$T_X(f_n)=\langle X, f_n\rangle \to \langle X, f\rangle~\forall{\omega}\in\Omega$であるから,$g=\langle X, f\rangle$ a.e. となる.従って$T_X(f_n)\to T_X(f)$ in $L^1(A, P)$である.$T_X$は閉作用素であり,閉グラフ定理から連続である. 以上より \[ E_X:\mathcal{H}\ni f\mapsto \int_A \langle X, f\rangle \mathrm{d}P \in\mathbf{R}, \] は連続かつ線形な写像である.リースの表現定理から$\exists{x}_A\in\mathcal{H}$で \[ \int_A \langle X,f\rangle_{\mathcal{H}} \mathrm{d}P =\langle x_A, f\rangle, \] を満たすものが存在する.

*4:内積となるためには埋め込みが単射になる必要があります.

再生核ヒルベルト空間(演習その1)

演習問題のまとめその1(第1章)です.

[9]

$K$を$E\times E$上の正定値関数とする.$E$上の関数$f$が$K$を再生核に持つRKHS$\mathcal{H}_K$の元であることと,ある正の定数$\lambda$が存在して$K(s,t)-\lambda f(s)\overline{f(t)}$が正定値関数であることは同値であることを示せ.

[証明]
($\Rightarrow$)
$0<\lambda <1/||f||$を固定する.
任意の$(a_1,\ldots, a_n)\in\mathbf{C}^n$と$(t_1,\ldots,t_n)\in E^n$に対して再生核の性質とコーシー・シュワルツの不等式から
\begin{eqnarray*}
&&\sum_{i=1}^n\sum_{j=1}^n a_i\overline{a_j}(K(t_i,t_j)-\lambda f(t_i)\overline{f(t_j)})\\
&=&\sum_{i,j} a_i\overline{a_j}\langle K(\cdot, t_j), K(\cdot, t_i)\rangle -\sum_{i,j}(\sqrt{\lambda}a_if(t_i))\overline{(\sqrt{\lambda}a_jf(t_j))}\\
&=&\bigg\|\sum_i \overline{a_i}K(\cdot, t_i)\bigg\|^2-\bigg|\sum_i \sqrt{\lambda}a_if(t_i) \bigg|^2\\
&=&\bigg\|\sum_i \overline{a_i}K(\cdot, t_i)\bigg\|^2-\bigg|\sum_i \sqrt{\lambda}a_i\langle f, K(\cdot, t_i)\rangle \bigg|^2\\
&=&\bigg\|\sum_i \overline{a_i}K(\cdot, t_i)\bigg\|^2-\bigg|\langle f, \sum_i \sqrt{\lambda}\overline{a_i}K(\cdot, t_i)\rangle \bigg|^2\\
&\geq & \bigg\|\sum_i \overline{a_i}K(\cdot, t_i)\bigg\|^2-\lambda \|f\|\bigg\|\sum_i \overline{a_i}K(\cdot, t_i)\bigg\|^2\\
&\geq &0.
\end{eqnarray*}

($\Leftarrow$)
$K_f(s,t)\stackrel{\rm def}{=}f(s)\overline{f(t)}$とする.明らかに$K_f$は正定値であり,$K_f$を再生核に持つRKHSは$\{\alpha f|\alpha \in \mathbf{C}\}$に等しい.何故ならば$\mathcal{H}_0\stackrel{\rm def}{=}\{\sum_i a_iK_f(\cdot, t_i)|a_i\in\mathbf{C},t_i\in E\}=\{af|a\in \mathbf{C}\}$とすると$\mathcal{H}_0$は完備である.実際$\{u_n=\sum_{i=1}^{k_n}a_i^{(n)} f\stackrel{\rm def}{=}\alpha_n f\}\subset\mathcal{H}_0$をコーシー列とすると,$\|u_m-u_n\|_{\mathcal{H}_0}=|\sum_{j=1}^{k_m}a_j^{(m)}-\sum_{i=1}^{k_n}a_i^{(n)}|=|\alpha_m-\alpha_n|$.$\{u_m\}$はコーシー列であり$\mathbf{C}$は完備なので,$\alpha_n \to\exists{\alpha} \in \mathbf{C}$.$u\stackrel{\rm def}{=}\alpha f\in\mathcal{H}_0$とすると$||u_n-u||_{\mathcal{H}_0}=|\alpha_n-\alpha|\to 0$より$\mathcal{H}_0$は完備.これより$\bar{\mathcal{H}}_0=\mathcal{H}_{K_f}=\{\alpha f|\alpha\in\mathbf{C}\}$.

ここで仮定から$\frac{1}{\lambda}K-K_f$は正定値であるから,Aronszajn(1950)の結果(p30, THEOREM12)から$\mathcal{H}_{K_f}\subset \mathcal{H}_{K}$.よって$f\in\mathcal{H}_K$.(終)


[10]
$l^2(\mathbf{N})$(数列空間)と$L^2(0,1)$を考えることによって,RKHSという性質ヒルベルト空間の同型に関して不変ではないことを示せ.

[証明]
$l^2(\mathbf{N})$と$L^2(0,1)$は共に可分なヒルベルト空間であるから同型である(互いの完全正規直交系を$\{e_i\}$と$\{f_i\}$とすれば$e_i\mapsto f_i$が同型を与える).

$l^2(\mathbf{N})$はRKHSである.実際$K(i,j)\stackrel{\rm def}{=}\delta_{ij}$(クロネッカーのデルタ)とすれば再生核の定義を満たす.

一方$L^2(0,1)$はRKHSでない.これを背理法で示す.$L^2(0,1)$がRKHSとすると,再生核$K$が存在する.再生核の定義から
\[
\forall{\varphi}\in L^2(0,1),\quad \int_0^1 \varphi(x){\overline{K(x, y)}}\mathrm{d}x=\varphi(y).
\]
この積分作用素$\mathcal{K}\varphi(y)\stackrel{\rm def}{=}\int_0^1 \varphi(x){\overline{K(x, y)}}\mathrm{d}x$はコンパクト作用素であり(関数解析の本参照),再生核の定義から恒等作用素ということになる.このことから$L^2(0,1)$は有限次元でなければならない.何故ならばコンパクト作用素でかつ恒等作用素が存在することから,$L^2(0,1)$の単位閉球はコンパクト,つまり$L^2(0,1)$は局所コンパクトである.局所コンパクトなバナッハ空間は有限次元であるから$L^2(0,1)$は有限次元である.これは$L^2(0,1)$が無限次元であることに矛盾する.よって$L^2(0,1)$はRKHSでない.(終)


[11]
[任意のヒルベルト空間はあるRKHSに等長同型である]

$\mathcal{E}$をヒルベルト空間とし,内積を$K(\alpha, \beta)\stackrel{\rm def}{=}\langle \alpha, \beta\rangle_{\mathcal{E}}$とする.このとき明らかに$K$は$\mathcal{E}\times\mathcal{E}$上の正定値な関数である.従って$K$を再生核に持つRKHS$\mathcal{H}_K \subset \mathbf{C}^{\mathcal{E}}$が存在する.このとき次の二つが成立する:

  1. $\mathcal{H}_K=\{K(\cdot,\beta)|\beta\in\mathcal{E}\}$
  2. $\mathcal{E}$と$\mathcal{H}_K$は等長同型である.

[証明]
(1)
$\mathcal{H}_0\stackrel{\rm def}{=}\bigg\{\sum_{i=1}^n a_iK(\cdot,\beta_i) \bigg| a_i\in\mathbf{C},\beta_i\in\mathcal{E}, n\in\mathbf{N}\bigg\}$とすると,
RKHSの性質から$\bar{\mathcal{H}}_0$($\mathcal{H}_0$の$\mathcal{H}_K$での閉包)は$\mathcal{H}_K$に等しい.

内積の線形性から
\begin{eqnarray*}
\sum_{i=1}^n a_iK(\cdot,\beta_i)=K(\cdot, \sum_{i=1}^n\bar{a}_i\beta_i),
\end{eqnarray*}
より$\mathcal{H}_0\subset \{K(\cdot,\beta)|\beta\in\mathcal{E}\}$.

つぎに$\{K(\cdot,\beta)|\beta\in\mathcal{E}\}$が閉であることを示す.点列$\{u_n\}\subset \{K(\cdot,\beta)|\beta\in\mathcal{E}\}$で$u_n\to u$ in $\mathcal{H}_K$とする.このときある点列$\beta_n\in \mathcal{E}$で$u_n=K(\cdot, \beta_n)$となるものが存在する.内積の線形性と再生核の定義から

\begin{eqnarray*}
||u_m-u_n||_{\mathcal{H}_K}&=&||K(\cdot, \beta_m)-K(\cdot, \beta_n)||_{\mathcal{H}_K}\\
&=&||K(\cdot, \beta_m-\beta_n)||_{\mathcal{H}_K}\\
&=&[K(\beta_m-\beta_n,\beta_m-\beta_n)]^{1/2}\\
&=&||\beta_m-\beta_n||_{\mathcal{E}}.
\end{eqnarray*}

$\{u_n\}$は$\mathcal{H}_K$内のコーシー列であるから,$\{\beta_n\}$は$\mathcal{E}$内のコーシー列である.$\mathcal{E}$は完備であるから$\{\beta_m\}$は収束列であり,$\exists{\beta}\in\mathcal{E}$ s.t. $\beta_n\to\beta$ in $\mathcal{E}$.

任意の$\varepsilon>0$に対し,ある$n_0$が存在して$n\geq n_0$ならば
\begin{eqnarray*}
||u_n-u||_{\mathcal{H}_K}<\varepsilon/2, \quad ||\beta_n-\beta||_{\mathcal{E}}<\varepsilon/2,
\end{eqnarray*}
となる.よってこの$n$に対し
\begin{eqnarray*}
||u-K(\cdot, \beta)||_{\mathcal{H}_K}&\leq &||u-u_n||_{\mathcal{H}_K}+||u_n-K(\cdot, \beta)||_{\mathcal{H}_K}\\
&=&||u-u_n||_{\mathcal{H}_K}+||K(\cdot, \beta_n-\beta)||_{\mathcal{H}_K}\\
&=&||u-u_n||_{\mathcal{H}_K}+||\beta_n-\beta||_{\mathcal{E}}\\
&<& \varepsilon.
\end{eqnarray*}
よって$u=K(\cdot,\beta)$.以上から$\mathcal{H}_K=\bar{\mathcal{H}}_0\subset \{K(\cdot,\beta)|\beta\in\mathcal{E}\}$.
再生核の定義から明らかに$\mathcal{H}_K\supset \{K(\cdot,\beta)|\beta\in\mathcal{E}\}$なので
\[
\mathcal{H}_K=\{K(\cdot,\beta)|\beta\in\mathcal{E}\}.
\]

(2)
(1)より$\Psi:\mathcal{E}\to \mathcal{H}_K$を
\[
\Psi(\beta)=K(\cdot, \beta) \quad (\beta\in \mathcal{E})
\]
と定める.このとき
\[
||\beta||_{\mathcal{E}}=[K(\beta,\beta)]^{1/2}=||K(\cdot,\beta)||_{\mathcal{H}_K}=||\Psi(\beta)||_{\mathcal{H}_K}.
\]
明らかに$\Psi$は線形で全射なので,$\Psi$は$\mathcal{E}$と$\mathcal{H}_K$の間の等長同型を与える.(終)

この問題より任意のヒルベルト空間はあるRKHSに等長同型である.


[12]
[再生核の性質から再生核ヒルベルト空間の元の性質が分かる例]

$K$が$E\times E$の対角線上で有界ならば,$K$を再生核に持つRKHS$\mathcal{H}_K$の任意の元は$E$上有界な関数であることを示せ.また$\mathcal{H}_K$の点列がノルムの意味で収束すれば,一様ノルムの意味でも収束する(つまり一様収束)ことを示せ.

[証明]
$K$が$E\times E$の対角線上で有界なので
\[
\sup_{t\in E} K(t,t)\stackrel{\rm def}{=}M<\infty.
\]
$\forall{\varphi}\in\mathcal{H}_K$に対し
\[
|\varphi(t)|=|\langle \varphi, K(\cdot,t)\rangle|\leq ||\varphi||[K(t,t)]^{1/2}\leq \sqrt{M}||\varphi||,
\]
なので
\[
\sup_{t\in E} |\varphi(t)|<\infty.
\]
よって有界.また$\{\varphi_n\}\subset\mathcal{H}_K$を$\varphi\in \mathcal{H}_K$に収束する収束列とすると,上と同様にして
\[
\sup_{t}|\varphi_n(t)-\varphi(t)|\leq \sqrt{M}||\varphi_n-\varphi||
\]
なので,$\varphi_n\to \varphi$ uniformly on $E$となる.(終)

再生核ヒルベルト空間(まとめその1)

Reproducing Kernel Hilbert Spaces in Probability and Statistics

Reproducing Kernel Hilbert Spaces in Probability and Statistics

の勉強(Exercise)のまとめその1です.


定義(再生核ヒルベルト空間)
$E$を集合としその上の複素数値関数の全体を$\mathbf{C}^E$とする.
\begin{eqnarray*}
K:E\times E&\to &\mathbf{C}\\
(s,t)&\mapsto &K(s,t)
\end{eqnarray*}
ヒルベルト空間$\mathcal{H}\subset \mathbf{C}^E$の再生核(reproducing kernel)であるとは
\begin{eqnarray*}
&{\rm i)}&~ \forall{t}\in E,\quad K(\cdot,t)\in\mathcal{H}\\
&{\rm ii)}&~ \forall{t}\in E,\quad \forall{\varphi}\in\mathcal{H}, \quad \langle \varphi, K(\cdot,t)\rangle_{\mathcal{H}}=\varphi(t)
\end{eqnarray*}
を満たすことをいう.また再生核$K$をもつヒルベルト空間$\mathcal{H}_K$を再生核ヒルベルト空間(Reproducing Kernel Hilbert Space, RKHS)という.


RKHSの性質
RKHSの定義自体はシンプルですが,色々な数学的な性質が導けます.

[性質1]
$u_n \to u$ in $\mathcal{H}_K$ $\Rightarrow$ $u_n(t)\to u(t)\quad (\forall{t}\in E)$
(ノルムの意味での収束が各点収束を意味します)

[証明]
$|u_n(t)-u(t)|=|\langle u_n-u, K(\cdot, t)\rangle_{\mathcal{H}_K}|\leq ||u_n-u||_{\mathcal{H}_K}||K(\cdot, t)||_{\mathcal{H}_K}$,
で$||K(\cdot,t)||_{\mathcal{H}_K}=[K(t,t)]^{1/2}<\infty$なので各点収束する.(終)


[性質2]
部分空間$\mathcal{H}_0\subset \mathcal{H}_K$を
\begin{eqnarray*}
\mathcal{H}_0\stackrel{\rm def}{=}\bigg\{ \sum_{i=1}^n a_iK(\cdot, t_i)\in \mathcal{H}_K \bigg| (a_1,\ldots,a_n)\in\mathbf{C}^n, (t_1,\ldots, t_n)\in E^n, n\in \mathbf{N}\bigg\}
\end{eqnarray*}
とする.つまり$\mathcal{H}_0$を$\{K(\cdot, t),t\in E\}$の張る線形空間とする.

このとき$\mathcal{H}_0$は$\mathcal{H}_K$の中で稠密である.つまり$\forall{\varphi}\in\mathcal{H}_K$は$K(\cdot, t)$の線形和で近似できる.


[性質3]
再生核$K$は正定値な関数である.逆に正定値な関数$K$があれば,$K$を再生核に持つRKHS$\mathcal{H}_K$が一意に存在する.
つまり
\begin{eqnarray*}
K({\rm positive~ definite}) \leftrightarrow \mathcal{H}_K
\end{eqnarray*}
は一対一対応である.


これら以外にも$K$と$\mathcal{H}_K$の関係がいろいろと存在します.つまり$K$が何らかの条件を満たすと$\mathcal{H}_K$の元である$f$についてもある程度の性質がわかります(たとえば$K$の連続性及び有界性と$f$の連続性など).