clear all
clc
format short
 X = [25 25 28 32 32 32 38 42 48 51 51 58 62 65];
Y = [180 195 186 180 210 197 239 183 204 221 243 208 222 269];

% X = [45 50 55 60 65 70 75]
% Y = [24.2 25 23.3 22 21.5 20.6 19.8]

%X = [3 5 7 8 10 11 12 12 13 15 15 16 18 19 20];
 %Y = [7 20 20 15 25 17 20 35 26 25 35 32 44 37 45];

% X = [60 70 80 85 95];
% Y = [70 65 70 95 85];

sum_x = sum(X)
sum_y = sum(Y)

xsqr = X.*X
xy = X.*Y
ysqr = Y.*Y

sum_xsqr = sum(xsqr)
sum_xy = sum(xy)
sum_ysqr = sum(ysqr)

n = length(X);

S_xx = sum_xsqr - (n*((sum_x/n)^2))
S_xy = sum_xy - (n*(sum_x/n)*(sum_y/n))
S_yy = sum_ysqr - (n*((sum_y/n)^2))

B = S_xy/S_xx
A = (sum_y/n) - B*(sum_x/n)

% Step II
% Build intervals of confidence

SS_R = (S_xx * S_yy - (S_xy^2)) / S_xx
t_cr = 1.77;

left = B - (sqrt(SS_R/((n-2)*S_xx)) * t_cr)
right = B + (sqrt(SS_R/((n-2)*S_xx)) * t_cr)

beta = 0;
if( beta >= left)
    disp('stop')
else
    disp('yes, beta =/= 0, go on')
end

for k = 1:n
    y_cap(k) = A+B*X(k);
    res(k) = Y(k) - y_cap(k);
end
y_cap
res
res_sqr = res.*res;
sum_res_sqr = sum(res_sqr)
sum_res = sum(res)
mean_appx =  (sum(res))/n;
if (mean_appx >= 0.1)
    disp('stop'), disp(mean_appx)
else
    disp('go on'), disp(mean_appx)
end


% Step III
% Calculate Proportion

i_uno = (1-0.375)/(n+0.25);
prop = ((1:n)-1)/(n+0.25) + i_uno
z_score(1:n) = NORMINV(prop(1:n),0,1)
 ordered_res = sort(res)

 sum_z = sum((z_score.*z_score));

mid_flag = 0;
for k = 1:length(prop)
    if(prop(k) >= 0.50 & ~mid_flag)
        mid = k;
        mid_flag = 1;
    end
    if(mid_flag)
        break;
    end
    
end
mid
%mid = find(prop == 0.500);

if(mod(n,2))
    t(1:mid-1) = abs(z_score(1:mid-1)).*(abs(ordered_res(1:mid-1)) + abs(ordered_res(n:-1:mid+1)));
else
    t(1:mid-1) = abs(z_score(1:mid-1)).*(abs(ordered_res(1:mid-1)) + abs(ordered_res(n:-1:mid)));
end

top = sum(t)
rb = sum_res/n;
zb = sum(z_score)/n;
bleft = sqrt(sum_res_sqr - n*(rb^2));
bright = sqrt(sum_z - n*(zb^2));

top/(bleft*bright)