*Post by Serigne Lo*Hi Bill,

Thanks again for this function. It was very useful

however I

have one thing I don't understand well about the

convergence

of the mean. Here is waht I did

P=[1 4]; % set the parameters P and M

M=[1 1];

CorrMatrix=[1 .7;.7 1]; % the covariance matrix

n=1000000; % sample size

% simulate bivariate gamma distribution

R=corrgammarnd(P, M, CorrMatrix, n, 'sim');

mean(R)

corr(R)

MEAM = 0.8863 1.7742

CORR = 1.0000 0.6879

0.6879 1.0000

The theoritical means are P./M ie. [1 4] instead of

[0.8863

1.7742]. Though the sample size is enough large

1'000'000. Could you please explain me why the the

simulated

data didn't converge to the true mean.

I think the correlation is fine there is no big

difference

between 0.6879 and 0.7.

Thanks a lot for your help

Regards

SL

forum.org>...

*Post by Bill August**Post by Serigne Lo*Hi all

I would like to generate Bivariate Exponential

distribution

or Bivariate Gamma distribution with Matlab. Can

someone

help me to do this.

THANKS for your help

Regards

SLo

function R = corrgammarnd(P, M, CorrMatrix, n,

method)

*Post by Bill August*% generate correlected gamma variables

% P: mean value

% M: fading index

% CorrMatrix

% n: number of sample

% method: sim or decomposition

% R

% Bill August

% 31/03/2008

if nargin < 5

method = 'decomposition' ;

else

end

switch method

case 'sim'

N = length(P); % Number of variates

var = P.^2 ./ M; % Variance of variates

covar = sqrt(diag(var)) * CorrMatrix *

sqrt(diag(var)); % covariance matrix of variates.

*Post by Bill August*l = M ./ P; % lambda(i) = m(i)/P(i)

