% 
% Written by:
% -- 
% John L. Weatherwax                2009-04-21
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

close all; drawnow; clear functions; 
rehash; 
clc; 

xpy = [ 10; 40; 40 ];  % the vector X' y 
M   = [ 10, 0, 0 ;     % the matrix X' X 
        0, 20, 0; 
        0, 0, 40 ]; 

beta = M \ xpy;        % estimate of beta

bpxpy = beta' * xpy;   % beta * X' * y 

sum_y_square = 165; n = 10; y_bar = 1; p = 2; 

SSR = bpxpy - n*y_bar^2;                         % sum squared regression
SSE = 165 - 2 * bpxpy + beta' * M * beta;        % sum squared error 
SSTO = sum_y_square - n * y_bar^2;               % sum squared total 

R2 = SSR / SSTO;                                 % the R^2 statstic

MSR = SSR / p; 
MSE = SSE / (n-p-1); 

F_ratio = MSR / MSE; 
fprintf('F_ratio = %10.6f\n', F_ratio); 

MI = inv(M)

dM = diag(MI); 

s2 = SSE/ ( n - p - 1 );                         % an estimate of \sigma^2 (the error variance)

beta_hat_se = sqrt(s2) * sqrt( dM )              % an estimate of the standard error of each \beta_i estimate

fprintf('the t statistic of the null hypothsis beta_i=0\n');
t = beta ./ beta_hat_se                          

alpha=0.05; 
t_star = tinv( 1 - 0.5 * alpha, n-p-1 ); 
fprintf('t_star = %10.6f\n', t_star ); 

f_star = finv( 1 - alpha, p, n-p-1 ); 
fprintf('f_star = %10.6f\n', f_star );