3
$\begingroup$

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

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:

enter image description here

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:

enter image description here

  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<Y) = P(V<0) = \Phi(0)$ which is: (take care here with non-standard Mma notation)

enter image description here

  1. Conditional joint pdf:

The conditional pdf $f\big((x,y) \; \big| \; X<Y\big) = \frac{f(x,y)}{P(X<Y)}$ is then fcon:

enter image description here

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<Y\big)$ ... let $Z = \frac{X+Y}{2}$ and $W = X$. Then the joint conditional pdf of $(Z,W)$, say $g(z,w)$, is obtained with:

enter image description here

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:

enter image description here

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

enter image description here

... 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.

enter image description here

$\endgroup$
  • 1
    $\begingroup$ +1 Very nice!!! $\endgroup$ – JimB Apr 21 at 4:35
  • $\begingroup$ @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. $\endgroup$ – user120911 Apr 21 at 20:51
  • $\begingroup$ 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? $\endgroup$ – user120911 Apr 21 at 21:42
  • $\begingroup$ Also, can mathStatica find the expectation of your solution? Mathematica cannot, it would appear. $\endgroup$ – user120911 Apr 21 at 22:32
  • $\begingroup$ You may have to resort to NIntegrate to find moments, given numerical parameter values. $\endgroup$ – wolfies Apr 22 at 14:41
2
$\begingroup$

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"}]

pdf and cdf

$\endgroup$
  • $\begingroup$ How long did int1 take to compute on your system? $\endgroup$ – user120911 Apr 20 at 9:08
  • 1
    $\begingroup$ 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. $\endgroup$ – JimB Apr 20 at 14:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.