「大人の教養・知識・気付き」を伸ばすブログ

一流の大人（ビジネスマン、政治家、リーダー…）として知っておきたい、教養・社会動向を意外なところから取り上げ学ぶことで“気付く力”を伸ばすブログです。データ分析・語学に力点を置いています。　→現在、コンサルタントの雛になるべく、少しずつ勉強中です(※2024年1月21日改訂)。

MENU

統計的機械学習の数理100問（16/20）

機械学習 統計学 Python R

　いい加減時代の潮流に乗ろうということで機械学習を学びたいと思う。またRはともかくとしてPythonは未経験であるため、丁度良い書籍として

統計的機械学習の数理100問 with R (with R) (機械学習の数理100問シリーズ)

統計的機械学習の数理100問 with R (with R) (機械学習の数理100問シリーズ)

Amazon

統計的機械学習の数理100問 with Python (機械学習の数理100問シリーズ)

統計的機械学習の数理100問 with Python (機械学習の数理100問シリーズ)

Amazon

を用いることにする。

前回

power-of-awareness.com

6.　非線形回帰

  • 平滑化スプライン、局所回帰、一般化加法モデルを扱う。

6.5　平滑化スプライン

　観測した標本(x_1,y_1),\cdots,(x_i,y_i),\cdots,(x_n,y_n),i=1,2,\cdots,n,x_1\lt x_2,\cdots\lt x_nに対して


\begin{aligned} L(f)=\displaystyle{\sum_{i=1}^n(y_i-f(x_i))^2+\lambda\int_{-\infty}^{\infty}\{f^{\prime\prime}(x)\}^2dx} \end{aligned}


を最小にするようなf:\mathbb{R}\rightarrow\mathbb{R}を求めたい。ここで\lambda\geq0は所与とする。この第2項は曲率が大きい、すなわちfが複雑であることに対するペナルティと解釈できる。
　まずデータx_1,\cdots,x_nを区切り点とする自然なスプラインによってそのような最適なfが実現されることを示す。



平滑化スプラインの存在性　x_1,\cdots,x_nを区切り点とする自然なスプラインの中に


\begin{aligned} L(f)=\displaystyle{\sum_{i=1}^n(y_i-f(x_i))^2+\lambda\int_{-\infty}^{\infty}\{f^{\prime\prime}(x)\}^2dx} \end{aligned}


を最小にするようなf:\mathbb{R}\rightarrow\mathbb{R}が存在する。

