Takezawa, Introduction to Nonparametric Regression


Takezawa, Introduction to Nonparametric Regression　の一部を実装しました。

```{r}

library(dplyr)

#smoothing spline

data=data.frame(num=1:20,y=c(167,167.5,168.4,172,155.3,151.4,163,174,168,160.4,164.7,171,162.6,164.8,163.3,167.6,169.2,168,167.4,172))

S=unit=matrix(0,nrow=nrow(data),ncol=nrow(data))

S[1,]=c(c(1,-2,1),rep(0,nrow(data)-3))
S[2,]=c(c(-2,5,-4,1),rep(0,nrow(data)-4))
S[nrow(data),]=c(rep(0,nrow(data)-3),c(1,-2,1))
S[nrow(data)-1,]=c(rep(0,nrow(data)-4),c(1,-4,5,-2))
for(j in 3:(nrow(data)-2)){
k=j-3
vec=c(rep(0,k),c(1,-4,6,-4,1),rep(0,nrow(data)-k-5))
S[j,]=vec
}
diag(unit)<-1
lambda=3

Y=solve(unit+lambda*S)%*%t(t(data$y))

result=data.frame(num=1:nrow(data),y=data$y,predict=Y)

plot(result$num,result$y,xlim=c(1,nrow(data)),ylim=c(min(result$y,result$predict),max(result$y,result$predict)),xlab="sample",ylab="values",type="o",col=2)
par(new=T)
plot(result$num,result$predict,xlim=c(1,nrow(data)),ylim=c(min(result$y,result$predict),max(result$y,result$predict)),xlab="sample",ylab="values",type="o",col=3)


```

局所線形回帰(Introduction to Nonparametric Regression, Takezawa Kunio, p.49)を実装してみました。Gaussianを重みとした平均２乗誤差を最小化する方法です。各(x,y)値ごとに傾きと切片を求める方法です。

```{r}
library(dplyr)
data=data.frame(num=1:20,y=c(167,167.5,168.4,172,155.3,151.4,163,174,168,160.4,164.7,171,162.6,164.8,163.3,167.6,169.2,168,167.4,172),x=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82))

a=b=array(0,dim=c(1,nrow(data)))
h=1

for(j in 1:nrow(data)){
 val_a1=(sum(data$x*exp(-((data$x-data$x[j])^2)/(2*h^2)))*sum(data$y*exp(-((data$x-data$x[j])^2)/(2*h^2)))-sum(data$x*data$y*exp(-((data$x-data$x[j])^2)/(2*h^2)))*sum(exp(-((data$x-data$x[j])^2)/(2*h^2))) )
 val_a2=sum(data$x*exp(-((data$x-data$x[j])^2)/(2*h^2)))^2-sum(data$x^2*exp(-((data$x-data$x[j])^2)/(2*h^2)))*sum(exp(-((data$x-data$x[j])^2)/(2*h^2)))
 
 a[1,j]=val_a1/val_a2
 
 val_b1=sum(data$x*data$y*exp(-((data$x-data$x[j])^2)/(2*h^2)))*sum(data$x*exp(-((data$x-data$x[j])^2)/(2*h^2)))-sum(data$x^2*exp(-((data$x-data$x[j])^2)/(2*h^2)))*sum(data$y*exp(-((data$x-data$x[j])^2)/(2*h^2)))

val_b2=sum(data$x*exp(-((data$x-data$x[j])^2)/(2*h^2)))^2-sum(exp(-((data$x-data$x[j])^2)/(2*h^2)))*sum(data$x^2*exp(-((data$x-data$x[j])^2)/(2*h^2))) 

b[1,j]=val_b1/val_b2
 
 }
X=a*data$x+b
result=data.frame(num=1:nrow(data),y=data$y,predict=t(X))

plot(result$num,result$y,xlim=c(1,nrow(data)),ylim=c(min(result$y,result$predict),max(result$y,result$predict)),xlab="sample number",ylab="values",type="o",col=2)
par(new=T)
plot(result$num,result$predict,xlim=c(1,nrow(data)),ylim=c(min(result$y,result$predict),max(result$y,result$predict)),xlab="sample number",ylab="values",type="o",col=3)




```

