function regression

global x ynoisy

% First do a linear regression example
% Assume y = a1+a2*x+a3*x2+a4*cos(x)+a5*sin(x)+a6*cos(2*x)+a7*sin(2*x)

fprintf('Linear regression example solved using linear least squares . . .\n\n')

coefs = [3; 0.5; -0.2; -35; 55; -10; 4]
a1 = coefs(1,1);
a2 = coefs(2,1);
a3 = coefs(3,1);
a4 = coefs(4,1);
a5 = coefs(5,1);
a6 = coefs(6,1);
a7 = coefs(7,1);

x = linspace(0,100,101)';
x2 = x.*x;
x3 = x2.*x;
cosx = cos(x);
sinx = sin(x);
cos2x = cos(2*x);
sin2x = sin(2*x);

y = a1+a2*x+a3*x2+a4*cosx+a5*sinx+a6*cos2x+a7*sin2x;

%Add noise to original data
noisemag = 20; % Try bigger and smaller values to see what happens
ynoisy = y+noisemag*(2*rand(101,1)-1);

plot(x,y,'k-',x,ynoisy,'kx')
legend('Original','Noisy')
title('Linear regression sample points')
pause

A = [ones(101,1) x x2 cosx sinx cos2x sin2x];
b = ynoisy;
coefs = A\b
a1 = coefs(1,1);
a2 = coefs(2,1);
a3 = coefs(3,1);
a4 = coefs(4,1);
a5 = coefs(5,1);
a6 = coefs(6,1);
a7 = coefs(7,1);

% Compute fitted value of y
yfit = a1+a2*x+a3*x2+a4*cosx+a5*sinx+a6*cos2x+a7*sin2x;
plot(x,yfit,'k-',x,ynoisy,'kx')
legend('Fitted','Noisy')
title('Linear regression fit')
pause

% Next do a nonlinear regression example three ways
% 1) Taylor series with iteration
% 2) Optimization
% 3) Transformation of variables
% Assume y = a1*exp(a2*x)

fprintf('Nonlinear regression example solved using three different approaches . . .\n\n')

coefs = [5; 0.5]
a1 = coefs(1,1);
a2 = coefs(2,1);

x = linspace(0,10,101)';
y = a1*exp(a2*x);

%Add noise to original data
noisemag = 20; % Try bigger and smaller values to see what happens
ynoisy = y+noisemag*(1-rand(101,1));

plot(x,y,'k-',x,ynoisy,'kx')
legend('Original','Noisy')
title('Nonlinear regression sample points')
pause

% Approach 1 - Taylor series with iteration
% f(a1new,a2new) = f(a1old,a2old) + df/da1(a1old,a2old)*(a1new-a1old)
%                                 + df/da2(a1old,a2old)*(a2new-a2old)
% where df/da1 = exp(a2*x)
%       df/da2 = a1*x*exp(a2*x)
% Steps
% 1) Guess a1old and a2old
% 2) Set f(a1new,a2new) in first order Taylor series approximation
%    and solve 101 linear equations in 2 unknowns for a1new and a2new
%    using linear least squares
% 3) Set a1new and a2new equal to a1old and a2old, respectively, and
%    iterate until convergence tolerance is met

fprintf('Approach 1: Taylor series with iteration\n')

a1old = 1;
a2old = 1;
tol = 1e-4;
ea = 100;
iter = 0;

while ea > tol
    % Calculate first derivative terms
    dfda1 = exp(a2old*x);
    dfda2 = a1old*x.*exp(a2old*x);
    
    % Form A matrix for redundant equations
    A = [dfda1 dfda2];
    
    % Form b vector for redundant equations
    yfit = a1old*exp(a2old*x);
    b = ynoisy-yfit+dfda1*a1old+dfda2*a2old;
    
    % Solve linear least squares problem for a1new and a2new
    coefs = A\b;
    a1new = coefs(1,1);
    a2new = coefs(2,1);
    
    % Calculate approximate percent error
    a1error = abs((a1old-a1new)/a1new)*100;
    a2error = abs((a2old-a2new)/a2new)*100;
    ea = max(a1error,a2error);
    
    % Set new values to old values
    a1old = a1new;
    a2old = a2new;
    iter = iter+1;
end

iter
coefs

% Compute fitted value of y
a1 = coefs(1,1);
a2 = coefs(2,1);
yfit = a1*exp(a2*x);
plot(x,yfit,'k-',x,ynoisy,'kx')
legend('Fitted','Noisy')
title('Nonlinear regression fit using Taylor series with iteration')
pause

% Approach 2 - Nonlinear least squares optimization

fprintf('Approach 2: Nonlinear least squares optimization\n')

options = optimset('Display','iter');
coefs0 = [1; 1];
coefs = lsqnonlin(@CostFunc,coefs0,[],[],options)

% Compute fitted value of y
a1 = coefs(1,1);
a2 = coefs(2,1);
yfit = a1*exp(a2*x);
plot(x,yfit,'k-',x,ynoisy,'kx')
legend('Fitted','Noisy')
title('Nonlinear regression fit using nonlinear least squares optimization')
pause

% Approach 3 - Transformation of variables

fprintf('Approach 3: Transformation of variables\n')

lny = log(ynoisy);
lna1 = log(a1);

A = [ones(101,1) x];
b = lny;
coefs = A\b;
coefs(1,1) = exp(coefs(1,1))

% Compute fitted value of y
a1 = coefs(1,1);
a2 = coefs(2,1);
yfit = a1*exp(a2*x);
plot(x,yfit,'k-',x,ynoisy,'kx')
legend('Fitted','Noisy')
title('Nonlinear regression fit using transformation of variablesn')

function errors = CostFunc(coefs)

global x ynoisy

a1 = coefs(1,1);
a2 = coefs(2,1);
yfit = a1*exp(a2*x);

errors = ynoisy-yfit;