(\because　命題中のL(f), すなわち


\begin{aligned} L(f)=\displaystyle{\sum_{i=1}^n(y_i-f(x_i))^2+\lambda\int_{-\infty}^{\infty}\{f^{\prime\prime}(x)\}^2dx} \end{aligned}


を最小にするような任意の関数をf(x)とし、g(x)をx_1,\cdots,x_nを区切り点とする自然なスプライン関数だとする。さらに


\begin{aligned} r(x)=f(x)-g(x) \end{aligned}


とおく。まずg(x)の次元がnであるからh_1(x),\cdots,h_n(x)を基底として


\begin{aligned} g(x)=\displaystyle{\sum_{i=1}^{n}\gamma_i h_i(x)} \end{aligned}


として、


\begin{aligned} g(x_1)=f(x_1),\cdots,g(x_n)=f(x_n) \end{aligned}


を満たすように\gamma_1,\cdots,\gamma_nを決めることができる。すなわち


\begin{aligned} \begin{bmatrix}h_1(x_1)&\cdots&h_n(x_1)\\ \vdots&\ddots&\vdots\\ h_1(x_n)&\cdots&h_n(x_n) \end{bmatrix}\begin{bmatrix} \gamma_1\\ \vdots\\ \gamma_n \end{bmatrix}=\begin{bmatrix} f(x_1)\\ \vdots\\ f(x_n) \end{bmatrix} \end{aligned}


を解くことで\gamma_1,\cdots,\gamma_nが得られる。
　このときr(x_1)=\cdots=r(x_n)=0となり、g(x)がx\leq x_1,x_n\leq xで直線で、その他の区間で3次多項式であって、g^{\prime\prime\prime}(x)が各区間[x_i,x_{i+1}]で一定(=\gamma_iとおく。)であることから部分積分により


\begin{aligned} \displaystyle{\int_{x_1}^{x_n}g^{\prime\prime}(x)r^{\prime\prime}(x)dx}=\left[g^{\prime\prime}(x)r^{\prime}(x)\right]_{x_1}^{x_n}-\displaystyle{\int_{x_1}^{x_n}g^{\prime\prime\prime}(x)r^{\prime}(x)dx}=-\displaystyle{\sum_{i=1}^{n-1}\gamma_i\left[r(x) \right]_{x_i}^{x_{i+1}} }=0 \end{aligned}


が成立する。したがって


\begin{aligned} \displaystyle{\int_{-\infty}^{\infty}\{f^{\prime\prime}(x)\}^2 dx}&\geq\displaystyle{\int_{x_1}^{x_n}\{g^{\prime\prime}(x)+r^{\prime\prime}(x)\}^2 dx}\\ &\geq\displaystyle{\int_{x_1}^{x_n}\{g^{\prime\prime}(x)\}^2dx}+\displaystyle{\int_{x_1}^{x_n}\{r^{\prime\prime}(x)\}^2dx}+2\displaystyle{\int_{x_1}^{x_n}g^{\prime\prime}(x)r^{\prime\prime}(x)dx}\\ &\geq\displaystyle{\int_{x_1}^{x_n}\{g^{\prime\prime}(x)\}^2 dx} \end{aligned}


が得られる。すなわちL(f)を最小にする任意のfのそれぞれに対して


\begin{aligned} L(f)&=\displaystyle{\sum_{i=1}^{n}(y_i-f(x_i))^2}+\lambda\displaystyle{\int_{-\infty}^{\infty}\{f^{\prime\prime}(x)\}^2dx}\\ &\geq\displaystyle{\sum_{i=1}^{n}(y_i-g(x_i))^2}+\lambda\displaystyle{\int_{-\infty}^{\infty}\{g^{\prime\prime}(x)\}^2dx}=L(g) \end{aligned}


となるような自然なスプラインgが存在する。　\blacksquare)


　次にそのような自然なスプラインf(x)=\displaystyle{\sum_{i=1}^{n}\gamma_i h_i(x) }の係数\gamma_1,\cdots,\gamma_nを求める。


\begin{aligned} g_{i,j}=\displaystyle{\int_{-\infty}^{\infty}h^{\prime\prime}_i(x)h^{\prime\prime}_j(x)dx} \end{aligned}


を(i,j)成分に持つ行列をGと書くと、L(g)の第2項は


\begin{aligned} \lambda\displaystyle{\int_{-\infty}^{\infty}\{g^{\prime\prime}(x)\}^2dx }&=\displaystyle{\int_{-\infty}^{\infty}\sum_{i=1}^{n}\gamma_i h_i^{\prime\prime}(x)\sum_{j=1}^{n}\gamma_jh_j^{\prime\prime}(x)dx }\\ &=\lambda\displaystyle{\sum_{i=1}^{n}\sum_{j=1}^{n}\gamma_i\gamma_j\int_{-\infty}^{\infty}h_i^{\prime\prime}(x)h_j^{\prime\prime}(x)dx }\\ &=\lambda{}^{t}\boldsymbol{\gamma}G\boldsymbol{\gamma} \end{aligned}


となる。ここで\boldsymbol{\gamma}={}^{t}(\gamma_1,\cdots,\gamma_n)である。したがってL(g)を\boldsymbol{\gamma}で微分すると、\mathrm{Lidge}回帰の場合と同様に、


