Home > wafo > wstats > wgamrnd.m

# wgamrnd

## PURPOSE

Random matrices from a Gamma distribution.

## SYNOPSIS

R = wgamrnd(a,b,varargin);

## DESCRIPTION

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

## 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:
 dist2drnd Random points from a bivariate DIST2D distribution mdist2drnd Random points from a bivariate MDIST2D distribution wchi2rnd Random matrices from a Chi squared distribution.

## SOURCE CODE

```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 %
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```

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