randraw.m

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

                    %  1.   y = randraw('anglit', [], [1 1e5]);
                    %  2.   y = randraw('anglit', [], 1, 1e5);
                    %  3.   y = randraw('anglit', [], 1e5 );
                    %  4.   y = randraw('anglit', [10 3], [1e5 1] );
                    %  5.   randraw('anglit');
                    %
                    % END anglit HELP

                    checkParamsNum(funcName, 'Anglit', 'anglit', distribParams, [0, 2]);
                    if numel(distribParams)==2
                         t  = distribParams(1);
                         s  = distribParams(2);
                         validateParam(funcName, 'Anglit', 'anglit', '[t, s]', 's', s, {'> 0'});
                    else
                         t = 0;
                         s = pi/4;
                    end

                    out = t + s * (4/pi*asin(sqrt(rand(sampleSize)))-1);

               case {'arcsin'}
                    % START arcsin HELP
                    % THE ARC-SINE DISTRIBUTION
                    %
                    % pdf(y) = 1 ./ (pi*sqrt(y.*(1-y))); 0<y<1;
                    % cdf(y) = 2*asin(sqrt(y))/pi; 0<y<1;
                    %
                    % Mean = 0.5;
                    % Variance = 0.125;
                    %
                    % PARAMETERS:
                    %  None
                    %
                    % SUPPORT:
                    %  y,    0<y<1
                    %
                    % CLASS:
                    %   Continuous symmetric distributions
                    % NOTES:
                    %  The arc-sine distribution is a special case of the beta distribution
                    %  with both parameters equal to 1/2. The generalized arc-sine distribution
                    %  is the special case of the beta distribution where the two parameters sum
                    %  to 1 but are not necessarily equal to 1/2.
                    %
                    % USAGE:
                    %   randraw('arcsin', [], sampleSize) - generate sampleSize number
                    %         of variates from the Arc-sine distribution;
                    %   randraw('arcsin') - help for the Arc-sine distribution;
                    %
                    % EXAMPLES:
                    %  1.   y = randraw('arcsin', [], [1 1e5]);
                    %  2.   y = randraw('arcsin', [], 1, 1e5);
                    %  3.   y = randraw('arcsin', [], 1e5 );
                    %  4.   randraw('arcsin');
                    %  SEE ALSO:
                    %    U distribution
                    % END arcsin HELP

                    checkParamsNum(funcName, 'Arcsin', 'arcsin', distribParams, 0);
                    out = sin( rand(sampleSize)*pi/2 ).^2;

               case {'bernoulli', 'bern'}
                    % START bernoulli HELP START bern HELP
                    % THE BERNOULLI DISTRIBUTION
                    %
                    % pdf(y) = p.^y .* (1-p).^(1-y);
                    % cdf(y) = (y==0)*(1-p) + (y==1)*1;
                    %
                    % PARAMETERS:
                    %    p is a probability of success; (0<p<1)
                    %
                    % SUPPORT:
                    %     y = [0 1];
                    %
                    % CLASS:
                    %   Discrete distributions
                    %
                    % USAGE:
                    %   randraw('bern', p, sampleSize) - generate sampleSize number
                    %         of variates from the Bernoulli distribution with probability of success p
                    %   randraw('bern') - help for the Bernoulli distribution;
                    %
                    % EXAMPLES:
                    %  1.   y = randraw('bern', 0.5, [1 1e5]);
                    %  2.   y = randraw('bern', 0.1, 1, 1e5);
                    %  3.   y = randraw('bern', 0.9, 1e5 );
                    %  4.   randraw('bern');
                    % END bernoulli HELP END bern HELP

                    checkParamsNum(funcName, 'Bernoulli', 'bernoulli', distribParams, 1);
                    validateParam(funcName, 'Bernoulli', 'bernoulli', 'p', 'p', distribParams(1), {'>=0','<=1'});

                    out = double( rand(  sampleSize  ) < distribParams );

               case {'beta', 'powerfunction', 'powerfunc'}
                    % START beta HELP
                    %  THE BETA DISTRIBUTION
                    %
                    %  ( sometimes: Power Function distribution )
                    %
                    % Standard form of the Beta distribution:
                    %  pdf(y) = y.^(a-1).*(1-y).^(b-1) / beta(a, b);
                    %  cdf(y) = betainc(y,a,b), if (y>=0 & y<=1); 0, if x<0; 1, if x>1
                    %
                    %  Mean = a/(a+b);
                    %  Variance = (a*b)/((a+b)^2*(a+b+1));
                    %
                    % General form of the Beta distribution:
                    %  pdf(y) = (y-m).^(a-1).*(n-y).^(b-1) / (beta(a, b)*(n-m)^(a+b-1));
                    %  cdf(y) = betainc((y-m)/(n-m),a,b), if (y>=m & y<=n); 0, if x<m; 1, if x>n
                    %
                    %  Mean = (n*a + m*b)/(a+b);
                    %  Variance = (a*b)*(n-m)^2/((a+b)^2*(a+b+1));
                    %
                    % PARAMETERS:
                    %   a>0 - shape parameter
                    %   b>0 - shape parameter
                    %   m - location
                    %   n -scale (upper bound); n>=m
                    %
                    % SUPPORT:
                    %   y,   0<=y<=1 - standard beta distribution
                    %    or
                    %   y,   m<=y<=n - generalized beta distribution
                    %
                    % CLASS:
                    %   Continuous skewed distributions
                    %
                    % USAGE:
                    %   randraw('beta', [a, b], sampleSize) - generate sampleSize number
                    %         of variates from standard beta distribution with shape parameters
                    %         'a' and 'b'
                    %   randraw('beta', [m, n, a, b], sampleSize) - generate sampleSize number
                    %         of variates from generalized beta distribution on the interval [m, n]
                    %         with shape parameters 'a' and 'b';
                    %   randraw('beta') - help for the Beta distribution;
                    % EXAMPLES:
                    %  1.   y = randraw('beta', [0.2 0.9], [1 1e5]);
                    %  2.   y = randraw('beta', [0.6 3.2], 1, 1e5);
                    %  3.   y = randraw('beta', [-10 20 3.1 6.2], 1e5 );
                    %  4.   y = randraw('beta', [3 4 5.3 0.7], [1e5 1] );
                    %  5.   randraw('beta');
                    % END beta HELP

                    % Refernce:
                    %      Dagpunar, John.
                    %      Principles of Random Variate Generation.
                    %      Oxford University Press, 1988.
                    %
                    %  max_ab < 0.5            Joehnk's algorithm
                    %  1 < min_ab              Cheng's algortihm BB
                    %  min_ab <= 1 <= max_ab   Atkinson's switching algorithm
                    %  0.5<= max_ab < 1        Atkinson's switching algorithm


                    checkParamsNum(funcName, 'Beta', 'beta', distribParams, [2, 4]);

                    if numel(distribParams) == 2
                         a = distribParams(1);
                         b = distribParams(2);
                         m = 0;
                         n = 1;
                         validateParam(funcName, 'Beta', 'beta', '[a, b]', 'a', a, {'> 0'});
                         validateParam(funcName, 'Beta', 'beta', '[a, b]', 'b', b, {'> 0'});
                    else
                         m = distribParams(1);
                         n = distribParams(2);
                         a = distribParams(3);
                         b = distribParams(4);
                         validateParam(funcName, 'Beta', 'beta', '[m, n, a, b]', 'n-m', n-m, {'>= 0'});
                         validateParam(funcName, 'Beta', 'beta', '[m, n, a, b]', 'a', a, {'> 0'});
                         validateParam(funcName, 'Beta', 'beta', '[m, n, a, b]', 'b', b, {'> 0'});

                    end

                    sampleSizeIn = sampleSize;
                    sampleSize = [ 1, prod( sampleSizeIn ) ];

                    max_ab = max( a, b );
                    min_ab = min( a, b );
                    if max_ab < 0.5
                         %  Use log(u1^a) and log(u2^b), rather than a and b, to avoid
                         %  underflow for very small a or b.

                         loga = log(rand( sampleSize ))/a;
                         logb = log(rand( sampleSize ))/b;
                         logsum = (loga>logb).*(loga + log(1+ exp(logb-loga))) + ...
                              (loga<=logb).*(logb + log(1+ exp(loga-logb)));
                         out = exp(loga - logsum);

                         indxs = find( logsum > 0);

                         while ~isempty( indxs )
                              indxsSize = size( indxs );
                              loga = log(rand( indxsSize ))/a;
                              logb = log(rand( indxsSize ))/b;
                              logsum = (loga>logb).*(loga + log(1+ exp(logb-loga))) + ...
                                   (loga<=logb).*(logb + log(1+ exp(loga-logb)));

                              l = (logsum <= 0);
                              out( indxs( l ) ) = exp(loga(l) - logsum(l));
                              indxs = indxs( ~l );
                         end

                    elseif min_ab > 1
                         % Algorithm BB

                         sum_ab = a + b;
                         lambda = sqrt((sum_ab-2)/(2*a*b-sum_ab));
                         c = min_ab+1/lambda;

                         u1 = rand( sampleSize );
                         u2 = rand( sampleSize );
                         v = lambda*log(u1./(1-u1));
                         z = u1.*u1.*u2;
                         clear('u1'); clear('u2');
                         w = min_ab*exp(v);
                         r = c*v-1.38629436112;
                         clear('v');
                         s = min_ab+r-w;
                         if a == min_ab
                              out = w./(max_ab+w);
                         else
                              out = max_ab./(max_ab+w);
                         end

                         t = log(z);
                         indxs = find( (s+2.609438 < 5*z) & (r+sum_ab*log(sum_ab./(max_ab+w)) < t) );

                         clear('v');
                         clear('z');
                         clear('w');
                         clear('r');
                         while ~isempty( indxs )
                              indxsSize = size( indxs );

                              u1 = rand( indxsSize );
                              u2 = rand( indxsSize );
                              v = lambda*log(u1./(1-u1));
                              z = u1.*u1.*u2;
                              clear('u1'); clear('u2');
                              w = min_ab*exp(v);
                              r = c*v-1.38629436112;
                              clear('v');
                              s = min_ab+r-w;
                              t = log(z);

                              l = (s+2.609438 >= 5*z) | (r+sum_ab*log(sum_ab./(max_ab+w)) >= t);
                              if a == min_ab
                                   out( indxs( l ) ) = w(l)./(max_ab+w(l));
                              else
                                   out( indxs( l ) ) = max_ab./(max_ab+w(l));
                              end
                              indxs = indxs( ~l );

                         end

                    elseif min_ab < 1 & max_ab > 1
                         %  Atkinson's switching method

                         t = (1-min_ab)/(1+max_ab - min_ab);
                         r = max_ab*t/(max_ab*t + min_ab*(1-t)^max_ab);

                         u1 = rand( sampleSize );
                         w = zeros( sampleSize );
                         l = u1 < r;
                         w( l ) =  t*(u1( l )/r).^(1/min_ab);
                         l = ~l;
                         w( l ) =  1- (1-t)*((1-u1(l))/(1-r)).^(1/max_ab);
                         if a == min_ab
                              out = w;
                         else
                              out = 1 - w;
                         end
                         u2 = rand( sampleSize );

                         indxs1 = find(u1 < r);
                         indxs2 = find(u1 >= r);
                         clear('u1');
                         indxs = [ indxs1( log(u2(indxs1)) >= (max_ab-1)*log(1-w(indxs1)) ), ...
                              indxs2( log(u2(indxs2)) >= (min_ab-1)*log(w(indxs2)/t) ) ];

                         clear('u1');
                         clear('u2');
                         while ~isempty( indxs )
                              indxsSize = size( indxs );
                              u1 = rand( indxsSize );
                              w  = zeros( indxsSize );
                              l = u1 < r;
                              w( l ) =  t*(u1( l )/r).^(1/min_ab);
                              l = ~l;
                              w( l ) =  1- (1-t)*((1-u1(l))/(1-r)).^(1/max_ab);

                              u2 = rand( indxsSize );

                              indxs1 = find(u1 < r);
                              indxs2 = find(u1 >= r);
                              clear('u1');
                              l = logical( zeros( indxsSize ) );
                              l( [ indxs1( log(u2(indxs1)) < (max_ab-1)*log(1-w(indxs1)) ), ...