# TransformedDistribution with Conditioned

Is the following attempt beyond Mathematica 11?

Z = TransformedDistribution[ (A + B)/2 \[Conditioned] A < B, {A \[Distributed] NormalDistribution[mA , sA], B \[Distributed] NormalDistribution[mB , sB]}]


When I try to get Mathematica to show me the PDF of Z, it doesn't work. I tried:

 PDF[Z, y]

• Would that not work: Z = TransformedDistribution[ (A + B)/2 , {A [Distributed] NormalDistribution[mA , sA], B [Distributed] NormalDistribution[mB , sB]},Assumptions -> A < B ] ? – amator2357 Apr 19 at 21:10
• @amator2357 The documentation says that the Assumptions are for parameters rather than the random variables. – JimB Apr 19 at 21:19
• OK, yeah, I thought that might have been the case, haven't read through the documentation properly, thanks for the heads up @JimB – amator2357 Apr 19 at 21:21
• @amator2357. Join the club. – JimB Apr 19 at 21:23
• So, my code is fine but Mathematica 11 stumbles? – user120911 Apr 20 at 1:42

It is possible to derive an exact solution to this problem.

Given: $$X$$ and $$Y$$ are independent random variables where $$X \sim N(\mu_1, \sigma_1^2)$$ and $$Y \sim N(\mu_2, \sigma_2^2)$$, with parameter conditions:

Problem: Find the pdf of $$\frac{X+Y}{2} \; \big| \; X < Y$$

1. Joint pdf of $$(X,Y)$$:

By independence, the joint pdf of $$(X,Y)$$, say $$f(x,y)$$ is simply the product of the individual pdf's:

1. Let $$V = X - Y$$. Then $$V \sim N(\mu_1 - \mu_2, \sigma_1^2 + \sigma_2^2)$$ with cdf $$\Phi(v)$$.

