
disp('Bootstrap results for regressions');
disp('Compare 95% confidence intervals and P-values for the slope')
disp('  for least-squares and Theil`s regression with boostrap');
disp('  P-values and bootstrap percentile confidence intervals for');
disp('     (i) bootraps of residuals');
disp('    (ii) bootraps of observations');

disp(' ');
disp('Initializing the uniform random-number generator at');
Startseed=12345678;
fprintf('    %d\n',Startseed);
rand('state',Startseed);

disp('Initializing the normal random-number generator at');
Startseed=12345678;
fprintf('    %d\n',Startseed);
randn('state',Startseed);
fprintf('(In MATLAB, rand() and randn() must both be seeded.)\n');

disp(' ');
xx=rand(1,5);
fprintf('Five U(0,1)s:  %g   %g   %g   %g   %g\n', xx);
xx=randn(1,5);
fprintf('Five N(0,1)s:  %g   %g   %g   %g   %g\n', xx);

disp(' ');
beta_true=5;  mu_true=10;  sigma=100;  xmax=20;
nn=30;  tt30=2.04841;  % 5% critical value for t(30)
disp('Generating (X,Y) for a normal regression');
fprintf('  with  Y = %g X + %g + %g*N(0,1), n=%d pairs\n', ...
  beta_true, mu_true, sigma, nn);

xx = xmax*rand(nn,1);
yy = beta_true*xx + mu_true + sigma*randn(nn,1);
%% Round to two places beyond decimal point
xx = floor(100*xx+0.50)/100;
yy = floor(100*yy+0.50)/100;

