randraw.m

能够产生多达50种分布的随机变量(Nakagami,Rayleigh,Poisson...).

function varargout = randraw(distribName, distribParams, varargin)
%
%   EFFICIENT RANDOM VARIATES GENERATOR
%
% See alphabetical list of the supported distributions below (over 50 distributions)
%
% 1)  randraw
%           presents general help.
% 2)  randraw( distribName )
%           presents help for the specific distribution defined
%           by usage string distribName (see table below).
% 3)  Y = randraw( distribName, distribParams, sampleSize );
%           returns array Y of size = sampleSize of random variates from distribName
%           distribution with parameters distribParams
%
%               ALPHABETICAL LIST OF THE SUPPORTED DISTRIBUTIONS:
%  ____________________________________________________________________
% |      DISTRIBUTION NAME                    |   USAGE STRING         |
% |___________________________________________|________________________|
% |        Alpha                              |    'alpha'             |
% |        Anglit                             |    'anglit'            |
% |        Antilognormal                      |    'lognorm'           |
% |        Arcsin                             |    'arcsin'            |
% |        Bernoulli                          |    'bern'              |
% |        Bessel                             |    'bessel'            |
% |        Beta                               |    'beta'              |
% |        Binomial                           |    'binom'             |
% |        Bradford                           |    'bradford'          |
% |        Burr                               |    'burr'              |
% |        Cauchy                             |    'cauchy'            |
% |        Chi                                |    'chi'               |
% |        Chi-Square (Non-Central)           |    'chisqnc'           |
% |        Chi-Square (Central)               |    'chisq'             |
% |        Cobb-Douglas                       |    'lognorm'           |
% |        Cosine                             |    'cosine'            |
% |        Double-Exponential                 |    'laplace'           |
% |        Erlang                             |    'erlang'            |
% |        Exponential                        |    'exp'               |
% |        Extreme-Value                      |    'extrval'           |
% |        F (Central)                        |    'f'                 |
% |        F (Non-Central)                    |    'fnc'               |
% |        Fisher-Tippett                     |    'extrval'           |
% |        Fisk                               |    'fisk'              |
% |        Frechet                            |    'frechet'           |
% |        Furry                              |    'furry'             |
% |        Gamma                              |    'gamma'             |
% |        Generalized Inverse Gaussian       |    'gig'               |
% |        Generalized Hyperbolic             |    'gh'                |
% |        Geometric                          |    'geom'              |
% |        Gompertz                           |    'gompertz'          |
% |        Gumbel                             |    'gumbel'            |
% |        Half-Cosine                        |    'hcos'              |
% |        Hyperbolic Secant                  |    'hsec'              |
% |        Hypergeometric                     |    'hypergeom'         |
% |        Inverse Gaussian                   |    'ig'                |
% |        Laplace                            |    'laplace'           |
% |        Logistic                           |    'logistic'          |
% |        Lognormal                          |    'lognorm'           |
% |        Lomax                              |    'lomax'             |
% |        Lorentz                            |    'lorentz'           |
% |        Maxwell                            |    'maxwell'           |
% |        Nakagami                           |    'nakagami'          |
% |        Negative Binomial                  |    'negbinom'          |
% |        Normal                             |    'norm'              |
% |        Normal-Inverse-Gaussian (NIG)      |    'nig'               |
% |        Pareto                             |    'pareto'            |
% |        Pareto2                            |    'pareto2'           |
% |        Pascal                             |    'pascal'            |
% |        Planck                             |    'planck'            |
% |        Poisson                            |    'po'                |
% |        Quadratic                          |    'quadr'             |
% |        Rademacher                         |    'rademacher'        |
% |        Rayleigh                           |    'rayl'              |
% |        Rice                               |    'rice'              |
% |        Semicircle                         |    'semicirc'          |
% |        Skellam                            |    'skellam'           |
% |        Student's-t                        |    't'                 |
% |        Triangular                         |    'tri'               |
% |        Truncated Normal                   |    'normaltrunc'       |
% |        Tukey-Lambda                       |    'tukeylambda'       |
% |        U-shape                            |    'u'                 |
% |        Uniform (continuous)               |    'uniform'           |
% |        Von Mises                          |    'vonmises'          |
% |        Wald                               |    'wald'              |
% |        Weibull                            |    'weibull'           |
% |        Wigner Semicircle                  |    'wigner'            |
% |        Yule                               |    'yule'              |
% |        Zeta                               |    'zeta'              |
% |        Zipf                               |    'zipf'              |
% |___________________________________________|________________________|

