% Example of computing mean / variance for a simple linear model
%
% Model:  y(t) = q_1 + q_2 t + q_3 t^2

 n = 100;                  % number of times

 t = linspace(0,1,n)';     % take measurements at these times
 X = [t.^0 t.^1 t.^2];     % construct design matrix

 p = size(X,2); 

 q0 = [1 2 3]';            % true solution

 sig0sq = .01;             % true variance of noise

 e = sqrt(sig0sq)*randn(n,1); 

 Y = X*q0 + e;             % noise-polluted "observations"

 qhat = (X'*X)\(X'*Y)      % least squares solution

 Rhat = Y - X*qhat;        % residual

 sig0sqhat = norm(Rhat)^2/(n-p);    % estimate of variance

 fprintf('   Estimate of variance (sigma_0_hat)^2:  %10.7f\n', sig0sqhat)
 fprintf('          True variance     (sigma_0)^2:  %10.7f\n', sig0sq)
 
