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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

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

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