fprintf('\nData (n=%d)  (X then Y):\n',nn);
%% The following assumes THAT nn IS EVEN!
n2=nn/2; 
xx1=xx(1:n2);  yy1=yy(1:n2);  xx2=xx(n2+1:nn);  yy2=yy(n2+1:nn);
fprintf('    %2d. %6.2f  %8.2f        %2d. %6.2f  %8.2f\n', ...
 [ (1:n2)'  xx1  yy1  (n2+1:nn)'  xx2  yy2]'); %'

% Are there any ties?
[midrx  xties] = getmidranks(xx);
xg=0;
for i=1:nn;  
  if xties(i)>0 xg=i;  else break;  end;  end;
[midry  yties] = getmidranks(yy);
yg=0;
for i=1:nn;  
  if yties(i)>0 yg=i;  else break;  end;  end;
fprintf('\nWithin-sample tie groups:  X: %d   Y: %d\n', xg, yg);
if xg==0 && yg==0  disp('NO TIES in Xs or Ys');
elseif xg>0  nt=0;
  for i=1:nn;  for j=i+1:nn;
    if xx(i)==xx(j)
      fprintf('  X(%d)=X(%d)=%g changed to  %g  %g\n', ...
        i,j, xx(i), xx(i)-0.004, xx(j)+0.004);
      xx(i)=xx(i)-0.004;  xx(j)=xx(j)+0.004;  nt=nt+1;  end;
    end;  end;
  fprintf ('  %d X tie(s) corrected.\n',nt);
  [midrx xties] = getmidranks(xx);
  if xties>0  fprintf('\nTOO MANY TIES!\n');  return;  end;
  end;

         

%% LEAST-SQUARES REGRESSION

ymean=sum(yy)/nn;  xmean=sum(xx)/nn;
yyc=yy-ymean;  xxc=xx-xmean;
ssxx = sum(xxc.*xxc);
beta = sum(yyc.*xxc)/ssxx;
mu = ymean - beta*xmean;
fprintf('\nLEAST-SQUARES REGRESSION:  Y = %g X + %g\n', beta,mu);

resid = yy - beta*xx - mu;
disp(' ');
disp('Normal-theory test of H_0:beta=0  (Student-t test)');
mse = sum(resid.*resid)/(nn-2);
sigbeta=sqrt(mse/ssxx);
tt = beta/sigbeta;  df=nn-2;
fprintf('  beta_hat=%6.4f   T=%6.4f   P=%6.4f  (two-sided, df=%d)\n', ...
  beta, tt, pvtdist(tt,df), df);
fprintf('  95%% Conf.Int. for beta:  (%6.4f   %6.4f   %6.4f)\n', ...
  beta-tt30*sigbeta, beta, beta+tt30*sigbeta);

%% Bootstrap procedures
nboot=10000;

% Offsets for 95% percentile credible intervals
blow=floor(nboot*0.025);  if blow==0 blow=1;  end;
bhigh=ceil(nboot*0.975);
bboot=zeros(nboot,1);

fprintf('\nFor nb=%d bootstrap replications:\n',nboot);

fprintf('(LEAST SQUARES) BOOTSTRAP ON RESIDUALS (beta=%g)\n',beta);
resid = yy - beta*xx - mu;  %% This repeats this definition
for ns=1:nboot;
  %% Bootstrap sample of residuals
  residstar=resid(ceil(nn*rand(nn,1)));
  %% Observed values with new residuals
  yystar = beta*xx + mu + residstar;
  %% Since X_i are the same, xxstarc = xxc
  bboot(ns) = sum(yystar.*xxc)/ssxx;
  end;
bboot=sort(bboot);
fprintf('  95%% Cred.Interval  (%6.4f   %6.4f   %6.4f)\n', ...
  bboot(blow), median(bboot), bboot(bhigh));

%% Find bootstrap bias
mb=0;
for b=1:nboot;  if bboot(b)<=beta mb=mb+1;  end;  end;
fprintf('  Bootstrap bias: %d/%d = %6.4f\n', mb,nboot, mb/nboot);

%% Find bootstrap P-values for H_0:beta=0 assumed observed beta>0
fprintf ('  Bootstrap P-value for H_0:beta=0:\n');
ns=0;
for b=1:nboot;  if bboot(b)<=0 ns=ns+1;  end;  end;
pp=ns/nboot;  pwid=1.960*sqrt(pp*(1-pp)/nboot);
fprintf ('    One-sided P:  P=%d/%d=%7.5f\n', ns,nboot, pp);
fprintf ('    Two-sided P:  P=%7.5f   95%% CI: (%7.5f, %7.5f)\n', ...
   2*pp, 2*(pp-pwid), 2*(pp+pwid));


fprintf('(LEAST SQUARES) BOOTSTRAP ON OBSERVATIONS (beta=%g)\n',beta);
ww = [xx yy];
for ns=1:nboot;
  %% Bootstrap X,Y pairs --- Xs now change as well
  wstar = ww(ceil(nn*rand(nn,1)),:);
  xstar=wstar(:,1);  ystar=wstar(:,2);
  xstarc=xstar-sum(xstar)/nn;
  bboot(ns) = sum(ystar.*xstarc)/sum(xstarc.*xstarc);
  end;
bboot=sort(bboot);
fprintf('  95%% Cred.Interval  (%6.4f   %6.4f   %6.4f)\n', ...
  bboot(blow), median(bboot), bboot(bhigh));

mb=0;
for b=1:nboot;  if bboot(b)<=beta mb=mb+1;  end;  end;
fprintf('  Bootstrap bias: %d/%d = %6.4f\n', mb,nboot, mb/nboot);

fprintf ('  Bootstrap P-value for H_0:beta=0:\n');
ns=0;
for b=1:nboot;  if bboot(b)<=0 ns=ns+1;  end;  end;
pp=ns/nboot;  pwid=1.960*sqrt(pp*(1-pp)/nboot);
fprintf ('    One-sided P:  P=%d/%d=%7.5f\n', ns,nboot, pp);
fprintf ('    Two-sided P:  P=%7.5f   95%% CI: (%7.5f, %7.5f)\n', ...
   2*pp, 2*(pp-pwid), 2*(pp+pwid));



%% THEIL REGRESSION

%% Set  Delta Y/(Delta X)  array
sratio=zeros(1,nn*nn+3);  r=0;
for i=1:nn;  for j=i+1:nn;
  xdiff=xx(j)-xx(i);
  % Since X_i have no ties, this always happens
  if abs(xdiff)>1e-8 r=r+1;  sratio(r)=(yy(j)-yy(i))/xdiff;  end;
  end;  end;
% Theil`s estimate of beta
% There are no X ties, so nslo=nn*(nn-1)/2
nslo=r;
sratio=sratio(1:nslo);   %% DON`T COUNT unused values in sratio()
beta=median(sratio);

% Estimate mu as the median of Walsh averages of Y(i)-beta*X(i)
resid = yy - beta*xx;
wsums=zeros(1,nn*nn+3);  nw=0;
for i=1:nn;  for j=i:nn;
  nw=nw+1;  wsums(nw)=(resid(i)+resid(j))/2;  end;  end;
wsums=wsums(1:nw);    % DON`T COUNT unused values
mu=median(wsums);

fprintf('\n\nTHEIL REGRESSION:  Y = %g X + %g\n', beta,mu);

%% Find Theil`s 95% large-sample approx. CI for beta
varC = nn*(nn-1)*(2*nn+5)/18;
stdC = sqrt(varC);
cca = floor(1.960*stdC);  %% Smallest value consistent
%% X_i have no ties, so there are n*(n-1)/2 slopes
mm=nn*(nn-1)/2;
mlow  = (mm - cca)/2;        if mlow<1 mlow=1; end;
mhigh = ((mm + cca)/2) + 1;  if mhigh>mm  mhigh=mm;  end;
sratio=sort(sratio);
betamed=median(sratio);
disp('  Theil`s 95% Confidence interval (with median of S-ratios):');
fprintf('    (%6.4f    %6.4f    %6.4f)\n', ...
  sratio(mlow), betamed,  sratio(mhigh))
fprintf('  with offsets  %d  %g  %d  in sratio array  (n=%d)\n', ...
  mlow, mm/2, mhigh, mm);

disp(' ');
disp('Theil`s test of H_0:beta=0  (same as Kendall`s test)');
disp('(NOT assuming that X(i) are increasing)');
cc=0;
for i=1:nn;  for j=i+1:nn;
  zz = (yy(j)-yy(i))*(xx(j)-xx(i));
  if (zz>0) cc=cc+1;
  elseif zz<0 cc=cc-1;  end;  end;  end;
zz = cc/stdC;   %% See above
fprintf('  C=%d   Z=%6.4f   P=%6.4f (Assuming no ties)\n', ...
  cc, zz, normpval(zz));


%% Bootstrap procedures
nboot=10000;
% For 95% credible intervals
% Repeat in case nboot has changed
blow=floor(nboot*0.025);  if blow==0 blow=1;  end;
bhigh=ceil(nboot*0.975);
bboot=zeros(nboot,1);

fprintf('\nFor nb=%d bootstrap replications:\n',nboot);
fprintf('(THEIL) BOOTSTRAP ON RESIDUALS (beta=%g)\n',beta);
resid = yy - beta*xx - mu;
for ns=1:nboot;
  residstar=resid(ceil(nn*rand(nn,1)));
  ystar = beta*xx + mu + residstar;
  % Compute betastar
  nd=0;
  for i=1:nn;  for j=i+1:nn;
    xdiff=xx(j)-xx(i);
    if abs(xdiff)>1e-8
      nd=nd+1;  sratio(nd)=(ystar(j)-ystar(i))/xdiff;  end;
  end;  end;
  % Theil`s estimate of beta
  bboot(ns) = median(sratio(1:nd));
  end;
bboot=sort(bboot);
betamed=median(bboot);
fprintf('  95%% Cred.Interval  (%6.4f   %6.4f   %6.4f)\n', ...
  bboot(blow), betamed, bboot(bhigh));

mb=0;
for b=1:nboot;  if bboot(b)<=beta  mb=mb+1;  end;  end;
fprintf('  Bootstrap bias: %d/%d = %6.4f\n', mb,nboot, mb/nboot);

ns=0;
for b=1:nboot;  if bboot(b)<=0 ns=ns+1;  end;  end;
fprintf ('  Bootstrap P-value for H_0:beta=0:\n');
pp=ns/nboot;  pwid=1.960*sqrt(pp*(1-pp)/nboot);
fprintf ('    One-sided P:  P=%d/%d=%7.5f\n', ns,nboot, pp);
fprintf ('    Two-sided P:  P=%7.5f   95%% CI: (%7.5f, %7.5f)\n', ...
   2*pp, 2*(pp-pwid), 2*(pp+pwid));


fprintf('(THEIL) BOOTSTRAP ON OBSERVATIONS (beta=%g)\n',beta);
ww = [xx yy];
for ns=1:nboot;
  %% Bootstrap rows of X,Y data
  wstar = ww(ceil(nn*rand(nn,1)),:);
  xstar=wstar(:,1);  ystar=wstar(:,2);
  % Compute betastar
  nd=0;
  for i=1:nn;  for j=i+1:nn;
    xdiff=xstar(j)-xstar(i);
    if abs(xdiff)>1e-8
      nd=nd+1;  sratio(nd)=(ystar(j)-ystar(i))/xdiff;  end;
    end;  end;
  if nd==0 ns=ns-1;  continue;  end;
  % Theil`s estimate of beta
  bboot(ns) = median(sratio(1:nd));
  end;
bboot=sort(bboot);
betamed=median(bboot);
fprintf('  95%% Cred.Interval  (%6.4f   %6.4f   %6.4f)\n', ...
  bboot(blow), betamed, bboot(bhigh));

mb=0;
for b=1:nboot;  if bboot(b)<=beta  mb=mb+1;  end;  end;
fprintf('  Bootstrap bias: %d/%d = %6.4f\n', mb,nboot, mb/nboot);

ns=0;
for b=1:nboot;  if bboot(b)<=0 ns=ns+1;  end;  end;
fprintf ('  Bootstrap P-value for H_0:beta=0:\n');
pp=ns/nboot;  pwid=1.960*sqrt(pp*(1-pp)/nboot);
fprintf ('    One-sided P:  P=%d/%d=%7.5f\n', ns,nboot, pp);
fprintf ('    Two-sided P:  P=%7.5f   95%% CI: (%7.5f, %7.5f)\n', ...
   2*pp, 2*(pp-pwid), 2*(pp+pwid));



disp(' ');