%  Version 2.0 - August 2007
%         1) New distributions support: Nakagami and Rician !
%         2) Small typo corrections in comments
%  Version 1.8 - February 2007
%         GIG distribution (thanks to Mr. Demetris Lamnisos)
%           Computational exceptions in the reparameterized GIG generation were fixed
%  Version 1.7 - December 2006
%         GIG distribution (thanks to Dr. Junbin Gao)
%           Computational exceptions in the reparameterized GIG generation were fixed
%  Version 1.6 - September 2006
%         Exception handling: BINOMIAL distribution - special case for n*p~=0
%         Geometric distibution: additional note in help section
%  Version 1.5 - December 2005
%        'true' and 'false' functions were replased by ones and zeros to support Matlab releases
%         below 6.5
%  Version 1.4 - September 2005 -
%      Bugs fix:
%        1) GAMMA distribution (thanks to Earl Lawrence):
%             special case for a<1
%        2) GIG distribution (thanks to Panagiotis Braimakis):
%            typo in help
%            code adjustment to overcome possible computational overflows
%        3) CHI SQUARE distribution
%            typo in help
%  Version 1.3 - July 2005 -
%      Bug fix:
%         Typo in GIG distribution generation:
%         should be 'out' instead of 'x' in lines 1852 and 1858
%  Version 1.2 - May 2005  -
%      Bugs fix:
%        1) Poisson distribution did not work for lambda < 21.4. Typo ( ti instead of t )
%        2) GIG distribution:  support to chi=0 or psi=0 cases
%        3) Beta distribution: column sampleSize
%        4) Cauchy distribution: typo in example
%        5) Chi distribution:   typo in example
%        6) Non-central F distribution:  number of input parameters
%        7) INVERSE GAUSSIAN (IG) distribution: typo in example
%
%  Version 1.1 - April 2005 -  Bug fix:   Generation from binomial distribution using only 'binomial'
%                                   usage string was changed to 'binom' ( 'binomial' works too ).
%  Version 1.0 - March 2005 -  Initial version
%  Alex Bar Guy  &  Alexander Podgaetsky
%    [e-mail]

% These programs are distributed in the hope that they will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

% Any comments and suggestions please send to:
%    [e-mail]

% Reference links:
%   1) http://mathworld.wolfram.com/topics/StatisticalDistributions.html
%   2) http://en.wikipedia.org/wiki/Category:Probability_distributions
%   3) http://www.brighton-webs.co.uk/index.asp
%   4) http://www.jstatsoft.org/v11/i03/v11i03.pdf
%   5) http://www.quantlet.com/mdstat/scripts/csa/html/node236.html

funcName = mfilename;

if nargin == 0
     help(funcName);
     return;
elseif nargin == 1
     runMode = 'distribHelp';
elseif nargin == 2
     runMode = 'genRun';
     sampleSize = [1 1];
else
     runMode = 'genRun';
     sampleSize = [varargin{1:end}];
end

distribNameInner = lower( distribName( ~isspace( distribName ) ) );

if strcmp(runMode, 'distribHelp')
     fid = fopen( [ funcName '.m' ], 'r' );
     printHelpFlag = 0;
     while 1
          tline = fgetl( fid );
          if ~ischar( tline )
               fprintf( '\n Unknown distribution name ''%s''.\n', distribName );
               break;
          end
          if ~isempty( strfind( tline, [ 'END ', distribNameInner,' HELP' ] ) )
               printHelpFlag = 0;
               break;
          end
          if printHelpFlag
               startPosition = strfind( tline, ' % ' ) + 3;
               printLine = tline( startPosition : end );
               if ~strcmp( funcName, 'randraw' )
                    indxs = strfind( printLine, 'randraw' );
                    while ~isempty( indxs )
                         headLine = printLine( 1:indxs(1)-1 );
                         tailLine = printLine( indxs(1)+7:end );
                         printLine = [ headLine, funcName, tailLine ];
                         indxs = strfind( printLine, 'randraw' );
                    end
               end
               pause(0.02);
               fprintf( '\n%s', printLine );
          end
          if ~isempty( strfind( tline, [ 'START ', distribNameInner,' HELP' ] ) )
               printHelpFlag = 1;
          end
     end
     fprintf( '\n\n' );
     fclose( fid );
     if nargout > 0
          varargout{1} = [];
     end
     return;