\begin{aligned} \ -X{}^{t}(\boldsymbol{y}-X\boldsymbol{\gamma})+\lambda G\boldsymbol{\gamma}=\boldsymbol{0} \end{aligned}

となり、解は


\begin{aligned} \hat{\boldsymbol{\gamma}}=({}^{t}XX+\lambda G)^{-1}{}^{t}X\boldsymbol{y} \end{aligned}


で与えられる。

6.6　局所回帰

　\mathrm{Nadaraya}-\mathrm{Watson}推定量と局所線形回帰を導入する。
　\mathcal{X}を集合とし、以下の条件を満たすk:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}を(狭義の)カーネルと呼ぶ。

  1. 任意のn\geq1とx_1,\cdots,x_nに対してK_{i,j}=k(x_i,x_j)を元とするK\in\mathcal{X}^{n\times n}が非負定値である。
  2. 任意のx,y\in\mathcal{X}についてk(x,y)=k(y,x)

　たとえば\mathcal{X}をベクトル空間とすれば、その内積はカーネルとなる。カーネルは集合\mathcal{X}の元の間の類似度を表す(類似していればk(x,y)は大きくなる。)。

　

\begin{aligned} K_{\lambda}(x,y)&=D\left(\displaystyle{\frac{|x-y|}{\lambda}}\right),\\ D(t)&=\begin{cases} \displaystyle{\frac{3}{4}(1-t^2)},|t|\leq1,\\ 0,\mathrm{otherwise} \end{cases} \end{aligned}

によって定義されるk:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}は正定値性を満足しないが、\mathrm{Epanechnikov}カーネルと呼ぶ。

6.6.1　Nadaraya-Watson推定量

　\mathrm{Nadaraya}-\mathrm{Watson}推定量は\mathcal{X}をある集合、k:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}をカーネル関数として、観測データ(x_1,y_1),\cdots,(x_n,y_n)\in\mathcal{X}\times\mathcal{X}から、


\begin{aligned} \hat{f}(x)=\displaystyle{\frac{\displaystyle{\sum_{i=1}^{n}}K(x,x_i)y_i }{\displaystyle{\sum_{j=1}^{n}}K(x,x_j) }} \end{aligned}


により構成された関数である。観測データとは別のx_{*}\in\mathcal{X}が与えられた際に、y_1,\cdots,y_nに


\begin{aligned} \displaystyle{\frac{K(x_{*},x_1) }{\displaystyle{\sum_{j=1}^{n}K(x_{*},x_j)} } },\cdots, \displaystyle{\frac{K(x_{*},x_n) }{\displaystyle{\sum_{j=1}^{n}K(x_{*},x_j)} } } \end{aligned}

のウェイト付けした値を\hat{f}(x_{*})の値として返す。k(u,v)はu,v\in\mathcal{X}が類似していれば大きいことを仮定しているため、x_{*}と類似度の高いx_iに対してのy_iの重みが大きくなる。

6.6.2　局所線形回帰

　各点の近傍を直線で近似して観測点をつい準する方法を局所線形回帰という。
　k:\mathbb{R}^p\times\mathbb{R}^p\rightarrow\mathbb{R}をカーネルとして各点x\in\mathbb{R}で


\begin{aligned} \displaystyle{\sum_{i=1}^{n} k(x,x_i)(y_i-{}^{t}(1,x_i)\boldsymbol{\beta})^2} \end{aligned}


を最小にするような\boldsymbol{\beta}\in\mathbb{R}^{p+1}を求める。
　行列を用いると、本式は


\begin{aligned} {}^{t}(\boldsymbol{y}-X\boldsymbol{\beta}(x))\begin{bmatrix} k(x,x_1)&\cdots&0\\ \vdots&\ddots&\vdots\\ 0&\cdots&k(x,x_n) \end{bmatrix}(\boldsymbol{y}-X\boldsymbol{\beta}(x)) \end{aligned}


