%% MATLAB Examples 2 (covering Statistics Lectures 3 and 4)

%% Example 1: Fit a linearized regression model

% generate some random data
x = -20:20;
y = 1.5*x.^2 - 60*x + 10 + 30*randn(size(x));

% given these data, we will now try to fit a model
% to the data.  we assume we know what model to use,
% namely, y = ax^2 + bx + c where a, b, and c are
% free parameters.  notice that this was the model
% we used to generate the data.

% since our model is a linearized model, we can express
% it using simple matrix notation.

% construct the regressor matrix
X = [x(:).^2 x(:) ones(length(x),1)];

% estimate the free parameters of the model
% using ordinary least-squares (OLS) estimation
w = inv(X'*X)*X'*y(:);

% what are the estimated weights?
w
%%

% what is the model fit?
modelfit = X*w;

% what is the squared error achieved by this model fit?
squarederror = sum((y(:)-modelfit).^2);

% visualize
figure;
hold on;
scatter(x,y,'r.');
ax = axis;
xx = linspace(ax(1),ax(2),100)';
yy = [xx.^2 xx ones(length(xx),1)]*w;
plot(xx,yy,'k-');
xlabel('x');
ylabel('y');
title(sprintf('squared error = %.1f',squarederror));
%%

%% Example 2: Fit a parametric nonlinear model

% generate some random data
x = 0:.05:3;
y = 5*x.^3.5 + 5*randn(size(x));

% given these data, we will now try to fit a model
% to the data.  we assume we know what model to use,
% namely, y = ax^n where a and n are free parameters.
% notice that this was the model we used to generate the data.

% define optimization options
options = optimset('Display','iter','FunValCheck','on', ...
                   'MaxFunEvals',Inf,'MaxIter',Inf, ...
                   'TolFun',1e-6,'TolX',1e-6);

% define bounds for the parameters
%              a    n
paramslb = [-Inf    0];  % lower bound
paramsub = [ Inf  Inf];  % upper bound

% define the initial seed
%            a     n
params0 = [  1     1];

% estimate the free parameters of the model using nonlinear optimization
  % hint: pp is a vector of parameters like [a n],
  %       xx is a vector of x-values like [.1 .4 .5],
  %       and what modelfun does is to output the predicted y-values.
modelfun = @(pp,xx) pp(1)*xx.^pp(2);
[params,resnorm,residual,exitflag,output] = lsqcurvefit(modelfun,params0,x,y,paramslb,paramsub,options);
%%

% what are the estimated parameters?
params
%%

% what is the model fit?
modelfit = modelfun(params,x);

% what is the squared error achieved by this model fit?
squarederror = sum((y(:)-modelfit(:)).^2);
% (note that the 'resnorm' output also provides the squared error.)

% visualize
figure;
hold on;
scatter(x,y,'r.');
ax = axis;
xx = linspace(ax(1),ax(2),100)';
yy = modelfun(params,xx);
plot(xx,yy,'k-');
xlabel('x');
ylabel('y');
title(sprintf('squared error = %.1f',squarederror));
%%

%% Example 3: Another optimization example

% generate some random data
x = rand(1,1000).^2;

% we want to determine the value that minimizes the sum of the 
% absolute differences between the value and the data points

% define some stuff
options = optimset('Display','iter','FunValCheck','on', ...
                   'MaxFunEvals',Inf,'MaxIter',Inf, ...
                   'TolFun',1e-6,'TolX',1e-6);
paramslb = [];  % [] means no bounds
paramsub = [];
params0 = [0];

% perform the optimization
  % hint: here we have to apply a square-root transformation so that when 
  % lsqnonlin squares the output of costfun, we will be back to absolute error.
costfun = @(pp) sqrt(abs(x-pp));
[params,resnorm,residual,exitflag,output] = lsqnonlin(costfun,params0,paramslb,paramsub,options);
%%

% what is the solution?
params
%%

% what is the median of the data?
median(x)
%%

% notice that the solution is nearly identical to the median of the data.
% (the slop comes from the tolerance used in the optimization and
% from the interpolation used in the calculation of the median.)
