Expressing a joint distribution for use in ContourPlot

Question

I have

$$X_1=T + W$$ $$X_2=2 T + 3 W$$

where

$$W \sim N(3,5) \,\text{;}\,\, T \sim N(1,2)$$ and $$P(X=X_1)=P(X=X_2)=\tfrac{1}{2}$$

which I have awkwardly expressed as

X2 = TransformedDistribution[
  2 T + 3 W, {W \[Distributed] NormalDistribution[3, Sqrt[5]], 
   T \[Distributed] NormalDistribution[1, Sqrt[2]]}]

X1 = TransformedDistribution[
  T + W, {W \[Distributed] NormalDistribution[3, Sqrt[5]], 
   T \[Distributed] NormalDistribution[1, Sqrt[2]]}]

XX = TransformedDistribution[
  i a + (1 - i) b, {a \[Distributed] X1, b \[Distributed] X2, 
   i \[Distributed] BernoulliDistribution[0.5]}]

This serves fine for generating RandomVariates and for calculating Mean, Variance, etc.; but I would like to plot the PDF of $T$ vs $X$ and am struggling with how to correctly construct the joint PDF to do this.

How should I go about constructing these distributions so that the joint PDF, $f_{X,T}(x,t)$, can be plotted using, for example with ContourPlot or Plot3D?

Answer 1

As a warm up, let's find the PDF of XX. Now if we just do PDF[XX, x], Mathematica will seemingly spin forever, which one might think is caused by Integrate. Let's see if that's correct.

First let's Print each time we Integrate:

Block[{Integrate},
  Integrate[e__] /; (Print[{e}]; False) := Null;
  PDF[XX, x]
]

Ok so there is one integral, and notice we can simplify it by realizing x3 is either 0 or 1 as seen in the Assumptions option above. Making our own custom integrator for this case we get an answer for the PDF of XX:

With[{X = \[FormalX]3},
  Internal`InheritedBlock[{Integrate},
    Unprotect[Integrate];
    Integrate[e__] /; !FreeQ[{e}, X] := 
      Piecewise[{{Integrate @@ ({e} /. X->0), X==0}, {Integrate @@ ({e} /. X->1), X==1}}];

    PDF[XX, x]
  ]
]
(* 0.0273995 E^(-(1/106) (-11+x)^2)+0.075393 E^(-(1/14) (-4+x)^2) *)

Plot[0.0273995 E^(-(1/106) (-11+x)^2)+0.075393 E^(-(1/14) (-4+x)^2), {x, -10, 30}]

Now I'm a little fuzzy on how to compute a joint PDF but perhaps it can be computed in a similar fashion. (Let me know if this is the wrong way to find the joint pdf!):

With[{X=\[FormalX]3},
  Internal`InheritedBlock[{Integrate},
    Unprotect[Integrate];
    Integrate[e__] /; !FreeQ[{e}, X] := 
      Piecewise[{{Integrate @@ ({e} /. X->0), X==0}, {Integrate @@ ({e} /. X->1), X==1}}];

    jcdf = Probability[
      X <= x && T <= t, 
      {X\[Distributed]XX, T\[Distributed]NormalDistribution[1,Sqrt[2]]}
    ];

    jpdf = D[jcdf, x, t]
  ]
]
(* 0.0705237 E^(-0.25 (-1.+t)^2) (0.109598 E^(-0.00943396 (-11.+x)^2)+0.301572 E^(-0.0714286 (-4.+x)^2)) *)

Plot3D[jpdf, {x, -10, 30}, {t, -10, 10}, PlotRange -> All, 
  PlotPoints -> 50, AxesLabel -> {x, t}]

Edit:

I feel this is wrong because Probability thinks X and T are independent... and I'm pretty sure they're dependent, is that correct?