```{r}

#kringing

library(dplyr)

library(MASS)

library(stringr)

#Kriging p.238

#example data

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22)+50)

X=data[,str_count(colnames(data), c("x"))>0 | str_count(colnames(data), c("y"))>0]

covariance=cov(X)

#estimate parameters

d0=t(t(covariance[1,2:ncol(covariance)]))

C=covariance[-1,-1]

Y=apply(X,2,mean);Y=Y[2:length(Y)]

alpha=solve(C)%*%d0

beta0_1=c()

beta0_2=c()

for(i in 1:ncol(C)){
for(j in 1:nrow(C)){
  
  beta0_1=c(beta0_1,solve(C)[i,j]*Y[j]) 
  beta0_2=c(beta0_2,solve(C)[i,j])
  
}  
}

beta0=sum(beta0_1)/sum(beta0_2)

beta0+sum(alpha*(Y-beta0))

E_general=c()

for(i in 1:ncol(C)){
for(j in 1:nrow(C)){
  
E_general=c(E_general,(beta0-Y[i])*solve(C)[i,j]*(beta0-Y[j]))  
  
}
}

alpha0=beta0*(1-sum(alpha))

Y0_hat=alpha0+sum(alpha*Y)

epsilon0_hat=Y_hat-beta0

values1=c();

values2=c();

for(j in 1:ncol(C)){
for(i in 1:nrow(C)){
  
values1=c(values1,solve(C)[i,j]*d0[j]) 
values2=c(values2,solve(C)[i,j])
  
}  
}

W_hat=solve(C)%*%(d0+(1-sum(values1))/sum(values2))

sum(W_hat*Y)

```

```{r}

library(dplyr)

library(MASS)

library(stringr)

#Kriging p.248

#example data

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22)+50)

X=data[,str_count(colnames(data), c("x"))>0 | str_count(colnames(data), c("y"))>0]

Y=X[,1];X=X[,2]

g1=function(x){
  
return(log(abs(x)))  
  
}

g2=function(x){
  
return(log(x^2))  
  
}

p=2

G=cbind(g1(X),g2(X));

G=G[2:nrow(G),]

beta=c(1/2,1/2)

ms=G%*%beta

C=array(0,dim=c(length(X),length(X)))

diag(C)=1;C[C!=1]=1/2

d0=t(t(C[1,2:ncol(C)]))

C=C[-1,-1]

alpha=solve(C)%*%d0

Y=Y[2:length(Y)]

beta0_1=c()

beta0_2=c()

for(i in 1:ncol(C)){
for(j in 1:nrow(C)){
  
  beta0_1=c(beta0_1,solve(C)[i,j]*Y[j]) 
  beta0_2=c(beta0_2,solve(C)[i,j])
  
}  
}

beta0=sum(beta0_1)/sum(beta0_2)

Y0_hat=beta0+sum(alpha*(Y-beta0))

epsilon0_hat=Y0_hat-beta0

epsilon_hat=Y-ms

m0=Y0_hat-sum(alpha*(Y-beta0))

alpha0=m0-sum(alpha*ms)



```
```{r}

#p.180 natural spline and smoothing spline

library(dplyr)

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22)+50)

Y=data$y;X=data$x1

g1=function(x){
  
return(exp(-abs(x)/10))  
  
}

g2=function(x){
  
return(exp(-3*abs(x)/1000))  
  
}

G=cbind(g1(X),g2(X))

#一変数関数の積分
func=function(x){
  h=0.001
  y<-(g1(x+h)-2*g1(x)+g1(x-h))*(g2(x+h)-2*g2(x)+g2(x-h))/(h^4)
  return(y)
}

k12=integrate(func,-Inf,Inf)$value

func=function(x){
  h=0.001
  y<-(g1(x+h)-2*g1(x)+g1(x-h))*(g1(x+h)-2*g1(x)+g1(x-h))/(h^4)
  return(y)
}

k11=integrate(func,-Inf,Inf)$value

func=function(x){
  h=0.001
  y<-(g2(x+h)-2*g2(x)+g2(x-h))*(g2(x+h)-2*g2(x)+g2(x-h))/(h^4)
  return(y)
}

k22=integrate(func,-Inf,Inf)$value

K=array(0,dim=c(2,2))

K[1,1]=k11;K[1,2]=K[2,1]=k12;K[2,2]=k22

lambda=0

beta=solve(t(G)%*%G+lambda*K)%*%t(G)%*%Y

H=G%*%solve(t(G)%*%G+lambda*K)%*%t(G)

predict=G%*%beta



plot(c(1:length(Y)),Y,col=2,type="o",xlim=c(1,length(Y)),ylim=c(min(c(Y,predict)),max(c(Y,predict))),xlab="index",ylab="values")

par(new=T)

plot(c(1:length(Y)),predict,col=3,type="o",xlim=c(1,length(Y)),ylim=c(min(c(Y,predict)),max(c(Y,predict))),xlab="index",ylab="values",main=paste0("cor:",floor(cor(Y,predict)*1000)/1000))

```

