Discussion:
Simulation of Bivariate Exponential / Gamma distribution
(too old to reply)
Serigne Lo
2008-08-15 03:09:02 UTC
Permalink
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
Bill August
2008-08-15 07:30:41 UTC
Permalink
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)
% Usage:
% generate correlected gamma variables
% Input:
% P: mean value
% M: fading index
% CorrMatrix
% n: number of sample
% method: sim or decomposition
% 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 = 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
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
Serigne Lo
2008-08-15 09:12:01 UTC
Permalink
Thanks a lot Bill for the function,
I'll test it over the we and will come back to you later.
Cheers
Serigne
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)
% 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)
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);
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.
Post by Bill August
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
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
Post by Bill August
b(:,i,j) = 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);
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)) ) ;
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
Serigne Lo
2008-08-26 06:12:03 UTC
Permalink
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)

and here are the results :

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
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)
% 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)
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);
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.
Post by Bill August
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
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
Post by Bill August
b(:,i,j) = 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);
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)) ) ;
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
Bill August
2008-08-28 17:06:51 UTC
Permalink
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
Post by Bill August
b(:,i,j) =
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
Peter Perkins
2008-08-15 13:47:46 UTC
Permalink
Post by Serigne Lo
Hi all
I would like to generate Bivariate Exponential distribution
or Bivariate Gamma distribution with Matlab.
My understanding of these is that

1) there are various constructions that work in the bivariate case
2) these become unmanageable very quickly for the N-D multivariate case
3) they can be fairly restrictive in the marginal distributions and in the correlation structure

If you have access to the Statistics Toolbox, you might consider using copulas. For example:

u = copularnd('Gumbel',3,1000);
g = gaminv(u,repmat([1,2],size(u,1),1));
scatterhist(g(:,1),g(:,2),[25,25]);

creates 1000 bivariate vectors, with a gamma(1,1) distribution for the first column, a gamma(2,1) distribution for the second, and a Spearman's rank correlation between those components of about 0.85.

Hope this helps.
Serigne Lo
2008-08-16 06:52:01 UTC
Permalink
Hi Peter,
Thanks a lot for your answer.
Unfortunately I my Matlab Statistical Toolbox is not
complete. I don't have "scatterhist" function
Regards
Serigne
Post by Peter Perkins
Post by Serigne Lo
Hi all
I would like to generate Bivariate Exponential distribution
or Bivariate Gamma distribution with Matlab.
My understanding of these is that
1) there are various constructions that work in the
bivariate case
Post by Peter Perkins
2) these become unmanageable very quickly for the N-D
multivariate case
Post by Peter Perkins
3) they can be fairly restrictive in the marginal
distributions and in the correlation structure
Post by Peter Perkins
If you have access to the Statistics Toolbox, you might
u = copularnd('Gumbel',3,1000);
g = gaminv(u,repmat([1,2],size(u,1),1));
scatterhist(g(:,1),g(:,2),[25,25]);
creates 1000 bivariate vectors, with a gamma(1,1)
distribution for the first column, a gamma(2,1) distribution
for the second, and a Spearman's rank correlation between
those components of about 0.85.
Post by Peter Perkins
Hope this helps.
a***@gmail.com
2018-03-01 17:43:10 UTC
Permalink
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
Hi.... can some one help me what is the expression for bi-variate exponential distribution given two random variables having exponential distributions with means u and v?
Loading...