と書ける。X\in\mathbb{R}^{n\times(p+1)}の最も左の1列の成分はすべて1である。この式を\boldsymbol{\beta}で微分すると


\begin{aligned} W=\mathrm{diag}(k(x,x_1),\cdots,k(x,x_n)) \end{aligned}


とおけば、


\begin{aligned} \displaystyle{\frac{\partial}{\partial\boldsymbol{\beta}}}\left({}^{t}(\boldsymbol{y}-X\boldsymbol{\beta}(x))\begin{bmatrix}k(x,x_1)&\cdots&0\\\vdots&\ddots&\vdots\\0&\cdots&k(x,x_n)\end{bmatrix}(\boldsymbol{y}-X\boldsymbol{\beta}(x))\right)=-2{}^{t}XW(\boldsymbol{y}-X\boldsymbol{\beta}(x)) \end{aligned}


を得、これを0とおいて解くことで、


\begin{aligned} \hat{\boldsymbol{\beta}}(x)=({}^{t}XWX)^{-1}{}^{t}XW\boldsymbol{y} \end{aligned}


である。

6.7　一般化加法モデル

　基底となる関数が有限個ならば、線形回帰と同様の最小二乗法により係数を推定することができる。平滑化スプラインのように基底の個数が多い場合、逆行列を求めることが困難である。また局所回帰のように有限個の基底で表現できない場合もある。そうした場合、バックフィッティングという方法を適用する。

バックフィッティング
　関数f(x)を複数の関数f_1(x),f_2(x),\cdots,f_p(x)の和で表現する場合、まずf_1(x),f_2(x),\cdots,f_p(x)=0としてから、各j=1,2,\cdots,pに対して残差


\begin{aligned} r_j(x)=f(x)-\displaystyle{\sum_{k\neq j}f_k(x) } \end{aligned}


をf_j(x)で説明する操作を繰り返す。

6.8　Rによるシミュレーション

6.8.1　平滑化スプライン
library("ggplot2")

########################
### 平滑化スプライン ###
########################

d <- function(j, x, knots) {
  K <- length(knots)
  return((max((x - knots[j]) ^ 3, 0) - max((x - knots[K]) ^ 3, 0))
         / (knots[K] - knots[j]))
}

h <- function(j, x, knots) {
  K <- length(knots)
  if (j == 1) {
    return(1)
  } else if (j == 2) {
    return(x)
  } else {
    return(d(j - 2, x, knots) - d(K - 1, x, knots))
  }
}

G <- function(x) {  # xの各値が昇順になっていることを仮定している
  n <- length(x)
  g <- matrix(0, nrow = n, ncol = n)
  for (i in 3:n) {
    for (j in i:n) {
      g[i, j] <- (12 * (x[n] - x[n-1]) * (x[n-1] - x[j-2]) * (x[n-1] - x[i-2])
                  / (x[n] - x[i-2]) / (x[n] - x[j-2])
                  + (12 * x[n-1] + 6 * x[j-2] - 18 * x[i-2])
                  * (x[n-1] - x[j-2]) ^ 2 / (x[n] - x[i-2]) / (x[n] - x[j-2]))
      g[j, i] <- g[i, j]
    }
  }
  return(g)
}


### テストデータ
n <- 100
x <- runif(n, -5, 5)
y <- x + sin(x) * 2 + rnorm(n)

idx <- order(x)
x <- x[idx]
y <- y[idx]
X <- matrix(nrow = n, ncol = n)
X[, 1] <- 1

# 行列Xの生成
for (j in 2:n){
  for(i in 1:n){
    X[i, j] <- h(j, x[i], x)
  }
}

GG <- G(vc_x)

# テストシミュレーションを実施
vc_lambda <- c(1, 30, 80)

