1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269
| function [r,type,coefs] = pearsrnd(mu,sigma,skew,kurt,varargin)
%PEARSRND Random arrays from the Pearson system of distributions.
% R = PEARSRND(MU,SIGMA,SKEW,KURT,M,N) returns an M-by-N matrix of random
% numbers drawn from the distribution in the Pearson system with mean MU,
% standard deviation SIGMA, skewness SKEW, and kurtosis KURT. MU, SIGMA,
% SKEW, and KURT must be scalars.
%
% Note: Because R is a random sample, its sample moments, especially the
% skewness and kurtosis, will typically differ somewhat from the specified
% distribution moments.
%
% Some combinations of moments are not valid for any random variable, and in
% particular, the kurtosis must be greater than the square of the skewness
% plus 1. The kurtosis of the normal distribution is defined to be 3.
%
% R = PEARSRND(MU,SIGMA,SKEW,KURT) returns a scalar value.
% R = PEARSRND(MU,SIGMA,SKEW,KURT,M,N,...) or
% R = PEARSRND(MU,SIGMA,SKEW,KURT,[M,N,...]) returns an M-by-N-by-... array.
%
% [R,TYPE] = PEARSRND(...) returns the type of the specified distribution
% within the Pearson system. Type is a scalar integer from 0 to 7. Set M
% and N to zero to identify the distribution type without generating any
% random values.
%
% The seven distribution types in the Pearson system correspond to the
% following distributions:
%
% Type 0: Normal distribution
% Type 1: Four-parameter beta
% Type 2: Symmetric four-parameter beta
% Type 3: Three-parameter gamma
% Type 4: Not related to any standard distribution. Density proportional
% to (1+((x-a)/b)^2)^(-c) * exp(-d*arctan((x-a)/b)).
% Type 5: Inverse gamma location-scale
% Type 6: F location-scale
% Type 7: Student's t location-scale
%
% [R,TYPE,C] = PEARSRND(...) returns the coefficients of the quadratic
% polynomial that defines the distribution via the differential equation
% d(log(p(x)))/dx = (a + x) / (c(0) + c(1)*x + c(2)*x^2).
%
% Examples
% % Generate random values from the standard normal distribution
% r = pearsrnd(0,1,0,3,100,1); % equivalent to randn(100,1)
%
% % Determine the distribution type
% [r,type] = pearsrnd(0,1,1,4,0,0); % returns [] for r
%
% See also RANDOM, JOHNSRND.
% PEARSRND uses transformations of various standard random variates for types
% 0-III and types V-VII, and a rejection algorithm for type IV.
% References:
% [1] Johnson, N.L., S. Kotz, and N. Balakrishnan (1994) Continuous
% Univariate Distributions, Volume 1, Wiley-Interscience.
% [2] Devroye, L. (1986) Non-Uniform Random Variate Generation,
% Springer-Verlag.
% Copyright 2005 The MathWorks, Inc.
% $Revision: 1.1.6.2 $ $Date: 2005/12/12 23:33:54 $
if nargin < 4
error('stats:pearsrnd:TooFewInputs','Requires at least four input arguments.');
elseif ~(isscalar(mu) && isscalar(sigma) && isscalar(skew) && isscalar(kurt))
error('stats:pearsrnd:NonScalarInputs','MU, SIGMA, SKEW, and KURT must be scalars.');
end
[err, sizeOut] = statsizechk(4,mu,sigma,skew,kurt,varargin{:});
if err > 0
error('stats:pearsrnd:InputSizeMismatch','Size information is inconsistent.');
end
outClass = superiorfloat(mu,sigma,skew,kurt);
beta1 = skew.^2;
beta2 = kurt;
% Return NaN for illegal parameter values.
if (sigma < 0) || (beta2 <= beta1 + 1)
r = NaN(sizeOut,outClass);
type = NaN;
coefs = NaN(1,3,outClass);
return
end
% Classify the distribution and find the roots of c0 + c1*x + c2*x^2
c0 = (4*beta2 - 3*beta1); % ./ (10*beta2 - 12*beta1 - 18);
c1 = skew .* (beta2 + 3); % ./ (10*beta2 - 12*beta1 - 18);
c2 = (2*beta2 - 3*beta1 - 6); % ./ (10*beta2 - 12*beta1 - 18);
if c1 == 0 % symmetric dist'ns
if beta2 == 3
type = 0;
else
if beta2 < 3
type = 2;
elseif beta2 > 3
type = 7;
end
a1 = -sqrt(abs(c0./c2));
a2 = -a1; % symmetric roots
end
elseif c2 == 0 % kurt = 3 + 1.5*skew^2
type = 3;
a1 = -c0 ./ c1; % single root
else
kappa = c1.^2 ./ (4*c0.*c2);
if kappa < 0
type = 1;
elseif kappa < 1-eps
type = 4;
elseif kappa <= 1+eps
type = 5;
else
type = 6;
end
% Solve the quadratic for general roots a1 and a2 and sort by their real parts
tmp = -(c1 + sign(c1).*sqrt(c1.^2 - 4*c0.*c2)) ./ 2;
a1 = tmp ./ c2;
a2 = c0 ./ tmp;
if (real(a1) > real(a2)), tmp = a1; a1 = a2; a2 = tmp; end
end
denom = (10*beta2 - 12*beta1 - 18);
if abs(denom) > sqrt(realmin)
c0 = c0 ./ denom;
c1 = c1 ./ denom;
c2 = c2 ./ denom;
coefs = [c0 c1 c2];
else
type = 1; % this should have happened already anyway
% beta2 = 1.8 + 1.2*beta1, and c0, c1, and c2 -> Inf. But a1 and a2 are
% still finite.
coefs = Inf(1,3,outClass);
end
% generate standard (zero mean, unit variance) values
switch type
case 0
% normal: standard support (-Inf,Inf)
m1 = zeros(outClass);
m2 = ones(outClass);
r = normrnd(m1,m2,sizeOut);
case 1
% four-parameter beta: standard support (a1,a2)
if abs(denom) > sqrt(realmin)
m1 = (c1 + a1) ./ (c2 .* (a2 - a1));
m2 = -(c1 + a2) ./ (c2 .* (a2 - a1));
else
% c1 and c2 -> Inf, but c1/c2 has finite limit
m1 = c1 ./ (c2 .* (a2 - a1));
m2 = -c1 ./ (c2 .* (a2 - a1));
end
r = a1 + (a2 - a1) .* betarnd(m1+1,m2+1,sizeOut);
case 2
% symmetric four-parameter beta: standard support (-a1,a1)
m = (c1 + a1) ./ (c2 .* 2*abs(a1));
r = a1 + 2*abs(a1) .* betarnd(m+1,m+1,sizeOut);
case 3
% three-parameter gamma: standard support (a1,Inf) or (-Inf,a1)
m = (c0./c1 - c1) ./ c1;
r = c1 .* gamrnd(m+1,1,sizeOut) + a1;
case 4
% Pearson IV is not a transformation of a standard distribution: density
% proportional to (1+((x-lambda)/a)^2)^(-m) * exp(-nu*arctan((x-lambda)/a)),
% standard support (-Inf,Inf)
m = 1 ./ (2*c2);
nu = 2.*c1.*(1 - m) ./ sqrt((4.*c0.*c2 - c1.^2));
b = 2*(m-1);
a = sqrt(b.^2 .* (b-1) ./ (b.^2 + nu.^2)); % gives unit variance
lambda = a.*nu ./ b; % gives zero mean
r = pearson4rnd(m,nu,a,lambda,sizeOut);
case 5
% inverse gamma location-scale: standard support (-C1,Inf) or (-Inf,-C1)
C1 = c1 ./ (2*c2);
r = -((c1 - C1) ./ c2) ./ gamrnd(1./c2 - 1,1,sizeOut) - C1;
case 6
% F location-scale: standard support (a2,Inf) or (-Inf,a1)
m1 = (a1 + c1) ./ (c2.*(a2 - a1));
m2 = -(a2 + c1) ./ (c2.*(a2 - a1));
% a1 and a2 have the same sign, and they've been sorted so a1 < a2
if a2 < 0
nu1 = 2*(m2 + 1);
nu2 = -2*(m1 + m2 + 1);
r = a2 + (a2 - a1) .* (nu1./nu2) .* frnd(nu1,nu2,sizeOut);
else % 0 < a1
nu1 = 2*(m1 + 1);
nu2 = -2*(m1 + m2 + 1);
r = a1 + (a1 - a2) .* (nu1./nu2) .* frnd(nu1,nu2,sizeOut);
end
case 7
% t location-scale: standard support (-Inf,Inf)
nu = 1./c2 - 1;
r = sqrt(c0 ./ (1-c2)) .* trnd(nu,sizeOut);
end
% scale and shift
r = r.*sigma + mu;
function r = pearson4rnd(m,nu,a,lambda,sizeOut)
% PEARSON4RND Generate Pearson type 4 random variates.
%
% Based on the exponential rejection method for log-concave densities in
% Devroye, Section VII.2. Valid only when m>1, which is if called by PEARSRND.
%
% References:
% [1] Devroye, L. (1986) Non-Uniform Random Variate Generation,
% Springer-Verlag. Also available in PDF format on-line at
% <a href="http://cgm.cs.mcgill.ca/~luc/rnbookindex.html" target="_blank">http://cgm.cs.mcgill.ca/~luc/rnbookindex.html</a>.
% [2] Heinrich, J. (2004) "A Guide to the Pearson Type IV Distribution",
% CDF/MEMO/STATISTICS/PUBLIC/6820, available on-line at
% <a href="http://www-cdf.fnal.gov/publications/cdf6820_pearson4.pdf" target="_blank">http://www-cdf.fnal.gov/publications...0_pearson4.pdf</a>.
K = (1./HypGeo(m,nu/2)).*exp(gammaln(m) - gammaln(m-.5)) ./ (sqrt(pi)*a);
% Generate y = arctan(x) with density g(y) = K*cos(y)^(2m-1)*exp(-nu*y)
b = 2*(m-1);
M = atan(-nu./b); % mode of y = arctan(x)
cosM = a ./ sqrt(b-1);
loggM = b.*log(cosM) - nu.*M; % log(g(mode)) + log(K)
invgM = exp(-loggM) ./ K; % 1/g(mode)
outClass = superiorfloat(m,nu,a,lambda);
r = zeros(sizeOut,outClass);
j = 1:numel(r);
while length(j) > 0
U = 4*rand(size(j)); % dist'd Unif([0,4])
S = (U>2); % use this to get a random +1/-1
U(S) = U(S) - 2; % now dist'd Unif([0,2])
negEstar = log(max(U,1)-(U>1)); % zero for U<=1, dist'd Exp(1) for U>1
X = min(U,1) - negEstar; % U or 1+Estar
Z = log(rand(size(j))) + negEstar; % -E or -E-Estar
X = M + (2*S-1).*X.*invgM;
k = (abs(X) < pi/2) & (Z <= b.*log(abs(cos(X))) - nu.*X - loggM);
r(j(k)) = X(k);
j(k) = [];
end
% Transform, scale, and shift to standard Pearson type IV
r = a.*tan(r) + lambda;
function F = HypGeo(x,y)
% HYPGEO A special case of the hypergeometric function.
%
% Returns F(-iy,iy,x,1) = abs(gamma(x)/gamma(x+iy))^2, where F is the complex
% hypergeometric function. Based on methods described in Heinrich, J. (2004) "A
% Guide to the Pearson Type IV Distribution", CDF/MEMO/STATISTICS/PUBLIC/6820.
% For small x, compute (1+(y/x)^2)*...*(1+(y/(x+n))^2) which scales F(-iy,iy,x,1)
% to F(-iy,iy,x+n,1), which we can compute quickly if x+n is large.
if x < 100
xstep = x:1:99;
r = prod(1 + (y./xstep).^2);
x = xstep(end) + 1;
else
r = 1;
end
% Compute F(-iy,iy,x+n,1), then multiply by r to get F(-iy,iy,x,1)
s = ones(class(y)); p = ones(class(y)); f = zeros(class(y));
while p > eps(s)
p = p .* (y.^2 + f.^2) ./ (x.*(f+1));
x = x + 1;
f = f + 1;
s = s + p;
end
F = r.*s; |
Partager