From owner-octave-sources at bevo dot che dot wisc dot edu  Mon Jun 10 08:24:12 1996
Subject: prob.m
From: Eyal Doron <doron at mickey dot mpi-hd dot mpg dot de>
To: octave-sources at bevo dot che dot wisc dot edu (Octave Sources)
Date: Mon, 10 Jun 1996 15:22:20 +0100 (MET DST)

Hi,
   Here is one of my favorites. 

Eyal Doron

========================prob.m===============================

function [Nout,Xout,Wid]=prob(Y,n,St1,St2,St3,St4);
% function [N,X,Wid]=prob(Y,n,mode,color,holdon,plotit) plots a 
% probability density.
%
% Y       - data vector.
% n       - number of bins, default 10.
% mode    - 'lin'   (default) gives a normal probability plot;
%           'log'   (or 'logy') gives a semi-logarithmic plot;
%           'logx'  gives a log-x bin spacing and plot
%           'logxy' gives a log-log bin spacing and plot;
%           'log10',... - plot the log10 rather than using log scale.
% color   - color and line style to plot with, in the usual syntax.
% holdon  - 'hold' causes the plot to be overlaid on top of an existing
%           plot. The hold state is preserved.
% plotit  - 'plot' plot results even if output arguments are specified.
% -------------------------------------------------------------
% N,X,Wid - Height, position and width of the bins. If specified, result
%           is not plotted, only returned, unless plotting is explicitly
%           requested (see above).

Max_Str=4;

mode='lin'; color=''; hold_on=0; Do_Plot_Also=0;
if nargin<1
  help prob
  return
elseif nargin==1 
  n=10; 
else
  n_str=min(nargin-2,Max_Str);
  if isstr(n)
    n_str=min(nargin-1,Max_Str);
    for l=n_str:-1:2
      eval(sprintf('St%g=St%g;',l,l-1));
    end
    St1=n;
    n=10;
  end
  for l=1:n_str
    eval(sprintf('St=St%g;',l));
    if strcmp(St,'hold')
      hold_on=1;
    elseif strcmp(St,'plot')
      Do_Plot_Also=1;
    elseif any(St=='l')
      mode=St;
    else
      color=St;
    end
  end
end


log_X=strcmp(mode,'logx') | strcmp(mode,'logxy') | strcmp(mode,'logx10') |...
      strcmp(mode,'logxy10');
log_Y=strcmp(mode,'log') | strcmp(mode,'logy') | strcmp(mode,'logxy') |...
      strcmp(mode,'logy10') | strcmp(mode,'logxy10');
Use10=any(mode=='1') & any(mode=='0');

Y=Y(:);
if any(imag(Y)~=0)
  Y=abs(Y);
end
if log_X 
  Y=Y(find(Y>0));
end
% if log_X & min(Y)<=0,  error('Can''t use log x axis with <=0 data!'); end

% Do the work

if log_X
  [N,X]=hist(log10(Y),n);
else
  [N,X]=hist(Y,n);
end
delta=X(2)-X(1);
edges=[X(1)-delta/2;X(:)+delta/2];
if log_X, edges=(10).^edges; end
barwidth=edges(2:length(edges))-edges(1:length(edges)-1);
N=N./barwidth.'/length(Y);

%                          Output only
if nargout>0
  Nout=N;
  Xout=X;
  Wid=barwidth.';
  if Use10
%    if log_X, Xout=log10(Xout); end
    if log_Y, Nout=log10(Nout); end
  else
    if log_X, Xout=10.^Xout; end
  end
  if ~Do_Plot_Also, return; end
end

% Prepare plot data

[xb,yb]=bar(X,N);
if log_X & ~Use10, xb=(10).^xb; end
if log_Y
  yzero=find(yb==0);
  if ~isempty(yzero)
    ymin=min(yb(find(yb~=0)));
    yb(yzero)=ymin/2*ones(size(yzero));
  end
  if Use10, yb=log10(yb); end
end

plot_str1='plot(xb,yb';
if ~Use10
  if log_X & log_Y, plot_str1='loglog(xb,yb';
  elseif log_X, plot_str1='semilogx(xb,yb';
  elseif log_Y, plot_str1='semilogy(xb,yb';
  end
end
plot_str2=');';
if ~isempty(color), plot_str2=[',''',color,''');']; end

plot_str=[plot_str1,plot_str2];

OldHold=ishold;
if hold_on, hold on; end
eval(plot_str);
if hold_on & ~OldHold, hold off; end