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

 trials = 1000;            % # experiments to run
 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

 Qhat = zeros(p,trials);
 Sig0sqhat = zeros(trials,1); 

 for k=1:trials

    e = sqrt(sig0sq)*randn(n,1);  % generate noise with variance sig0sq

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

    Qhat(:,k) = (X'*X)\(X'*Y);    % least squares solution
    Rhat = Y - X*Qhat(:,k);       % residual

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

 end

 qhatmean = sum(Qhat,2)/trials;   % average all q_hat solutions 
 fprintf(' Mean of q_hat over %d trials: \n', trials)
 disp(qhatmean)


 fprintf(' Mean of variance estimates (sigma_0_hat)^2:  %10.7f\n', mean(Sig0sqhat))
 fprintf('              True variance     (sigma_0)^2:  %10.7f\n', sig0sq)
 
