# Expressing a joint distribution for use in ContourPlot

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?

-
I don't know how probability distributions work in Mathematica. Does Mathematica know that the T and W in the definitions of X1 and X2 are the same? Otherwise you won't get the right correlations. For example, consider y = TransformedDistribution[x, x \[Distributed] NormalDistribution[]]; z = TransformedDistribution[x, x \[Distributed] NormalDistribution[]]; are y and z independent? –  Rahul Apr 23 '14 at 22:42
@RahulNarain: Does that matter here, where I'm just looking at $X$, which is either $X_1$ or $X_2$? –  raxacoricofallapatorius Apr 23 '14 at 23:00
Hmm... probably not, you're right. –  Rahul Apr 23 '14 at 23:05
@RahulNarain: Good catch though, for the general case (i not limited to $\{0, 1\}$). –  raxacoricofallapatorius Apr 23 '14 at 23:36

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},
InternalInheritedBlock[{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},
InternalInheritedBlock[{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?

-
My Mathematica formulation may not be the correct way to capture the mathematical formulation (which is definitive) at the start of the question. It works for some applications (see my comment) but not generally, and perhaps not here. –  raxacoricofallapatorius Apr 23 '14 at 23:38