%MATLAB SCRIPT FILE s322eg7.m  */
%
%    EXAMPLE 7 FOR STATISTICS 322
%    POSTED 9/17/96
%    EXPLORES Q-Q AND P-P PLOTS

ipar = 1 ;     %  1,2

autop = 0 ;    %  1 TO SAVE PLOTS AS POSTSCRIPT FILES


disp('Running MATLAB script s322eg7.m') ;


if ipar == 1 ;    %  THEN DO PART 1: Q-Q PLOTS

    npgrid = 20 ;
    pdel = 1 / (npgrid + 1) ;
  pgrid = linspace(pdel,1-pdel,npgrid) ;
          %  EQUALLY SPACED GRID IN SENSE OF MEAN OF UNIF. ORDER STATS.

  for ieg = 1:4 ;

    if ieg == 1 ;
      subplot(2,2,1) ;
        qx = norminv(pgrid,0,1) ;
          %  VECTOR OF N(0,1) QUANTILES
        qy = norminv(pgrid,5,2) ;
          %  VECTOR OF N(5,4) QUANTILES
          %  CAREFUL:  THAT LAST ARGUMENT IS S.D. NOT VAR!
      plot(qx,qy,'+',qx,qy,'-') ;
          %  USE TWO VERSIONS TO GET PLUS SIGNS CONNECTED BY LINES
        title('Q-Q Plot of N(0,1) vs. N(5,4)') ;
        xlabel('N(0,1)') ;
        ylabel('N(5,4)') ;
        axis('square') ;
         %  THIS CONTRAINS THE PLOT TO LIE IN A SQUARE REGION,
          %  VS. THE USUAL "FILLING UP AS MUCH SPACE AS POSSIBLE"
        axis('equal') ;
          %  THIS MAKES LENGTHS ON THE AXES BE IN THE SAME UNITS

          %  THE ABOVE TWO COMMANDS ARE USEFUL WHEN "VERTICAL LENGTH"
          %  SHOULD HAVE THE SAME VISUAL CONTENT AS "HORIZONTAL LENTH"
          %  BE CAREFUL, THEY CAN SURPRISE YOU SOMETIMES...

    elseif ieg == 2 ;
      subplot(2,2,2) ;
      qx = norminv(pgrid,0,1) ;
      qy = norminv(pgrid,-1,2/3) ;
      plot(qx,qy,'+',qx,qy,'-') ;
        title('Q-Q Plot of N(0,1) vs. N(-1,4/9)') ;
        xlabel('N(0,1)') ;
        ylabel('N(-1,4/9)') ;
        axis('square') ;
        axis('equal') ;

    elseif ieg == 3 ;
      subplot(2,2,3) ;
        qy = tan(pi * (pgrid - 1/2)) ;
          %  RECALL CAUCHY CDF IS 1/2 + ARCTAN(X)/PI
          q3c = tan(pi * (.75 - 1/2)) ;
          %  THIRD QUARTILE OF CAUCHY
          q3n = norminv(.75,0,1) ;
          %  THIRD QUARTILE OF N(0,1)
        qx = norminv(pgrid,0,q3c/q3n) ;
      plot(qx,qy,'+',qx,qy,'-') ;
        title('Q-Q Plot of N(0,q3c/qcn) vs. Cauchy') ;
        xlabel('N(0,q3c/q3n)') ;
        ylabel('Cauchy') ;
        axis('square') ;
        axis('equal') ;

    elseif ieg == 4 ;
      subplot(2,2,4) ;
      qx = norminv(pgrid,1/2,sqrt(1/12)) ;
      qy = pgrid ;
          %  UNIFORM(0,1) C.D.F. 
      plot(qx,qy,'+',qx,qy,'-') ;
        title('Q-Q Plot of N(1/2,1/12) vs. Uniform') ;
        xlabel('N(1/2,1/12)') ;
        ylabel('Uniform') ;
        axis('square') ;
        axis('equal') ;

    end ;

  end ;

  if autop == 1 ;
    print s322eg7a.ps -dps ;
  end ;
          %  PUTTING THIS INSIDE THIS if BLOCK ALLOWS ONE TO CHOOSE
          %  WHETHER OR NOT TO SAVE TO FILE (USEFUL WHILE DEBUGGING).


