#step wise regression
#Regression and the Moore-Penrose Pseudoinverse p.46
f1=function(x){
return(x^2+x^3+x^4)
}
f2=function(x){
return(x)
}
f3=function(x){
return(x^3)
}
#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),sign=0)
tau=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82)
z=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)
n=length(z)
C=cbind(c(f1(tau)),c(f2(tau)))
C_plus=array(0,dim=c(ncol(C),n))
for(j in 1:ncol(C)){
if(j==1){
C1_plus=t(C[,j])/as.numeric(t(C[,j])%*%C[,j])
C_plus[j,]=C1_plus
}else{
C1_plus=t(C_plus[(j-1),])
k2=t((diag(1,n)-t(t(C[,(j-1)]))%*%C1_plus)%*%C[,j])/(t((diag(1,n)-t(t(C[,(j-1)]))%*%C1_plus)%*%C[,j])^2)
C2_plus=C1_plus%*%(diag(1,n)-C[,j]%*%k2)/k2
C_plus[j,]=c(C2_plus)
}
}
x_hat=C_plus%*%z
predict=C%*%x_hat
plot(c(1:length(z)),z,type="o",col=2,xlim=c(1,length(z)),ylim=c(min(c(z,predict)),max(c(z,predict))),xlab="index",ylab="values")
par(new=T)
plot(c(1:length(z)),predict,type="o",col=3,xlim=c(1,length(z)),ylim=c(min(c(z,predict)),max(c(z,predict))),xlab="index",ylab="values")
More than 3 years have passed since last update.
Register as a new user and use Qiita more conveniently
- You get articles that match your needs
- You can efficiently read back useful information
- You can use dark theme
List of users who liked
00