function [beta,r,s2]=regre2new(Y,X,alfa,x0)
%
% Y = var. risposta
% X = var. predittrici
% alfa = liv.sign.
% x0 = nuova osservazione
%
% beta = stima con minimi quadrati
% r    = residuo con stima di beta
% s2   = stima corretta della varianza
%

format short e;
[n,p]=size(X);
X=[ones(n,1),X];
p1=p+1;
fprintf('Stima per beta \n')
beta=(X'*X)\(X'*Y);
disp(beta)
r=Y-X*beta;
normres2=r'*r;
fprintf('\n')
fprintf('Norma quadro del residuo , ||r||^2: %d \n', normres2)
fprintf('\n')
s2=r'*r/(n-p1);
fprintf('Stima della varianza sigma^2 : %d \n', s2)
fprintf('\n')
%fprintf('regione di confidenza\n')
 param=finv(1-alfa,p+1,n-p-1); %F_{p+1,n-p-1}(alfa)
 c2=p1*s2*param;
 varbeta=s2*inv(X'*X);

fprintf('T2-intervalli di confidenza\n')
fprintf('        i       beta_min       beta_max \n')
for k=1:p1,
disp([k,beta(k)-sqrt(varbeta(k,k))*sqrt(p1*param),...
      beta(k)+sqrt(varbeta(k,k))*sqrt(p1*param)])
%[beta(2)-sqrt(varbeta(2,2))*sqrt(p1*param),...
%     beta(2)+sqrt(varbeta(2,2))*sqrt(p1*param)]
%[beta(3)-sqrt(varbeta(3,3))*sqrt(p1*param),...
%      beta(3)+sqrt(varbeta(3,3))*sqrt(p1*param)]
end

R2=1-(r'*r)/norm(Y-mean(Y)*ones(n,1))^2;
fprintf('\n')
fprintf('Coeff. determinazione R^2: %d \n',R2)
fprintf('\n')


yhat=X*beta;
subplot(1,2,1)
plot(yhat,r,'*'); xlabel('y');ylabel('r')
subplot(1,2,2)
histfit(r)
%fprintf('intervallo di stima\n')
figure(2)
for k=1:p
subplot(1,p,k)
plot(X(:,k+1),r,'*')
end

x0=[1;x0];
yhat=x0'*beta;
fprintf('Stima per nuova osservazione : %d \n',yhat)
fprintf('\n')
tt=tinv(1-0.05/2,n-p-1);
fprintf('Intervallo di confidenza per nuova osservazione\n')
disp([x0'*beta-tt*sqrt(s2*(1+x0'*((X'*X)\x0))), x0'*beta+tt*sqrt(s2*(1+x0'*((X'*X)\x0)))])

