How create joint distribution of two dependent variables?

Hi,

I have defined two random variables, A and B, which are correlated. I have defined their respective (marginal) pdfs using lognpdf. I want to define their joint pdf. How do I do this?

Or, equivalently, if I define the pdfs of A and the joint for (A,B), can I back out the pdf for B?

mu_A = 0;

mu_B = 0;

sigma_A = sqrt(1);

sigma_B = sqrt(1);

sigma_AB = 0.2;

pdfA = @(a) lognpdf(a,par.mu_A,par.sigma_A);

pdfB = @(b) lognpdf(lb,par.mu_B,par.sigma_B);

Answer 1

[Jeff Miller]

I don't think your two questions are equivalent.

If you have the two marginal pdfs and the correlation, there are lots of ways to define the joint distribution--it is not uniquely determined by the marginals and correlation, except in the special case of the bivariate normal. So, you need to choose one of the possible joint distributions depending on your situation. For more information, read about "copulas".

On the other hand, if you have the joint distribution, you can compute everything from it. For example, you can integrate it across one variable to get the marginal distribution for the other variable. Similarly, you can compute the correlation. So, the joint distribution does uniquely determine everything else. Just not vice versa.

[ASKH]

Hi Jeff,

Thanks for that advice! One thing I could do, then, is define

pdfA = @(a) lognpdf(a,par.mu_A,par.sigma_A);

and

jointpdfAB = @(a,b) ...

and then integrate jointpdfAB over b. But how do I create a joint pdf of A and B? I.e. how do I generate a multivariate log normal distribution...?

[Jeff Miller]

As it says here, there are infinite ways to do that. None of them are very easy, though. Google "lognormal copulas" and you will see that this is not a simple problem. Here is some ugly code that might help you get started:

aMu = 1; bMu = 2;    % Means & SDs of the normals,
aSD = 1; bSD = 1.5;  %  not of the lognormals
rho = 0.7;           % Correlation of the normals, not the lognormals.

amin = 0.12;  % Ranges on the lognormal values, not the normals.
amax = 59;
bmin = 0.07;
bmax = 760;

global mu
global sigma
global intconstant

mu = [aMu, bMu];
sigma = [aSD^2, aSD*bSD*rho; ...
         aSD*bSD*rho, bSD^2];
intconstant = 1;  % A default value that will be recomputed immediately.
intconstant = 1 / integral2(@mvnlnpdf2,amin,amax,bmin,bmax)

% Make sure the integral of the pdf is 1.
newint = integral2(@mvnlnpdf2,amin,amax,bmin,bmax)

% Plot the pdf to see what it looks like:
[X,Y] = meshgrid(linspace(amin,amax,50),linspace(bmin,bmax,50));
Z = mvnlnpdf2(X,Y);
figure; surf(X,Y,Z)

function pdf = mvnlnpdf2(a,b)
    global mu
    global sigma
    global intconstant
    pdf = zeros(size(a));
    lna = log(a);
    lnb = log(b);
    for i=1:numel(a)
        pdf(i) = mvnpdf([lna(i),lnb(i)],mu,sigma) * intconstant;
    end
end

[ASKH]

Thanks Jeff!