clear all
close all

 disp('Cyclosporine Example')
 lw = 3; 
 set(0, 'DefaultAxesFontSize', 14);
 fs = 14;
 msize = 10;

%Cyclosporine Clinical Trials
%
n = 37; %subjects in cyclosporine group
%three months 
X1 = 22;     p1hat = X1/n;   q1hat = 1-p1hat;
% six monthe
X2 = 7;       p2hat = X2/n;   q2hat = 1- p2hat;
%===============================
%Wald Intervals @ 95% confidence
W3 = [p1hat  - norminv(0.975) * sqrt( p1hat * q1hat / n), ...
    p1hat  + norminv(0.975) * sqrt( p1hat * q1hat / n)];
W6 = [p2hat  - norminv(0.975) * sqrt( p2hat * q2hat / n), ...
    p2hat  + norminv(0.975) * sqrt( p2hat * q2hat / n)];
 %===========  
W3
W6
% plot
    figure(1)
    plot( [1,1],W3, 'k-', 'linewidth', 10)  
    hold on
  
%==================================    
% Wilson Intervals  
% three months
    p1hats = (X1+2)/(n+4);   q1hats = 1-p1hats;
% six monthe
     p2hats = (X2+2)/(n+4);   q2hats = 1- p2hats;
   Wi3 = [p1hats  - norminv(0.975) * sqrt( p1hats * q1hats / (n+4)), ...
    p1hats  + norminv(0.975) * sqrt( p1hats * q1hats / (n+4))];
   Wi6 = [p2hats  - norminv(0.975) * sqrt( p2hats * q2hats / (n+4)), ...
    p2hats  + norminv(0.975) * sqrt( p2hats * q2hats / (n+4))];
Wi3
Wi6
% plot
    plot( [2,2],Wi3, 'r-', 'linewidth', 10)  

%==========================
%Wilson Score- 
z=norminv(0.975);

Wis3 = [  1/(1 + z^2/n) *  (p1hat  + z^2/(2 * n)  -  z * sqrt( p1hat * q1hat / n  + z^2/(4 * n^2))), ...
    1/(1 + z^2/n) * (p1hat  + z^2/(2 * n)  + z * sqrt( p1hat * q1hat / n  + z^2/(4 * n^2)))];
Wis6 = [  1/(1 + z^2/n) *  (p2hat  + z^2/(2 * n)  -  z * sqrt( p2hat * q2hat / n  + z^2/(4 * n^2))), ...
    1/(1 + z^2/n) * (p2hat  + z^2/(2 * n)  + z * sqrt( p2hat * q2hat / n  + z^2/(4 * n^2)))];
% plot
    plot( [3,3],Wis3, 'g-', 'linewidth', 10)  
  
%=========================
% Clopper-Pearson
       Fs = finv(0.975, 2*(n-X1 + 1), 2*X1);
       Fss = finv(0.975, 2*(X1+1), 2*(n-X1));
      CP3 = [ X1/(X1 + (n-X1+1).*Fs),  (X1+1).*Fss./(n - X1 + (X1+1).*Fss)];

       Fs = finv(0.975, 2*(n-X2 + 1), 2*X2);
       Fss = finv(0.975, 2*(X2+1), 2*(n-X2));
      CP6 = [ X2/(X2 + (n-X2+1).*Fs),  (X2+1).*Fss./(n - X2 + (X2+1).*Fss)];
      CP3
      CP6
% plot
    plot( [4,4],CP3, 'm-', 'linewidth', 10)  
 
 
%==========================================
% Anscombe ARCSIN intervals
%
p1h = (X1 + 3/8)/(n + 3/4);  p2h = (X2 + 3/8)/(n + 3/4); 

AA3 = [ (sin( asin( sqrt(p1h)) - norminv(0.975)/(2 * sqrt(n)) ))^2 , ...
           (sin( asin( sqrt(p1h)) + norminv(0.975)/(2 * sqrt(n)) ))^2];
AA6 = [ (sin( asin( sqrt(p2h)) - norminv(0.975)/(2 * sqrt(n)) ))^2 , ...
           (sin( asin( sqrt(p2h)) + norminv(0.975)/(2 * sqrt(n)) ))^2];
     AA3
     AA6
% plot
    plot( [5,5], AA3, '-', 'linewidth', 10)  
   
  axis([0 17 0 1])
  text(0.8,0.8,'3 month follow-up')
  % 
   plot( [6,6],W6, 'k-', 'linewidth', 10) 
   plot( [7,7],Wi6,  'r-', 'linewidth', 10)  
   plot( [8,8],Wis6,  'g-', 'linewidth', 10)  
   plot( [9,9],CP6,  'm-', 'linewidth', 10)  
   plot( [10,10],AA6,  '-', 'linewidth', 10)  
   
  text(5.8,0.4,'6 month follow-up')
  legend('1:Wald_3','2:Wilson_3','3:WScore_3','4:Clopper_3','5:ArcSin_3','6:Wald_6',...
      '7:Wilson_6','8:WScore_6','9:Clopper_6','10:ArcSin_6')
  hold off
 % print -depsc 'C:/Springer/Estim/Estimeps/Cycloint.eps'


