いい加減時代の潮流に乗ろうということで機械学習を学びたいと思う。またRはともかくとしてPythonは未経験であるため、丁度良い書籍として
を用いることにする。
6. 非線形回帰
- 平滑化スプライン、局所回帰、一般化加法モデルを扱う。
6.5 平滑化スプライン
観測した標本に対して
を最小にするようなを求めたい。ここでは所与とする。この第2項は曲率が大きい、すなわちが複雑であることに対するペナルティと解釈できる。
まずデータを区切り点とする自然なスプラインによってそのような最適なが実現されることを示す。
平滑化スプラインの存在性 を区切り点とする自然なスプラインの中に
を最小にするようなが存在する。
を最小にするような任意の関数をとし、をを区切り点とする自然なスプライン関数だとする。さらに
とおく。まずの次元がであるからを基底として
として、
を満たすようにを決めることができる。すなわち
を解くことでが得られる。
このときとなり、がで直線で、その他の区間で3次多項式であって、が各区間で一定(とおく。)であることから部分積分により
が成立する。したがって
が得られる。すなわちを最小にする任意ののそれぞれに対して
となるような自然なスプラインが存在する。 )
次にそのような自然なスプラインの係数を求める。
を成分に持つ行列をと書くと、の第2項は
となる。ここでである。したがってをで微分すると、回帰の場合と同様に、
となり、解は
で与えられる。
6.6 局所回帰
-推定量と局所線形回帰を導入する。
を集合とし、以下の条件を満たすを(狭義の)カーネルと呼ぶ。
- 任意のとに対してを元とするが非負定値である。
- 任意のについて
たとえばをベクトル空間とすれば、その内積はカーネルとなる。カーネルは集合の元の間の類似度を表す(類似していればは大きくなる。)。
によって定義されるは正定値性を満足しないが、カーネルと呼ぶ。
6.7 一般化加法モデル
基底となる関数が有限個ならば、線形回帰と同様の最小二乗法により係数を推定することができる。平滑化スプラインのように基底の個数が多い場合、逆行列を求めることが困難である。また局所回帰のように有限個の基底で表現できない場合もある。そうした場合、バックフィッティングという方法を適用する。
バックフィッティング
関数を複数の関数の和で表現する場合、まずとしてから、各に対して残差
をで説明する操作を繰り返す。
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("局所線形回帰")