📄 randraw.m
字号:
indxs2( log(u2(indxs2)) < (min_ab-1)*log(w(indxs2)/t) ) ] ) = 1; clear('u1'); clear('u2'); if a == min_ab out( indxs(l) ) = w(l); else out( indxs(l) ) = 1 - w(l); end indxs = indxs( ~l ); end else % Atkinson's Algorithm if min_ab == 1 t = 0.5; r = 0.5; else t = 1/(1+sqrt(max_ab*(1-max_ab)/(min_ab*(1-min_ab)))); r = max_ab*t / (max_ab*t + min_ab*(1-t)); end u1 = rand( sampleSize ); out = zeros( sampleSize ); w = zeros( sampleSize ); l1 = u1 < r; w(l1) = t*(u1(l1)/r).^(1/min_ab); l2 = u1 >= r; w(l2) = 1 - (1-t)*((1-u1(l2))/(1-r)).^(1/max_ab); if a == min_ab out = w; else out = 1 - w; end u2 = rand( sampleSize ); indxs1 = find(l1); indxs2 = find(l2); indxs = [ indxs1( log(u2(l1)) >= (max_ab -1)*log((1-w(l1))/(1-t)) ), ... indxs2( log(u2(l2)) >= (min_ab -1) * log(w(l2)/t) ) ]; clear('u2'); while ~isempty( indxs ) indxsSize = size( indxs ); u1 = rand( indxsSize ); w = zeros( indxsSize ); l1 = u1 < r; w(l1) = t*(u1(l1)/r).^(1/min_ab); l2 = u1 >= r; w(l2) = 1 - (1-t)*((1-u1(l2))/(1-r)).^(1/max_ab); u2 = rand( indxsSize ); indxs1 = find(l1); indxs2 = find(l2); clear('u1'); l = logical( zeros( indxsSize ) ); l( [ indxs1(log(u2(l1)) < (max_ab -1)*log((1-w(l1))/(1-t))), ... indxs2(log(u2(l2))< (min_ab -1) * log(w(l2)/t)) ] ) = 1; if a == min_ab out(indxs(l)) = w(l); else out(indxs(l)) = 1 - w(l); end indxs = indxs( ~l ); end end out = m + (n-m) * out; reshape( out, sampleSizeIn ); case {'bessel'} % START bessel HELP % THE BESSEL DISTRIBUTION % % Bessel distribution arises in the theory of stochastic processes. % Bessel(nu,a) is a discrete distribution on the non-negative integers with % parameters nu > -1 and a > 0. % % pdf(y) = (a/2).^(2*y+nu) ./ (besseli(nu,a).*factorial(y).*gamma(y+nu+1)); % % PARAMETERS: % nu > -1, a > 0 % SUPPORT: % y = 0, 1, 2, 3, ... % CLASS: % Discrete distributions % % USAGE: % randraw('bessel', [nu, a], sampleSize) - generate sampleSize number % of variates from the Bessel distribution with parameters % 'nu' and 'a' % randraw('bessel') - help for the Bessel distribution; % EXAMPLES: % 1. y = randraw('bessel', [2 0.9], [1 1e5]); % 2. y = randraw('bessel', [0.6 3.2], 1, 1e5); % 3. y = randraw('bessel', [-0.2 8.1], 1e5 ); % 4. y = randraw('bessel', [4 5.3], [1e5 1] ); % 5. randraw('bessel'); % END bessel HELP % Method: % % We implemented Condensed Table-Lookup method suggested in % George Marsaglia, "Fast Generation Of Discrete Random Variables," % Journal of Statistical Software, July 2004, Volume 11, Issue 4 % % Reference: % L. Devroye, "Simulating Bessel random variables," % Statistics and Probability Letters, vol. 57, pp. 249-257, 2002. % checkParamsNum(funcName, 'Bessel', 'bessel', distribParams, [2]); nu = distribParams(1); a = distribParams(2); validateParam(funcName, 'Bessel', 'bessel', '[nu, a]', 'nu', nu, {'> -1'}); validateParam(funcName, 'Bessel', 'bessel', '[nu, a]', 'a', a, {'> 0'}); % mu = 0.5*a*besseli(nu+1,a)/besseli(nu,a); % chi2 = mu + 0.25*a^2*besseli(nu+1,a)/besseli(nu,a)*... % (besseli(nu+2,a)/besseli(nu+1,a)-besseli(nu+1,a)/besseli(nu,a)); besseliNuA = besseli(nu, a); proceed = 1; if ~isfinite( besseliNuA ) warnStr{1} = [upper(funcName), ' - Bessel Variates Generation: ']; warnStr{2} = ['besseli(', num2str(nu), ', ' num2str(a), ') returns Inf.']; warnStr{3} = ['Unable to proceed, return zeros ...']; warning('%s\n %s\n %s',warnStr{1},warnStr{2},warnStr{3}); %warning([upper(funcName), ' - Bessel Variates Generation: besseli(', num2str(nu), ', ' num2str(a), ') returns Inf. Unable to proceed, return zeros ...']); out = zeros( sampleSize ); proceed = 0; end if besseliNuA == 0 warnStr{1} = [upper(funcName), ' - Bessel Variates Generation: ']; warnStr{2} = ['besseli(', num2str(nu), ', ' num2str(a), ') returns 0.']; warnStr{3} = ['Unable to proceed, return zeros ...']; warning('%s\n %s\n %s',warnStr{1},warnStr{2},warnStr{3}); %warning([upper(funcName), '- Bessel Variates Generation: besseli(', num2str(nu), ', ' num2str(a), ') returns 0. Unable to proceed, return zeros ...']); out = zeros( sampleSize ); proceed = 0; end if proceed p0 = exp( nu*log(a/2) - gammaln(nu+1) ) / besseliNuA; if p0 >= 5e-10 t = p0; aa = (a/2)^2; nu1 = nu+1; i = 1; while t*2147483648 > 1 t = t * aa/((i)*(i+nu)); i = i + 1; end sizeP = i-1; offset = 0; P = round( 2^30*p0*cumprod([1, aa./((1:sizeP-1).*((1:sizeP-1)+nu))] ) ); else % if p0 >= 5e-10 m = floor(0.5*(sqrt(a^2+nu^2)-nu)); pm = exp( (2*m+nu)*log(a/2) - log(besseliNuA) - ... gammaln(m+1) - gammaln(m+nu+1) ); aa = (a/2)^2; t = pm; i = m + 1; while t * 2147483648 > 1 t = t * aa/((i)*(i+nu)); i = i + 1; end last = i-2; t = pm; j = -1; for i = m-1:-1:0 t = t * (i+1)*(i+1+nu)/aa; if t*2147483648 < 1 j=i; break; end end offset = j+1; sizeP = last-offset+1; P = zeros(1, sizeP); P(m-offset+1:last-offset+1) = ... round( 2^30*pm*cumprod([1, aa./(((m+1):last).*(((m+1):last)+nu))] ) ); P(m-offset:-1:1) = ... round( 2^30*pm*cumprod((m:-1:(offset+1)).*((m:-1:(offset+1))+nu)/aa) ); end % if p0 >= 5e-10, else ... out = randFrom5Tbls( P, offset, sampleSize); end % if proceed case {'binom', 'binomial'} % START binom HELP START binomial HELP % THE BINOMIAL DISTRIBUTION % % pdf(y) = nchoosek(n,y)*p^y*(1-p)^(n-y) = ... % exp( gammaln(n+1) - gammaln(n-y+1) - gammaln(y+1) + ... % y*log(p) + (n-y)*log(1-p) ); 0<p<1, n>1 % % Mean = n*p; % Variance = n*p*(1-p); % Mode = floor( (n+1)*p ); % % PARAMETERS: % p - probability of success in a single trial; (0<p<1) % n - total number of trials; (n= 1, 2, 3, 4, ...) % % SUPPORT: % y - number of success, y = 0, 1, 2, 3 ... % % CLASS: % Discrete distributions % % NOTES: % Constructive definition: % We consider a random experiment with n independent trials; in each trial % a certain random event A can occur (the urn model with replacement is % a special case of such an experiment). Let % p = probability of A in a single trial; % n = total number of trials; % y = number of successes (= number of trials where A occurs). % % USAGE: % randraw('binom', [n, p], sampleSize) - generate sampleSize number % of variates from the Binomial distribution with total number of trials % 'n' and probability of success in a single trial 'p' % randraw('binom') - help for the Binomial distribution; % % EXAMPLES: % 1. y = randraw('binom', [10 0.9], [1 1e5]); % 2. y = randraw('binom', [100 0.15], 1, 1e5); % 3. y = randraw('binom', [5 0.5], 1e5 ); % 4. y = randraw('binom', [1000 0.02], [1e5 1] ); % 5. randraw('binom'); % END binom HELP END binomial HELP % Method: % % We implemented Condensed Table-Lookup method suggested in % George Marsaglia, "Fast Generation Of Discrete Random Variables," % Journal of Statistical Software, July 2004, Volume 11, Issue 4 checkParamsNum(funcName, 'Binomial', 'binomial', distribParams, [2]); n = distribParams(1); p = distribParams(2); validateParam(funcName, 'Binomial', 'binomial', '[n, p]', 'n', n, {'> 0','==integer'}); validateParam(funcName, 'Binomial', 'binomial', '[n, p]', 'p', p, {'> 0','< 1'}); if n*p > 1e3 & n > 1e3 out = round( n*p + sqrt(n*p*(1-p))*randn( sampleSize ) ); elseif p<1e-4 & n*p > 1 & n*p < 100 out = feval(funcName,'poisson',n*p, sampleSize); else % if n large and p near 1, generate j=Binom(n,1-p), return n-j switchFlag = 0; if n>100 & p>0.99 p = 1-p; switchFlag = 1; end mode = floor( (n+1)*p ); q = 1 - p; h = p/q; pmode = exp( gammaln(n+1) - gammaln(n-mode+1) - gammaln(mode+1) + ... mode*log(p) + (n-mode)*log(1-p) ); if( 1-pmode < 5e-10 ) % "Fast Generation of Discrete Random Variables", p.3 citation: % "For ... discrete variates with an infinite number of probabilities, we select only % those for which, for a sample of size 2^31(10^9.33), the expected number of occurences exceeds % 0.5. The other probabilities are assumed zero. For those unusual situations where occurrences % with probability less than 5 ⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -