mvnrnd_fast.asv

function [r] = mvnrnd_fast(mu,sigma)
%MVNRND Random vectors from the multivariate normal distribution.
%   R = MVNRND(MU,SIGMA) returns an N-by-D matrix R of random vectors
%   chosen from the multivariate normal distribution with mean vector MU,
%   and covariance matrix SIGMA.  MU is an N-by-D matrix, and MVNRND
%   generates each row of R using the corresponding row of MU.  SIGMA is a
%   D-by-D symmetric positive semi-definite matrix, or a D-by-D-by-N array.
%   If SIGMA is an array, MVNRND generates each row of R using the
%   corresponding page of SIGMA, i.e., MVNRND computes R(I,:) using MU(I,:)
%   and SIGMA(:,:,I).  If the covariance matrix is diagonal, containing
%   variances along the diagonal and zero covariances off the diagonal,
%   SIGMA may also be specified as a 1-by-D matrix or a 1-by-D-by-N array,
%   containing just the diagonal. If MU is a 1-by-D vector, MVNRND
%   replicates it to match the trailing dimension of SIGMA.
%
%   R = MVNRND(MU,SIGMA,N) returns a N-by-D matrix R of random vectors
%   chosen from the multivariate normal distribution with 1-by-D mean
%   vector MU, and D-by-D covariance matrix SIGMA.
%
%   Example:
%
%      mu = [1 -1]; Sigma = [.9 .4; .4 .3];
%      r = mvnrnd(mu, Sigma, 500);
%      plot(r(:,1),r(:,2),'.');
%
%   See also MVTRND, MVNPDF, MVNCDF, NORMRND.

%   R = MVNRND(MU,SIGMA,N,T) supplies the Cholesky factor T of
%   SIGMA, so that SIGMA(:,:,J) == T(:,:,J)'*T(:,:,J) if SIGMA is a 3D array or SIGMA
%   == T'*T if SIGMA is a matrix.  No error checking is done on T.
%
%   [R,T] = MVNRND(...) returns the Cholesky factor T, so it can be
%   re-used to make later calls more efficient, although there are greater
%   efficiency gains when SIGMA can be specified as a diagonal instead.

%   Copyright 1993-2020 The MathWorks, Inc.

[n,d] = size(mu);
sz = size(sigma);
sz(1) = sz(2);
sigmaIsDiag = true;

% Special case: if mu is a column vector, then use sigma to try
% to interpret it as a row vector.
if d == 1 && sz(1) == n
    mu = mu';
    [n,d] = size(mu);
end

% Get size of data.
if nargin < 3 || isempty(cases)
    nocases = true; % cases not supplied
else
    nocases = false; % cases was supplied
    if n == cases
        % mu is ok
    elseif n == 1 % mu is a single row, make cases copies
        n = cases;
        mu = repmat(mu,n,1);
    else
        error(message('stats:mvnrnd:InputSizeMismatchMu'));
    end
end

outtype = internal.stats.dominantType(mu, sigma); % single if mu or sigma is

% Single covariance matrix
if ismatrix(sigma)
    % Make sure sigma is the right size
    
    if sz(1) ~= sz(2)
        error(message('stats:mvnrnd:BadCovariance2DSize'));
    elseif ~isequal(sz, [d d])
        error(message('stats:mvnrnd:InputSizeMismatch2DSigma'));
    end
    
        % Just the diagonal of sigma has been specified.
        t = sqrt(sigma);
        r = randn(n,d,'like',outtype).*t + mu;
        r(:,t==0) = mu(:,t==0); % force exact mean when variance is 0
    else
        % Factor sigma using a function that will perform a Cholesky-like
        % factorization as long as the sigma matrix is positive
        % semi-definite (can have perfect correlation). Cholesky requires a
        % positive definite matrix.  sigma == T'*T
        [T,err] = cholcov(sigma);
        if isnan(err)
            error(message('stats:mvnrnd:BadCovariance2DSym'));
        elseif err ~= 0
            error(message('stats:mvnrnd:BadCovariance2DPos'));
        end
        r = randn(n,size(T,1),'like',outtype) * T + mu;
        t = diag(sigma);
        r(:,t==0) = mu(:,t==0); % force exact mean when variance is 0
    end
    
% Multiple covariance matrices
elseif ndims(sigma) == 3
    % mu is a single row and cases not given, rep mu out to match sigma
    if n == 1 && nocases % already know size(sigma,3) > 1
        n = sz(3);
        mu = repmat(mu,n,1);
    end
    
    % Make sure sigma is the right size
    if sz(1) ~= sz(2) % Sigma is 3-D
        error(message('stats:mvnrnd:BadCovariance3DSize'));
    elseif (sz(1) ~= d) || (sz(2) ~= d) % Sigma is 3-D
        error(message('stats:mvnrnd:InputSizeMismatch3DSigma'));
    elseif sz(3) ~= n
        error(message('stats:mvnrnd:InputSizeMismatchSigmaDimension'));
    end
  
    if nargin < 4
        if nargout > 1
            T = zeros(sz,'like',outtype);
        end
        if sigmaIsDiag
            sigma = reshape(sigma,sz(2),sz(3))';
            if any(any(sigma<0))
                error(message('stats:mvnrnd:BadDiagSigma'));
            end
            R = sqrt(sigma);
            r = randn(n,d,'like',outtype).*R + mu;
            if nargout > 1
                for i=1:n
                    t = R(i,:);
                    T(:,:,i) = diag(t);
                    r(i,t==0) = mu(i,t==0); % force exact mean
                end
            end
        else
            r = zeros(n,d,'like',outtype);
            for i = 1:n
                [R,err] = cholcov(sigma(:,:,i));                
                if isnan(err)
                    error(message('stats:mvnrnd:BadCovariance3DSym'));
                elseif err ~= 0
                    error(message('stats:mvnrnd:BadCovariance3DPos'));
                end
                Rrows = size(R,1);
                r(i,:) = randn(1,Rrows,'like',outtype) * R + mu(i,:);
                t = diag(sigma(:,:,i));
                r(i,t==0) = mu(i,t==0); % force exact mean
                if nargout > 1
                    T(1:Rrows,:,i) = R;
                end
            end
        end
    else
        % T specified
        r = zeros(n,d,'like',outtype);
        for i = 1:n
            r(i,:) = randn(1,d,'like',outtype) * T(:,:,i) + mu(i,:);
        end
    end
end