Home > wafo > wstats > wggamrnd.m

wggamrnd

PURPOSE ^

Random matrices from a Generalized Gamma distribution.

SYNOPSIS ^

R = wggamrnd(a,b,c,varargin);

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group


Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003