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,'.')

## CROSS-REFERENCE INFORMATION

This function calls:
 comnsize Check if all input arguments are either scalar or of common size. error Display message and abort function. nan Not-a-Number.
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 %
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:
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