WGGAMRND Random matrices from a Generalized Gamma distribution. CALL: R = wggamrnd(a,b,c,sz); a,b,c = parameters (see wggampdf) 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). Example: R=wggamrnd(1,2,4,1,100); plot(R,'.') See also wggaminv
Check if all input arguments are either scalar or of common size. | |
Display message and abort function. | |
Not-a-Number. |
001 function R = wggamrnd(a,b,c,varargin); 002 %WGGAMRND Random matrices from a Generalized Gamma distribution. 003 % 004 % CALL: R = wggamrnd(a,b,c,sz); 005 % 006 % a,b,c = parameters (see wggampdf) 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 % Example: 011 % R=wggamrnd(1,2,4,1,100); 012 % plot(R,'.') 013 % 014 % See also wggaminv 015 016 % Reference: 017 % rejection sampling: there is no reference 018 % inversion method: 019 % Cohen & Whittle, (1988) "Parameter Estimation in Reliability 020 % and Life Span Models", p. 220 ff, Marcel Dekker. 021 022 % Tested on; Matlab 5.3 023 % History: 024 % added ms 09.08.2000 025 % revised pab 23.10.2000 026 % % - added default b and c 027 % - added comnsize, nargchk 028 % - added greater flexibility on the sizing of R 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 if nargin<2|isempty(c), c=1;end 041 042 if nargin<3, 043 [errorcode a b c] = comnsize(a,b,c); 044 else 045 [errorcode a b c] = comnsize(a,b,c,zeros(varargin{:})); 046 end 047 if errorcode > 0 048 error('a, b and c must be of common size or scalar.'); 049 end 050 051 % 052 %R=wggaminv(rand(size(a)),a,b,c); % slow 053 %return 054 055 R = zeros(size(a)); 056 057 ok = ((a>0) & (b>0) & c>0); 058 k = find(ok); 059 if any(k), 060 ak=a(k); 061 y0 = log(ak)-1./sqrt(ak); 062 c0 = ak - exp(y0); 063 c1 =(ak.*y0 - exp(y0)); 064 c2 = abs((y0-log(ak))); 065 066 accept=k; omit=[]; 067 while ~isempty(accept) 068 ak(omit)=[]; c0(omit) =[]; 069 c1(omit)=[]; c2(omit)=[]; 070 sz = size(ak); 071 la = log(ak); 072 y = log(rand(sz)).*sign(rand(sz)-0.5)./c0 + la; 073 074 f = ak.*y-exp(y) - c1; 075 g = c0.*(c2 - abs(y-la)); 076 077 omit = find((log(rand(sz)) + g) <= f); 078 if ~isempty(omit) 079 R(accept(omit)) = exp(y(omit)); 080 accept(omit)=[]; 081 end 082 end % while 083 R(k) = R(k).^(1./c(k)).*b(k); 084 end 085 086 k1=find(~ok); 087 if any(k1); 088 tmp=NaN; 089 R(k1)=tmp(ones(size(k1))); 090 end 091 return 092 093 094
Comments or corrections to the WAFO group