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

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')';

## CROSS-REFERENCE INFORMATION

This function calls:
 genchol Generalized Cholesky factorization chol Cholesky factorization. error Display message and abort function. full Convert sparse matrix to full matrix. issparse True for sparse matrix. lower Convert string to lowercase. nnz Number of nonzero matrix elements. num2str Convert number to string. (Fast version) sqrtm Matrix square root. svd Singular value decomposition. svds Find a few singular values and vectors.
This function is called by:
 cov2csdat generates conditionally simulated values

## 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
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 %
031
032 % cutoff = cut off value for eigenvalues
033
034 % History
035 % Revised pab 22.11.2003
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