octave_packages/statistics-1.1.3/normplot.m

   1 ## Author: Paul Kienzle <pkienzle@users.sf.net>
   2 ## This program is granted to the public domain.
   3
   4 ## -*- texinfo -*-
   5 ## @deftypefn {Function File} normplot (@var{X})
   6 ##
   7 ## Produce a normal probability plot for each column of @var{X}.
   8 ##
   9 ## The line joing the 1st and 3rd quantile is drawn on the
  10 ## graph.  If the underlying distribution is normal, the
  11 ## points will cluster around this line.
  12 ##
  13 ## Note that this function sets the title, xlabel, ylabel,
  14 ## axis, grid, tics and hold properties of the graph.  These
  15 ## need to be cleared before subsequent graphs using 'clf'.
  16 ## @end deftypefn
  17
  18 function normplot(X)
  19   if nargin!=1, print_usage; end
  20   if (rows(X) == 1), X=X(:); end
  21
  22   # plot labels
  23   title "Normal Probability Plot"
  24   ylabel "% Probability"
  25   xlabel "Data"
  26
  27   # plot grid
  28   t = [0.00001;0.0001;0.001;0.01;0.1;0.3;1;2;5;10;25;50;
  29       75;90;95;98;99;99.7;99.9;99.99;99.999;99.9999;99.99999];
  30   tics ('y',normal_inv(t/100),num2str(t));
  31   grid on
  32
  33   # Transform data
  34   n = rows(X);
  35   if n<2, error("normplot requires a vector"); end
  36   q = normal_inv([1:n]'/(n+1));
  37   Y = sort(X);
  38
  39   # Set view range with a bit of space around data
  40   miny = min(Y(:)); minq = min(q(1),normal_inv(0.05));
  41   maxy = max(Y(:)); maxq = max(q(end),normal_inv(0.95));
  42   yspace = (maxy-miny)*0.05; qspace = (q(end)-q(1))*0.05;
  43   axis ([miny-yspace, maxy+yspace, minq-qspace, maxq+qspace]);
  44
  45   # Find the line joining the first to the third quartile for each column
  46   q1 = ceil(n/4);
  47   q3 = n-q1+1;
  48   m = (q(q3)-q(q1))./(Y(q3,:)-Y(q1,:));
  49   p = [ m; q(q1)-m.*Y(q1,:) ];
  50
  51   # Plot the lines one at a time.  Plot the lines before overlaying the
  52   # normals so that the default label is 'line n'.
  53   if columns(Y)==1,
  54     leg = "+;;";
  55   else
  56     leg = "%d+;Column %d;";
  57   endif
  58
  59   for i=1:columns(Y)
  60     plot(Y(:,i),q,sprintf(leg,i,i)); hold on;
  61
  62     # estimate the mean and standard deviation by linear regression
  63     # [v,dv] = wpolyfit(q,Y(:,i),1)
  64   end
  65
  66   # Overlay the estimated normal lines.
  67   for i=1:columns(Y)
  68     # Use the end points and one point guaranteed to be in the view since
  69     # gnuplot skips any lines whose points are all outside the view.
  70     pts = [Y(1,i);Y(q1,i);Y(end,i)];
  71     plot(pts, polyval(p(:,i),pts), [num2str(i),";;"]);
  72   end
  73   hold off;
  74 end
  75