🚀 go-pugleaf

Hoping to start a trend, here is an unsolicited matlab m-file.
RANDG generates samples from an arbitrary probability distribution
function, specified symbolically.  Suggested improvements are welcome.

#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create the files:
#	randg.m
# This archive created: Fri Jan 22 16:56:20 1993
export PATH; PATH=/bin:$PATH
echo shar: extracting "'randg.m'" '(3479 characters)'
if test -f 'randg.m'
then
	echo shar: will not over-write existing file "'randg.m'"
else
sed 's/^XX//' << \SHAR_EOF > 'randg.m'
XXfunction x=randg(PDF_,BND_,ICBND_,N_,opts_)
XX
XX%RANDG
XX%	Generates samples from arbitrary probability distributions.
XX%	X=randg(ICDF,N) uses the inverse cumulative distribution
XX%	function ICDF with argument 'y' and result 'x'.  For example,
XX%	ICDF='x=-log(1-y)' is the inverse CDF for 'y=exp(-x)' (x>0),
XX%	which has CDF 'y=1-exp(-x)'.  N is the number of samples.
XX%
XX%	X=randg(PDF,BND,ICBND,N) works when the CDF is not easily inverted.
XX%	PDF is a string like 'yp=exp(-x/2)' that describes the
XX%	probability density function 'yp' as a function of 'x'.  BND is
XX%	another string that defines a function that is everywhere
XX%	not less than PDF, for example 'yb=exp(-x/4)' for x>0.  ICBND 
XX%	is the inverse of the CDF of BND, in this case the inverse of
XX%	'y=1-exp(-x/4)' which is 'x=-4*log(1-y)', defined for 0<=y<=1.
XX%	Note that the names of the variables must be the same as shown
XX%	here (x, y, yp, and yb).  Intermediate variables may have any name
XX%	that does not end with '_' and which is not a built-in M-file
XX%	variable (such as 'nargin').
XX%
XX%	X=randg(PDF,BND,ICBND,N,'plot') will plot the algorithm's progress.
XX%
XX%	Here is a complicated example of the second method, taken from
XX%	Leon-Garcia "Probability and Random Processes for Electrical
XX%	Engineering", p. 174.  It generates samples X from a gamma pdf with
XX%	parameter alpha, 0<alpha<1.
XX%
XX%	PDF = 'alpha=0.33;yp=x.^(alpha-1).*exp(-x)/gamma(alpha);'
XX%
XX%	BND = [	'alpha=0.33;                          ';
XX%		'Ip=find(x>1); In=find(x<=1);         ';
XX%		'yb=x;                                ';
XX%		'yb(In)=x(In).^(alpha-1)/gamma(alpha);';
XX%		'yb(Ip)=exp(-x(Ip))/gamma(alpha);     ']
XX%	ICBND= ['alpha=0.33;                             '
XX%		'e=exp(1);                               '
XX%		'In=find(y<(e/(alpha+e)));               '
XX%		'Ip=find(y>=(e/(alpha+e)));              '
XX%		'x=y;                                    '
XX%		'x(In)=((alpha+e)*y(In)/e).^(1/alpha);   '
XX%		'x(Ip)=-log((alpha+e)*(1-y(Ip))/alpha/e);']
XX
XX%	Michael Maurer, 5 May 1991
XX%	Copyright (c) 1993 by Michael Maurer
XX
XXdist_=rand('dist');
XXrand('uniform');
XXif nargin==2,
XX   ICDF_=PDF_;
XX   N_=BND_;
XX   y=rand(N_,1);
XX   for l_=1:height(ICDF_),
XX      eval(ICDF_(l_,:));				% x=ICDF_(y)
XX   end
XXelse
XX   xf_=zeros(N_,1);
XX   fnd_=0;
XX   eff_=0.99;				% initial guess to efficiency
XX   while fnd_<N_,
XX      Ntry_=fix((N_-fnd_)/eff_+10);		% slightly overestimate misses
XX      y=rand(Ntry_,1);
XX      for l_=1:height(ICBND_),
XX	 eval(ICBND_(l_,:));			% x=ICBND_(y)
XX      end
XX      for l_=1:height(BND_),
XX	 eval(BND_(l_,:));			% yb=BND_(x)
XX      end
XX      for l_=1:height(PDF_),
XX	 eval(PDF_(l_,:));			% yp=PDF_(x)
XX      end
XX      t_=rand(Ntry_,1);
XX      Ifnd_=find(t_.*yb <= yp);
XX      Nfnd_=length(Ifnd_);
XX      Nfnd_=min([N_-fnd_ Nfnd_]);
XX      xf_(fnd_+1:fnd_+Nfnd_)=x(Ifnd_(1:Nfnd_));
XX      fnd_=fnd_+Nfnd_;
XX      eff_=Nfnd_/Ntry_;
XX      
XX      if nargin>4,				% nifty graphical debugging
XX	 disp2('Begin loop with ', N_-(fnd_-Nfnd_), ' samples needed.');
XX	 disp2('Efficiency this iteration was ', eff_);
XX	 if (fnd_-Nfnd_+length(Ifnd_)>N_),
XX	    disp2('Computed ', fnd_-Nfnd_+length(Ifnd_)-N_, ' extra samples'); 
XX	 end;
XX	 hold off;
XX	 axis([min(x) max(x) 0 min(4,max(yp))]);
XX	 [xx_,ii_]=sort(x);
XX	 plot(xx_,yb(ii_));
XX	 hold on;
XX	 plot(xx_,yp(ii_));
XX	 tt_=t_.*yb;
XX	 plot(x(Ifnd_),tt_(Ifnd_),'.')
XX	 Infnd_=find(tt_>yp);
XX	 plot(x(Infnd_),tt_(Infnd_),'o')
XX	 hold off;
XX	 if fnd_<N_,
XX	    disp('Hit a key for next iteration')
XX	    pause
XX	 end
XX      end
XX   end
XX   x=xf_;
XXend
XXrand(dist_)
SHAR_EOF
if test 3479 -ne "`wc -c < 'randg.m'`"
then
	echo shar: error transmitting "'randg.m'" '(should have been 3479 characters)'
fi
fi # end of overwriting check
#	End of shell archive
exit 0
--
______________________________________________________________________
Michael Maurer          maurer@nova.stanford.edu        (415) 723-1024