Let constant $$c = P(X which is: (take care here with non-standard Mma notation)

1. Conditional joint pdf:

The conditional pdf $$f\big((x,y) \; \big| \; X is then fcon:

where all the dependence is captured within the fcon statement using the Boole statement, and we can enter the 'domain' as a rectangular structure on the real line, i.e.

domain[fcon] = domain[f]

1. Transformation $$Z = \frac{X+Y}{2}$$

Given the conditional joint pdf $$f\big((x,y) \; \big| \; X ... let $$Z = \frac{X+Y}{2}$$ and $$W = X$$. Then the joint conditional pdf of $$(Z,W)$$, say $$g(z,w)$$, is obtained with:

where I am using the Transform function from the mathStatica package for Mathematica, and the domain can again be entered as a rectangular set as:

Then, the marginal pdf of $$Z = \frac{X+Y}{2}$$ is:

... which is the exact solution. All done.

# Monte Carlo check

The following plot compares:

• the exact symbolic pdf derived above (red dashed curve)

• ... to the Monte Carlo simulated pdf (squiggly blue curve)

... here when: $$\mu_1 = -1, \mu_2 = 4, \sigma_1 = 1, \sigma_2 = 12$$

Looks fine.

• +1 Very nice!!! – JimB Apr 21 at 4:35
• @Wolfies, how do I place the two normal distributions and your solution on the same axis? I am trying to use Manipulate to examine your solution together with the normals. I will gladly open another question for this, if needed. – user120911 Apr 21 at 20:51
• When I place your solution and the normals on different axis, the scale of the normals and your solution seem off. I suspect they may not be directly comparable. Is that correct? – user120911 Apr 21 at 21:42
• Also, can mathStatica find the expectation of your solution? Mathematica cannot, it would appear. – user120911 Apr 21 at 22:32
• You may have to resort to NIntegrate to find moments, given numerical parameter values. – wolfies Apr 22 at 14:41

This is not the answer you desired but here is an approach to get the cdf and pdf using numerical integration.

Proportion of the time that x1 < x2 given than x1 and x2 are independent (this took 80 seconds):

int0 = Integrate[PDF[NormalDistribution[μ1, σ1], x1] PDF[NormalDistribution[μ2, σ2], x2],
{x2, -∞, ∞}, {x1, -∞, x2}, Assumptions -> {σ1 > 0, σ2 > 0}]

(* Integrate[(E^(-((x2 - μ2)^2/(2 σ2^2))) (1 + Erf[(x2 - μ1)/(Sqrt[2] σ1)]))/(2 Sqrt[2 π] σ2),
{x2, -∞, ∞}, Assumptions -> {σ1 > 0, σ2 > 0}] *)


Proportion of the time that (x1+x2)/2 < t and x1 < x2 given that x1 and x2 are independent (this took 100 minutes):

int1 = Integrate[PDF[NormalDistribution[μ1, σ1], x1] PDF[NormalDistribution[μ2, σ2], x2],
{x2, -∞, ∞}, {x1, -∞, Min[x2, 2 t - x2]}, Assumptions -> {σ1 > 0, σ2 > 0}]

(* Integrate[(E^(-((x2-μ2)^2/(2 σ2^2)))Erfc[(μ1 - Min[2 t - x2, x2])/(Sqrt[2] σ1)])/(2 Sqrt[2 π] σ2),
{x2, -∞, ∞}, Assumptions -> {σ1 > 0, σ2 > 0}] *)


So we can define a cdf using numerical integration:

cdf[t_, μ1_, μ2_, σ1_, σ2_] := NIntegrate[(
E^(-((x2 - μ2)^2/(2 σ2^2))) Erfc[(μ1 - Min[2 t - x2, x2])/(Sqrt[2] σ1)])/(2 Sqrt[2 π] σ2),
{x2, -∞, ∞}]/
NIntegrate[(E^(-((x2 - μ2)^2/(2 σ2^2))) (1 + Erf[(x2 - μ1)/(Sqrt[2] σ1)]))/(2 Sqrt[2 π] σ2),
{x2, -∞, ∞}]


For the pdf we can differentiate the part of the cdf that depends on t (this took 52 minutes):

FullSimplify[D[Integrate[(E^(-((x2 - μ2)^2/(2 σ2^2))) Erfc[(μ1 - Min[2 t - x2, x2])/(Sqrt[2] σ1)])/
(2 Sqrt[2 π] σ2), {x2, -∞, ∞}, Assumptions -> {σ1 > 0, σ2 > 0}], t]]

(* (E^(-((-2 t + μ1 + μ2)^2/(2 (σ1^2 + σ2^2)))) Erfc[((t - μ2) σ1^2 + (-t + μ1) σ2^2)/
(Sqrt[2] σ1 σ2 Sqrt[σ1^2 + σ2^2])])/(Sqrt[2 π] Sqrt[σ1^2 + σ2^2]) *)


We can now define a pdf function:

pdf[t_, μ1_, μ2_, σ1_, σ2_] := ((E^(-((-2 t + μ1 + μ2)^2/(2 (σ1^2 + σ2^2))))
Erfc[((t - μ2) σ1^2 + (-t + μ1) σ2^2)/(Sqrt[2] σ1 σ2 Sqrt[σ1^2 + 2^2])])/
Sqrt[2 π] Sqrt[σ1^2 + σ2^2]))/
NIntegrate[(E^(-((x2 - μ2)^2/(2 σ2^2))) (1 + Erf[(x2 - μ1)/(Sqrt[2] σ1)]))/
(2 Sqrt[2 π] σ2), {x2, -∞, ∞}]


Here's a test example:

(* Generate a random sample *)
n = 1000000;
SeedRandom[12345];
z = RandomVariate[BinormalDistribution[{0, -1}, {1, 6}, 0], n];
z = Select[z, #[[1]] < #[[2]] &];
z = Total[#]/2 & /@ z;

(* Plot pdf's and cdf's *)
skd = SmoothKernelDistribution[z];
Plot[{PDF[skd, t], pdf[t, 0, -1, 1, 6]}, {t, Min[z], Max[z]},
PlotStyle -> {{LightGray, Thickness[0.03]}, {Red, Thickness[0.001]}},
PlotLegends -> {"Simulations", "Numerical integration"}]

Plot[{CDF[skd, t], cdf[t, 0, -1, 1, 6]}, {t, Min[z], Max[z]},
PlotStyle -> {{LightGray, Thickness[0.03]}, {Red, Thickness[0.001]}},
PlotLegends -> {"Simulations", "Numerical integration"}]


• How long did int1 take to compute on your system? – user120911 Apr 20 at 9:08
• A long, long time. I'll run it again later today and list the timing. I also let some code run to compute the pdf. I'll add that later today. – JimB Apr 20 at 14:15