R <-> Matlab, generating correlating normal variates

[fuzzyswat]

R <-> Matlab, generating correlated normal variates

HI, I need to generate correlated normal variates, with certain mean and sd.
I don't have Statistics Toolbox, so I have to code it myself.
I have a running R code from:
http://www.unt.edu/benchmarks/archives/2001/december01/rss.htm
but, cannot translate it in Matlab.

R code:

# n=sample size; p=number of columns
# u=mean of columns; s=standard deviation
# of columns; S=correlation matrix

mvrnorm <- function(n, p, u, s, S) {
 
	Z <- matrix(rnorm(n*p), p, n)
	t(u + s*t(chol(S))%*% Z)
}

# Sample from Bivariate Normal - R
z<-mvrnorm(n=100, p=2, u=c(100,100), s=c(15,15), S=matrix(c(1, .4, .4, 1), ncol=2, nrow=2,byrow=T))

it works fine, but....

MATLAB version:

%   function [X,Y] = corrNormal(Npar, Nobs, mu1, mu2, sd1, sd2, rho) 
% Input:
%   Npar - number of correlated parameters
%   Nobs - number of observations
%   rho - correlation between two parameters
% Output:
%   M - correlated normal distributed numbers
 
function M = corrNormal(Npar, Nobs, mu1, mu2, sd1, sd2, rho) 
 
mu = [mu1,mu2];
sd = [sd1,sd2];
 
% create uncorrelated observations
X = randn(Npar, Nobs);
 
% compute cholesky factor
S = [1 rho; rho 1];
cF = chol(S);
 
% R code: Y = t(u + s*t(chol(S))%*% Z)
M = mu + sd*cF'*X;
 
end

# Sample from Bivariate Normal - MATLAB

M = corrNormal(2, 100, 10, 10, .5, .5, 0.7);

produces an error

Error using  + 
Matrix dimensions must agree.

Error in corrNormal (line 26)
Y = mu + sd*cF'*X;

The error point is the product:
sd*cF'

it produces in MATLAB a vector 1x2, whereas R produces a matrix 2x2.

Any comments on the vector matrix product differences would be appreciated,
MJ

PS: I don't have the Statistics Toolbox. M

 

[fuzzyswat]

Re: R <-> Matlab, generating correlated normal variates

I found the solution. Thanks to the author of
http://comisef.wikidot.com/tutorial:correlation

one needs to extend the algorithm a bit, because it generates 'correlated normal random numbers' with mean=0 and sd=1.
what you have to do is to multiply the result with desired SDs and add the desired MEANs, so simple as that...
MJ

clear all; clc
close all;

Npar = 2; % number of correlated parameters
Nobs = 100; % number of observations
mu1 = 8; mu2 = 0.13; % the means and sd's
sd1 = .2; sd2 = .2;
rho = 0.8; % correlation between two parameters

% create uncorrelated observations
X = randn(Nobs, Npar);

% check X
figure
plotmatrix(X); corrcoef(X)

% set correlation matrix
M = ones(Npar, Npar) * rho;
M(1:(Npar + 1):(Npar * Npar)) = 1;

% compute cholesky factor
cF = chol(M);

% induce correlation, check
Y = X * cF
figure
plotmatrix(Y); corrcoef(Y)

% transform with means and sd's, check
Y1 = mu1.*ones(Nobs,1) + Y(:,1)*sd1
Y2 = mu2.*ones(Nobs,1) + Y(:,2)*sd2
YY = [Y1,Y2]

figure
plotmatrix(YY); corrcoef(YY)