www.gusucode.com > matlab非局部均值工具箱 > matlab非局部均值工具箱/matlab非局部均值工具箱/toolbox_nlmeans/toolbox/stabrnd.m
function [x] = stabrnd(alpha, beta, c, delta, m, n) % STABRND.M % Stable Random Number Generator (McCulloch 12/18/96) % % x = stabrnd(alpha, beta, c, delta, m, n); % % Returns m x n matrix of iid stable random numbers with % characteristic exponent alpha in [.1,2], skewness parameter % beta in [-1,1], scale c > 0, and location parameter delta. % Based on the method of J.M. Chambers, C.L. Mallows and B.W. % Stuck, "A Method for Simulating Stable Random Variables," % JASA 71 (1976): 340-4. % Encoded in MATLAB by J. Huston McCulloch, Ohio State % University Econ. Dept. (mcculloch.2@osu.edu). This 12/18/96 % version uses 2*m*n calls to RAND, and does not rely on % the STATISTICS toolbox. % The CMS method is applied in such a way that x will have the % log characteristic function % log E exp(ixt) = i*delta*t + psi(c*t), % where % psi(t) = -abs(t)^alpha*(1-i*beta*sign(t)*tan(pi*alpha/2)) % for alpha ~= 1, % = -abs(t)*(1+i*beta*(2/pi)*sign(t)*log(abs(t))), % for alpha = 1. % With this parameterization, the stable cdf S(x; alpha, beta, % c, delta) equals S((x-delta)/c; alpha, beta, 1, 0). See my % "On the parametrization of the afocal stable distributions," % _Bull. London Math. Soc._ 28 (1996): 651-55, for details. % When alpha = 2, the distribution is Gaussian with mean delta % and variance 2*c^2, and beta has no effect. % When alpha > 1, the mean is delta for all beta. When alpha % <= 1, the mean is undefined. % When beta = 0, the distribution is symmetrical and delta is % the median for all alpha. When alpha = 1 and beta = 0, the % distribution is Cauchy (arctangent) with median delta. % When the submitted alpha is > 2 or < .1, or beta is outside % [-1,1], an error message is generated and x is returned as a % matrix of NaNs. % Alpha < .1 is not allowed here because of the non-negligible % probability of overflows. % % If you're only interested in the symmetric cases, you may just % set beta = 0 and skip the following considerations: % When beta > 0 (< 0), the distribution is skewed to the right % (left). % When alpha < 1, delta, as defined above, is the unique fractile % that is invariant under averaging of iid contributions. I % call such a fractile a "focus of stability." This, like the % mean, is a natural location parameter. % When alpha = 1, either every fractile is a focus of stability, % as in the beta = 0 Cauchy case, or else there is no focus of % stability at all, as is the case for beta ~=0. In the latter % cases, which I call "afocal," delta is just an arbitrary % fractile that has a simple relation to the c.f. % When alpha > 1 and beta > 0, med(x) must lie very far below % the mean as alpha approaches 1 from above. Furthermore, as % alpha approaches 1 from below, med(x) must lie very far above % the focus of stability when beta > 0. If beta ~= 0, there % is therefore a discontinuity in the distribution as a function % of alpha as alpha passes 1, when delta is held constant. % CMS, following an insight of Vladimir Zolotarev, remove this % discontinuity by subtracting % beta*c*tan(pi*alpha/2) % (equivalent to their -tan(alpha*phi0)) from x for alpha ~=1 % in their program RSTAB, a.k.a. RNSTA in IMSL (formerly GGSTA). % The result is a random number whose distribution is a contin- % uous function of alpha, but whose location parameter (which I % call zeta) is a shifted version of delta that has no known % interpretation other than computational convenience. % The present program restores the more meaningful "delta" % parameterization by using the CMS (4.1), but with % beta*c*tan(pi*alpha/2) added back in (ie with their initial % tan(alpha*phi0) deleted). RNSTA therefore gives different % results than the present program when beta ~= 0. However, % the present beta is equivalent to the CMS beta' (BPRIME). % Rather than using the CMS D2 and exp2 functions to compensate % for the ill-condition of the CMS (4.1) when alpha is very % near 1, the present program merely fudges these cases by % computing x from their (2.4) and adjusting for % beta*c*tan(pi*alpha/2) when alpha is within 1.e-8 of 1. % This should make no difference for simulation results with % samples of size less than approximately 10^8, and then % only when the desired alpha is within 1.e-8 of 1, but not % equal to 1. % The frequently used Gaussian and symmetric cases are coded % separately so as to speed up execution. % % Additional references: % V.M. Zolotarev, _One Dimensional Stable Laws_, Amer. Math. % Soc., 1986. % G. Samorodnitsky and M.S. Taqqu, _Stable Non-Gaussian Random % Processes_, Chapman & Hill, 1994. % A. Janicki and A. Weron, _Simulaton and Chaotic Behavior of % Alpha-Stable Stochastic Processes_, Dekker, 1994. % J.H. McCulloch, "Financial Applications of Stable Distributons," % _Handbook of Statistics_ Vol. 14, forthcoming early 1997. % Errortraps: if alpha < .1 | alpha > 2 disp('Alpha must be in [.1,2] for function STABRND.') alpha x = NaN * zeros(m,n); return end if abs(beta) > 1 disp('Beta must be in [-1,1] for function STABRND.') beta x = NaN * zeros(m,n); return end % Generate exponential w and uniform phi: w = -log(rand(m,n)); phi = (rand(m,n)-.5)*pi; % Gaussian case (Box-Muller): if alpha == 2 x = (2*sqrt(w) .* sin(phi)); x = delta + c*x; return end % Symmetrical cases: if beta == 0 if alpha == 1 % Cauchy case x = tan(phi); else x = ((cos((1-alpha)*phi) ./ w) .^ (1/alpha - 1) ... .* sin(alpha * phi) ./ cos(phi) .^ (1/alpha)); end % General cases: else cosphi = cos(phi); if abs(alpha-1) > 1.e-8 zeta = beta * tan(pi*alpha/2); aphi = alpha * phi; a1phi = (1 - alpha) * phi; x = ((sin(aphi) + zeta * cos(aphi)) ./ cosphi) ... .* ((cos(a1phi) + zeta * sin(a1phi)) ... ./ (w .* cosphi)) .^ ((1-alpha)/alpha); else bphi = (pi/2) + beta * phi; x = (2/pi) * (bphi .* tan(phi) - beta * log((pi/2) * w ... .* cosphi ./ bphi)); if alpha ~= 1 x = x + beta * tan(pi * alpha/2); end end end % Finale: x = delta + c * x; return % End of STABRND.M