第6章 Gauss過程と関数データ解析

6.1 回帰

例83

## (m, k) の定義
m <- function(x) {
  return(0)
}

k <- function(x, y) {
  return(exp(-(x-y)^2/2))
}

## 関数 gp.sample の定義
gp.sample <- function(x, m, k) {
  n <- length(x)
  m.x <- m(x)
  k.xx <- matrix(0, n, n)
  for (i in 1:n)
    for (j in 1:n)
      k.xx[i, j] <- k(x[i], x[j])
  R <- t(chol(k.xx))
  u <- rnorm(n)
  return(as.vector(R %*% u + m.x))
}
## 乱数を発生して、共分散行列を生成して k.xx と比較
x <- seq(-2, 2, 1)
n <- length(x)
r <- 100
z <- matrix(0, r, n)
for (i in 1:r)
  z[i, ] <- gp.sample(x, m, k)
k.xx <- matrix(0, n, n)
for (i in 1:n)
  for (j in 1:n)
    k.xx[i, j] <- k(x[i], x[j])

cov(z)
##             [,1]        [,2]       [,3]        [,4]        [,5]
## [1,]  1.08880903  0.67243959 0.09594867 -0.17497804 -0.21306343
## [2,]  0.67243959  1.02513038 0.61171018  0.05903026 -0.07896771
## [3,]  0.09594867  0.61171018 0.91013896  0.51992286  0.19617466
## [4,] -0.17497804  0.05903026 0.51992286  0.97136699  0.75986754
## [5,] -0.21306343 -0.07896771 0.19617466  0.75986754  1.23720531
k.xx
##              [,1]      [,2]      [,3]      [,4]         [,5]
## [1,] 1.0000000000 0.6065307 0.1353353 0.0111090 0.0003354626
## [2,] 0.6065306597 1.0000000 0.6065307 0.1353353 0.0111089965
## [3,] 0.1353352832 0.6065307 1.0000000 0.6065307 0.1353352832
## [4,] 0.0111089965 0.1353353 0.6065307 1.0000000 0.6065306597
## [5,] 0.0003354626 0.0111090 0.1353353 0.6065307 1.0000000000

例84

## (m, k) の定義
m <- function(x) {
  return(x[, 1] - x[, 2])
}

k <- function(x, y) {
  return(exp(-sum((x-y)^2)/2))
}

## 関数 gp.sample の定義
gp.sample <- function(x, m, k) {
  n <- nrow(x)
  m.x <- m(x)
  k.xx <- matrix(0, n, n)
  for (i in 1:n)
    for (j in 1:n)
      k.xx[i, j] <- k(x[i, ], x[j, ])
  R <- t(chol(k.xx))
  u <- rnorm(n)
  return(R %*% u + m.x)
}

## 乱数を発生して、共分散行列を生成して k.xx と比較
n <- 5
x <- matrix(rnorm(n*2), n, n)
r <- 100
z <- matrix(0, r, n)
for (i in 1:r)
  z[i, ] <- gp.sample(x, m, k)
k.xx <- matrix(0, n, n)
for (i in 1:n)
  for (j in 1:n)
    k.xx[i, j] <- k(x[i], x[j])

cov(z)
##              [,1]        [,2]       [,3]        [,4]         [,5]
## [1,]  1.101240270 -0.01143859 0.05126216 -0.08733825 -0.007069143
## [2,] -0.011438593  1.00318703 0.70327875  0.64024212  0.017702691
## [3,]  0.051262157  0.70327875 1.25712058  0.15333950  0.016796672
## [4,] -0.087338248  0.64024212 0.15333950  1.02948949  0.105466695
## [5,] -0.007069143  0.01770269 0.01679667  0.10546669  1.000941026
k.xx
##              [,1]        [,2]         [,3]       [,4]      [,5]
## [1,] 1.0000000000 0.002558694 0.0001935038 0.01523558 0.1160772
## [2,] 0.0025586944 1.000000000 0.7933973164 0.85385171 0.3861066
## [3,] 0.0001935038 0.793397316 1.0000000000 0.46214732 0.1198318
## [4,] 0.0152355814 0.853851711 0.4621473195 1.00000000 0.7159672
## [5,] 0.1160772396 0.386106646 0.1198317618 0.71596721 1.0000000
gp.1 <- function(x.pred) {
  h <- array(dim = n)
  for (i in 1:n)
    h[i] <- k(x.pred, x[i])
  R <- solve(K + sigma.2 * diag(n))               # O(n^3) の計算
  mm <- mu(x.pred) + t(h) %*% R %*% (y - mu(x))
  ss <- k(x.pred, x.pred) - t(h) %*% R %*% h
  return(list(mm = mm, ss = ss))
}

gp.2 <- function(x.pred) {
  h <- array(dim = n)
  for (i in 1:n)
    h[i] <- k(x.pred, x[i])
  L <- chol(K + sigma.2 * diag(n))                # O(n^3/3) の計算
  alpha <- solve(L, solve(t(L), y - mu(x)))       # O(n^2) の計算
  mm <- mu(x.pred) + sum(t(h) * alpha)
  gamma <- solve(t(L), h)                         # O(n^2) の計算
  ss <- k(x.pred, x.pred) - sum(gamma^2)
  return(list(mm = mm, ss = ss))
}

例85

sigma.2 <- 0.2

k <- function(x, y) {                    # 共分散関数
  return(exp(-(x-y) ^ 2 / 2 / sigma.2))
}

mu <- function(x) {                      # 平均関数
  return(x)
}

# データ生成
n <- 1000
x <- runif(n)*6 - 3
y <- sin(x/2) + rnorm(n)

K <- array(dim = c(n, n))
for (i in 1:n)
  for (j in 1:n)
    K[i, j] <- k(x[i], x[j])

## 実行時間を測定
library(tictoc)

tic()
gp.1(0)
## $mm
##            [,1]
## [1,] -0.1127373
## 
## $ss
##             [,1]
## [1,] 0.003663641
toc()
## 0.73 sec elapsed
tic()
gp.2(0)
## $mm
## [1] -0.1127373
## 
## $ss
## [1] 0.003663641
toc()
## 0.72 sec elapsed
## 平均の前後で 3 sigma の幅も記載
u.seq <- seq(-3, 3, 0.1)
v.seq <- NULL
w.seq <- NULL
for (u in u.seq) {
  res <- gp.1(u)
  v.seq <- c(v.seq, res$mm)
  w.seq <- c(w.seq, sqrt(res$ss))
}
plot(u.seq, v.seq, xlim = c(-3, 3), ylim = c(-3, 3), type = "l")
lines(u.seq, v.seq + 3*w.seq, col = "blue")
lines(u.seq, v.seq - 3*w.seq, col = "blue")
points(x, y)

## サンプルを変えて 5 回
plot(0, xlim = c(-3, 3), ylim = c(-3, 3), type = "n")
n <- 100
for (h in 1:5) {
  x <- runif(n)*6 - 3
  y <- sin(pi*x/2) + rnorm(n)
  sigma2 <- 0.2
  K <- array(dim = c(n, n))
  for (i in 1:n)
    for (j in 1:n)
      K[i, j] <- k(x[i], x[j])
  u.seq <- seq(-3, 3, 0.1)
  v.seq <- NULL
  for (u in u.seq) {
    res <- gp.1(u)
    v.seq <- c(v.seq, res$mm)
  }
  lines(u.seq, v.seq, col = h+1)
}

6.2 分類

例86

## Iris データ
df <- iris
x <- df[1:100, 1:4]
y <- c(rep(1, 50), rep(-1, 50))
n <- length(y)

## 4 個の共変量でカーネルを計算
k <- function(x, y) {
  return(exp(sum(-(x-y)^2)/2))
}

K <- matrix(0, n, n)
for (i in 1:n)
  for (j in 1:n)
    K[i, j] <- k(x[i, ], x[j, ])
