
R=10000;
nVec=[round(exp(6)),round(exp(9)),round(exp(12)),round(exp(15)),round(exp(18))]';
nN=length(nVec);
storeHBar=zeros(R,nN);
storeHBar2=zeros(R,nN);
for nidx=1:nN
    n=nVec(nidx)
    for r=1:R
        X=randn(n,1);
        psi=1./(pi.*(1+(X.^2)));
        phi=exp(-(X.^2)./2)./sqrt(2*pi);
        laplace=exp(-sqrt((X.^2)))./2; % faster than abs()
        H=psi./phi;
        S=H'*ones(n,1); 
        hBar=S/n;
        storeHBar(r,nidx)=hBar;
        H2=laplace./phi;
        S2=H2'*ones(n,1);
        hBar2=S2/n;
        storeHBar2(r,nidx)=hBar2;
    end
end
indExceed=(storeHBar>ones(size(storeHBar)));
probExceed=mean(indExceed,1)
indExceed2=(storeHBar2>ones(size(storeHBar2)));
probExceed2=mean(indExceed2,1)

clear H
clear H2
clear laplace
clear phi
clear psi
clear X


fineGrid=(6:1:18)';
figure()
plot(log(nVec),probExceed,'-blacko','LineWidth',1.5,'MarkerSize',10,'MarkerFaceColor','black')
grid on
xlabel('$\log n$', 'Interpreter', 'Latex', 'FontSize', 14)
ylabel('P$(\hat{\theta}>0)$', 'Interpreter', 'Latex', 'FontSize', 14)
hold on
plot(fineGrid,(0.5./fineGrid),'--red','LineWidth',1.5)
legend('Simulated probabilities (Cauchy)', '$0.5/\log n$', 'Interpreter', 'Latex', 'Location','northeast', 'FontSize', 14)
xlim([5,19])
ylim([0.01, 0.1])


fineGrid2=(6:1:18)';
figure()
plot(log(nVec),probExceed2,'-blacko','LineWidth',1.5,'MarkerSize',10,'MarkerFaceColor','black')
grid on
xlabel('$\log n$', 'Interpreter', 'Latex', 'FontSize', 14)
ylabel('P$(\hat{\theta}>0)$', 'Interpreter', 'Latex', 'FontSize', 14)
hold on
plot(fineGrid2,1./sqrt(2.*fineGrid2),'--red','LineWidth',1.5)
legend('Simulated probabilities (Laplace)', '$1/\sqrt{2\log n}$', 'Interpreter', 'Latex', 'Location','northeast', 'FontSize', 14)
xlim([5,19])
ylim([0.14, 0.3])



