
graphics_toolkit("qt")

close all

randn("seed",1)

# skutecny model

N=100
x=linspace(0,5,N)';
ymod=5*x.^2-x.^3+1;

# zasumeni dat
sigma=ymod/5;
y=ymod+randn(size(ymod)).*sigma;

figure(1)
plot(x,y,"o")


# bazove fce

baze=[x.^0 x.^1 x.^2 x.^3];


b=y./sigma;

# fit funkcemi 1,x

M=2
baze1=baze(:,1:M);
A=baze1./sigma;

pseudoinv=(A'*A)^(-1)*A';
a=pseudoinv*b;
yfit=baze1*a;

chi2=sum(((y-yfit)./sigma).^2)
pval=gammainc((N-M)/2,chi2/2)



figure(2)
plot(x,y,"o",x,yfit,"b-","linewidth",2)

#fit funkcemi 1,x,x^2

M=3
baze2=baze(:,1:M);
A=baze2./sigma;

pseudoinv=(A'*A)^(-1)*A';
a=pseudoinv*b;
yfit=baze2*a;

chi2=sum(((y-yfit)./sigma).^2)
pval=gammainc((N-M)/2,chi2/2)


figure(3)
plot(x,y,"o",x,yfit,"b-","linewidth",2)



#fit funkcemi 1,x,x^2,x^4

M=4
baze3=baze(:,1:M);
A=baze3./sigma;

pseudoinv=(A'*A)^(-1)*A';
a=pseudoinv*b;
yfit=baze3*a;

chi2=sum(((y-yfit)./sigma).^2)
pval=gammainc((N-M)/2,chi2/2)



h=figure(4)
plot(x,y,"o",x,yfit,"b-","linewidth",1.5,x,ymod,"r--","linewidth",1.5)
axis("tight")
xlabel("x")
ylabel("y")
print -dpdf  -mono "-S320,250" fig1.pdf


h=figure(5)
xx=linspace(0,160,100);
plot(xx,gammainc((N-M)/2,xx/2),"b-","linewidth",2,
[chi2, chi2],[0, 1],"r--","linewidth",1.5)
axis("tight")
set(gca,'xtick',[0,40,80,120,160]);
set(gca,'ytick',[0,0.2,0.4,0.6,0.8,1]);
xlabel("chi2")
ylabel("P")
print -dpdf  -mono "-S320,250" fig2.pdf


# vysledne koeficienty

a

# chyby

kovar=(A'*A)^(-1);
chyby=sqrt(diag(kovar))





