% Test example for polynomial regression
x = 0.1:0.1:1;
Data = x.^3 - 0.1*x  + 1;
%DataC = Data+0.02*randn(1,10);

xP = 0:0.01:1.2;

close all
figure(1);
%subplot(2,2,1);
plot(x,DataC,'+');
set(gca,'XLim',[0,1.4]);
title('y=x^3-0.1x+1 + 0.1 N(0,1)');
set(gcf,'Position',[100 100 350 250]);
set(gcf,'PaperPositionMode','auto');
print -depsc2 f0


ind = 2;
for n=1:5
	Pn = polyfit(x,DataC,n);
	devn = DataC-polyval(Pn,x);
	Yn = polyval(Pn,xP);
	
	%	subplot(2,2,ind);
	figure
	ind = ind+1;
	plot(x,DataC,'+');
	hold on;
	plot(xP,Yn);
	title(['Order=' num2str(n) ', MSE=' num2str(devn*devn'/10)]);
	set(gcf,'Position',[100 100 350 250]);
	set(gcf,'PaperPositionMode','auto');
	eval(['print -depsc2 f' num2str(n)]);
	close all
	mse(n) = devn*devn'/10;
end

for o=1:5
	D = DataC;
	for n=1:o
		D = [D; x.^n];
	end
	[th,tll]=arx([D;ones(size(x))]);
	marll(o) = tll;
end
emar = exp(marll);
figure;
subplot(1,2,1);
bar(emar/sum(emar));
title('Probability of the correct order');

subplot(1,2,2);
plot(-marll/100);
hold on
plot(1:5,mse*10,'r');
title('Bayes factor and MSE')

set(gcf,'Position',[100 100 600 250]);
set(gcf,'PaperPositionMode','auto');

print -depsc2 prob