Best way to represent a joint distribution

Question

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.

Answer 1

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

Answer 2

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.