end

if length(sampleSize) == 1
     sampleSize = [ sampleSize, 1 ];
end

if strcmp(runMode, 'genRun')
     runExample = 0;
     plotFlag = 0;

     dbclear if warning;
     out = [];
     if prod(sampleSize) > 0
          switch lower( distribNameInner )
               case {'alpha'}
                    % START alpha HELP
                    % THE ALPHA DISTRIBUTION
                    %
                    % pdf(y) = b*normpdf(a-b./y) ./ (y.^2*normcdf(a)); y>0; a>0; b>0;
                    % cdf(y) = normcdf(a-b./y)/normcdf(a); y>0; a>0; b>0;
                    %   where normpdf(x) = 1/sqrt(2*pi) * exp(-1/2*x.^2); is the standard normal PDF
                    %         normcdf(x) = 0.5*(1+erf(y/sqrt(2))); is the standard normal CDF
                    %
                    % PARAMETERS:
                    %   a - shape parameter (a>0)
                    %   b - shape parameter (b>0)
                    %
                    % SUPPORT:
                    %   y,  y>0
                    %
                    % CLASS:
                    %   Continuous skewed distributions
                    %
                    % USAGE:
                    %   randraw('alpha', [], sampleSize) - generate sampleSize number
                    %         of variates from Alpha distribution with shape parameters a and b;
                    %   randraw('alpha') - help for Alpha distribution;
                    %
                    % EXAMPLES:
                    %  1.   y = randraw('alpha', [1 2], [1 1e5]);
                    %  2.   y = randraw('alpha', [2 3], 1, 1e5);
                    %  3.   y = randraw('alpha', [10 50], 1e5 );
                    %  4.   y = randraw('alpha', [20.5 30.5], [1e5 1] );
                    %  5.   randraw('alpha');
                    % END alpha HELP

                    % References:
                    % 1. doc erf

                    checkParamsNum(funcName, 'Alpha', 'alpha', distribParams, [2]);
                    a  = distribParams(1);
                    b  = distribParams(2);
                    validateParam(funcName, 'Alpha', 'alpha', '[a, b]', 'a', a, {'> 0'});
                    validateParam(funcName, 'Alpha', 'alpha', '[a, b]', 'b', b, {'> 0'});

                    out = b ./ ( a - norminv(normcdf(a)*rand(sampleSize)) );

               case {'anglit'}
                    % START anglit HELP
                    % THE ANGLIT DISTRIBUTION
                    %
                    % Standard form of anglit distribution:
                    %   pdf(y) = sin(2*y+pi/2);  -pi/4 <= y <= pi/4;
                    %   cdf(y) = sin(y+pi/4).^2; -pi/4 <= y <= pi/4;
                    %
                    %   Mean = Median = Mode = 0;
                    %   Variance = (pi/4)^2 - 0.5;
                    %
                    % General form of anglit distribution:
                    %   pdf(y) = sin(pi/2*(y-t)/s+pi/2);  t-s <= y <= t+s; s>0
                    %   cdf(y) = sin(pi/4*(y-t)/s+pi/4).^2;  t-s <= y <= t+s; s>0
                    %
                    %   Mean = Median = Mode = t;
                    %   Variance = ???????;
                    %
                    % PARAMETERS:
                    %   t - location
                    %   s -scale; s>0
                    %
                    % SUPPORT:
                    %   y,   -pi/4 <= y <= pi.4   - standard Anglit distribution
                    %    or
                    %   y,   t-s <= y <= t+s  - generalized Anglit distribution
                    %
                    % CLASS:
                    %   Continuous distributions
                    %
                    % USAGE:
                    %   randraw('anglit', [], sampleSize) - generate sampleSize number
                    %         of variates from standard Anglit distribution;
                    %   randraw('anglit', [t, s], sampleSize) - generate sampleSize number
                    %         of variates from generalized Anglit distribution
                    %         with location parameter 't' and scale parameter 's';
                    %   randraw('anglit') - help for Anglit distribution;
                    %
                    % EXAMPLES: