
disp('Water Diffusion across Human Chorioamnion Tissue');
disp('  ');
disp('Data from textbook Table 4.1 p110 (2nd edn)');
disp('  X: At term (approx 39weeks, Control)    Y: 12-26 weeks');
disp('  ');
disp('Wilcoxon Rank-Sum tests assumes Var(X)=Var(Y): Is this reasonable?');
disp('(See also Section 5.2 p158 in text.)');

% The data (same as in WcxRankSum.m):

xx = [ 0.80  0.83  1.89  1.04  1.45  1.38  1.91  1.64  0.73  1.46 ];

yy = [ 1.15  0.88  0.90  0.74  1.21 ];

mm=length(xx);    % Sample sizes
nn=length(yy);

% Display the two samples
% (There is no easy way to do this in MATLAB other than by using fprintf)
disp('  ');
disp('The data:');
disp('  ');
fprintf('X (m=%2d):',mm);
for i=1:mm;  fprintf('  %4.2f',xx(i));  end;  fprintf('\n');

disp('  ');
fprintf('Y (n=%2d):',nn);
for i=1:nn;  fprintf('  %4.2f',yy(i));  end;  fprintf('\n');

disp('  ');
disp('Let  Gamma=Var(X)/Var(Y)');

disp('  ');
disp('Observed VarX, VarY, Gamma=VarX/VarY, StdevX, StdevY:');
Mx=sum(xx)/mm;  xc=xx-Mx;  VarX=sum(xc.*xc)/(mm-1);  StdevX=sqrt(VarX);
My=sum(yy)/nn;  yc=yy-My;  VarY=sum(yc.*yc)/(nn-1);  StdevY=sqrt(VarY);

ObsGamma = VarX/VarY;

% Display VarX, VarY, Gamma=VarX/VarY, StdevX, StdevY
disp(' ');
disp([VarX VarY ObsGamma StdevX StdevY]);

disp('  ');
disp('Use Miller Jackknife/Bootstrap procedure to test  H_0:Gamma=1');
disp('  (H_0:Var(X)=Var(Y)) and find a 95% Conf.Int. for Gamma');
disp(' ');
disp('BOOTSTRAP TEST OF H_0:Gamma=1 (Var(X)=Var(Y)) and');
disp('  BOOTSTRAP 95% Confidence Interval for Gamma=Var(X)/Var(Y):');

disp(' ');
disp('Use a `stratified bootstrap` to generate an array');
disp('  of bootstrap values gamma_star = Var(Xstar)/Var(Ystar)');

% Initialize random-number seed
disp(' ');
disp('Starting random-number seed:');
disp(' ');
Startseed=1234567;
disp(Startseed);

rand('state',Startseed);

disp('Three example pairs of Xstar,Ystar bootstrap samples with Gamma=Var(X)/Var(Y):');
for Example=1:3;
  % Bootstrap resample Xstar from X
  xb=xx(ceil(mm*rand(1,mm)));
  for i=1:mm;  fprintf('  %4.2f',xb(i));  end;  fprintf('\n');
  % Bootstrap resample Ystar from Y
  yb=yy(ceil(nn*rand(1,nn)));
  for i=1:nn;  fprintf('  %4.2f',yb(i));  end;  fprintf('\n');
    xm=sum(xb)/mm;  xc=xb-xm;  xvar=sum(xc.*xc)/(mm-1);
    ym=sum(yb)/nn;  yc=yb-ym;  yvar=sum(yc.*yc)/(nn-1);
  % GammaStar = Var(Xstar)/Var(Ystar)
  fprintf('  Gamma=%g\n', xvar/yvar);
end;  % Of three trials

disp(' ');
disp('Number of bootstrap replications:');
disp(' ');
Nboot=10000;
disp(Nboot);

% Fill the bootstrap array
disp('Generate gamma-star bootstrap array:');
% Array for gamma-star bootstrap replication values
gamboot=zeros(1,Nboot);

for ns=1:Nboot;
  % Bootstrap resample for X
  xb=xx(ceil(mm*rand(1,mm)));
  xm=sum(xb)/mm;  xc=xb-xm;  xvar=sum(xc.*xc)/(mm-1);
  % Bootstrap resample for Y
  % Pick another resample if yvar<1e-8
  yvar=0;  % Fool the `while` command so that it always does at least one
  while yvar<1e-8
    yb=yy(ceil(nn*rand(1,nn)));
    ym=sum(yb)/nn;  yc=yb-ym;  yvar=sum(yc.*yc)/(nn-1);
  end;
  % gamma-star value
  gamboot(ns)=xvar/yvar;  end;

% Sort the bootstrap array
gamboot=sort(gamboot);

disp(' ');
disp('The first 5 sorted gamma-star values are:');
for i=1:5;  fprintf('  %.4f',gamboot(i));  end;
fprintf('\n');

disp(' ');
disp('The last 5 sorted gamma-star values are:');
for i=Nboot-4:Nboot;  fprintf('  %.4f',gamboot(i));  end;
fprintf('\n');

