Home > wafo > wstats > wmnormrnd.m

wmnormrnd

PURPOSE ^

Random vectors from a multivariate Normal distribution

SYNOPSIS ^

r = wmnormrnd(mu,sa,cases,method,cutoff);

DESCRIPTION ^

 WMNORMRND Random vectors from a multivariate Normal distribution 
  
  CALL:  r = wmnormrnd(mu,S,n,method) 
  
        r = matrix of random numbers from the multivariate normal 
            distribution with mean mu and covariance matrix S. 
        n = number of sample vectors 
   method = 'svd'  Singular value decomp.  (stable, quite fast) (default) 
            'chol'  Cholesky decomposition (fast, but unstable)  
            'sqrtm'  sqrtm                 (stable and slow)  
            'genchol' 
    
    S must be a symmetric, semi-positive definite matrix with size equal 
    to the length of MU. METHOD used for calculating the square root of S  
    is either svd, cholesky or sqrtm. (cholesky is fastest but least accurate.) 
    When cholesky is chosen and S is not positive definite, the svd-method  
    is used instead. 
  
  Example 
    mu = [0 5]; S = [1 0.45; 0.45 0.25]; 
    r = wmnormrnd(mu,S,100)'; 
    plot(r(:,1),r(:,2),'.') 
  
    d = 40; rho = 2*rand(1,d)-1; 
    mu = zeros(0,d); 
    S = (rho.'*rho-diag(rho.^2))+eye(d); 
    r = wmnormrnd(mu,S,100,'genchol')';  
  
  See also  chol, svd, sqrtm

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function r = wmnormrnd(mu,sa,cases,method,cutoff); 
002 %WMNORMRND Random vectors from a multivariate Normal distribution 
003 % 
004 % CALL:  r = wmnormrnd(mu,S,n,method) 
005 % 
006 %       r = matrix of random numbers from the multivariate normal 
007 %           distribution with mean mu and covariance matrix S. 
008 %       n = number of sample vectors 
009 %  method = 'svd'  Singular value decomp.  (stable, quite fast) (default) 
010 %           'chol'  Cholesky decomposition (fast, but unstable)  
011 %           'sqrtm'  sqrtm                 (stable and slow)  
012 %           'genchol' 
013 %   
014 %   S must be a symmetric, semi-positive definite matrix with size equal 
015 %   to the length of MU. METHOD used for calculating the square root of S  
016 %   is either svd, cholesky or sqrtm. (cholesky is fastest but least accurate.) 
017 %   When cholesky is chosen and S is not positive definite, the svd-method  
018 %   is used instead. 
019 % 
020 % Example 
021 %   mu = [0 5]; S = [1 0.45; 0.45 0.25]; 
022 %   r = wmnormrnd(mu,S,100)'; 
023 %   plot(r(:,1),r(:,2),'.') 
024 % 
025 %   d = 40; rho = 2*rand(1,d)-1; 
026 %   mu = zeros(0,d); 
027 %   S = (rho.'*rho-diag(rho.^2))+eye(d); 
028 %   r = wmnormrnd(mu,S,100,'genchol')';  
029 % 
030 % See also  chol, svd, sqrtm 
031  
032 % cutoff = cut off value for eigenvalues 
033  
034 % History 
035 % Revised pab 22.11.2003 
036 % -added option genchol   
037 % revised jr 19.09.2001 
038 % - Changed *name* in function call.  
039 % - Fixed a bug: the spurious default option 'swd' -> 'svd' 
040 % revised jr 18.06.2001 
041 % - Changed default from 'cholesky' to 'svd' 
042 % revised pab 17.01.2000  
043 % - updated documentation 
044 %  by Per A. Brodtkorb 17.09.98,20.08.98 
045  
046 [m n] = size(sa); 
047 if m ~= n 
048   error('S must be square'); 
049 end 
050 if  nargin<1 | isempty(mu); 
051   mu=zeros(m,1); 
052 end 
053  
054 musiz = size(mu); 
055 rows = max(musiz); 
056 if prod(musiz) ~= rows 
057    error('Mu must be a vector.'); 
058 end 
059 mu=mu(:); % make sure it is a column vector 
060  
061  
062 if m ~= rows 
063   error('The length of mu must equal the number of rows in S.'); 
064 end 
065 if nargin<5|isempty(cutoff),  
066   cutoff=0; 
067 else 
068   cutoff=abs(cutoff); % make sure cutoff is positive 
069 end 
070  
071 if nargin<4|isempty(method),  
072   method='svd'; % default method 
073 end 
074  
075 switch lower(method(1:2)), % 
076   case 'sq', % SQRTM slow but stable 
077     T=sqrtm(full(sa));%squareroot of S 
078      
079   case 'sv', % SVD stable quite fast 
080     if issparse(sa) % approximate method 
081       nz=nnz(sa); 
082       K=floor(min([10 3*sqrt(n) n/4 4*nz/n])); 
083       %[U S V]=svds(sa,floor(p/2)); 
084       [U S V]=svds(sa,K,'L'); 
085     else % exact method 
086       [U S V]=svd(full(sa),0); 
087     end 
088     L = n; 
089     if cutoff>0 
090       L=sum((diag(S)>=cutoff)); 
091       disp(['cutting off ' num2str(L) ' eigenvalues']) 
092     end 
093      
094     T=(U(:,1:L)*sqrt(S(1:L,1:L))*V(:,1:L)'); %squareroot of S 
095      
096   case 'ch', % CHOL not stable , but fast. 
097     [T, p] = chol(sa); 
098     if p ~= 0 
099       disp('S is not a positive definite matrix.'); 
100       disp('Cholesky factorization failed; trying svd instead.'); 
101       [U S V]=svd(full(sa),0); 
102       T=(U*sqrt(S)*V'); %squareroot of S 
103     end 
104  case 'ge', % GENCHOL stable 
105   [T,p,rank1] = genchol(sa,eps); 
106   r = T*randn(rows,cases)  + mu(p,ones(1,cases)); 
107   [tmp, ix] = sort(p); 
108   r = r(ix,:); 
109   return 
110  otherwise, error('unknown option') 
111 end 
112  
113  
114 r = T*randn(rows,cases)  + mu(:,ones(1,cases)); 
115  
116

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