
clear all
close all
%
disp('Coverage Probability by Arcsin Confidence Interval')
  lw = 2; 
  set(0, 'DefaultAxesFontSize', 16);
  fs = 15;
  msize = 10;
% end of preabmle--------------------------------------

n = 20;  %sample
alpha = 0.05; %confidence 1-alpha
M=100000;
covs=[];
for p=0.05:0.01:0.95
    q=1-p;
        X = binornd(n, p, [1 M]);
        pstar= (X+3/8)/(n + 3/4);  
    cov =  p <  (sin(asin(sqrt(pstar)) + norminv(1-alpha/2)/(2 * sqrt(n)) )).^2 & ...
                p >  (sin(asin(sqrt(pstar))  - norminv(1-alpha/2)/(2 * sqrt(n)) )).^2;
    covs=[covs sum(cov)/M];
end
figure(1)
plot(0.05:0.01:0.95, covs, 'linewidth', 2)
hold on
plot(0.05:0.01:0.95, covs, 'o')
plot(0.00:0.01:1, (1-alpha)*ones(1,101), 'r-','linewidth',3)
 xlabel('p'); ylabel('coverage');
% print -depsc 'C:\Springer\Estim\Estimeps\arcsincover.eps'