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
Generalized Cholesky factorization | |
Cholesky factorization. | |
Display message and abort function. | |
Convert sparse matrix to full matrix. | |
True for sparse matrix. | |
Convert string to lowercase. | |
Number of nonzero matrix elements. | |
Convert number to string. (Fast version) | |
Matrix square root. | |
Singular value decomposition. | |
Find a few singular values and vectors. |
generates conditionally simulated values |
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
Comments or corrections to the WAFO group