WGAMRND Random matrices from a Gamma distribution. CALL: R = wgamrnd(a,b,sz); a,b = parameters, a,b>0 sz = size(R) (Default common size of a and b) sz can be a comma separated list or a vector giving the size of R (see zeros for options). The random numbers are generated by rejection sampling. Example: R = wgamrnd(2,1,1,100); plot(R,'.') See also wgaminv
Check if all input arguments are either scalar or of common size. | |
Display message and abort function. | |
Not-a-Number. |
Random points from a bivariate DIST2D distribution | |
Random points from a bivariate MDIST2D distribution | |
Random matrices from a Chi squared distribution. |
001 function R = wgamrnd(a,b,varargin); 002 %WGAMRND Random matrices from a Gamma distribution. 003 % 004 % CALL: R = wgamrnd(a,b,sz); 005 % 006 % a,b = parameters, a,b>0 007 % sz = size(R) (Default common size of a and b) 008 % sz can be a comma separated list or a vector 009 % giving the size of R (see zeros for options). 010 % 011 % The random numbers are generated by rejection sampling. 012 % 013 % Example: 014 % R = wgamrnd(2,1,1,100); 015 % plot(R,'.') 016 % 017 % See also wgaminv 018 019 020 % Tested on: matlab 5.3 021 % History: 022 % based on rgamma.m (stixbox) by (c) Anders Holtsberg 10-05-2000. 023 % revised pab 23.10.2000 024 % - added default b 025 % - added comnsize, nargchk 026 % - added greater flexibility on the sizing of R 027 % - rejection sampling method modified to give any size of a and b 028 % the recursive call in rgamma is replaced with a while loop. 029 % - replaced inversion method with a modified version of 030 % rgamma from stixbox (Anders Holtsberg) 031 % The algorithm is a rejection method. The logarithm of the gamma 032 % variable is simulated by dominating it with a double exponential. 033 % The proof is easy since the log density is convex! 034 % 035 % Reference: There is no reference! Send me (Anders Holtsberg) an email 036 % if you can't figure it out. 037 038 error(nargchk(1,inf,nargin)) 039 if nargin<2|isempty(b), b=1; end 040 041 if nargin<3, 042 [errorcode a b] = comnsize(a,b); 043 else 044 [errorcode a b] = comnsize(a,b,zeros(varargin{:})); 045 end 046 047 if errorcode > 0, 048 error('a and b must be of common size or scalar.'); 049 end 050 051 %R = wgaminv(rand(size(a)),a,b); % slower 052 %return 053 R = zeros(size(a)); 054 055 ok = ((a>0) & (b>0)); 056 k = find(ok); 057 if any(k), 058 ak=a(k); 059 y0 = log(ak)-1./sqrt(ak); 060 c = ak - exp(y0); 061 c1 =(ak.*y0 - exp(y0)); 062 c2 = abs((y0-log(ak))); 063 064 accept=k; omit=[]; 065 while ~isempty(accept) 066 ak(omit)=[]; c(omit) =[]; 067 c1(omit)=[]; c2(omit)=[]; 068 sz = size(ak); 069 la = log(ak); 070 y = log(rand(sz)).*sign(rand(sz)-0.5)./c + la; 071 072 f = ak.*y-exp(y) - c1; 073 g = c.*(c2 - abs(y-la)); 074 075 omit = find((log(rand(sz)) + g) <= f); 076 if ~isempty(omit) 077 R(accept(omit)) = exp(y(omit)); 078 accept(omit)=[]; 079 end 080 end % while 081 R(k) = R(k).*b(k); 082 end 083 084 k1=find(~ok); 085 if any(k1); 086 tmp=NaN; 087 R(k1)=tmp(ones(size(k1))); 088 end 089 return 090 091 092 093
Comments or corrections to the WAFO group