function nppdemo( option )
% SHOWNORM  Shows the histograms and normal prob plots for actual normal rvs.
% USE:  Type  nppdemo  at the MATLAB prompt

% Peter Dunn
% 17 October 2000

if nargin==0,

   % Create figure
   hOverFit = figure('Name','Normal Random Vars Demonstration', ...
                     'tag','tagShowNorm');

   % Set up canvas
   hNPPlotist = axes( 'Position', [0.1 0.08 0.32 0.85], 'tag','tagNPPlot' );
   set(gca, 'FontSize',12);
   title('Normal Probability Plot');
   hHist = axes( 'Position', [0.43 0.08 0.32 0.85], 'tag','tagHist' );
   set(gca, 'FontSize',12);
   title('Histogram');
   set( gca,'YTickLabel',[]);

   % Control Buttons
   uicontrol('Style','Frame','Units','normalized',...
             'Position',[0.78 0.02 0.20 0.20] );

   uicontrol('Style','PushButton','Units','normalized',...
             'Position',[0.80 0.04 0.16 0.07], 'String','Quit',...
             'Callback','delete(gcf);');
   uicontrol('Style','PushButton','Units','normalized',...
             'Position',[0.80 0.13 0.16 0.07], 'String','Plot',...
             'tag','tagPauseBtn','Callback','nppdemo(1);');

   % Parameter Buttons

   uicontrol('Style','Frame','Units','normalized',...
             'Position',[0.78 0.24 0.20 0.38] );
              
   uicontrol('Style','Text','Units','normalized',...
             'Position',[0.80 0.53 0.16 0.07], 'String','Mean:');
   uicontrol('Style','Edit','Units','normalized',...
             'Position',[0.80 0.50 0.16 0.06], 'String','1',...
             'tag','tagMean','Callback','nppdemo(20);');
   uicontrol('Style','Text','Units','normalized',...
             'Position',[0.80 0.41 0.16 0.07], 'String','Variance:');
   uicontrol('Style','Edit','Units','normalized',...
             'Position',[0.80 0.38 0.16 0.06], 'String','1',...
             'tag','tagVar','Callback','nppdemo(21);');
   uicontrol('Style','Text','Units','normalized',...
             'Position',[0.80 0.29 0.16 0.07], 'String','Number Pts:');
   uicontrol('Style','Edit','Units','normalized',...
             'Position',[0.80 0.26 0.16 0.06], 'String','100',...
             'tag','tagNumPts','Callback','nppdemo(22);');

   % Plot Type Buttons
   uicontrol('Style','Frame','Units','normalized',...
             'Position',[0.78 0.64 0.20 0.34] );

   uicontrol('Style','Radio','Units','normalized',...
             'Position',[0.80 0.90 0.16 0.06], 'String','Normal', ...
             'Value',1,...
             'tag','tagShowNormal','Callback','nppdemo(10);');
   uicontrol('Style','Radio','Units','normalized',...
             'Position',[0.80 0.84 0.16 0.06], 'String','Right Skew',...
             'tag','tagRightSkew','Callback','nppdemo(11);');
   uicontrol('Style','Radio','Units','normalized',...
             'Position',[0.80 0.78 0.16 0.06], 'String','Left Skew',...
             'tag','tagLeftSkew','Callback','nppdemo(12);');
   uicontrol('Style','Radio','Units','normalized',...
             'Position',[0.80 0.72 0.16 0.06], 'String','Long Tail',...
             'tag','tagLongTail','Callback','nppdemo(13);');
   uicontrol('Style','Radio','Units','normalized',...
             'Position',[0.80 0.66 0.16 0.06], 'String','Short Tail',...
             'tag','tagShortTail','Callback','nppdemo(14);');
   
