source("http://macosa.dima.unige.it/r.R")    # If I have not already loaded the library
---------- ---------- ---------- ---------- ---------- ---------- ---------- ----------
S 12 Regression. Splines. Curve fitting.
 
Linear, quadratic and cubic regression:  the linear function with graph for a given
point (p,q) and the generic linear, quadratic and cubic functions that "best"
approximate some experimental data (x,y).  Commands:
 
regression(x,y,p,q), regression1(x,y), regression2(x,y), regression3(x,y)
 
If the points are 2, 3 or 4 regression1, regression2 and regression3 find exactly the
curve.   We also see how splines can be calculated.
 
For a more complete method of studyng linear regression see here.

x = c(-0.2, 1.3, 1.7,  3.2,  4.1,  4.6)
y = c(9.1, 17.5, 20.7, 25.4, 35.7, 51.2)
Plane(-1,6,-10,70);  POINT(x,y,1)
regression(x,y, 0,0); POINT(0,0,"blue")
# 9.829 * x
f = function(x) 9.829 * x; graph(f,-1,6, "red")
regression1(x,y)
# 7.628 * x + 7.911
g = function(x) 7.628 * x + 7.911; graph(g,-1,6, "blue")
regression2(x,y)
# 1.37 * x^2 + 1.38 * x + 11.2 
h = function(x) 1.37 * x^2 + 1.38 * x + 11.2; graph(h,-1,6, "violet")
regression3(x,y)
# 1.21 * x^3 + -6.3 * x^2 + 11.9 * x + 11.6 
k = function(x) 1.21 * x^3 + -6.3 * x^2 + 11.9 * x + 11.6; graph(k,-1,6, "green")
POINT(0,0,"blue"); POINT(x,y,1)
 
          
 
# On the right is the cubic spline that interpolates the points (see):
 
Plane(-1,6,-10,70);  POINT(x,y,1); L=spline(x,y,1000)
 
# I calculated and put in L the 1000 points of the spline passing through the
# points x,y
 
polyline(L$x,L$y,"red")                              #  I traced it in red
 
# The "extended" (blue) curve that passes by the same points, lowered by 10.
 
POINT(x,y-10,1); L1=spline(x,y-10,1000, xmin=-1,xmax=5)
polyline(L1$x,L1$y, "blue")
#
# The curve that goes exactly for (-1,-10),(0,-1),(1,-2),(2,5)
 
Plane(-2,3,-20,20)
POINT(-1,-10, 1); POINT(0,-1,"red"); POINT(1,-2,"green"); POINT(2,5,"blue")
regression3( c(-1,0,1,2), c(-10,-1,-2,5) )
# 3 * x^3 + -5 * x^2 + 1 * x + -1
f = function(x) 3 * x^3 + -5 * x^2 + 1 * x + -1
graph(f,-2,3, 6)
POINT(-1,-10, 1); POINT(0,-1,"red"); POINT(1,-2,"green"); POINT(2,5,"blue")
                
#
# See here to see other methods of curve fitting.
  
# # How to proceed for a regression of the second order when I have a fixed point? # Without using specific methods, I can do the following: x = c(-0.2, 1.3, 1.7, 3.2, 4.1, 4.6) y = c(9.1, 17.5, 20.7, 25.4, 35.7, 51.2) Plane(-1,6,-10,70); POINT(x,y,"black") # Without a fixed point: regression2(x,y) # 1.37 * x^2 + 1.38 * x + 11.2 h = function(x) 1.37*x^2 + 1.38*x + 11.2; graph2(h,-1,6, "magenta") # If I have the fixed point (3,20) I can force the passage for a large amount of points # (for example 1e6) that are equal to (3.20)! x1=c(rep(3,1e6),x); y1=c(rep(20,1e6),y) regression2(x1,y1) # 2.95 * x^2 + -6.54 * x + 13 h1 = function(x) 2.95 * x^2 + -6.54 * x + 13 graph2(h1,-2,7, "brown"); POINT(3,20, "red") # # Another example of use of regression: # A car (initially) goes at a steady speed. The driver sees an obstacle and (after 1.5 # seconds) starts braking. Here are the photos taken every half second from the instant # in which he saw the obstacle. What was the initial speed of the car? # Times (in seconds) and spaces (in meters): s0 = c(1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5,6.0,6.5,7.0,7.5,8.0,8.5,9.0,9.5) m0 = c( 70, 92,114,133,151,168,183,198,209,221,230,238,244,250,253,255,257) # Since he has begun to brake: s = s0-1.5; m = m0-70; s; m # 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 # 0 22 44 63 81 98 113 128 139 151 160 168 174 180 183 185 187 # The graph (time,space): Plane(-1/2,8.5, -1,190); POINT(s,m, "blue") abovex("sec"); abovey("m") # I want to impose that s(0)=0. I proceed as in the last example: s1 = c(rep(0,1e6), s); m1 = c(rep(0,1e6), m) regression2(s1,m1) # -2.89 * x^2 + 46.5 * x + -2.89e-07 # I neglect -2.89e-07 F = function(x) -2.89 * x^2 + 46.5 * x graph2(F,0,8, "brown") # See the chart above # I find the derivative (I can find it "by hand" but ...) deriv(F,"x") # 46.5 - 2.89 * (2 * x) # The initial speed is 46.5 m/s. 46.5/1000*60*60 # 167.4 167 km/h Other examples of use