disp(' ');
disp('Median and Mean of bootstrapped values:');
Gmedian=median(gamboot);
Gmean=sum(gamboot)/Nboot;
disp([Gmedian Gmean]);

disp('Median bias (proportion of gamma-star <= Observed Gamma):');
mb=0;
for i=1:Nboot;  if (gamboot(i)<=ObsGamma) mb=mb+1;  end;  end;
fprintf('  %d/%d = %6.4f   ObsGamma=%g\n', mb,Nboot, mb/Nboot,ObsGamma);

disp(' ');
disp('Bootstrap (two-sided) P-value for H_0:gamma=1:');
disp('(Twice the proportion of gamma(star)<=1)');
mb=0;
for i=1:Nboot;  if (gamboot(i)<=1) mb=mb+1;  end;  end;
fprintf ('  2*%d/%d = %6.4f (2-sided)\n',mb,Nboot,2*mb/Nboot);

disp(' ');
disp('95% bootstrap percentile Conf.Int. bracketing Median');
gam_cilow=gamboot(floor(0.025*Nboot));
gam_cihigh=gamboot(ceil(0.975*Nboot));

disp([gam_cilow  Gmedian  gam_cihigh]);

disp(' ');
disp('ALTERNATIVE:  JACKKNIFE PROCEDURE (See Section 5.2 p158 in text)');
disp(' ');
disp('Define two sets of jackknife pseudo-values:');

% Define mm pseudo-values for X_i
xfull=log(VarX);  psx=zeros(1,mm);  xxd1=zeros(1,mm-1);
fprintf('Log variance of X sample: %.4f\n', xfull);
disp('Pseudo-values of log variance for X sample:');
disp('(These are Ps_i = mm*FullVarX - (mm-1)*Var(X_Delete_Xi) )');
for i=1:mm;
  % Put X-delete-one-value (delete X_i) in xxd1()
  r=0;  for j=1:mm;  if (j~=i) r=r+1;  xxd1(r)=xx(j);  end;  end;
  % X-delete-X_i has sample size mm-1, not mm
  xxm=sum(xxd1)/(mm-1);  xxc=xxd1-xxm;
  psx(i)=mm*xfull - (mm-1)*log(sum(xxc.*xxc)/(mm-2));  end;
% Display the pseudo-values
for i=1:mm;  fprintf('  %4.2f',psx(i));  end;  fprintf('\n');

% Now define nn pseudo-values for &_i
yfull=log(VarY);  psy=zeros(1,nn);  yyd1=zeros(1,nn-1);
disp(' ');
fprintf('Log variance of Y sample: %.4f\n', yfull);
disp('Pseudo-values of log variance for Y sample:');
for i=1:nn;
  r=0;  for j=1:nn;  if (j~=i) r=r+1;  yyd1(r)=yy(j);  end;  end;
  yym=sum(yyd1)/(nn-1);  yyc=yyd1-yym;
  psy(i)=nn*yfull - (nn-1)*log(sum(yyc.*yyc)/(nn-2));  end;
for i=1:nn;  fprintf('  %4.2f',psy(i));  end;  fprintf('\n');

disp(' ');
disp('Sample means and sample variances of the two sets');
disp('  of pseudo-values (of log variances):');
Jxmean=sum(psx)/mm;  psc=psx-Jxmean;  Jxvar=sum(psc.*psc)/(mm-1);
Jymean=sum(psy)/nn;  psc=psy-Jymean;  Jyvar=sum(psc.*psc)/(nn-1);

Jxmean_Jymean_Jxvar_Jyvar = [Jxmean Jymean Jxvar Jyvar]

disp('Unpooled variance estimator and standard deviation for Jxmean-Jymean:');
Jpoolvar = (Jxvar/mm) + (Jyvar/nn);
Jpool_Stdev = sqrt(Jpoolvar);
disp(' ');
disp([Jpoolvar Jpool_Stdev]);

disp('Miller 2-sample Z-test for H_0:Var(X)=Var(Y) (two-sided):');
Zval = (Jxmean-Jymean)/Jpool_Stdev;
disp(' ');
fprintf('  Z=%6.4f   P=%6.4f (two-sided)\n', Zval, normpval(Zval));

% `Bias-corrected` estimate of gamma=Var(X)/Var(y)
disp(' ');
disp('Bias-corrected gamma=Var(X)/Var(Y) estimator');

Jxdiff=Jxmean-Jymean;
% Recall that Jxmean etc are for log variances
Jgamma=exp(Jxdiff);

Exp_of_Jxmean_Minus_Jymean = Jgamma

disp('Jackknife 95% Conf.Int. for gamma bracketing Exp(Jxmean-Jymean):');
% Recall that Jxmean etc are for log variances
Jwid=1.960*Jpool_Stdev;

disp([exp(Jxdiff-Jwid)  Jgamma  exp(Jxdiff+Jwid)]);





