x = [41.53 18.73 2.99 30.34 12.33 117.52 73.02 223.63 4.00 26.78];
n = 10;
mu_hat = mean(x)
lambdahat = 1/mu_hat

%Parametrisk bootstrap
B = 1000;
lambdasim = zeros(1,B);
musim = zeros(1,B);
for b=1:B
    xstar = exprnd(1/lambdahat,1,n);
    lambdasim(b) = 1/mean(xstar);
    musim(b) = mean(xstar);
end
subplot(2,2,2)
xx = 0:0.005:0.12;
nn = histc(lambdasim,xx);
bar(xx,nn/(sum(nn)*0.005),1);
axis([0 0.12 0 60]);
xx = 0:10:150;
nn = histc(musim,xx);
subplot(2,2,1)
bar(xx,nn/(sum(nn)*0.005),1);
axis([0 150 0 60]);
parboot_bias = mean(lambdasim)-lambdahat
parboot_SE = sqrt(var(lambdasim))


%Ikke-parametrisk bootstrap
subplot(2,2,3)
cdfplot(x)

B = 10000;
lambdasim = zeros(1,B);
musim = zeros(1,B);
for b=1:B
    xstar = randsample(x,n,true);
    lambdasim(b) = 1/mean(xstar);
    musim(b) = mean(xstar);
end
subplot(2,2,4)
xx = 0:0.005:0.12;
nn = histc(lambdasim,xx);
bar(xx,nn/(sum(nn)*0.005),1);
axis([0 0.12 0 60]);
xx = 0:10:150;
nn = histc(musim,xx);
%bar(xx,nn/(sum(nn)*0.005),1);
%axis([0 150 0 60]);
nonparboot_bias = mean(lambdasim)-lambdahat
nonparboot_SE = sqrt(var(lambdasim))