NIPALS(regression)

#入門機械学習による異常検知　コロナ社　井出剛先生　p.176

#NIPALS method(regression)

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),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))

X=t(matrix(c(data$x1,data$x2),ncol=ncol(data)-2,nrow=nrow(data)))

Y=t(t(as.numeric(data$y)));N=length(Y)

cols=nrow(X);X2=X

P=array(0,dim=c(nrow(X),cols))

for(j in 1:cols){

P[,j]=X2%*%Y/sqrt(sum(X2%*%Y^2))  

d=t(X)%*%P[,j]/sqrt(sum((t(X)%*%P[,j])^2))

X2=X2-X2%*%d%*%t(d)

}

R=t(P)%*%X

Beta=solve(R%*%t(R))%*%R%*%Y

result=data.frame(num=1:nrow(data),y=data$y,predict=t(t(Beta)%*%R))

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)

cor(result$y,result$predict)