randist.m

很多matlab的源代码

function z = randist(x,ty,a)
% RANDIST Random numbers with distributions
%
%         Z=RANDIST(X,TY,M) computes random distributions of type TY and mean M.
%         Z has the same size as X
%         If X is a matrix, each column of X has the appropriate distribution
%         If X is a vector or array, X has the given distribution
%         If X is a scalar with X=N, X is set to a COLUMN vector of length N
%         TY can be: 'nor'mal, 'lap'lace, 'uni'form, 'exp'onential, 'gam'ma,
%	  'ray'leigh, or 'poi'sson. You may use just the first three letters. 
% 
%	  DEFAULTS: TY='nor'     M: M=0(Nor & Lap), M=0.5(uni), M=1(others)
%
%         RANDIST (with no arguments) invokes the following example:
%
%         % Generate 1000 exponentially distributed values with mean=3, 
%         % Plot their histogram and compare with expected distribution 
%         >>ze=randist(1000,'exp',3);meanz=mean(ze) 	%display mean
%         >>[x0,l] = hist(ze,20);			%find binwidth
%         >>bw = l(2)-l(1);da=1000*(bw/3)*exp(-l/3);	%expected
%         >>hist(ze,20),hold on,plot(l,da),hold off	%plot


% ADSP Toolbox: Version 2.0 
% For use with "Analog and Digital Signal Processing", 2nd Ed.
% Published by PWS Publishing Co.
%
% Ashok Ambardar, EE Dept. MTU, Houghton, MI 49931, USA
% http://www.ee.mtu/faculty/akambard.html
% e-mail: [e-mail]
% Copyright (c) 1998


%REF: Numerical recipes, Press et al, Cambridge, 1986
%NOTE: ADD variance for laplace, normal and uniform?? by finding zero mean 
%random variable z and then scaling to sqrt(var)*z + mean

if nargin==0,help randist,disp('Strike a key to see results of the example')
pause,ze=randist(1000,'exp',3);meanz=mean(ze),[x0,l]=hist(ze,20);bw=l(2)-l(1);
da=1000*(bw/3)*exp(-l/3);hist(ze,20),hold on,
plot(l,da,'--',l,da,'o'),hold off,return,end

if exist('version')==5,v=4;else,v=3;end,if prod(size(x))==1,x=1:x;x=x(:);end
if nargin<2,ty='normal';end,ty=ty(1:3);[n1,n2]=size(x);

%UNIFORM
if ty=='uni',if nargin<3,a=.5;end,if v==3,rand('uniform');end
a=a-.5;z=a+rand(n1,n2);return,end 

%NORMAL
if ty=='nor',if nargin<3,a=0;end,
if v==3,rand('normal');z=a+rand(n1,n2);else,z=a+randn(n1,n2);end
return,end 

if v==3,rand('uniform');end
%LAPLACE
if ty=='lap',if nargin<3,a=0;end
%z=zeros(n1,n2);for k=1:n2,l2=fix(n1/2);z1=rand(l2,1);z1=log(z1);
%z2=rand(n1-l2,1);z2=-log(z2);z(:,k)=a+[z1;z2];end,return,end
%%%NEW
z=zeros(n1,n2);for k=1:n2,l2=fix(n1/2);z1=rand(l2,1);z1=log(2*z1);
z2=rand(n1-l2,1);z2=-log(2*z2);z(:,k)=a+[z1;z2];end,return,end

if nargin<3,a=1;end
%EXPONENTIAL
if ty=='exp',z=rand(n1,n2);z=-a*log(z);return,end 

%RAYLEIGH
if ty=='ray',a=a*sqrt(2/pi);
%z=rand(n1,n2);z=sqrt(2*a*a*log(1 ./z));return,end  %%%CHECK THIS
z=rand(n1,n2);z=sqrt(2*a*a*log(1 ./(1-z)));return,end 

z=zeros(n1,n2);
%GAMMA
if ty=='gam',if a<6,z=1;for j=1:a,z=z.*rand(n1,n2);end,z=-log(z);return,end
for k=1:n1*n2,v1=2;v2=2;e=-1;
while (v1*v1+v2*v2) > 1 | z(k) <= 0 | rand > e
v1=2*rand-1;v2=2*rand-1;y=v2/v1;m=a-1;s=sqrt(2*m+1);z(k)=s*y+m;
e=(1+y*y)*exp(m*log(z(k)/m)-s*y);end,end,return,end

%POISSON
if ty=='poi',if a < 12,g=exp(-a);
for k=1:n1*n2,em=-1;t=1;while t > g,em=em+1;t=t*rand;end,z(k)=em;end
else,sq=sqrt(2*a);alxm=log(a);for k=1:n1*n2,em=-1;t=-1;g=a*alxm-gmln(a+1);
while em < 0 | rand > t,y=tan(pi*rand);em=sq*y+a;e1=fix(em);
t=0.9*(1+y*y)*exp(e1*alxm-gmln(e1+1)-g);end,z(k)=e1;end,end,return,end

error('Use exp,gam(ma),lap(lace),nor,poi(sson),ray(leigh) or uni')