% IAP 2006 Introduction to MATLAB
% Linear Algebra and Calculus

% EXAMPLE 2 Curve Fitting
% A. Fit a straight line, a quadratic, and a cubic equation through a set of points (x,y).
% A straight line is y = c1*x + c2
% A quadratic is y = c1*x^2 + c2*x +c3
% A cubic is y = c1*x^3 + c2*x^2 + c3*x + c4

% A1. Define point coordinates 
% You can use any set of points; for example, let's use a set of five points:
clear all
X = [2 3 5 6 8]'
Y = [75 86 82 70 40]'

% A2. Plot the points to see them
plot (X, Y, 'ro')

% A3. Fit a straight line, a quadratic, and a cubic line

% Two methods
% Method 1 - solve a system of linear equations
% A = [X^n ... X 1], then C = A \ Y
% For a straight line:
A1 = [X   ones(5,1)]
C1 = A1 \ Y
% For a quadratic:
A2 = [X.^2   X   ones(5,1)]
C2 =  A2 \ Y
% For a cubic:
A3 = [X.^3   X.^2   X   ones(5,1)]
C3 = A3 \ Y

% Method 2 - using the polyfit function
% Straight line
C1 = polyfit(X, Y, 1); C1 = C1'
% Quadratic
C2 = polyfit(X, Y, 2); C2 = C2'
% Cubic
C3 = polyfit(X, Y, 3); C3 = C3'

% A4. Plot the points and the curves (optional - a preview for Graphics)

% Plot the points as red circles
plot (X, Y, 'ro')
hold on
% plot the straight line as a blue line
x = [0 : 0.01 : 10];
y1 = polyval (C1, x);
plot (x, y1, 'b-')
hold on
% plot the quadratic as a green curve
y2 = polyval (C2, x);
plot (x, y2, 'g-')
hold on
% plot the cubic as a yellow curve
y3 = polyval (C3, x);
plot (x, y3, 'y-')
hold off

% A5. Save the computed curves (we shall use them later in Graphics)
save('x_coord.txt', 'x', '-ascii'); 
save('y_line.txt', 'y1', '-ascii'); 
save('y_quadratic.txt', 'y2', '-ascii'); 
save('y_cubic.txt', 'y3', '-ascii'); 

% B. Compute residuals to find the best fit

% B1. For a straight line fit
Y1 = polyval(C1, X)
% Alternatively:
Y1 = A1*C1
% Compute the sum of squared differences between actual and computed Y values
r1 = Y1 - Y
res1 = sum(r1.^2)

% B2. For a quadratic fit
Y2 = polyval(C2, X)
% Alternatively:
Y2 = A2*C2
% Compute the sum of squared differences between actual and computed Y values
r2 = Y2 - Y
res2 = sum(r2.^2)

% B3. For a cubic fit
Y3 = polyval(C3, X)
% Alternatively:
Y3 = A3*C3
% Compute the sum of squared differences between actual and computed Y values
r3 = Y3 - Y
res3 = sum(r3.^2)

% B4. Create a matrix of squared residual values
R = [1 res1; 2 res2; 3 res3]
% Which is the best fit?
% You can also type this command like this, for readability:
R = [1  res1
     2  res2
     3  res3]
% When you create matrices ONLY, MATLAB understands that pressing Return means
% the same as ; i.e. the end of a row - but ONLY while you are inside [].

% D. Do a polynomial fit of the 4th order. 
% This requires one MATLAB command only! Which one?
% WRITE IT HERE

% Is that a better fit?