Rで等方的ガウスカーネルを用いたNadaraya-Watsonモデルによる回帰

PRML図6.3と同様に、Nadaraya-Watsonモデル（カーネル回帰）により、以下を図示します。

• 回帰関数（条件付き期待値） $E[t|x]=\sum_{n=1}^Nk(x,x_n)t_n$
• 条件付き確率密度 $p(t|x)=\frac{1}{(2\pi\sigma^2)^{1/2}}\sum_{n=1}^{N}\exp(-\frac{(t-t_n)^2}{2\sigma^2})k(x,x_n)$
• 条件付き分散 $Var(t|x)=\sum_{n=1}^N(\sigma^2+t_n^2)k(x,x_n)-E[t|x]^2$

また、各観測点に関連付けられたカーネル関数 $k(x,x_n)$ も図示します。

なお、条件付き確率密度、条件付き分散は、演習6.18と同じく、分散を$\sigma^2$とする等方的ガウス分布 $f(x,t)$ によるカーネル $k(x,x_n)$ を用いて求めています。

f(x,t)=N((x,t)^T|\mathbf{0},\sigma^2\mathbf{I})=\frac{1}{2\pi\sigma^2}\exp(-\frac{x^2+t^2}{2\sigma^2}) \\
g(x)=\int f(x,t) dt \\
k(x,x_n)=\frac{g(x-x_n)}{\sum_m g(x-x_m)}

frame()
set.seed(0)
par(mfrow=c(3, 2))
par(mar=c(2, 2, 1, 0.1))
par(mgp=c(1, 0.2, 0))
xrange <- c(-0.1, 1.1)
yrange <- c(-1.5, 1.5)
xcell <- seq(xrange[1], xrange[2], .01)
ycell <- seq(yrange[1], yrange[2], .01)
colors <- rainbow(450)[256:1]
N <- 10
x <- seq(0, 1, 1 / (N - 1))
# y <- rnorm(length(x), sin(2 * pi * x), 0.2)
y <- c(0.4, 0.8, 1.1, 1.0, 0, 0.1, -1, -0.5, -0.6, 0.3)  

doplot <- function(sigma) {
	
	# g(x)
	g <- function(x) {
		1 / sqrt(2 * pi * sigma ^ 2) * exp(-x ^ 2 / (2 * sigma ^ 2))
	}
	
	# k(x, x_n)
	kernel <- Vectorize(function(xx, xn) {
		g(xx - xn) / sum(g(xx - x))
	})
	
	# E[t|x]
	estimate <- Vectorize(function(xx) {
		sum(kernel(xx, x) * y)
	})
	
	# p(t|x)
	p <- Vectorize(function(xx, yy) {
		1 / sqrt(2 * pi * sigma ^ 2) * sum(exp(-(yy - y) ^ 2 / (2 * sigma ^ 2)) * kernel(xx, x))
	})
	
	# Var(t|x)
	variance <- Vectorize(function(xx) {
		sum((sigma ^ 2 + y ^ 2) * kernel(xx, x)) - estimate(xx) ^ 2
	})

	# k(x, xn)の描画
	for (n in 1:N) {
		if (n > 1) {
			par(new=T)
		}
		xn <- x[n]
		curve(kernel(x, xn), xlim=c(-0.1, 1.1), ylim=c(-0.1, 1.1), 
			xlab=ifelse(n == 1, "x", ""), ylab="", axes=(n == 1), col=n)
	}
	title(bquote(paste(k(x, x[n]), ~ sigma, "=", .(sigma))))

	# p(t|x)の描画
	plot(x, y, xlim=xrange, ylim=yrange, xlab="x", ylab="t")
	image(xcell, ycell, outer(xcell, ycell, p), xlim=xrange, ylim=yrange, col=colors, axes=F, add=T)
	
	# 観測点の描画
	par(new=T)
	plot(x, y, xlim=xrange, ylim=yrange, xlab="", ylab="", axes=F)
	
	# E[t|x], Var(t|x)の描画
	par(new=T)
	curve(estimate, xlim=xrange, ylim=yrange, xlab="", ylab="", axes=F, col=2)
	hi <- function(xx){ estimate(xx) + sqrt(variance(xx)) * 2 }
	lo <- function(xx){ estimate(xx) - sqrt(variance(xx)) * 2 }
	par(new=T)
	curve(hi, xlim=xrange, ylim=yrange, xlab="", ylab="", axes=F, col=7)
	par(new=T)
	curve(lo, xlim=xrange, ylim=yrange, xlab="", ylab="", axes=F, col=7)

	legend("bottomleft", c(
		expression("E[t|x]"), 
		expression("E[t|x]" %+-% 2 * sqrt("Var[t|x]"))
		), col=c(2, 7), lty=1, bg="gray")
	title(bquote(paste("p(t|x)", ~ sigma, "=", .(sigma))))

	# p(t|x)のtに関する積分が1であることを、いくつかのxで確認
	for (xx in seq(0, 1, 0.1)) {
		cat(paste0(
			"integral(p(t|x=", xx ,")dt)=", 
			sum(sapply(seq(-10, 10, 0.01), function(yy) { p(xx, yy) * 0.01 })),
			"\n"
			))
	}
}

doplot(0.02)
doplot(0.06)
doplot(0.1)