%% Generate observations from Multivariate Normal distribution
function draw=multi_normrnd(mu, cov_matrix,sample)

[row,col]=size(mu);
if row>col
    mu=mu'; % make sure mu is a row vector
end
k=size(cov_matrix,1);
R=chol(cov_matrix); % upper-triangular matrix
draw=repmat(mu,[sample, 1])+randn(sample,k)*R;
