load('/home/ingelins/Documents/Fag/Statistikk/TMA4245V11/Data/StudDataV2011.txt')
%Finn indeks til dei med oppgitt høgde på både mor og far
f=find(1-isnan(StudDataV2011(:,3)) & (1-isnan(StudDataV2011(:,4))));
kjonn=StudDataV2011(f,1)-1;
%Hogde barn, jenter justert til gutar
hogdeBarn=StudDataV2011(f,2).*(1+0.08*kjonn);
%Gj.snitt hogda til foreldra. Mor*1.08
hogdeForeldre=0.5*StudDataV2011(f,3)+0.5*1.08*StudDataV2011(f,4);

plot(hogdeForeldre, hogdeBarn,'.')
xlabel('Høgde foreldre','FontSize',16)
ylabel('Høgde barn','FontSize',16)
%estimerer regresjonsparametra
[b,bint,r,rint,stats]=regress(hogdeBarn,[ones(length(f),1),hogdeForeldre]);
hold
plot([160,200],[b(1)+b(2)*160, b(1)+b(2)*200],'r','LineWidth',2)

%plotter hogdebarn=hogdeforeldre linje
plot([160,200],[160,200],'k','LineWidth',2)

a=b(1)
bb=b(2)
sigmae=stats(4)
sqrt(sigmae)



%Predikerer voksenhogde for eit barn med gj.snitt foreldrehogde 182.5
x0=182.5;
ypred=a+bb*x0


%Finn Pred.varians
Sxx=sum((hogdeForeldre-hogdeBarn).^2)
predvar=sigmae*(1+1/length(hogdeForeldre)+(x0-mean(hogdeBarn)).^2/Sxx)
sqrt(predvar)

%Finn 95% pred.intervall
ypred-1.96*sqrt(predvar)
ypred+1.96*sqrt(predvar)


%Finn KI for barn med gj.snitt foreldrehogde 182.5
muvar=sigmae*(1/length(hogdeForeldre)+(x0-mean(hogdeBarn)).^2/Sxx)
sqrt(muvar)
%Finn 95% KI
ypred-1.96*sqrt(muvar)
ypred+1.96*sqrt(muvar)


%Finn KI for mangen x0
x0=160:1:200;
ypred=a+bb*x0;
muvar=sigmae*(1/length(hogdeForeldre)+(x0-mean(hogdeBarn)).^2/Sxx);
KInedre=ypred-1.96*sqrt(muvar);
KIovre=ypred+1.96*sqrt(muvar);
plot(x0,KIovre,'r--', 'LineWidth',2)

plot(x0,KInedre,'r--','LineWidth',2)


%Finn PI for mangen x0
x0=160:1:200;
ypred=a+bb*x0;
predvar=sigmae*(1+1/length(hogdeForeldre)+(x0-mean(hogdeBarn)).^2/Sxx);
PInedre=ypred-1.96*sqrt(predvar);
PIovre=ypred+1.96*sqrt(predvar);
plot(x0,PIovre,':g','LineWidth',2)
plot(x0,PInedre,':g','LineWidth',2)


%Hogde far og hogde mor
far=StudDataV2011(f,3);
mor=StudDataV2011(f,4);
figure
plot(far,mor,'.')

xlabel('Høgde far','FontSize',16)
ylabel('Høgde mor','FontSize',16)

[b,bint,r,rint,stats]=regress(mor,[ones(length(f),1),far]);
hold
plot([150,200],[b(1)+b(2)*150, b(1)+b(2)*200],'r','LineWidth',2)