eps <- 0.00001
f <- rep(0, n)
g <- rep(0.1, n)
while (sum((f-g)^2) > eps) {
  g <- f     # 比較のため、更新前の値を保存する
  v <- exp(-y * f)
  u <- y * v / (1+v)
  w <- as.vector(v / (1+v) ^ 2)
  W <- diag(w)
  W.p <- diag(w ^ 0.5)
  W.m <- diag(w ^ (-0.5))
  L <- chol(diag(n) + W.p %*% K %*% W.p)
  L <- t(L)  # R 言語の chol 関数は転置した行列を出力する
  gamma <- W %*% f + u
  beta <- solve(L, W.p %*% K %*% gamma)
  alpha <- solve(t(L) %*% W.m, beta)
  f <- K %*% (gamma - alpha)
}
as.vector(f)
##   [1]  2.901760  2.666188  2.736000  2.596215  2.888826  2.422990  2.712837
##   [8]  2.896583  2.263840  2.722794  2.675787  2.804277  2.629917  2.129059
##  [15]  1.994737  1.725577  2.502404  2.894768  2.211715  2.757889  2.580703
##  [22]  2.788434  2.450147  2.598253  2.493633  2.589272  2.799560  2.854338
##  [29]  2.858034  2.682198  2.663180  2.652952  2.409809  2.157029  2.738196
##  [36]  2.777507  2.605493  2.848624  2.342636  2.882683  2.887406  1.561917
##  [43]  2.454169  2.649399  2.407117  2.633906  2.727124  2.673216  2.749571
##  [50]  2.884288 -1.870442 -2.537382 -1.932737 -2.531886 -2.579367 -2.785015
##  [57] -2.381783 -1.467521 -2.486519 -2.356974 -1.600772 -2.811324 -2.381371
##  [64] -2.734406 -2.330770 -2.341664 -2.614817 -2.748088 -2.372972 -2.628800
##  [71] -2.251908 -2.789266 -2.345693 -2.704238 -2.712489 -2.502956 -2.162480
##  [78] -2.015710 -2.840261 -2.194749 -2.474090 -2.342643 -2.755051 -2.189084
##  [85] -2.415174 -2.382286 -2.275106 -2.496755 -2.695880 -2.664357 -2.648632
##  [92] -2.771846 -2.807510 -1.546826 -2.814728 -2.756341 -2.834644 -2.849801
##  [99] -1.182284 -2.845863
pred <- function(z) {
  kk <- array(0, dim = n)
  for (i in 1:n)
    kk[i] <- k(z, x[i, ])

  # 平均
  mu <- sum(kk * as.vector(u))

  # 分散
  alpha <- solve(L, W.p %*% kk)
  sigma2 <- k(z, z) - sum(alpha^2)

  # 予測値
  m <- 1000
  b <- rnorm(m, mu, sigma2)
  pi <- sum((1 + exp(-b))^(-1)) / m

  return(pi)
}

例87

z <- array(0, dim = 4)
for (j in 1:4)
  z[j] <- mean(x[1:50, j])
pred(z)
## [1] 0.9464529
for (j in 1:4)
  z[j] <- mean(x[51:100, j])
pred(z)
## [1] 0.05268072

6.3 補助変数法

例88

sigma.2 <- 0.05  # 本来は推定すべき

k <- function(x, y) {                    # 共分散関数
  return(exp(-(x-y) ^ 2 / 2 / sigma.2))
}

mu <- function(x) {
  return(x)                              # 平均関数
}

# データ生成
n <- 200
x <- runif(n)*6 - 3
y <- sin(x/2) + rnorm(n)

eps <- 10 ^ (-6)

m <- 100
index <- sample(1:n, m, replace = FALSE)
z <- x[index]
m.x <- 0
m.z <- 0
K.zz <- array(dim = c(m, m))
for (i in 1:m)
  for (j in 1:m)
    K.zz[i, j] <- k(z[i], z[j])
K.xz <- array(dim = c(n, m))
for (i in 1:n)
  for (j in 1:m)
    K.xz[i, j] <- k(x[i], z[j])
K.zz.inv <- solve(K.zz + diag(rep(10 ^ eps, m)))
lambda <- array(dim = n)
for (i in 1:n)
  lambda[i] <- k(x[i], x[i]) - K.xz[i, 1:m] %*% K.zz.inv %*% K.xz[i, 1:m]
Lambda.0.inv <- diag(1 / (lambda + sigma.2))
Q <- K.zz + t(K.xz) %*% Lambda.0.inv %*% K.xz  # Q の計算は O(n^3) を要求しない
Q.inv <- solve(Q + diag(rep(eps, m)))
muu <- Q.inv %*% t(K.xz) %*% Lambda.0.inv %*% (y - m.x)
dif <- K.zz.inv - Q.inv
K <- array(dim = c(n, n))
for (i in 1:n)
  for (j in 1:n)
    K[i, j] <- k(x[i], x[j])
