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")
regression**3**( 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**