# Best way to represent a joint distribution

The motivation for this question is producing a distribution that produces the gender and age of an individual when the distribution of ages depends on gender.

Suppose I want a distribution which, when RandomVariate is applied to it, produces two values. The first value it produces is either 0 or 1 and the second value it produces is some real value between 0 and 65. But the second value is correlated with the first. By way of example, it might be the case that when the first value is 0, the second value tends to be lower than the when the first value is 1. I want the end product to be a distribution, so that one could apply things such as CDF to it.

Obviously, the gender distribution can be modeled as a BernoulliDistribution. And, as it happens I know how to model the age distribution conditioned on the person being a male and the age distribution conditioned on the person being a female. But how do I put them together? I've looked at CopulaDistribution and, while conceptually it might be appropriate, I don't see the kind of kernel I would need.

I have the feeling I am being stupid and that there is some obvious and elegant representation of this situation, but, at the moment, it escapes me. Help appreciated.

-
If the mixed age distribution is the only thing you're after you could perhaps use MixtureDistribution. but since you say you want both values that's probably not what you want. – Sjoerd C. de Vries Jan 24 '14 at 7:40
But for one thing, this would be a great candidate for the ProbabilityDistribution function, which lets you build a distribution from a specified distribution function. Thing is, combining discrete and continuous variables is not allowed (by ProbabilityDistribution, not in general). – Ian Jan 24 '14 at 14:43

I believe TransformedDistribution can do the job. I'll do the conditional part first and finish with producing the number pairs asked for.

The trick is to let the transformation function contain the conditional using Boole:

dist[μm_, σm_, μf_, σf_, p_] =
TransformedDistribution[
Boole[mf == 0] f + Boole[mf == 1]  m,
{
m \[Distributed] NormalDistribution[μm, σm],
f \[Distributed] NormalDistribution[μf, σf],
mf \[Distributed] BernoulliDistribution[p]
}
]


Let's try it:

Mean[dist[μm, σm, μf, σf, p]]
(* μf - p μf + p μm *)

StandardDeviation[dist[μm, σm, μf, σf, p]]
(* Sqrt[p μf^2 - p^2 μf^2 - 2 p μf μm + 2 p^2 μf μm +
p μm^2 - p^2 μm^2 + σf^2 - p σf^2 + p σm^2] *)

RandomVariate[dist[6, 1, 0, 1/2, 1/4], 10000] // Histogram


It also works for more complicated functions (if you're lucky):

dist2[am_, bm_, af_, bf_, p_] :=
TransformedDistribution[
Boole[mf == 0] f + Boole[mf == 1]  m,
{
m \[Distributed] GompertzMakehamDistribution[am, bm],
f \[Distributed] GompertzMakehamDistribution[af, bf],
mf \[Distributed] BernoulliDistribution[p]
}
]

Mean[dist2[am, bm, af, bf, p]]
(* (
am E^bf MeijerG[{{}, {1}}, {{0, 0}, {}}, bf] -
am E^bf p MeijerG[{{}, {1}}, {{0, 0}, {}}, bf] +
af E^bm p MeijerG[{{}, {1}}, {{0, 0}, {}}, bm])/(af am) *)


The above can be easily extended to generate the numbers pairs you need:

dist[μm_, σm_, μf_, σf_, p_] :=
TransformedDistribution[
{mf, Boole[mf == 0] f + Boole[mf == 1]  m},
{
m \[Distributed] NormalDistribution[μm, σm],
f \[Distributed] NormalDistribution[μf, σf],
mf \[Distributed] BernoulliDistribution[p]
}
]

RandomVariate[dist[65, 5, 75, 5, 1/4], 100000] // DensityHistogram


-
This will be handy. +1 – ciao Jan 24 '14 at 23:37
Thanks. This looks very good. I had not thought of the idea of having a list of values inside a TransformedDistribution. – Seth Chandler Jan 28 '14 at 1:06
@seth glad to be of help, though I'm not really sure this matches your requirements. You mentioned that a CopulaDistribution is not what you're looking for but the above solution is actually equivalent to a copula with a product kernel. – Sjoerd C. de Vries Jan 28 '14 at 7:07
@Sjoerd C. de Vries For my purposes, this does exactly what I need – Seth Chandler Jan 28 '14 at 16:04

I don't think there's any direct way to do this currently with the MM statistics/probability features - it would be nice if WRI extended things like ProductDistribution et. al. to allow conditioning/other logic. In any case, here's one way you could accomplish the end result:

dist = EmpiricalDistribution[{#,
If[# == 0, RandomReal[50], RandomReal[65]]} & /@
RandomVariate[BernoulliDistribution[.45], 100]]

RandomVariate[dist]

(*  {0, 18.2127}  *)


So in this trivialized example, I created a distribution where 45% are gender 1, 55% gender 0, and the 'age' is tacked on based on the gender. In your case, you'd replace the RandomReal with RandomVariate on your actual age vs gender distribution.

You can use pretty much all of the MM statistical/probability functions on the distribution:

Probability[gender == 1 && age > 40, {gender, age} \[Distributed] dist]

(*  0.08  *)


You might also sniff around at ParameterMixtureDistribution, it could be of use here.

-