R <- solve(K + sigma.2 * diag(n))              # O(n^3) の計算が必要

gp.ind <- function(x.pred) {
  h <- array(dim = m)
  for (i in 1:m)
    h[i] <- k(x.pred, z[i])
  mm <- mu(x.pred) + h %*% muu
  ss <- k(x.pred, x.pred) - h %*% dif %*% h
  return(list(mm = mm, ss = ss))
}                                       # 補助変数法を用いる

gp.1 <- function(x.pred) {
  h <- array(dim = n)
  for (i in 1:n)
    h[i] <- k(x.pred, x[i])
  mm <- mu(x.pred) + t(h) %*% R %*% (y - mu(x))
  ss <- k(x.pred, x.pred) - t(h) %*% R %*% h
  return(list(mm = mm, ss = ss))
}                                       # 補助変数法を用いない

x.seq <- seq(-2, 2, 0.1)
mmv <- NULL
ssv <- NULL
for (u in x.seq) {
  mmv <- c(mmv, gp.ind(u)$mm)
  ssv <- c(ssv, gp.ind(u)$ss)
}
plot(0, xlim = c(-2, 2), ylim = c(min(mmv), max(mmv)), type = "n")
lines(x.seq, mmv, col = "red")
lines(x.seq, mmv + 3*sqrt(ssv), lty = 3, col = "red")
lines(x.seq, mmv - 3*sqrt(ssv), lty = 3, col = "red")

x.seq <- seq(-2, 2, 0.1)
mmv <- NULL
ssv <- NULL
for (u in x.seq) {
  mmv <- c(mmv, gp.1(u)$mm)
  ssv <- c(ssv, gp.1(u)$ss)
}

lines(x.seq, mmv, col = "blue")
lines(x.seq, mmv + 3*sqrt(ssv), lty = 3, col = "blue")
lines(x.seq, mmv - 3*sqrt(ssv), lty = 3, col = "blue")
points(x, y)

6.4 Karhunen-Loeve展開

例89

lambda <- function(j) {
  return(4 / ((2*j-1)*pi) ^ 2)            # 固有値
}

ee <- function(j, x) {
  return(sqrt(2) * sin((2*j-1)*pi/2*x))   # 固有関数の定義
}

n <- 10
m <- 7

f <- function(z, x) {                     # ガウス過程の定義
  n <- length(z)
  S <- 0
  for (i in 1:n)
    S <- S + z[i] * ee(i, x) * sqrt(lambda(i))
  return(S)
}

plot(0, xlim = c(-3, 3), ylim = c(-2, 2), type = "n",
     xlab = "x", ylab = "f(omega, x)")
for (j in 1:m) {
  z <- rnorm(n)
  x.seq <- seq(-3, 3, 0.001)
  y.seq <- NULL
  for (x in x.seq)
    y.seq <- c(y.seq, f(z, x))
  lines(x.seq, y.seq, col = j)
}
title("Brown Motion")

matern <- function(nu, l, r) {
  p <- nu - 1/2
  S <- 0
  for (i in 0:p)
    S <- S + gamma(p+i+1) / gamma(i+1) / gamma(p-i+1) * (sqrt(8*nu)*r/l) ^ (p-i)
  S <- S * gamma(p+2) / gamma(2*p+1) * exp(-sqrt(2*nu)*r/l)
  return(S)
}

例90

m <- 10
l <- 0.1
for (i in 1:1)
  curve(matern(i-1/2, l, x), 0, 0.5, ylim = c(0, 10), col = i+1)
for (i in 2:m)
  curve(matern(i-1/2, l, x), 0, 0.5, ylim = c(0, 10),
        ann = FALSE, add = TRUE, col = i+1)
legend("topright", legend = paste("nu = ", (1:m)+0.5), lwd = 2, col = 1:m)
title("Matern Kernel (l = 1)")

例91

rand.100 <- function(Sigma) {
  L <- t(chol(Sigma))       # 共分散行列を Cholesky 分解  
  u <- rnorm(100)
  y <- as.vector(L %*% u)   # 平均 0 共分散行列の乱数を 1 組生成
  return(y)
}

x <- seq(0, 1, length = 100)
z <- abs(outer(x, x, "-"))  # compute distance matrix, d_{i,j} = |x_i - x_j|
l <- 0.1