alpha = (1 ./ l') * l; % alpha(i,j) =

l(j) / l(i)

*Post by Bill August*a = zeros(size(covar));

Mg = zeros(1,N);

Mg(1) = M(1);

a(:,1) = (l(1) * l .* covar(1,:))'; %

a(i,1) =

lambda(1)*lambda(i)*C(1,i)

*Post by Bill August*for i = 2:N

for j = 2:i

s = 0;

for k = 1: j-1

s = s + a(i,k) * a(j,k) /

Mg(k);

*Post by Bill August*end;

a(i,j) = l(i) * l(j) * covar(i,j) -

s;

*Post by Bill August*s = 0;

for k = 1: i-1

s = s + a(i,k);

end;

Mg(i) = M(i) - s;

end;

end;

Pg = Mg * diag(1 ./ l); % Pg(i) =

Mg(i)./l(i)

*Post by Bill August*g = zeros(n,N);

for j = 1:N;

% Each column of G contains samples

of one

variate.

*Post by Bill August*g(:,j) = gamrnd(Mg(j),Pg(j)/Mg(j),[n

1]);

*Post by Bill August*end;

b = zeros(n,N,N);

s = zeros(n,1);

Gamma = zeros(n,N);

Gamma(:,1) = g(:,1);

for i = 2:N

s = zeros(n,1);

for j = 1:i-1

% b, Beta RVs container, is a 3D

array

with each (vertical) vector at

*Post by Bill August*% coordinations i,j (indicating

length and

width offset) contains samples

*Post by Bill August*% of a variate in its height. Note

that b

is computed for 1 <= j < i < N

betarnd(a(i,j),Mg(j)-a(i,j),[n

1]); % b[ij] samples

*Post by Bill August*if j < i

s = s + alpha(i,j) * b(:,i,j)

.* g(:,j);

*Post by Bill August*end;

end;

Gamma(:,i) = g(:,i) + s;

end;

R = sqrt(Gamma);

case 'decomposition'

% get the number of variable

VarNumber = length(P) ;

% calculate the std

Sigma = sqrt(P.^2 ./ M) ;

% transform the correlation matrix to

covariance

matrix

*Post by Bill August*CovMatrix = corr2cov(Sigma, CorrMatrix) ;

% Cholesky decomposition

L = chol(CovMatrix,'lower') ;

% calcualte Mw

Mw = zeros(1, VarNumber) ;

Mw(1) = M(1) ;

for i = 2 : length(P)

s = sum( L(i, 1:i-1) .* sqrt(Mw(1:i-1))

) ;

*Post by Bill August*Mw(i) = (s - P(i))^2 / (L(i, i)^2) ;

end

% calculate Pw

Pw = sqrt(Mw) ;

%

theta = Pw ./ Mw ;

k = Mw ;

W = zeros(VarNumber, n) ;

for num = 1 : VarNumber

W(num, :) = gamrnd(k(num), theta(num),

1, n) ;

*Post by Bill August*end

%

R = L*W ;

R = R' ;

otherwise

end

Hi Serigne,

I'm very sorry that i do not keep watching this thread.

I've checked the function and it is not the up-to-date version.:) Here it is...sorry again...

function R = corrgammarnd(P, M, CorrMatrix, n, method)

% Usage:

% generate correlected gamma variables

% Input:

% P

% M

% CorrMatrix

% n

% Output:

% R

% Author:

% Bill August

% Date:

% 31/03/2008

if nargin < 5

method = 'decomposition' ;

else

end

switch method

case 'sim'

N = length(P); % Number of variates

var = P.^2 ./ M; % Variance of variates

covar = sqrt(diag(var)) * CorrMatrix * sqrt(diag(var)); % covariance matrix of variates.

l = M ./ P; % lambda(i) = m(i)/P(i)

alpha = (1 ./ l') * l; % alpha(i,j) = l(j) / l(i)

a = zeros(size(covar));

Mg = zeros(1,N);

Mg(1) = M(1);

a(:,1) = (l(1) * l .* covar(1,:))'; % a(i,1) = lambda(1)*lambda(i)*C(1,i)

for i = 2:N

for j = 2:i

s = 0;

for k = 1: j-1

s = s + a(i,k) * a(j,k) / Mg(k);

end;

a(i,j) = l(i) * l(j) * covar(i,j) - s;

s = 0;

for k = 1: i-1

s = s + a(i,k);

end;

Mg(i) = M(i) - s;

end;

end;

Pg = Mg * diag(1 ./ l); % Pg(i) = Mg(i)./l(i)

g = zeros(n,N);

for j = 1:N;

% Each column of G contains samples of one variate.

g(:,j) = gamrnd(Mg(j),Pg(j)/Mg(j),[n 1]);

end;

b = zeros(n,N,N);

s = zeros(n,1);

Gamma = zeros(n,N);

Gamma(:,1) = g(:,1);

for i = 2:N

s = zeros(n,1);

for j = 1:i-1

% b, Beta RVs container, is a 3D array with each (vertical) vector at

% coordinations i,j (indicating length and width offset) contains samples

% of a variate in its height. Note that b is computed for 1 <= j < i < N

b(:,i,j) = betarnd(a(i,j),Mg(j)-a(i,j),[n 1]); % b[ij] samples

if j < i

s = s + alpha(i,j) * b(:,i,j) .* g(:,j);

end;

end;

Gamma(:,i) = g(:,i) + s;

end;

R = Gamma ;

case 'decomposition'

% get the number of variable

VarNumber = length(P) ;

% calculate the std

Sigma = sqrt(P.^2 ./ M) ;

% transform the correlation matrix to covariance matrix

CovMatrix = corr2cov(Sigma, CorrMatrix) ;

% Cholesky decomposition

L = chol(CovMatrix,'lower') ;

% calcualte Mw

Mw = zeros(1, VarNumber) ;

Mw(1) = M(1) ;

for i = 2 : length(P)

s = sum( L(i, 1:i-1) .* sqrt(Mw(1:i-1)) ) ;

Mw(i) = (s - P(i))^2 / (L(i, i)^2) ;

end

% calculate Pw

Pw = sqrt(Mw) ;

%

theta = Pw ./ Mw ;

k = Mw ;

W = zeros(VarNumber, n) ;

for num = 1 : VarNumber

W(num, :) = gamrnd(k(num), theta(num), 1, n) ;

end

%

R = L*W ;

R = R' ;

otherwise

end