for (i in 1:3) {
  lambda <- vc_lambda[i]
  gamma <- solve(t(X) %*% X + lambda * GG) %*% t(X) %*% y
  
  g <- function(u) {
    S <- gamma[1]
    for (j in 2:n)
      S <- S + gamma[j] * h(j, u, x)
    return(S)
  }
  
  u_seq <- seq(-8, 8, 0.02)

  if(i==1){
    df_smoothed_spline <- data.frame(x = u_seq,
                                     y = apply(data.frame(u_seq),1,g),
                                     lambda = paste0("λ=",lambda))
  }else{
    df_smoothed_spline <- rbind(df_smoothed_spline,
                                data.frame(x = u_seq,
                                           y = apply(data.frame(u_seq),1,g),
                                           lambda = paste0("λ=",lambda)))
  }  
}

df_smoothed_spline[,"lambda"] <- as.factor(df_smoothed_spline[,"lambda"])


g1 <- ggplot(data = data.frame(a = as.numeric(x),b = as.numeric(y)), mapping = aes(x = a, y = b)) + geom_point()
g1 <- g1 + geom_line(mapping = aes(x = x, y=y, col = lambda, group = lambda),data = df_smoothed_spline)
g1 <- g1 + labs(x = "X", y = "Y")
g1 <- g1 + theme_classic()
g1 <- g1 + theme(plot.title = element_text(hjust = 0.5),
                 axis.title.x = element_text(size = 10),
                 axis.title.y = element_text(size = 10),
                 legend.title = element_text(size = 14),
                 legend.text = element_text(size = 10),legend.position = "none")
g1 <- g1 + geom_hline(yintercept = 0, linetype =3)
plot(g1)
6.8.2　Nadaraya-Watson推定量
library("ggplot2")

# データ生成
n <- 250
x <- 2 * rnorm(n, 0 ,1)
y <- sin(2 * pi * x) + rnorm(n, 0 ,1) / 4

# 関数定義D
D <- function(t){
  return(max(0.75 * (1 - t^2), 0))
}

# 関数定義K
K <- function(x, y, lambda) {
  return(D(abs(x - y) / lambda))
}

# 関数定義f
f <- function(z, lambda) {
  S <- 0
  t <- 0
  for (i in 1:n) {
    S <- S + K(x[i], z, lambda) * y[i]
    t <- t + K(x[i], z, lambda)
  }
  return(S / t)
}


# lambda = 0.05, 0.25の曲線を図示
df_sample <- data.frame(x = x, y = y)

for(i in 1:2){
  if(i==1){
    df_spline <- data.frame(x = seq(-3,3, 0.1),
                            col = "λ=0.05",
                            y = apply(data.frame(seq(-3,3, 0.1)),1,function(nm){return(f(nm,0.05))}))
  }else{
    df_spline <- rbind(df_spline,
                       data.frame(x = seq(-3,3, 0.1),
                                  col = "λ=0.25",
                                  y = apply(data.frame(seq(-3,3, 0.1)),1,function(nm){return(f(nm,0.25))})))
  }
}


# 最適なlambdaの値を計算
m <- n / 10
lambda_seq <- seq(0.05, 1, 0.01)
SS_min <- Inf

for (lambda in lambda_seq) {
  SS <- 0
  for (k in 1:10) {
    test <- ((k - 1) * m + 1):(k * m)
    train <- setdiff(1:n, test)
    for (j in test) {
      u <- 0
      v <- 0
      for (i in train) {
        kk <- K(x[i], x[j], lambda)
        u <- u + kk * y[i]
        v <- v + kk
      }
      if (v == 0) {
        d_min <- Inf
        
        for (i in train) {
          d <- abs(x[j] - x[i])
          if (d < d_min) {
            d_min <- d
            index <- i
          }
        }
        z <- y[index]
      } else {
        z <- u / v
      }
      SS <- SS + (y[j] - z) ^ 2
    }
  }
  if (SS < SS_min) {
    SS_min <- SS
    lambda_best <- lambda
  }
}


