例59
Hermite <- function(j) { # R 言語は添字が 1 から
if (j == 0)
return(1)
a <- rep(0, j + 2)
b <- rep(0, j + 2)
a[1] <- 1
for (i in 1:j) {
b[1] <- -a[2]
for (k in 1:(i + 1))
b[k + 1] <- 2 * a[k] - (k + 1) * a[k + 2]
a <- b
}
return(b[1:(j + 1)]) # Hermite 多項式の係数を出力
}
Hermite(2) # 2 次の Hermite 多項式
## [1] -2 0 4
Hermite(3) # 3 次の Hermite 多項式
## [1] 0 -12 0 8
Hermite(4) # 4 次の Hermite 多項式
## [1] 12 0 -48 0 16
H <- function(j, x) {
coef <- Hermite(j)
S <- 0
for (i in 0:j)
S <- S + coef[i + 1] * x ^ i
return(S)
}
cc <- sqrt(5) / 4
a <- 1 / 4
phi <- function(j, x) {
return(exp(-(cc-a) * x ^ 2) * H(j, sqrt(2*cc) * x))
}
curve(phi(0, x), -2, 2, ylim = c(-2, 8), col = 1, ylab = "phi")
for (i in 1:3)
curve(phi(i, x), -2, 2, ylim = c(-2, 8), add = TRUE, ann = FALSE, col = i + 1)
legend("topright", legend = paste("j = ", 0:3), lwd = 1, col = 1:4)
title("Gauss カーネルの固有関数")

例62
# カーネルの定義
sigma <- 1
k <- function(x, y) {
return(exp(-(x - y) ^ 2 / sigma ^ 2))
}
# m サンプルの発生とグラム行列の設定
m <- 300
x <- rnorm(m) - 2 * rnorm(m) ^ 2 + 3 * rnorm(m) ^ 3
K <- matrix(0, m, m)
for (i in 1:m)
for (j in 1:m)
K[i, j] <- k(x[i], x[j])
# 固有値・固有ベクトルの計算
eig <- eigen(K)
lam.m <- eig$values
lam <- lam.m / m
U <- eig$vector
alpha <- array(0, dim = c(m, m))
for (i in 1:m)
alpha[, i] <- U[, i] * sqrt(m) / lam.m[i]
# グラフの表示
F <- function(y, i) {
S <- 0
for(j in 1:m)
S <- S + alpha[j, i] * k(x[j], y)
return(S)
}
i <- 1 # i の値を変えて実行する
G <- function(y) {
return(F(y, i))
}
plot(G, xlim = c(-2, 2))
title("Eigen Values and their Eigen Functions")
