B-spline regression（多変量解析入門　小西貞則著）

B-splineを「多変量解析入門　小西貞則著」をもとに作成してみました。



```{r}

#spline regression(正則化項)

library(dplyr);library(pforeach)

#fixed_data

lambda=-0.5

data=read.csv("~/spline_sample_data.csv")

sample=0

if(sample==1){

#sample_data

n=5;N=30;

data=pforeach(j=1:n,.c=rbind)({
  
Data=array(0,dim=c(N,2))

mean=sample(1:9,1);sd=0.5

Data[,1]=j*rnorm(N,mean=mean,sd=sd)

Data[,2]=rnorm(N,mean=mean,sd=sd)

colnames(Data)=c("y","x")

Data=data.frame(num=j,Data)

})

}

#estimate

mean=c();sd=c()

for(j in 1:n){
  
data_sub=data %>% dplyr::filter(num==j)

mean=c(mean,mean(data_sub$x));sd=c(sd,sd(data_sub$x))
  
}

#B

B=array(0,dim=c(nrow(data),n))


for(i in 1:ncol(B)){  
 
gaussian=function(x){
  
z<-exp(-((x-mean[i])^2)/(2*sd[i]^2))

return(z)
  
}

B[,i]=gaussian(data$x)


}  

B=cbind(array(1,dim=c(nrow(data),1)),B)

w=solve(t(B)%*%B)%*%t(B)%*%data$y

data=data %>% dplyr::mutate(predict=0)

for(j in 1:nrow(data)){
  
  
weighted_gaussian=function(x){
  
z=w[1];

for(i in 1:n){
  
z=z+w[i+1]*exp(-((x-mean[i])^2)/(2*sd[i]^2))  
  
}

return(z)  
  
}
    
data$predict[j]=weighted_gaussian(data$x[j]) 
  
}

#正則化

w=solve(t(B)%*%B+lambda*(sum((data$y-data$predict)^2)/nrow(data))* diag(1,n+1))%*%t(B)%*%data$y

data=data %>% dplyr::select(-predict) %>% dplyr::mutate(predict=0)

for(j in 1:nrow(data)){
  
  
weighted_gaussian=function(x){
  
z=w[1];

for(i in 1:n){
  
z=z+w[i+1]*exp(-((x-mean[i])^2)/(2*sd[i]^2))  
  
}

return(z)  
  
}
    
data$predict[j]=weighted_gaussian(data$x[j]) 
  
}

plot(data$x,data$y,type="p",col=2,xlim=c(min(data$x),max(data$x)),ylim=c(min(data$y,data$predict),max(data$y,data$predict)),xlab="x",ylab="y",main=paste0("spline:lambda=",lambda))

par(new=T)

plot(data$x,data$predict,type="p",col=3,xlim=c(min(data$x),max(data$x)),ylim=c(min(data$y,data$predict),max(data$y,data$predict)),xlab="x",ylab="y",main=paste0("spline:lambda=",lambda))




H=B%*%(solve(t(B)%*%B+lambda*(1/nrow(data))*sum((data$y-data$predict)^2)*diag(1,n+1)))%*%t(B)

effective_degrees_of_freedom=length(diag(H))

AIC=nrow(data)*(log(2*pi)+1)+n*log((1/nrow(data))*sum((data$y-data$predict)^2))+effective_degrees_of_freedom

AIC

```

```{r}

#spline regression(1次元)

library(dplyr);

#sample_data

data(anscombe)

X=anscombe[,colnames(anscombe) %in% c("x1")]

Y=anscombe$y4

num=c(1,2,3,2,1,3,2,2,3,2,2)

clusters=max(num)

data=data.frame(num,X,Y)

#estimate

mean=c()

sd=c()

for(j in 1:clusters){

data_sub=data %>% dplyr::filter(num==j)%>%select(-num,-Y)

mean=c(mean,mean(data_sub$X))

sd=c(sd,sum((data_sub-mean(data_sub$X))^2)/nrow(data_sub))

}

#B

B=array(0,dim=c(nrow(data),clusters))

for(i in 1:nrow(B)){

for(j in 1:clusters){

ave=mean[j]  

B[i,j]=exp(-sum((X[i]-ave)^2)/(2*sd[j]^2))

}  

}

B=cbind(rep(1,nrow(B)),B)

w=solve(t(B)%*%B)%*%t(B)%*%Y


#関数定義

weighted_gaussian=function(x){

z=w[1];

for(j in 1:clusters){

aves=mean[j]  

z=z+w[j+1]*exp(-sum((x-aves)^2)/(2*sd[j]^2))

}

return(z)  

}


#予測値の計算

predict=c()

for(j in 1:nrow(data)){

predict=c(predict,weighted_gaussian(X[j])) 

}

#結果のプロット

plot(c(1:length(Y)),Y,type="o",col=2,xlim=c(1,length(Y)),ylim=c(min(Y,predict),max(Y,predict)),xlab="x",ylab="y",main="spline")

par(new=T)

plot(c(1:length(Y)),predict,type="o",col=3,xlim=c(1,length(Y)),ylim=c(min(Y,predict),max(Y,predict)),xlab="x",ylab="y",main="spline")





```


```{r}

#spline regression(多次元)

library(dplyr);

#sample_data

data(anscombe)

X=anscombe[,colnames(anscombe) %in% c("x1","x2","x3")]

Y=anscombe$y4

num=c(1,2,3,2,1,3,2,2,3,2,2)

clusters=max(num)

data=data.frame(num,X,Y)

#estimate

mean=array(0,dim=c(clusters,ncol(X)))

sd=c()

for(j in 1:clusters){

data_sub=data %>% dplyr::filter(num==j)%>%select(-num,-Y)

mean[j,]=c(apply(data_sub,2,mean))

sd=c(sd,sum((data_sub-apply(data_sub,2,mean))^2)/nrow(data_sub))

}

#B

B=array(0,dim=c(nrow(data),clusters))

for(i in 1:nrow(B)){

for(j in 1:clusters){

ave=mean[j,]  

B[i,j]=exp(-sum((X[i,]-ave)^2)/(2*sd[j]^2))

}  

}

B=cbind(rep(1,nrow(B)),B)

w=solve(t(B)%*%B)%*%t(B)%*%Y


#関数定義

weighted_gaussian=function(x){

z=w[1];

for(j in 1:clusters){

aves=mean[j,]  

z=z+w[j+1]*exp(-sum((x-aves)^2)/(2*sd[j]^2))

}

return(z)  

}


#予測値の計算

predict=c()

for(j in 1:nrow(data)){

predict=c(predict,weighted_gaussian(X[j,])) 

}

#結果のプロット

plot(c(1:length(Y)),Y,type="p",col=2,xlim=c(1,length(Y)),ylim=c(min(Y,predict),max(Y,predict)),xlab="x",ylab="y",main="spline")

par(new=T)

plot(c(1:length(Y)),predict,type="p",col=3,xlim=c(1,length(Y)),ylim=c(min(Y,predict),max(Y,predict)),xlab="x",ylab="y",main="spline")






```