```{r}

#local polynomial regression p.140

library(dplyr)

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22)+50)

Y=data$y;X=data$x1

P=2;x_star=0

x=array(0,dim=c(length(X),P))

for(j in 1:P){
  
x[,j]=(X-x_star)^j
    
}

x=cbind(rep(1,length(X)),x)

W=array(0,dim=c(length(X),length(X)))

weight=function(z){
  
h=100;h=h^2  
  
return((1-(z-x_star)^2/h)^3)  
  
}

diag(W)=weight(X)

alpha=solve(t(x)%*%W%*%x)%*%t(x)%*%W%*%Y

predict=x%*%alpha

plot(c(1:length(Y)),Y,col=2,type="o",xlim=c(1,length(Y)),ylim=c(min(c(Y,predict)),max(c(Y,predict))),xlab="index",ylab="values")

par(new=T)

plot(c(1:length(Y)),predict,col=3,type="o",xlim=c(1,length(Y)),ylim=c(min(c(Y,predict)),max(c(Y,predict))),xlab="index",ylab="values",main=paste0("cor:",floor(cor(Y,predict)*1000)/1000))

```

```{r}

#p.176

library(dplyr)

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22)+50)

Y=data$y;X=data$x2

R=array(0,dim=c(length(X),length(X)))

for(j in 1:length(X)){
for(i in 1:length(X)){
  
R[i,j]=abs(X[i]-X[j])^3/12  
  
}  
}

Q=t(cbind(rep(1,length(X)),X))

lambda=0.01

t=solve(rbind(cbind(R+lambda*diag(1,length(X)),t(Q)),cbind(Q,array(0,dim=c(nrow(Q),nrow(Q))))))[1:length(X),1:length(X)]

H=(cbind(R,t(Q))%*%solve(rbind(cbind(R+lambda*diag(1,length(X)),t(Q)),cbind(Q,array(0,dim=c(nrow(Q),nrow(Q)))))))[1:length(X),1:length(X)]

sita=t%*%Y

inv_H=svd(H)$v%*%diag(svd(H)$d)%*%t(svd(H)$u)

L=t(t%*%inv_H)%*%R%*%t%*%inv_H

predict=c(solve(diag(1,length(X))+lambda*L)%*%Y)

plot(c(1:length(Y)),Y,col=2,type="o",xlim=c(1,length(Y)),ylim=c(min(c(Y,predict)),max(c(Y,predict))),xlab="index",ylab="values")

par(new=T)

plot(c(1:length(Y)),predict,col=3,type="o",xlim=c(1,length(Y)),ylim=c(min(c(Y,predict)),max(c(Y,predict))),xlab="index",ylab="values",main=paste0("cor:",floor(cor(Y,predict)*1000)/1000))

params=c(solve(rbind(cbind(R+lambda*diag(1,length(X)),t(Q)),cbind(Q,array(0,dim=c(nrow(Q),nrow(Q))))))%*%c(Y,rep(0,nrow(Q))))

C=params[(length(sita)+1):length(params)]

E_ss=t(Y-R%*%sita-t(Q)%*%C)%*%(Y-R%*%sita-t(Q)%*%C)+lambda*t(sita)%*%R%*%sita

```

```{r}

#p.182 

library(dplyr)

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22)+50)

Y=data$y;X=data$x2

D=array(0,dim=c(length(X)-2,length(X)))

for(j in 1:(length(X)-2)){
  
vec=c(1,-2,1)  
  
D[j,j:(j+2)]=vec  
  
}

lambda=0.2

predict=solve(diag(1,length(X))+lambda*t(D)%*%D)%*%Y

plot(c(1:length(Y)),Y,col=2,type="o",xlim=c(1,length(Y)),ylim=c(min(c(Y,predict)),max(c(Y,predict))),xlab="index",ylab="values")

par(new=T)

plot(c(1:length(Y)),predict,col=3,type="o",xlim=c(1,length(Y)),ylim=c(min(c(Y,predict)),max(c(Y,predict))),xlab="index",ylab="values",main=paste0("cor:",floor(cor(Y,predict)*1000)/1000))

```