#Simulate raw data for sine function
x1<-seq(-3*pi,3*pi,by=0.001)
y<-sin(x1)
data1<-as.data.frame(cbind(y,x1))

plot(data1$x1,data1$y,type = "l",col="orange",
     xlab = "x",ylab="y",lwd=5,cex.axis=1.5,cex.lab=2)

################################################################
#Taylor Series (Maclaurin Expansion) with just x
data1$Taylor<-data1$x1

# Regression fit with just x
fit1<-lm(y~x1,data=data1)
fit1<-step(fit1)
summary(fit1)

# Plot regression and Taylor Series with just x
plot(data1$x1,data1$y,type = "l",col="orange",xlab = "x",
     ylab="y",lwd=5,cex.axis=1.4,cex.lab=1.8)
points(data1$x1,fitted(fit1),col='blue',pch=21,cex=1)
points(data1$x,data1$Taylor,col='saddlebrown',pch=21,cex=1)
legend(6, -0.4, legend=c("Regression", "Taylor"),
       col=c("Blue", "saddlebrown"), lty=c(1,1), lwd=6,cex=0.8)

#plot residuals of fit1
plot(data1$x1,fit1$residuals,type = "l",col="red",xlab = "x",ylab="Residual",lwd=2)
#rm(fit1)

################################################################
#Add more polynomial terms to data i.e. till x^6
data1<-as.data.frame(cbind(y,poly(x1,degree = 6,raw=TRUE)))
colnames(data1)[-1] <- paste("x", colnames(data1)[-1], sep = "")

#Taylor Series (Maclaurin Expansion) with polynomial terms till x^6
data1$Taylor<-data1$x1-data1$x3/factorial(3)+data1$x5/factorial(5)

# Regression fit with polynomial terms till x^6
fit2<-lm(y~.,data=data1)
fit2<-step(fit2)
summary(fit2)

# Plot regression and Taylor Series with polynomial terms till x^6
plot(data1$x1,data1$y,type = "l",col="orange",xlab = "x",
     ylab="y",lwd=5,cex.axis=1.4,cex.lab=1.8)
points(data1$x1,fitted(fit2),col='blue',pch=21,cex=1)
points(data1$x1,data1$Taylor,col='saddlebrown',pch=21,cex=1)
legend(6, -0.4, legend=c("Regression", "Taylor"),
       col=c("Blue", "saddlebrown"), lty=c(1,1), lwd=6,cex=0.8)

#plot residuals of fit2
plot(data1$x1,fit2$residuals,type = "l",col="red",xlab = "x",ylab="Residual",lwd=2)
#rm(fit2)

################################################################
#Add more polynomial terms to data i.e. till x^10
data1<-as.data.frame(cbind(y,poly(x1,degree = 10,raw=TRUE)))
colnames(data1)[-1] <- paste("x", colnames(data1)[-1], sep = "")

#Taylor Series (Maclaurin Expansion) with polynomial terms till x^10
data1$Taylor<-data1$x1-data1$x3/factorial(3)+data1$x5/factorial(5)-data1$x7/factorial(7)+data1$x9/factorial(9)

# Regression fit with polynomial terms till x^10
fit3<-lm(y~.,data=data1)
fit3<-step(fit3)
summary(fit3)

# Plot regression and Taylor Series with polynomial terms till x^10
plot(data1$x1,data1$y,type = "l",col="orange",xlab = "x",
     ylab="y",lwd=5,cex.axis=1.4,cex.lab=1.8,
     ylim = c(min(fit3$fitted.values),max(fit3$fitted.values)))
points(data1$x1,fitted(fit3),col='blue',pch=21,cex=1)
points(data1$x1,data1$Taylor,col='saddlebrown',pch=21,cex=1)
legend(6, -0.4, legend=c("Regression", "Taylor"),
       col=c("Blue", "saddlebrown"), lty=c(1,1), lwd=6,cex=0.8)

#plot residuals of fit3
plot(data1$x1,fit3$residuals,type = "l",col="red",xlab = "x",ylab="Residual",lwd=2)
#rm(fit3)