df_spline <- rbind(df_spline,
                   data.frame(x = seq(-3,3, 0.1),
                              col = "λ=最適値",
                              y = apply(data.frame(seq(-3,3, 0.1)),1,function(nm){return(f(nm,lambda_best))}))
                   )

###
g1 <- ggplot(df_sample,aes(x = x, y = y)) + geom_point()
g1 <- g1 + geom_line(mapping = aes(x = x, y=y, col = col, group = col),data = df_spline)
g1 <- g1 + labs(x = "X", y = "Y", title = "Nadaraya-Watson推定量", color = "")
g1 <- g1 + theme_classic()
g1 <- g1 + theme(plot.title = element_text(hjust = 0.5),
                 axis.title.x = element_text(size = 10),
                 axis.title.y = element_text(size = 10),
                 legend.title = element_text(size = 14),
                 legend.text = element_text(size = 10),legend.position = "bottom")
g1 <- g1 + geom_hline(yintercept = 0, linetype =3)
plot(g1)
6.8.3　局所線形回帰
####################
### 局所線形回帰 ###
####################

library("ggplot2")

local <- function(x, y, z = x) {
  X <- cbind(rep(1, n), x)
  yy <- NULL
  beta.hat <- array(dim = 2)
  for (u in z) {
    w <- array(dim = n)
    for (i in 1:n)
      w[i] <- K(x[i], u, lambda = 1)
    W <- diag(w)
    beta.hat <- solve(t(X) %*% W %*% X) %*% t(X) %*% W %*% y
    yy <- c(yy, beta.hat[1] + beta.hat[2] * u)
  }
  return(yy)
}

# データ生成
n <- 300
x <- runif(n) * 2 * pi - pi
y <- sin(x) + rnorm(n, 0 ,1)


m <- 200
U <- seq(-pi, pi, pi / m)
V <- local(x, y, U)

df_sample <- data.frame(x = x, y = y)
df_smooth <- data.frame(x = U, y = V)

lines(U, V, col = "red", type = "l")


g1 <- ggplot(df_sample,aes(x = x, y = y)) + geom_point()
g1 <- g1 + geom_line(mapping = aes(x = x, y=y),data = df_smooth)
g1 <- g1 + labs(x = "X", y = "Y", title = paste0("局所線形回帰 (p = 1, N = ",n,")"))
g1 <- g1 + theme_classic()
g1 <- g1 + theme(plot.title = element_text(hjust = 0.5),
                 axis.title.x = element_text(size = 10),
                 axis.title.y = element_text(size = 10),
                 legend.title = element_text(size = 14),
                 legend.text = element_text(size = 10),legend.position = "bottom")
g1 <- g1 + geom_hline(yintercept = 0, linetype =3)
plot(g1)
6.8.4　一般化加法モデル
########################
### 一般化加法モデル ###
########################

library("ggplot2")

local <- function(x, y, z = x) {
  X <- cbind(rep(1, n), x)
  yy <- NULL
  beta_hat <- array(dim = 2)
  for (u in z) {
    w <- array(dim = n)
    for (i in 1:n)
      w[i] <- K(x[i], u, lambda = 1)
    W <- diag(w)
    beta_hat <- solve(t(X) %*% W %*% X) %*% t(X) %*% W %*% y
    yy <- c(yy, beta_hat[1] + beta_hat[2] * u)
  }
  return(yy)
}

poly <- function(x, y, z = x) {
  n <- length(x)
  m <- length(z)
  X <- cbind(rep(1, n), x, x^2, x^3)
  yy <- array(dim = n)
  beta_hat <- array(dim = 4)
  beta_hat <- solve(t(X) %*% X) %*% t(X) %*% y
  X <- cbind(rep(1, m), z, z^2, z^3)
  yy <- X %*% beta_hat
  return(yy)
}

