clear

% Enter Data - Hubble Data
x=[0.032 0.034 0.214 0.263 0.275 0.275 0.45 0.5 ... 
    0.5 0.63 0.8 0.9 0.9 0.9 0.9 1 1.1 1.1 1.4 ...
    1.7 2 2 2 2]; 
y=[170 290 -130 -70 -185 -220 200 290 270 200 ...
    300 -30 650 150 500 920 450 500 500 960 500 ...
    850 800 1090];
figure(1)

% Polynomial Fit
n = 1;                       % n = deg polynomial fit
[m,b]=polyfit(x,y,n);        % m = [slope intercept]

% Plot (xfit, yfit) for model
xfit = linspace(0,2,20);
yfit = polyval(m,xfit);
plot(x,y,'ok',xfit,yfit);
xlabel(['d (Mpc)'])
ylabel(['v (km/s)'])

% Residuals
ymodel = polyval(m,x);      % predicted values
resid  = y - ymodel;        % residuals
 
% R^2 Value (Rsq)
dev = y - mean(y);          % deviations - measure of spread
SST = sum(dev.^2);          % total variation accounted for
SSE = sum(resid.^2);        % variation NOT accounted for
Rsq = 1 - SSE/SST;          % percent of error explained

% F-Statistic
SSR = SST - SSE;            % variation of residual
dfe = length(x) - 1 - n;    % degrees of freedom
MSE = SSE/dfe;              % mean-square error of residuals
MSR = SSR/n;                % mean-square error for regression
f = MSR/MSE;                % f-statistic

% t-Statistics for coeffs from b in polyfit
R = b.R;                
covm = (R'*R)\eye(n+1);     % The covariance matrix
d = diag(covm)';            % ROW vector of the diagonal elements
MSE = (b.normr^2)/b.df;     % variance of the residuals
se = sqrt(MSE*d);           % the standard errors
t = m./se;                  % observed T-values 