################################################################
#Add more polynomial terms to data i.e. till x^15
data1<-as.data.frame(cbind(y,poly(x1,degree = 15,raw=TRUE)))
colnames(data1)[-1] <- paste("x", colnames(data1)[-1], sep = "")

#Taylor Series (Maclaurin Expansion) with polynomial terms till x^15
data1$Taylor<-c(data1$x1-data1$x3/factorial(3)+data1$x5/factorial(5)
                -data1$x7/factorial(7)+data1$x9/factorial(9)
                -data1$x11/factorial(11)+data1$x13/factorial(13)
                -data1$x15/factorial(15))

# Regression fit with polynomial terms till x^15
fit4<-lm(y~.,data=data1)
fit4<-step(fit4)
summary(fit4)

#Plot regression and Taylor Series  with polynomial terms till x^15
plot(data1$x1,data1$y,type = "l",col="orange",xlab = "x",
     ylab="y",lwd=5,cex.axis=1.4,cex.lab=1.8)
points(data1$x1,fitted(fit4),col='blue',pch=21,cex=1)
points(data1$x1,data1$Taylor,col='saddlebrown',pch=21,cex=1)
legend(6, -0.4, legend=c("Regression", "Taylor"),
       col=c("Blue", "saddlebrown"), lty=c(1,1), lwd=6,cex=0.8)

#plot residuals of fit4
plot(data1$x1,fit4$residuals,type = "l",col="red",xlab = "x",ylab="Residual",lwd=2)
#rm(fit4)

################################################################
#Add more polynomial terms to data i.e. till x^25
data1<-as.data.frame(cbind(y,poly(x1,degree = 25,raw=TRUE)))
colnames(data1)[-1] <- paste("x", colnames(data1)[-1], sep = "")

#Taylor Series (Maclaurin Expansion) with polynomial terms till x^21
data1$Taylor<-c(data1$x1-data1$x3/factorial(3)+data1$x5/factorial(5)
                -data1$x7/factorial(7)+data1$x9/factorial(9)
                -data1$x11/factorial(11)+data1$x13/factorial(13)
                -data1$x15/factorial(15)+data1$x17/factorial(17)
                -data1$x19/factorial(19)+data1$x21/factorial(21)
                -data1$x23/factorial(23)+data1$x25/factorial(25)
                  )

# Regression fit with polynomial terms till x^17
fit5<-lm(y~.,data=data1)
fit5<-step(fit5)
summary(fit5)

#Plot regression and Taylor Series  with polynomial terms till x^17 & x^25
plot(data1$x1,data1$y,type = "l",col="orange",xlab = "x",
     ylab="y",lwd=5,cex.axis=1.4,cex.lab=1.8)
points(data1$x1,fitted(fit5),col='blue',pch=21,cex=1)
points(data1$x1,data1$Taylor,col='saddlebrown',pch=21,cex=1)
legend(6, -0.4, legend=c("Regression", "Taylor"),
       col=c("Blue", "saddlebrown"), lty=c(1,1), lwd=6,cex=0.8)

#plot residuals of fit5
plot(data1$x1,fit5$residuals,type = "l",col="red",xlab = "x",ylab="Residual",lwd=2)
#rm(fit5)

x1<-seq(-10*pi,10*pi,by=0.01)
y<-sin(x1)
data2<-as.data.frame(cbind(y,x1))
data2<-as.data.frame(cbind(y,poly(x1,degree = 25,raw=TRUE)))
colnames(data2)[-1] <- paste("x", colnames(data2)[-1], sep = "")
#Plot regression and Taylor Series  with polynomial terms till x^17 & x^25
plot(data2$x1,data2$y,type = "l",col="orange",xlab = "x",
     ylab="y",lwd=5,cex.axis=1.4,cex.lab=1.8)
points(data2$x1,predict(fit5,data2),col='blue',pch=21,cex=1)
points(data1$x1,data1$Taylor,col='saddlebrown',pch=21,cex=1)
legend(6, -0.4, legend=c("Regression", "Taylor"),
       col=c("Blue", "saddlebrown"), lty=c(1,1), lwd=6,cex=0.8)

#plot residuals of fit5
plot(data1$x1,fit5$residuals,type = "l",col="red",xlab = "x",ylab="Residual",lwd=2)
#rm(fit5)