disp('Running MATLAB script file ScaSpaEg2.m') ;
%
%    Scale based Stationarity examples 2
%    Simulation study to check size for various simul types
%    This does not creat any graphical output,
%      instead storing summaries of number of significant pixels
%      in no signal cases
%
%    Caution:  this assumes an integer grid of times
%        which creates problems (especially tminew vs. nminew)
%        for non-integer times, as used in the internet traffic analysis
%        Should revise this to avoid that problem.


%addpath \matlab\stevesmo ;
%addpath \matlab\stevefun ;
addpath \matlabR12\work\SteveFun -end ;
addpath \matlabR12\work\SteveSmo -end ;


ipart = 2 ;    %  This only controls the output grid, 
               %  as in ScaStaEg1, where it meant:
               %  1 - Make several images (for different gaps)
               %  2 - Make movie through gaps

               
vn = [100 500 2500] ;               
               %  vector of smaple size
               

vinoise = [2 10] ;
                %  1 - no noise
                %  2 - lo Gaussian Noise, sig = 0.1
                %  3 - hi Gaussian Noise, sig = 1 
                %  4 - AR(1), 0.99, var = 0.1
                %  5 - AR(1), 0.99, var = 1
                %  6 - sin varying AR(1)
                %  7 - AR(1), 0.75, var = 0.1
                %  8 - N(0,sin varying), var = 0.1
                %  9 - sin varying mixture of heavy tail and bimodal
                %  10 - Cauchy, quartiles at +-0.1

                
visimul = [0 1 2 3 4] ;
                %  0 - no multiple comparison adjustment
                %  1 - adjust for multiple comparisons
                %  2 - Bonferroni adjustment
                %  3 - mixture
                %  4 - global simultaneous
               


vitest = [1 2 3 4 5 6 7 8] ;
               %  1 - means and t-test
               %  2 - Kolmogorov-Smirnov
               %  3 - variances and F-test
               %  4 - mean removed Kolmogorov-Smirnov
               %  5 - standardized Kolmogorov-Smirnov
               %  6 - median removed Kolmogorov-Smirnov
               %  7 - median and MAD st'dized Kolmogorov-Smirnov
               %  8 - Quartile st'dized Kolmogorov-Smirnov


               

%  Set common parameters
%  
nrepl = 1000 ;
          %  number of replications
ntest = length(vitest) ;
nsimul = length(visimul) ;
nwid = 15 ;
          %  number of widths to show at each scale
nloc = 100 ;
          %  number of locations to consider, for each scale
alpha = 0.05 ;
nminew = 5 ;
          %  minimum number in each window
          

randn('seed',21354979) ;
rand('seed',39045930) ;
          
          

for n = vn ;


  %  Set more parameters
  %
  left = 1 ;
  right = n ;
  vgap = (n / 100) * [0 5 20] + 1 ;
  ngap = length(vgap) ;
  mingap = min(vgap) ;

  
      
  %  loop through noise levels
  %
  for inoise = vinoise ;
    
    if inoise == 1 ;    %  don't add noise
      noisetit = 'no noise' ;
    elseif inoise == 2 ;    %  lo Gaussian Noise, sig = 0.1
      noisetit = 'N, s=0.1 noise' ;
    elseif inoise == 3 ;    %  hi Gaussian Noise, sig = 1
      noisetit = 'N, s=1 noise' ;
    elseif inoise == 4 ;    %  AR(1), 0.99, var = 0.1
      noisetit = '0.99 AR(1), v0.1' ;
    elseif inoise == 5 ;    %  AR(1), 0.99, var = 1
      noisetit = '0.99 AR(1), v1' ;
    elseif inoise == 6 ;    %  sin varying AR(1)
      noisetit = 'sin varying AR(1)' ;
    elseif inoise == 7 ;    %  AR(1), 0.75, var = 0.1
      noisetit = '0.75 AR(1), v0.1' ;
    elseif inoise == 8 ;    %  N(0,sin varying), var = 0.1
      noisetit = 'N, sin var.' ;
    elseif inoise == 9 ;    %  sin varying mixture of heavy tail and bimodal
      noisetit = 'sin mixed dist' ;
    elseif inoise == 10 ;    %  Cauchy, quartiles at +-0.1
      noisetit = 'Cauchy, Q13+-0.1' ;
    end;    %  of noise addition if-block


    
    
    %  Setup up cell arrays for storage
    %
    cellred = cell(ntest,nsimul) ;
    cellcolor = cell(ntest,nsimul) ;
    cellminpval = cell(ntest,nsimul) ;
    cellm = cell(1,nsimul) ;
    cellsalpha = cell(1,nsimul) ;

    isimulind = 0 ;
    for isimul = visimul ;
      isimulind = isimulind + 1 ;
 
      itestind = 0 ;
      for itest = vitest ;
        itestind = itestind + 1 ;
        
        cellred{itestind,isimulind) = zeros(ngap,nwid,nrepl) ;
        cellcolor{itestind,isimulind) = zeros(ngap,nwid,nrepl) ;
        cellminpval{itestind,isimulind) = ones(ngap,nwid,nrepl) ;
    
      end ;
      
      cellm{1,isimulind) = zeros(ngap,nwid) ;
      cellsalpha{1,isimulind) = zeros(ngap,nwid) ;
      
    end ;
    
    
    %  loop through replications of data
    %
    for irepl = 1:nrepl ;
  
    
    
      %  Generate this realization of the data
      %
      tgrid = linspace(left,right,n)' ;
      
      vsignal = zeros(n,1) ;
      
      %  Add noise
      %
      if inoise == 1 ;    %  don't add noise
        vx = vsignal;
                
      elseif inoise == 2 ;    %  lo Gaussian Noise, sig = 0.1
        sig = 0.1 ;
        vx = vsignal + sig * randn(n,1) ;
        
      elseif inoise == 3 ;    %  hi Gaussian Noise, sig = 1
        sig = 1 ;
        vx = vsignal + sig * randn(n,1) ;
        
      elseif inoise == 4 ;    %  AR(1), 0.99, var = 0.1
        phi = 0.99 ;
        sig0 = 0.1 ;
            %  marginal standard deviation
        sig = sig0 * sqrt(1 - phi^2) ;
            %  standard deviation of noise terms
        vnoise = sig * randn(n-1,1) ;
        varn = zeros(n,1) ;
        varn(1) = sig0 * randn(1,1) ;
        for i=2:n ;
          varn(i) = phi * varn(i-1) + vnoise(i-1) ;
        end ;
        vx = vsignal + varn ;
        
      elseif inoise == 5 ;    %  AR(1), 0.99, var = 1
        phi = 0.99 ;
        sig0 = 1 ;
            %  marginal standard deviation
        sig = sig0 * sqrt(1 - phi^2) ;
            %  standard deviation of noise terms
        vnoise = sig * randn(n-1,1) ;
        varn = zeros(n,1) ;
        varn(1) = sig0 * randn(1,1) ;
        for i=2:n ;
          varn(i) = phi * varn(i-1) + vnoise(i-1) ;
        end ;
        vx = vsignal + varn ;
        
      elseif inoise == 6 ;    %  sin varying AR(1)
        vphi = sqrt((1 - cos(2 * pi * tgrid / 5000)) / 2) ;
        sig0 = 0.1 ;
            %  marginal standard deviation
        vnoise = sig0 * randn(n-1,1) ;
        varn = zeros(n,1) ;
        varn(1) = sig0 * randn(1,1) ;
        for i=2:n ;
          varn(i) = vphi(i) * varn(i-1) + vnoise(i-1) ;
        end ;
        varn = 0.1 * varn / std(varn) ;
            %  put on same time scale as "lo noise"
        vx = vsignal + varn ;
                  
      elseif inoise == 7 ;    %  AR(1), 0.75, var = 0.1
        phi = 0.75 ;
        sig0 = 0.1 ;
            %  marginal standard deviation
        sig = sig0 * sqrt(1 - phi^2) ;
            %  standard deviation of noise terms
        vnoise = sig * randn(n-1,1) ;
        varn = zeros(n,1) ;
        varn(1) = sig0 * randn(1,1) ;
        for i=2:n ;
          varn(i) = phi * varn(i-1) + vnoise(i-1) ;
        end ;
        vx = vsignal + varn ;
        
      elseif inoise == 8 ;    %  N(0,sin varying), var = 0.1
        vsig = 0.2 * (1 - cos(2 * pi * tgrid / 500)) / 2 ;
        vx = vsignal + vsig .* randn(n,1) ;
      
      elseif inoise == 9 ;    %  sin varying mixture of heavy tail and bimodal
        vnoise = randn(n,1) ;
            %  N(0,1)
        vunif = rand(n,2) ;
            %  Uniform(0,1)
        vprobclass = (1 - cos(2 * pi * tgrid / 500)) / 2 ;
            %  Prob of choosing class 2
        classflag = vunif(:,1) > vprobclass ;
            %  1 for class 1
            %  0 for class 2
        mixflag = vunif(:,2) > 0.5 ;
            %  Uniform on {0,1}
        sig1 = 0.2 ;
        mu1 = sqrt(1 - sig1^2) ;
            %  ensures mean = 0, and var = 1
        sig21 = 0.2 ;
            %  smaller variance similar to above
        sig22 = sqrt(2 - sig21^2) ;
            %  ensures mean = 0, and var = 1
        vmu = classflag .* (mu1 * (2 * mixflag - 1)) ; 
            %  0's for class 2
            %  +- mu1 for class 1
        vsig = classflag .* (sig1) + ...
                  (1 - classflag) .* (mixflag .* sig21 + (1 - mixflag) .* sig22) ;
            %  sig1 for class 1
            %  mixture of sig21 and sig22 for class 2
        vx = vsignal + vmu + vsig .* vnoise ;

        elseif inoise == 10 ;    %  Cauchy, quartiles at +-0.1
        q3 = 0.1 ;
        vx = vsignal + q3 * tan(pi * randn(n,1)) ;
        
      end;    %  of noise addition if-block

    
    


      %  loop through simul levels
      %
      isimulind = 0 ;
      for isimul = visimul ;

        isimulind = isimulind + 1 ;
      
        if isimul == 0 ;   
         simultit = 'not simultaneous' ;      
        elseif isimul == 1 ; 
          simultit = 'simultaneous' ;      
        elseif isimul == 2 ; 
          simultit = 'Bonferroni' ;    
        elseif isimul == 3 ; 
          simultit = 'mixture' ;  
        elseif isimul == 4 ; 
          simultit = 'global' ;  
          %mglobal = 2 * nloc - 1 ;
          mglobal = 2 * min((right - left)/nminew,nloc) - 1 ;
          salphaglobal = 1 - (1 - alpha)^(1/mglobal) ; 
        end;
         
              
              
              
        %  loop through tests
        %
        itestind = 0 ;
        for itest = vitest ;
        
          itestind = itestind + 1 ;
        
          if itest == 1 ;    %  1 - means and t-test
            testtit = 'Mean t-test' ;      
          elseif itest == 2 ;    %  2 - Kolmogorov-Smirnov
            testtit = 'Kolm-Smir' ;      
          elseif itest == 3 ;    %  3 - variances and F-test
            testtit = 'Var. F-test' ;      
          elseif itest == 4 ;    %  4 - mean removed Kolmogorov-Smirnov
            testtit = 'Mean R. K-S' ;      
          elseif itest == 5 ;    %  5 - standardized Kolmogorov-Smirnov
            testtit = 'Std''zed K-S' ;  
          elseif itest == 6;    %  6 - median removed Kolmogorov-Smirnov
            testtit = 'Median R. K-S' ;    
          elseif itest == 7 ;    %  7 - median and MAD st'dized Kolmogorov-Smirnov
            testtit = 'MAD S''d K-S' ;  
          elseif itest == 8 ;    %  8 - Quartile st'dized Kolmogorov-Smirnov
            testtit = 'Q S''d K-S' ;  
          end ;    %  end of itest if-block
    
  
        
          %  loop through gaps
          %
          for igap = 1:ngap ;
          
            gap = vgap(igap) ;
            go2 = gap / 2 ;
                              %  gap over 2
            
            
            %  loop through wid's
            %
            if ipart == 1 ;    %  just making static plot, so use full size
            
              minwid = log10(nminew) ;
              maxwid = log10((right - left - gap) / 2) ;
              vwid = logspace(minwid,maxwid,nwid+1)' ;
              vwid = vwid(1:nwid) ;
                    %  cut off last one, since doesn't add to plot
              maxwid = log10(vwid(nwid)) ;
                    %  readjust maxwid
                  
              minloc = (left + nminew + go2) ;
              maxloc = (right - nminew - go2) ;
            
            elseif ipart == 2 ;    %  just making static plot, so use full size
            
              minwid = log10(nminew) ;
              maxwid = log10((right - left - mingap) / 2) ;
              vwid = logspace(minwid,maxwid,nwid+1)' ;
              vwid = vwid(1:nwid) ;
                    %  cut off last one, since doesn't add to plot
              maxwid = log10(vwid(nwid)) ;
                    %  readjust maxwid
                  
              minloc = (left + nminew + mingap / 2) ;
              maxloc = (right - nminew - mingap / 2) ;
            
            end ;
              
              
            vloc = linspace(minloc,maxloc,nloc) ;
            for iwid = 1:nwid ;

              wid = vwid(iwid) ;

              if iwid == 1 ;
                disp(['  Working on wid ' num2str(iwid) ' = ' ...
                               num2str(wid) ...
                               '    for ' simultit ', ' testtit ', ' ...
                               dattit ', ' noisetit ', gap ' num2str(gap)]) ;
              end ;
                               
              locflag = left + wid + go2 <= vloc ;
              locflag = locflag & (vloc <= right - wid - go2) ;
              
              vindloc = (1:nloc) ;
              vindloc = vindloc(locflag) ;
                  
  
              if isimul == 0 ;    %  don't bother about simultaneous inference
                m = 0 ;
                salpha = alpha ;
              elseif isimul == 1 ;
                %m = (right - left) / (2 * wid) ;
                m = (right - left) / wid ;
                m = min(length(vindloc),m+1) ;
                salpha = 1 - (1 - alpha)^(1/m) ;

              elseif isimul == 2 ;
                m = 0 ;
                salpha = alpha / length(vindloc) ;
       
              elseif isimul == 3 ;
                m = (right - left) / (2*wid) ;
                m = min(length(vindloc),m) ;
                salpha = 1 - (1 - (alpha/nwid))^(1/m) ;

              elseif isimul == 4 ;
                m = mglobal ;
                salpha = salphaglobal ;

              end ;
  
              
              vflagred = zeros(nloc,1) ;
              vflagcolor = zeros(nloc,1) ;
              vpval =  ones(nloc,1) ;
              for iloc = vindloc ;

                loc = vloc(iloc) ;
                  
                flagleft = (loc - go2 - wid <= tgrid) ;
                flagleft = flagleft & (tgrid <= loc - go2) ;
                    %  ones where data inside left group

                flagright = (loc + go2 <= tgrid) ;
                flagright = flagright & (tgrid <= loc + go2 + wid) ;
                    %  ones where data inside right group

                dataleft = vx(flagleft) ;
                dataright = vx(flagright) ;


                if itest == 1 ;    %  means and t-test

                  meanleft = mean(dataleft) ;
                  meanright = mean(dataright) ;
                  varleft = var(dataleft) ;
                  varright = var(dataright) ;
                  
                  if meanleft == meanright ;                      

                    h = 0 ;
                    p = 1 ;
                        %  no significant difference

                  else ;
                    
                    if varright == 0 ;
                      if varleft == 0 ;
                        h = 1 ;
                        p = 0 ;
                           %  significant difference

                      else ;    %  do 1 sample t-test
                        [h,p] = ttest(dataleft,meanright,salpha) ;
                           %  1 for significant difference found
                           %  0 otherwise
                      end ;

                    else ;

                      if varleft == 0 ;    %  do 1 sample t-test
                        [h,p] = ttest(dataright,meanleft,salpha) ;
                           %  1 for significant difference found
                           %  0 otherwise

                      else ;    %  do full two sample t-test

                        dfleft = length(dataleft) - 1;
                        dfright = length(dataright) - 1;
                        dfe  = dfleft + dfright;
                        msleft = dfleft * varleft ;
                        msright = dfright * varright ;
                        
                        difference = meanleft - meanright;
                        pooleds    = sqrt((msleft + msright) * (1/(dfleft + 1) + 1/(dfright + 1)) / dfe);
                        
                        ratio = difference / pooleds ;
                        p = 1 - tcdf(ratio,dfe) ;
                            %  one tail p-value
                        p = 2 * min(p, 1-p) ;
                            %  two tail p-value

                        h = 0 ;
                        if p <= salpha, 
                          h = 1; 
                        end ;
                            %  1 for significant difference found
                            %  0 otherwise
                        
                      end ;
                    
                    end ;
                      
                  end ;
                    
                
                      
                elseif itest == 2 ;    %  Kolmogorov-Smirnov

                  [h,p] = kstest2(dataright,dataleft,salpha) ;
                      %  1 for significant difference found
                      %  0 otherwise


                elseif itest == 3 ;    %  variances and F-test
                
                  varleft = var(dataleft) ;
                  varright = var(dataright) ;
              
                  if varright == 0 ;
                    if varleft == 0 ;
                      h = 0 ;
                      p = 1 ;
                         %  no significant difference

                    else ;    %  do 1 sample t-test
                      h = 1 ;
                      p = 0 ;
                         %  significant difference

                    end ;

                  else ;

                    if varleft == 0 ;    %  do 1 sample t-test
                      h = 1 ;
                      p = 0 ;
                         %  significant difference

                    else ;    %  do variance F-test

                      dfleft = length(dataleft) - 1;
                      dfright = length(dataright) - 1;
    
                      Fstat = varleft / varright ;
    
                      p = 1 - fcdf(Fstat,dfleft,dfright) ;
                          %  one tail p-value
                      p = 2 * min(p, 1-p) ;
                          %  two tail p-value

                      h = 0;
                      if p <= salpha, 
                        h = 1; 
                      end ;
                                  %  1 for significant difference found
                                  %  0 otherwise
                      
                    end ;
                
                  end ;
                
           
                elseif itest == 4 ;    %  mean removed Kolmogorov-Smirnov

                  dataleftmr = dataleft - mean(dataleft) ;
                  datarightmr = dataright - mean(dataright) ;
                  
                  [h,p] = kstest2(datarightmr,dataleftmr,salpha) ;
                      %  1 for significant difference found
                      %  0 otherwise
                      %  Note:  have not adjusted significance for removal of mean

                    
                elseif itest == 5 ;    %  Standardized Kolmogorov-Smirnov

                  dataleftst = dataleft - mean(dataleft) ;
                  datarightst = dataright - mean(dataright) ;

                  varleft = var(dataleftst) ;
                  varright = var(datarightst) ;

                  if varleft > 0 ;
                    dataleftst = dataleftst / varleft ;
                  end ;

                  if varright > 0 ;
                    datarightst = datarightst / varright ;
                  end ;
                  
                  [h,p] = kstest2(datarightst,dataleftst,salpha) ;
                      %  1 for significant difference found
                      %  0 otherwise
                      %  Note:  have not adjusted significance for standardization

                      
                elseif itest == 6 ;    %  median removed Kolmogorov-Smirnov

                  dataleftmr = dataleft - median(dataleft) ;
                  datarightmr = dataright - median(dataright) ;
                  countleft = length(dataleft);
                  countright = length(dataright);
                  ncoefficient = 2 * ...
                        countleft * countright /(countleft + countright);
                  
                  % The sqrt(2) comes from the fact that I have 2 Brownian
                  % bridges on the interval of length .5. The same applies
                  % for the sqrt in the calculation of p-value.

                  %              
                  % Calculate F1(x) and F2(x), the empirical (i.e., sample) CDFs.
                  %
                  
                  binEdges    =  [-inf ; sort([dataleftmr;datarightmr]) ; inf];
                  
                  binCounts1  =  histc (dataleftmr , binEdges);
                  binCounts2  =  histc (datarightmr , binEdges);
                  
                  sumCounts1  =  cumsum(binCounts1)./sum(binCounts1);
                  sumCounts2  =  cumsum(binCounts2)./sum(binCounts2);
                  
                  sampleCDF1  =  sumCounts1(1:end-1);
                  sampleCDF2  =  sumCounts2(1:end-1);
                  
                  
                  %
                  % Compute the test statistic of interest.
                  %
                  
                  deltaCDF  =  abs(sampleCDF1 - sampleCDF2);
                  KSstatistic   =  max(deltaCDF);
                  lambda =  (sqrt(ncoefficient) + 0.12 + 0.11/sqrt(ncoefficient))...
                                                * KSstatistic ;
                  
                  
                  %
                  %  Use the asymptotic Q-function to approximate the 2-sided P-value.
                  %
                  j       =  [1:101]';
                  p  =  sqrt(2 * sum((-1).^(j-1).*exp(-2*lambda*lambda*j.^2)));

                  h = 0;
                  if p <= salpha, 
                    h = 1; 
                  end ;
                              %  1 for significant difference found
                              %  0 otherwise
                              %  Note:  HAVE adjusted significance for removal of median

                      
                elseif itest == 7 ;    %  median and MAD st'dized Kolmogorov-Smirnov

                  dataleftst = dataleft - median(dataleft) ;
                  datarightst = dataright - median(dataright) ;

                  gpmadleft = gpmad(dataleftst) ;
                  gpmadright = gpmad(datarightst) ;

                  if gpmadleft > 0 ;
                    dataleftst = dataleftst / gpmadleft ;
                  end ;

                  if gpmadright > 0 ;
                    datarightst = datarightst / gpmadright ;
                  end ;
                  
                  [h,p] = kstest2(datarightst,dataleftst,salpha) ;
                      %  1 for significant difference found
                      %  0 otherwise
                      %  Note:  have not adjusted significance for standardization

                      
                elseif itest == 8 ;    %  Quartile st'dized Kolmogorov-Smirnov

                  medianleft = median(dataleft) ;
                  medianright = median(dataright) ;
                  
                  dataleftst = dataleft - medianleft ;
                  datarightst = dataright - medianright ;

                  Qoneleft =  prctile(dataleftst,25) ;
                  Qoneright =  prctile(datarightst,25) ;

                  Qthreeleft =  prctile(dataleftst,75) ;
                  Qthreeright =  prctile(datarightst,75) ;

                  flagoneleft = dataleftst < 0 ;
                  flagoneright = datarightst < 0 ;

                  if Qoneleft < 0 ;    %  then rescale guys on Q1 side
                    dataleftst(flagoneleft) = phiinv(0.25) * ...
                                        dataleftst(flagoneleft) / Qoneleft ;
                  end ;

                  if Qoneright < 0 ;    %  then rescale guys on Q1 side
                    datarightst(flagoneright) = phiinv(0.25) * ...
                                        datarightst(flagoneright) / Qoneright ;
                  end ;

                  if Qthreeleft > 0 ;    %  then rescale guys on Q3 side
                    dataleftst(~flagoneleft) = phiinv(0.75) * ...
                                        dataleftst(~flagoneleft) / Qthreeleft ;
                  end ;

                  if Qthreeright > 0 ;    %  then rescale guys on Q3 side
                    datarightst(~flagoneright) = phiinv(0.75) * ...
                                        datarightst(~flagoneright) / Qthreeright ;
                  end ;
                  
                  countleft    =  length(dataleftst);
                  countright     =  length(datarightst);
                  ncoefficient      =  4*...
                         countleft * countright /(countleft + countright);

                  % The 4 comes from the fact that I have 4 Brownian
                  % bridges on the subintervals of length .25. The same applies
                  % for the sqrt in the calculation of p-value.
        
                         
                  %
                  % Calculate F1(x) and F2(x), the empirical (i.e., sample) CDFs.
                  %
                  
                  binEdges    =  [-inf ; sort([dataleftst;datarightst]) ; inf];
                  
                  binCounts1  =  histc (dataleftst , binEdges);
                  binCounts2  =  histc (datarightst , binEdges);
                  
                  sumCounts1  =  cumsum(binCounts1)./sum(binCounts1);
                  sumCounts2  =  cumsum(binCounts2)./sum(binCounts2);
                  
                  sampleCDF1  =  sumCounts1(1:end-1);
                  sampleCDF2  =  sumCounts2(1:end-1);
                  
                  
                  %
                  % Compute the test statistic of interest.
                  %
                  
                  deltaCDF  =  abs(sampleCDF1 - sampleCDF2);
                  KSstatistic   =  max(deltaCDF);
                  lambda =  (sqrt(ncoefficient) + 0.12 + 0.11/sqrt(ncoefficient))...
                                                * KSstatistic ;
                  
            
                  %
                  %  Use the asymptotic Q-function to approximate the 2-sided P-value.
                  %
                  j       =  [1:101]';
                  p  =  (2 * sum((-1).^(j-1).*exp(-2*lambda*lambda*j.^2)))^(0.25) ;

                  h = 0;
                  if p <= salpha, 
                    h = 1; 
                  end ;
                              %  1 for significant difference found
                              %  0 otherwise
                              %  Note:  HAVE adjusted significance for standardization

                                
                end ;    %  end of itest if-block
                

                vflagred(iloc) = h ;
                vflagcolor(iloc) = 1 ;
                vpval(iloc) = p ;


              end ;    %  of iloc loop

              
              %  Get summaries
              numred = sum(vflagred) ;
              numcolor = sum(vflagcolor) ;
              minpval = min(vpval) ;

              
              %  Store results
              %
              mnumred = cellred{itestind,isimulind} ;
              mnumred(igap,iwid,irepl) = numred ;
              cellred{itestind,isimulind} = mnumred ;

              mnumcolor = cellcolor{itestind,isimulind} ;
              mnumcolor(igap,iwid,irepl) = numcolor ;
              cellcolor{itestind,isimulind} = mnumcolor ;

              mminpval = cellminpval{itestind,isimulind} ;
              mminpval(igap,iwid,irepl) = minpval ;
              cellminpval{itestind,isimulind} = mminpval ;

              mm = cellm{1,isimulind) ;
              mm(igap,iwid) = m ;
              cellm{1,isimulind) = mm ;
              
              mm = cellsalpha{1,isimulind) ;
              msalpha(igap,iwid) = salpha ;
              cellsalpha{1,isimulind) = msalpha ;
              
              
            end ;    %  of iwd loop

          
          end ;    %  of igap loop

  
        end ;    %  of itest loop

  
      end ;    %  of isimul if-block


    end ;    %  of irepl loop



    %  create .mat file
    outstr = ['\MatlabR12\work\ScaSta\SimOuts\ScaStaEg2' ...
                     'p' num2str(ipart) ...
                     'n' num2str(n) ...
                     'nois' num2str(inoise) ...
                                        ] ;

                                        
    eval(['save ' outstr 'cellred cellcolor cellminpval vitest visimul ' ...
                         'nloc alpha nminew nrepl nwid ngap;']) ;
    

    
        
  end ;    %  of inoise loop


end ;    %  of n loop