elseif nargin==1,

   if option==1, % Plot Data

      % Determine parameter values
      mn = str2num( get( findobj( 'tag','tagMean' ), 'String') );
      vr = str2num( get( findobj( 'tag','tagVar' ), 'String') );
      np = str2num( get( findobj( 'tag','tagNumPts' ), 'String') );

      % Determine what sort of data to plot
      pt1 = 0;
      pt2 = 0;
      pt3 = 0;
      pt4 = 0;
      pt5 = 0;
      pt1 = get( findobj('tag','tagShowNormal'), 'Value');
      pt2 = get( findobj('tag','tagRightSkew'), 'Value');
      pt3 = get( findobj('tag','tagLeftSkew'), 'Value');
      pt4 = get( findobj('tag','tagLongTail'), 'Value');
      pt5 = get( findobj('tag','tagShortTail'), 'Value');

      % Generate Data
      if pt1 == 1, % Normal Data
         r = ( randn( np, 1 ) * vr ) + mn;
      end;
      if pt2 == 1, % Right Skew
         % Could use gamma, but can generate 
         % exp easily without Stats Toolbox
         r = rand( np, 1);
         r = -mn * log( r );
      end;
      if pt3 == 1, % Left Skew
         % Could use beta or (K-gamma), but can generate 
         % (k-exp) easily without Stats Toolbox
         r = rand( np, 1);
         r = -mn * log( r );
         r = max(r) - r;
      end;
      if pt4 == 1, % Long Tail
         % Could use t, but easier to generate from the uniform
         % random nunmbers if we put df=2:

         r = rand( np, 1);
         r = sqrt(2) * (r - 0.5 ) ./ ( sqrt( r - r.^2 ) );
         r = r + mn;
      end;
      if pt5 == 1, % Short Tail
         r = rand( [np, 1] );
         a = mn - sqrt( 3*vr );
         b = mn + sqrt( 3*vr );
         r = ( b - a ) * r + a;

      end;

      % Histogram
      axes( findobj('tag','tagHist') );
      [n, x] = hist( r );
      barh( x, n, 1, 'g');
      set( gca, 'tag','tagHist', 'View',[180 270]);
      xlabel('Frequency');
      title('Histogram');
      ax1 = axis;
      set( gca,'YTickLabel',[]);

      % Normal Prob plot
      axes( findobj('tag','tagNPPlot') );
      limits = nprplot( r );
      rlim = max(limits) - min(limits);
      limits(1) = limits(1) - 0.1*rlim;
      limits(2) = limits(2) + 0.1*rlim;
      set( gca, 'tag','tagNPPlot');
      ax2 = axis;
      axis([ limits ax1(3:4) ] );

   elseif option==10, % Normal selected

      set( findobj('tag','tagRightSkew'),'Value',0);
      set( findobj('tag','tagLeftSkew'),'Value',0);
      set( findobj('tag','tagLongTail'),'Value',0);
      set( findobj('tag','tagShortTail'),'Value',0);

      set( findobj('tag','tagVar'), 'Enable', 'on');

   elseif option==11, % Right Skew selected

      set( findobj('tag','tagShowNormal'),'Value',0);
      set( findobj('tag','tagLeftSkew'),'Value',0);
      set( findobj('tag','tagLongTail'),'Value',0);
      set( findobj('tag','tagShortTail'),'Value',0);

      set( findobj('tag','tagVar'), 'Enable', 'off');

   elseif option==12, % Left Skew selected

      set( findobj('tag','tagShowNormal'),'Value',0);
      set( findobj('tag','tagRightSkew'),'Value',0);
      set( findobj('tag','tagLongTail'),'Value',0);
      set( findobj('tag','tagShortTail'),'Value',0);

      set( findobj('tag','tagVar'), 'Enable', 'off');

   elseif option==13, % Long Tail selected

      set( findobj('tag','tagShowNormal'),'Value',0);
      set( findobj('tag','tagRightSkew'),'Value',0);
      set( findobj('tag','tagLeftSkew'),'Value',0);
      set( findobj('tag','tagShortTail'),'Value',0);

      set( findobj('tag','tagVar'), 'Enable', 'off');

   elseif option==14, % Short Tail selected

      set( findobj('tag','tagShowNormal'),'Value',0);
      set( findobj('tag','tagRightSkew'),'Value',0);
      set( findobj('tag','tagLeftSkew'),'Value',0);
      set( findobj('tag','tagLongTail'),'Value',0);

      set( findobj('tag','tagVar'), 'Enable', 'on');

   elseif option == 20, % Mean entered

      mn = get( findobj( 'tag','tagMean'), 'String');
      mn = str2num(mn);

      if isempty(mn), % eg, text enetered
         mn = 0;
      end;

      mn = real( mn );

      if ( abs(mn) > 100 ),
         mn = sign(mn)*100;
      end;

      set( findobj('tag','tagMean'), 'String', num2str(mn) );

   elseif option == 21, % Variance entered

      vr = get( findobj( 'tag','tagVar'), 'String');
      vr = str2num(vr);

      if isempty(vr), % eg, text enetered
         vr = 1; 
      end;

      vr = real( vr );

      if ( vr<= 0 ),
         vr = 1;
      end;

      if ( vr > 100 ),
         vr = 100;
      end;

      set( findobj('tag','tagVar'), 'String', num2str(vr) );

   elseif option == 22, % Num Pts entered

      np = get( findobj( 'tag','tagNumPts'), 'String');
      np = str2num(np);

      if isempty(np), % eg, text enetered
         np = 10;
      end;

      np = real( round(np) );

      if ( np <= 5 ),
         np = 5;
      end;

      if ( np > 1000 ),
         np = 1000;
      end;

      set( findobj('tag','tagNumPts'), 'String', num2str(np) );

   end;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  SUBFUNCTIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function limits = nprplot( data )
% Plots the normal probability plot of  data

sy = sort( data );
nd = inorm( (1:length(sy))/(length(sy)+1) );

% Plot points
plot(nd,sy,'b+');
xlabel('Standard Normal');
ylabel('Data');
title('Normal Probability Plot');
 
% Add reference line
zx = [-3.5 0 3.5]; %fixed x-points
hold on;
cfplot = plot(zx, mean(sy) + zx*std(sy), 'r--');

limits = [ min(nd) max(nd) ];
hold off;

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5

function x = inorm(p,mu,sigma)
%INVNORM  Inverse cumulative normal ditsribution
%USE:  x=invnorm(p,mu,sigma)
%  where x  is the inverse of the cumulative normal distribution;
%        p  is the probability;
%        mu  is the mean (default is zero);
%        sigma  is the standard deviation (default is 1);
%

%Copyright 1996 Peter Dunn
% 11/11/1996

if nargin<3,
   sigma=1;
end;
if nargin==1,
   mu=0;
end;
 
%error checks
if ~(sigma>0),
   error('sigma  must be positive!');
end;
if size(p)~=size(mu) | size(p)~=size(sigma),
   error('Inputs arguments must be the same size!');
end;
%end error checks
 
x = mu + sqrt(2)*sigma.*erfinv(2*p-1);

return;