Home > wafo > wstats > winvgrnd.m

# winvgrnd

## PURPOSE

Random matrices from a Inverse Gaussian distribution.

## SYNOPSIS

R = winvgrnd(m0,l,varargin);

## DESCRIPTION

``` WINVGRND Random matrices from a Inverse Gaussian distribution.

CALL:  R = winvgrnd(m0,l,sz);

m0,l = parameters  (see winvgpdf)
sz = size(R)    (Default common size of m0 and l)
sz can be a comma separated list or a vector
giving the size of R (see zeros for options).

Examples:
R = winvgrnd(2,2,100,2);
R2 = winvgrnd(2,3,[100,2]);
wqqplot(R(:,1),R2(:,1))

## CROSS-REFERENCE INFORMATION

This function calls:
 comnsize Check if all input arguments are either scalar or of common size. wchi2rnd Random matrices from a Chi squared distribution. error Display message and abort function. nan Not-a-Number.
This function is called by:

## SOURCE CODE

```001 function R = winvgrnd(m0,l,varargin);
002 %WINVGRND Random matrices from a Inverse Gaussian distribution.
003 %
004 % CALL:  R = winvgrnd(m0,l,sz);
005 %
006 %      m0,l = parameters  (see winvgpdf)
007 %        sz = size(R)    (Default common size of m0 and l)
008 %             sz can be a comma separated list or a vector
009 %             giving the size of R (see zeros for options).
010 %
011 % Examples:
012 %   R = winvgrnd(2,2,100,2);
013 %   R2 = winvgrnd(2,3,[100,2]);
014 %   wqqplot(R(:,1),R2(:,1))
015 %
017
018 % Reference: Chhikara & Folks, "The Inverse Gaussian Distribution", p. 53
019
020 % Tested on; Matlab 5.3
021 % History:
022 % revised pab 24.10.2000
023 %  - added comnsize, nargchk
025
026 error(nargchk(2,inf,nargin))
027 if nargin==2,
028   [errorcode m0 l] = comnsize(m0,l);
029 else
030   [errorcode m0 l] = comnsize(m0,l,zeros(varargin{:}));
031 end
032 if errorcode > 0
033   error('m0 and l must be of common size or scalar.');
034 end
035 R=zeros(size(m0));
036 ok=((m0>0)&(l>0));
037 k=find(ok);
038 if any(k)
039   ksiz=size(k);
040   R1=rand(ksiz);
041   Y=wchi2rnd(1,ksiz);
042   X1=m0(k)./(2*l(k)).*(2*l(k)+m0(k).*Y-(4*l(k).*m0(k).*Y+m0(k).^2.*Y.^2).^(1/2));
043   X2=m0(k).^2./X1;
044   I=(R1<m0(k)./(m0(k)+X1));
045   R(k)=X1.*I+X2.*(1-I);
046 end
047 k1=find(~ok);
048 if any(k1)
049   tmp=NaN;
050   R(k1)=tmp(ones(size(k1)));
051 end
052
053
054```

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