Sigma_OU <- exp(-z / l)     # OU: matern(1/2, l, z) では遅い
y <- rand.100(Sigma_OU)

plot(x, y, type = "l", ylim = c(-3, 3))
for (i in 1:5) {
  y <- rand.100(Sigma_OU)
  lines(x, y, col = i+1)
}
title("OU process (nu = 1/2, l = 0.1)")

Sigma_M <- matern(3/2, l, z)   # Matern 
y <- rand.100(Sigma_M)
plot(x, y, type = "l", ylim = c(-3, 3))
for (i in 1:5) {
  y <- rand.100(Sigma_M)
  lines(x, y, col = i+1)
}
title("Matern process (nu = 3/2, l = 0.1)")

6.5 関数データ解析

例92

library(fda)
##  要求されたパッケージ splines をロード中です
##  要求されたパッケージ Matrix をロード中です
##  要求されたパッケージ fds をロード中です
##  要求されたパッケージ rainbow をロード中です
##  要求されたパッケージ MASS をロード中です
##  要求されたパッケージ pcaPP をロード中です
##  要求されたパッケージ RCurl をロード中です
## 
##  次のパッケージを付け加えます: 'fda'
##  以下のオブジェクトは 'package:graphics' からマスクされています: 
## 
##      matplot
g <- function(j, x) {        # 基底を p 個用意する
  if (j == 1)
    return(1 / sqrt(2*pi)) 
  if (j %% 2 == 0) {
    return(cos((j %/% 2) * x) / sqrt(pi))
  } else {
    return(sin((j %/% 2) * x) / sqrt(pi))
  }
}

beta <- function(x, y) {     # 関数の p 個の基底の前の係数を計算する
  X <- matrix(0, N, p)
  for (i in 1:N)
    for (j in 1:p)
      X[i, j] <- g(j, x[i])
  beta <- solve(t(X) %*% X + 0.0001 * diag(p)) %*% t(X) %*% y
  return(drop(beta))
}

N <- 365
n <- 35
m <- 5
p <- 100
df <- daily
C <- matrix(0, n, p)
for (i in 1:n) {
  x <- (1:N)*(2*pi/N) - pi
  y <- as.vector(df[[2]][, i])
  C[i, ] <- beta(x, y)
}
res <- prcomp(C)
B <- res$rotation
xx <- res$x
z <- function(i, m, x) {  # p 個の基底のうち、m 主成分で近似したもとの関数
  S <- 0
  for (j in 1:p)
    for (k in 1:m)
      for (r in 1:p)
        S <- S + C[i, j] * B[j, k] * B[r, k] * g(r, x)
  return(S)
}

x.seq <- seq(-pi, pi, 2 * pi / 100)
plot(0, xlim = c(-pi, pi), ylim = c(-15, 25), type = "n",
     xlab = "Days", ylab = "Temp (C)", main = "Reconstruction for each m")
lines(x, df[[2]][, 14], lwd = 2)
for (m in 2:6)
  lines(x.seq, z(14, m, x.seq), col = m, lty = 1)
legend("bottom", legend = c("Original", paste("m = ", 2:6)),
       lwd = c(2, rep(1, 5)), col = 1:6, ncol = 2)

lambda <- res$sdev ^ 2
ratio <- lambda / sum(lambda)
plot(1:5, ratio[1:5], xlab = "PC1 through PC5", ylab = "Ratio", type = "l",
     main = "Ratio")

h <- function(coef, x) {  # 係数を用いて関数を定義する
  S <- 0
  for (j in 1:p)
    S <- S + coef[j] * g(j, x)
  return(S)
}

plot(0, xlim = c(-pi, pi), ylim = c(-1, 1), type = "n")
for (j in 1:3)
  lines(x.seq, h(B[, j], x.seq), col = j)

index <- c(10, 12, 13, 14, 17, 24, 26, 27)
others <- setdiff(1:35, index)
first <- substring(df[[1]][index], 1, 1)
plot(0, xlim = c(-25, 35), ylim = c(-15, 10), type = "n",
     xlab = "PC1", ylab = "PC2", main = "Canadian Weather")
points(xx[others, 1], xx[others, 2], pch = 4)
points(xx[index, 1], xx[index, 2], pch = first, cex = 1, col = rainbow(8))
legend("bottom", legend = df[[1]][index], pch = first,
       col = rainbow(8), ncol = 2)