function x=randht(n, varargin) % RANDHT generates n observations distributed as some continous heavy- % tailed distribution. Options are power law, log-normal, stretched % exponential, power law with cutoff, and exponential. Can specify lower % cutoff, if desired. % % Example: % x = randht(10000,'powerlaw',alpha); % x = randht(10000,'xmin',xmin,'powerlaw',alpha); % x = randht(10000,'cutoff',alpha, lambda); % x = randht(10000,'exponential',lambda); % x = randht(10000,'lognormal',mu,sigma); % x = randht(10000,'stretched',lambda,beta); % % See also PLFIT, PLVAR, PLPVA % % Source: http://www.santafe.edu/~aaronc/powerlaws/ % Version 1.0 (2007 May) % Version 1.0.1 (2007 September) % Version 1.0.2 (2008 April) % Copyright (C) 2007 Aaron Clauset (Santa Fe Institute) % Distributed under GPL 2.0 % http://www.gnu.org/copyleft/gpl.html % RANDHT comes with ABSOLUTELY NO WARRANTY % % Notes: % type = ''; xmin = 1; alpha = 2.5; beta = 1; lambda = 1; mu = 1; sigma = 1; persistent rand_state; % parse command-line parameters; trap for bad input i=1; while i<=length(varargin), argok = 1; if ischar(varargin{i}), switch varargin{i}, case 'xmin', xmin = varargin{i+1}; i = i + 1; case 'powerlaw', type = 'PL'; alpha = varargin{i+1}; i = i + 1; case 'cutoff', type = 'PC'; alpha = varargin{i+1}; lambda = varargin{i+2}; i = i + 2; case 'exponential', type = 'EX'; lambda = varargin{i+1}; i = i + 1; case 'lognormal', type = 'LN'; mu = varargin{i+1}; sigma = varargin{i+2};i = i + 2; case 'stretched', type = 'ST'; lambda = varargin{i+1}; beta = varargin{i+2}; i = i + 2; otherwise, argok=0; end end if ~argok, disp(['(RANDHT) Ignoring invalid argument #' num2str(i+1)]); end i = i+1; end if (~isscalar(n) || n<1) fprintf('(RANDHT) Error: invalid ''n'' argument; using default.\n'); n = 10000; end; if (~isscalar(xmin) || xmin<1) fprintf('(RANDHT) Error: invalid ''xmin'' argument; using default.\n'); xmin = 1; end; if isempty(rand_state) rand_state = cputime; rand('twister',sum(100*clock)); end; switch type case 'EX', x = xmin - (1/lambda)*log(1-rand(n,1)); case 'LN', y = exp(mu+sigma*randn(10*n,1)); while true y(y0), r = randperm(length(y)); y(r(1:q)) = []; break; end; if (q<0), y = [y; exp(mu+sigma*randn(10*n,1))]; end; end; x = y; case 'ST', x = (xmin^beta - (1/lambda)*log(1-rand(n,1))).^(1/beta); case 'PC', x = []; y = xmin - (1/lambda)*log(1-rand(10*n,1)); while true y(rand(10*n,1)>=(y./xmin).^(-alpha)) = []; x = [x; y]; q = length(x)-n; if (q==0), break; end; if (q>0), r = randperm(length(x)); x(r(1:q)) = []; break; end; if (q<0), y = xmin - (1/lambda)*log(1-rand(10*n,1)); end; end; otherwise, x = xmin*(1-rand(n,1)).^(-1/(alpha-1)); end;