elseif ipar == 2 ;    %  THEN DO PART 2: P-P PLOTS

    npgrid = 2000 ;
    pdel = 1 / (npgrid + 1) ;
  pgrid = linspace(pdel,1-pdel,npgrid) ;
          %  NOW USE DENSER GRID FOR BETTER RESOLUTION

  for ieg = 1:4 ;

    if ieg == 1 ;
      subplot(2,2,1) ;
        px = pgrid ;
          %  VECTOR OF N(0,1) PROBABILITIES
        py = normcdf(norminv(pgrid,0,1),5,2) ;
          %  CORRESPONDING VECTOR OF N(5,4) PROBABILITIES
      plot(px,py,'-r') ;
          %  USE RED LINE TO SHOW UP WELL WHEN IT COINCIDES WITH THE AXIS
        title('P-P Plot of N(0,1) vs. N(5,4)') ;
        xlabel('N(0,1)') ;
        ylabel('N(5,4)') ;
        axis('square') ;
        axis([0 1 0 1]) ;
          %  NOW USE THIS TO CONTROL THE EXPLICIT NUMBERS OF THE AXES
          %  AGAIN BE CAREFUL.  IT IS EASY FOR PICS TO COMPLETELY
          %  DISAPPEAR FROM THE PAGE.

    elseif ieg == 2 ;
      subplot(2,2,2) ;
        px = pgrid ;
          %  VECTOR OF N(0,1) PROBABILITIES
        py = normcdf(norminv(pgrid,0,1),-1,2/3) ;
          %  CORRESPONDING VECTOR OF N(-1,4/9) PROBABILITIES
      plot(px,py,'-r') ;
        title('P-P Plot of N(0,1) vs. N(-1,4/9)') ;
        xlabel('N(0,1)') ;
        ylabel('N(-1,4/9)') ;
        axis('square') ;
        axis([0 1 0 1]) ;

    elseif ieg == 3 ;
      subplot(2,2,3) ;

          %  RECALL CAUCHY CDF IS 1/2 + ARCTAN(X)/PI
          q3c = tan(pi * (.75 - 1/2)) ;
          %  THIRD QUARTILE OF CAUCHY
          q3n = norminv(.75,0,1) ;
          %  THIRD QUARTILE OF N(0,1)
        qx = norminv(pgrid,0,q3c/q3n) ;

        py = pgrid ;
          %  VECTOR OF CAUCHY PROBABILITIES
          q3c = tan(pi * (.75 - 1/2)) ;
          q3n = norminv(.75,0,1) ;
        px = normcdf(tan(pi * (pgrid - 1/2)),0,q3c/q3n) ;
          %  CORRESPONDING N(0,q3c/q3n PROBABILITIES
          %  NOTE ALTERNATE FORM OF P-P PLOT IS USED HERE TO 
          %  MAKE RED CURVE COVER AXIS (COULD HAVE DONE SAME IN SUBPLOT 1)
      plot(px,py,'-r',px,py,':k') ;
          %  HERE IS A METHOD TO MAKE THE CURVE THAT COINCIDES WITH THE
          %  AXIS VISIBLE: FIRST PLOT THE CURVE IN RED, THEN PLOT IT AGAIN
          %  IN BLACK, DOTTED.
        title('P-P Plot of N(0,q3c/qcn) vs. Cauchy') ;
        xlabel('N(0,q3c/q3n)') ;
        ylabel('Cauchy') ;
        axis('square') ;
        axis([0 1 0 1]) ;

    elseif ieg == 4 ;
      subplot(2,2,4) ;
        px = pgrid ;
          %  VECTOR OF N(0,1) PROBABILITIES
        py = norminv(pgrid,1/2,sqrt(1/12)) ;
          %  CORRESPONDING VECTOR OF UNIF(0,1) PROBABILITIES
          flag = find(py < 0) ;
        py(flag) = zeros(1,length(flag)) ;
          flag = find(py > 1) ;
        py(flag) = ones(1,length(flag)) ;
          %  RECALL UNIF(0,1) C.D.F. IS x*1{0<x<1}. 
     plot(px,py,'-r') ;
        title('P-P Plot of N(1/2,1/12) vs. Uniform') ;
        xlabel('N(1/2,1/12)') ;
        ylabel('Uniform') ;
        axis('square') ;
        axis([0 1 0 1]) ;

    end ;

  end ;

  if autop == 1 ;
    print s322eg7b.ps -dps ;
  end ;
          %  PUTTING THIS INSIDE THIS if BLOCK ALLOWS ONE TO CHOOSE
          %  WHETHER OR NOT TO SAVE TO FILE (USEFUL WHILE DEBUGGING).


end ;

