clear;

%x = [ 1 2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19  20]';
%y = [ 2.6888   -0.4724    3.7286    5.6471    1.1164    4.3160    5.2080  -0.3875    0.0181    4.1423    2.3002    4.5800    4.9312    4.8238    6.0805    4.9372    6.0817    1.3951    3.8604    3.6866]';

x = [3  7 11 15 18 27 29 30 30 31 31 32 33 33 34 36 36 36 37 38 39 39 39 40 41 42 42 43 44 45 46 47 50]'
y = [5 11 21 16 16 28 27 25 35 30 40 32 34 32 34 37 38 34 36 38 37 36 45 39 41 40 44 37 44 46 46 49 51]'
%x = [0.0 0.0 0.0 0.0 0.0 2.3 2.3 2.3 2.3 2.3 4.6 4.6 4.6 4.6 4.6 9.2 9.2 9.2 9.2 9.2 18.4 18.4 18.4 18.4 18.4]'
%y = [10.9 10.6 10.8 9.8 9.0 11.0 11.3 9.9 9.2 10.1 10.6 10.4 8.8 11.1 8.4 9.7 7.8 9.0 8.2 2.3 2.9 2.2 3.4 5.4 8.2]'
n = length(x);
plot(x,y, '.');
title('Regresión lineal');
xlabel('x');
ylabel('y');
pause
d = [x, ones(n, 1)];

p = inv(d'*d)*d'*y;
m = p(1);
b = p(2);
yp = mean(y);

v1 = 1;
v2 = n-2;

SSR = (yp -m*x - b)'*(yp -m*x - b)
SSE = (y -m*x - b)'*(y -m*x - b);
SST = SSR + SSE;
s2 = SSE/v2;
F = SSR/s2;
alfa = 0.05;
Fc = finv(1-alfa, v1, v2);
fprintf('---------------------------------------------------------------------\n');
fprintf('Curva ajustada\n');
fprintf('y = %13.3f *x + %13.3f\n', m, b);
fprintf('---------------------------------------------------------------------\n');
fprintf('Analisis de Varianza\n');
fprintf('---------------------------------------------------------------------\n');
fprintf('Fuente de     Suma            Grados de      Cuadrado        F calc\n');
fprintf('Variacion     cuadrados       Libertad          medio         \n');
fprintf('---------------------------------------------------------------------\n');
fprintf('Regresion %13.2f  %13.f  %13.2f %13.2f\n', SSR, v1, SSR/v1, F);
fprintf('Error     %13.2f  %13.f  %13.2f \n', SSE, v2, s2);
fprintf('Total     %13.2f  %13.f  \n', SST, v1+v2);
fprintf('---------------------------------------------------------------------\n');
fprintf('Con (1-alfa) = %5.3f tenemos %13.2f > %13.2f\n', 1-alfa, F, Fc);


f = 0:0.05:8;
plot(f, fpdf(f, v1, v2));
title('Distribución F');
xlabel('f');
ylabel('h(f)');
[m, v] =fstat(6,24)

fa = finv(1-alfa, 6, 24)

hold on;
y = 0:0.1:0.8;
plot(ones(length(y),1)*fa, y, 'r');