# データ生成
n <- 300
x <- runif(n) * 2 * pi - pi
y <- sin(x) + rnorm(n,0,1)

y_1 <- 0
y_2 <- 0

for (k in 1:10) {
  y_1 <- poly(x, y - y_2)
  y_2 <- local(x, y - y_1)  # localは，前節で定義された関数を用いる
}

z <- seq(-2, 2, 0.1)
par(mfrow = c(1, 2))
plot(z, poly(x, y_1, z), type = "l", xlab = "x", ylab = "f(x)", 
     main = "多項式回帰（3次）", col = "green")
plot(z, local(x, y_2, z), type = "l", xlab = "x", ylab = "f(x)", 
     main = "局所線形回帰", col = "green")

6.9　Pythonによるシミュレーション

6.9.1　平滑化スプライン
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import japanize_matplotlib
import scipy
from scipy import stats
from numpy.random import randn

def d(j,x,knots):
    K=len(knots)
    return (np.maximum((x-knots[j])**3,0)-np.maximum((x-knots[K-1])**3,0))/(knots[K-1]-knots[j])

def h(j,x,knots):
    K=len(knots)
    if j==0:
        return 1
    elif j==1:
        return x
    else :
        return (d(j-2,x,knots)-d(K-2,x,knots))

def G(x):
    n = len(x)
    g = np.zeros((n, n))
    for i in range(2, n):
        for j in range(i, n):
            numer1 = ((x[n - 2] - x[j - 2]) ** 2) * (
                12 * x[n - 2] + 6 * x[j - 2] - 18 * x[i - 2]
            )
            numer2 = (
                12
                * (x[n - 2] - x[i - 2])
                * (x[n - 2] - x[j - 2])
                * (x[n - 1] - x[n - 2])
            )
            denom = (x[n - 1] - x[i - 2]) * (x[n - 1] - x[j - 2])
            g[i, j] = (numer1 + numer2) / denom
            g[j, i] = g[i, j]
    return g


### テストデータの生成
n = 50
a = -5;b = 5

x=(b-a)*np.random.rand(n)+ a #(-5,5)の一様分布
y=x+np.sin(x)*2+randn(n)
plt.scatter(x,y,c="black",s=10)

index=np.argsort(x)
x=x[index]
y=y[index]
X=np.zeros((n,n))
X[:,0]=1
for j in range(1,n):
    for i in range(n):
        X[i,j]=h(j,x[i],x)
GG=G(x)
lambda_set=[1,30,80]
col_set=["red","blue","green"]
plt.scatter(x,y,c="black",s=10)
plt.title("平滑化スプライン(n=" + str(n) + ")")
plt.xlabel("x")
plt.ylabel("g(x)")
plt.tick_params(labelleft=False)
for i in range(3):
    lam=lambda_set[i]
    gamma=np.linalg.inv(X.T@X+lam*GG)@X.T@y
    def g(u):
        S=gamma[0]
        for j in range(1,n):
            S=S+gamma[j]*h(j,u,x)
        return S
    u_seq=np.arange(-8,8,0.02)
    v_seq=[]
    for u in u_seq:
        v_seq.append(g(u))
    plt.plot(u_seq,v_seq,c=col_set[i],label="λ={}".format(lambda_set[i]))
plt.legend()
6.9.2　Nadaraya-Watson推定量
#############################
### Nadaraya-Watson推定量 ###
#############################

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import japanize_matplotlib
import scipy
from scipy import stats
from numpy.random import randn

n=250
x=2*randn(n)
y=np.sin(2*np.pi*x)+randn(n)/4
def D(t):
    return np.maximum(0.75*(1-t**2),0)
def K(x,y,lam):
    return D(np.abs(x-y)/lam)
def f(z,lam):
    S=0;T=0
    for i in range(n):
        S=S+K(x[i],z,lam)*y[i]
        T=T+K(x[i],z,lam)
    if T==0:
        return(0)
    else: 
        return S/T

