clear
format long

idx=(0:0.001:1)';
f=idx.^2;
niter=1000000;
ghat_true=idx-0.25;

ghat=0*idx;
ghat2=ghat;
R=0;
VAR=0;

for i=1:niter
    if (rem(i,1000)==0)
       i
    end
    x=rand(2,1);
    y=x.^2;
    a=sum(x);
    b=-prod(x);
    g=(a*idx+b);
    ghat=ghat+g/niter;
    ghat2=ghat2+g.^2/niter;
    R=R+mean((g-f).^2)/niter;
end

l=plot(idx,idx.^2,':',idx,ghat_true,idx,ghat,'--');
legend(l,'f(x)','ghat(x) analytic','ghat experimental',4)
title('f(x) and ghat(x)')
xlabel('x')
print -depsc ghatplot
bias=mean((ghat-f).^2);
bias_true=mean((ghat_true-f).^2);
VAR=mean(ghat2-ghat.^2);
[R,1/30]
[bias,1/80]
[VAR,1/48]