POISSRND Random matrices from Poisson distribution.
R = POISSRND(LAMBDA) returns a matrix of random numbers chosen
from the Poisson distribution with parameter LAMBDA.
The size of R is the size of LAMBDA. Alternatively,
R = POISSRND(LAMBDA,M,N) returns an M by N matrix.
Use R = POISSRND(LAMBDA,M,N,STATE) to initialise the random
number generator with STATE. Use R = POISSRND(LAMBDA,[],[],STATE)
when M and N are not given and POISSRND([],[],[],STATE) when
only the random number generator needs to be initialised.
POISSRND([],[],[],0) resets the generator to its initial state.
POISSRND([],[],[],J), for integer J, resets the generator to its J-th state.
POISSRND([],[],[],sum(100*clock)) resets it to a different state each time.
Note: POISSRND actually initialises RAND and RANDN.
Use R = POISSRND(LAMBDA,M,N,STATE,METHOD) to determine the method.
METHOD='normal' (default) uses a waiting time method.
METHOD='fast' uses a waiting time method at lambda<1000 and
normal distribution with mean and variance equal to lambda
at labda>=1000.
METHOD=value equal to 'fast' but with the border set at value
instead of 1000.
Use R = POISSRND(LAMBDA,[],[],[],METHOD) when M, N and STATE are
not given.
For more information see the M-file.
Addition 1 (12-3-1999):
This version of POISSRND has been improved in performance
by Mischa Tolsma, MMR-TN TU Delft. It removes elements of
which the result already has been found from the search matrix.
It's major speedup occurs when lambda contains a lot of
small elements with just a few large ones.
Speed up factors:
0.1 % large elements: 10
100 % large elements: 1.3 both at an array of 1000
Addition 2 (14-7-2000):
This version of POISSRND allows the use of normal distribution as an
approximation of the Poisson distribution at large values of lambda.
The Poisson distribution resembles a normal distribution with mean
lambda and variance lambda when lambda is large. But the standard
version of POISSRND needs a large amount of time for large values
of lambda.
Indication of speed up factor: 80 at lambda=10,000 and array of 100.
Addition 3 (7-12-2000):
The new version of POISSRND of release has been improved in speed
and has become faster starting from lambda = 15. These improvements
are therefore incorporated into this version.
For large lambda (between 15 and 1000), use the method of Ahrens
and Dieter as described in Knuth, Volume 2, 1998 edition.
note: this version calls the GAMRND and the BINORND functions.
References:
[1] L. Devroye, "Non-Uniform Random Variate Generation",
Springer-Verlag, 1986 page 504.
Based on poissrnd.
M.F.P. Tolsma, Signals, Systems and Control Group, Applied Physics, TU Delft
http://www.tn.tudelft.nl/mmr
copyright remains with author
$Revision: 1.1 $ $Date: 2004/05/12 16:22:33 $
imported from Matlab Central - File Exchange
$Id: poissonrnd.m,v 1.1 2004/05/12 16:22:33 dalai Exp $0001 function r = poissonrnd(lambda,m,n,state,method) 0002 %POISSRND Random matrices from Poisson distribution. 0003 % R = POISSRND(LAMBDA) returns a matrix of random numbers chosen 0004 % from the Poisson distribution with parameter LAMBDA. 0005 % 0006 % The size of R is the size of LAMBDA. Alternatively, 0007 % R = POISSRND(LAMBDA,M,N) returns an M by N matrix. 0008 % 0009 % Use R = POISSRND(LAMBDA,M,N,STATE) to initialise the random 0010 % number generator with STATE. Use R = POISSRND(LAMBDA,[],[],STATE) 0011 % when M and N are not given and POISSRND([],[],[],STATE) when 0012 % only the random number generator needs to be initialised. 0013 % POISSRND([],[],[],0) resets the generator to its initial state. 0014 % POISSRND([],[],[],J), for integer J, resets the generator to its J-th state. 0015 % POISSRND([],[],[],sum(100*clock)) resets it to a different state each time. 0016 % Note: POISSRND actually initialises RAND and RANDN. 0017 % 0018 % Use R = POISSRND(LAMBDA,M,N,STATE,METHOD) to determine the method. 0019 % METHOD='normal' (default) uses a waiting time method. 0020 % METHOD='fast' uses a waiting time method at lambda<1000 and 0021 % normal distribution with mean and variance equal to lambda 0022 % at labda>=1000. 0023 % METHOD=value equal to 'fast' but with the border set at value 0024 % instead of 1000. 0025 % Use R = POISSRND(LAMBDA,[],[],[],METHOD) when M, N and STATE are 0026 % not given. 0027 % 0028 % For more information see the M-file. 0029 % 0030 % Addition 1 (12-3-1999): 0031 % This version of POISSRND has been improved in performance 0032 % by Mischa Tolsma, MMR-TN TU Delft. It removes elements of 0033 % which the result already has been found from the search matrix. 0034 % It's major speedup occurs when lambda contains a lot of 0035 % small elements with just a few large ones. 0036 % 0037 % Speed up factors: 0038 % 0.1 % large elements: 10 0039 % 100 % large elements: 1.3 both at an array of 1000 0040 % 0041 % Addition 2 (14-7-2000): 0042 % This version of POISSRND allows the use of normal distribution as an 0043 % approximation of the Poisson distribution at large values of lambda. 0044 % The Poisson distribution resembles a normal distribution with mean 0045 % lambda and variance lambda when lambda is large. But the standard 0046 % version of POISSRND needs a large amount of time for large values 0047 % of lambda. 0048 % Indication of speed up factor: 80 at lambda=10,000 and array of 100. 0049 % 0050 % Addition 3 (7-12-2000): 0051 % The new version of POISSRND of release has been improved in speed 0052 % and has become faster starting from lambda = 15. These improvements 0053 % are therefore incorporated into this version. 0054 % For large lambda (between 15 and 1000), use the method of Ahrens 0055 % and Dieter as described in Knuth, Volume 2, 1998 edition. 0056 % note: this version calls the GAMRND and the BINORND functions. 0057 % 0058 % References: 0059 % [1] L. Devroye, "Non-Uniform Random Variate Generation", 0060 % Springer-Verlag, 1986 page 504. 0061 % 0062 % 0063 % Based on poissrnd. 0064 % M.F.P. Tolsma, Signals, Systems and Control Group, Applied Physics, TU Delft 0065 % http://www.tn.tudelft.nl/mmr 0066 % copyright remains with author 0067 % 0068 % $Revision: 1.1 $ $Date: 2004/05/12 16:22:33 $ 0069 % 0070 % imported from Matlab Central - File Exchange 0071 % $Id: poissonrnd.m,v 1.1 2004/05/12 16:22:33 dalai Exp $ 0072 0073 0074 %------------------------------Input argument verification 0075 if nargin < 1, 0076 error('Requires at least one input argument.'); 0077 end 0078 0079 if nargin > 3 0080 %random seed 0081 if ~isempty(state) 0082 rand('state',state) 0083 randn('state',state) 0084 if isempty(lambda) 0085 % F. WILLAME small change 0086 % change break to return to avoid Matlab complaining 0087 % break <-- original file 0088 return 0089 end 0090 end 0091 end 0092 0093 [rows columns] = size(lambda); 0094 0095 if nargin==2 0096 if ~isempty(m) 0097 if rows*columns>1 0098 error('Lambda should be a scalar.'); 0099 else 0100 rows=m; 0101 end 0102 end 0103 end 0104 0105 if nargin>2 0106 if ~isempty(n) 0107 if rows*columns>1 0108 error('Lambda should be a scalar.'); 0109 else 0110 if ~isempty(m) 0111 rows=m; 0112 end 0113 columns=n; 0114 end 0115 end 0116 end 0117 0118 if prod(size(lambda))==1 0119 lambda=lambda*ones(rows,columns); 0120 end 0121 0122 if nargin==5 0123 if ~isempty(findstr(method,'normal')) 0124 method=0; 0125 elseif ~isempty(findstr(method,'fast')) 0126 method=1000; 0127 elseif isempty(method) 0128 method=0; 0129 end 0130 else 0131 method=0; 0132 end 0133 0134 0135 %Initialize r to zero. 0136 0137 r = zeros(rows, columns); 0138 index = [1:rows]'*ones(1,columns)+rows*ones(rows,1)*[0:columns-1]; 0139 0140 % Return NaN if LAMBDA is not positive. 0141 if any(any(lambda <= 0)); 0142 if prod(size(lambda) == 1) 0143 r = NaN * ones(rows,columns); 0144 lambda=[]; 0145 else 0146 k = find(lambda <= 0); 0147 r(k) = NaN * ones(size(k)); 0148 0149 k = find(lambda > 0); 0150 lambda=lambda(k); % Remove elements for which a result has been calculated. 0151 index=index(k); 0152 end 0153 end 0154 0155 %values larger then 'method' are randomised using a discrete normal distribution. 0156 %This is approximatly equal to the poisson distribution. 0157 if method>0 0158 k = find(lambda>=method); % The lambda's larger than method 0159 if ~isempty(k) 0160 r(index(k))=round(lambda(k)+sqrt(lambda(k)).*randn(size(lambda(k)))); 0161 0162 k = find(lambda<method); % The lambda's smaller than method 0163 lambda=lambda(k); % Remove elements for which a result has been calculated. 0164 index=index(k); 0165 end 0166 end 0167 0168 %values larger (or equal to 15 are randomised using the method of Ahrens and Dieter 0169 k = find(lambda >= 15); 0170 if ~isempty(k) 0171 alpha = 7/8; 0172 lk=lambda(k); 0173 m = floor(alpha * lk); 0174 0175 % Generate m waiting times, all at once 0176 x = gamrnd(m,1); 0177 k2= find(x <= lk); 0178 0179 if ~isempty(k2) 0180 % If we did not overshoot, then the number of additional times 0181 % has a Poisson distribution with a smaller mean. 0182 r(index(k(k2))) = m(k2) + poissrnd(lk(k2)-x(k2)); 0183 end 0184 0185 k2= find(x > lk); 0186 if ~isempty(k2) 0187 % If we did overshoot, then the times up to m-1 are uniformly 0188 % distributed on the interval to x, so the count of times less 0189 % than lambda has a binomial distribution. 0190 r(index(k(k2))) = binornd(m(k2)-1, lk(k2)./x(k2)); 0191 end 0192 0193 k = find(lambda<15); 0194 lambda=lambda(k); 0195 index=index(k); 0196 end 0197 0198 kt=0; 0199 p=zeros(size(index)); 0200 while ~isempty(index) 0201 p = p - log(rand(size(p))); 0202 k = find(p >= lambda); % The r's larger than border 0203 r(index(k))=kt; 0204 0205 k = find(p < lambda); % The r's for which the border hasn't been reached. 0206 lambda=lambda(k); % Remove elements for which a result has been calculated. 0207 p=p(k); 0208 index=index(k); 0209 0210 kt=kt+1; 0211 end 0212