plt.scatter(x,y,c="black",s=10)
plt.xlim(-3,3)
plt.ylim(-2,3)
xx=np.arange(-3,3,0.1)
yy=[]
for zz in xx:
    yy.append(f(zz,0.05))
plt.plot(xx,yy,c="green",label="λ=0.05")
yy=[]
for zz in xx:
    yy.append(f(zz,0.25))
plt.plot(xx,yy,c="blue",label="λ=0.25")


# 最適値を取るλでの出力
m=int(n/10)
lambda_seq=np.arange(0.05,1,0.01)
SS_min=np.inf
for lam in lambda_seq:
    SS=0
    for k in range(10):
        test=list(range(k*m,(k+1)*m))
        train=list(set(range(n))-set(test))
        for j in test:
            u=0;v=0
            for i in train:
                kk=K(x[i],x[j],lam)
                u=u+kk*y[i]
                v=v+kk
            if v==0:
                d_min=np.inf
                for i in train:
                    dd=np.abs(x[j]-x[i])
                    if dd<d_min:
                        d_min=dd
                        index=i
                z=y[index]
            else:
                z=u/v
            SS=SS+(y[j]-z)**2
    if SS<SS_min:
        SS_min=SS
        lambda_best=lam
yy=[]
for zz in xx:
    yy.append(f(zz,lambda_best))
plt.plot(xx,yy,c="red",label="λ=λbest")
plt.title(" Nadaraya-Watson推定量")
plt.legend()
6.9.3　局所線形回帰
####################
### 局所線形回帰 ###
####################

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import japanize_matplotlib
import scipy
from scipy import stats
from numpy.random import randn

def local(x,y,z=None):
    if z is None:
        z=x
    n=len(y)
    x=x.reshape(-1,1)
    X=np.insert(x, 0, 1, axis=1)
    yy=[]
    for u in z:
        w=np.zeros(n)
        for i in range(n):
            w[i]=K(x[i],u,lam=1)
        W=np.diag(w)
        beta_hat=np.linalg.inv(X.T@W@X)@X.T@W@y
        yy.append(beta_hat[0]+beta_hat[1]*u)
    return yy

n=30
x=np.random.rand(n)*2*np.pi-np.pi
y=np.sin(x)+randn(n)
plt.scatter(x,y,s=15)
m=200
U=np.arange(-np.pi,np.pi,np.pi/m)
V=local(x,y,U)
plt.plot(U,V,c="red")
plt.title(" 局所線形回帰(p=1,N=30)")
6.9.4　一般化加法モデル
########################
### 一般化加法モデル ###
########################

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import japanize_matplotlib
import scipy
from scipy import stats
from numpy.random import randn

def poly(x,y,z=None):
    if z is None:
        z=x
    n=len(x)
    m=len(z)
    X=np.zeros((n,4))
    for i in range(n):
        X[i,0]=1; X[i,1]=x[i]; X[i,2]=x[i]**2; X[i,3]=x[i]**3
    beta_hat=np.linalg.inv(X.T@X)@X.T@y
    Z=np.zeros((m,4))
    for j in range(m):
        Z[j,0]=1; Z[j,1]=z[j]; Z[j,2]=z[j]**2; Z[j,3]=z[j]**3
    yy=Z@beta_hat
    return yy

n=30
x=np.random.rand(n)*2*np.pi-np.pi
x=x.reshape(-1,1)
y=np.sin(x)+randn(n)
y_1=0;y_2=0
for k in range(10):
    y_1=poly(x,y-y_2)
    y_2=local(x,y-y_1,z=x)
z=np.arange(-2,2,0.1)

plt.plot(z,poly(x,y_1,z))
plt.title("多項式回帰(3次)")

plt.plot(z,local(x,y_2,z))
plt.title("局所線形回帰")

次回

power-of-awareness.com

プライバシーポリシー お問い合わせ