7
$\begingroup$

I am new to Mathematica, and this question may reflect that.

I am trying to plot the probability density function of a random variable which equals the difference between two independent random variables that are Half Normal distributed.

Can someone help me out?

$\endgroup$
  • $\begingroup$ If I made the wrong changes to your question, please change it back. $\endgroup$ – JimB Jan 5 '17 at 22:46
  • $\begingroup$ Should be able to do PDF[TransformedDistribution[ x - y, {x \[Distributed] HalfNormalDistribution[\[Theta]], y \[Distributed] HalfNormalDistribution[\[Theta]]}], z], but it's taking a long time and returning unevaluated. $\endgroup$ – march Jan 5 '17 at 23:25
  • 1
    $\begingroup$ @march: PDF and CDF return unevaluated (even with values assigned to $\theta$) but Expectation and RandomVariate work on the TransformedDistribution. Maybe one needs to brute force the integration of the joint probability distribution or take a large random sample and use SmoothKernelDistribution. $\endgroup$ – JimB Jan 5 '17 at 23:30
  • $\begingroup$ @JimBaldwin. I wasn't aware of SmoothKernelDistribution. It ends up kind of nice: plotting the pdf of SmoothKernelDistribution@RandomVariate[dist /. \[Theta] -> 1, {1000000}] looks pretty good. I feel like an answer is in order there. $\endgroup$ – march Jan 5 '17 at 23:46
  • 1
    $\begingroup$ @JimBaldwin : "...Maybe one needs to brute force the integration of the joint probability distribution..." - that's gonna be ugly... $\endgroup$ – ciao Jan 6 '17 at 0:34
9
$\begingroup$

Hmmm. Maybe it's not so hard (with Mathematica) after all. If the two independent and identically distributed half-normal random variables are labeled $u$ and $v$, then the difference is $d=u-v$. If we can find the CDF of $d$, then we can differentiate to get the density of $d$.

Suppose $d_0 \ge 0$. Then

$$Pr(d\le d_0)=Pr(u-v\le d_0)=\int_0^\infty{\int_0^{d_0+u}(\frac{2 \theta e^{-\frac{\theta ^2 v^2}{\pi }}}{\pi })(\frac{2 \theta e^{-\frac{\theta ^2 u^2}{\pi }}}{\pi })}dv du$$

So we can get the density with

D[Integrate[PDF[HalfNormalDistribution[θ], u] PDF[HalfNormalDistribution[θ], v],
  {u, 0, ∞}, {v, 0, d + u}, Assumptions -> {θ > 0, d >= 0}], d]

(* (Sqrt[2] E^(-((d^2 θ^2)/(2 π))) θ Erfc[(d θ)/Sqrt[2 π]])/π *)

While the next step with $d_0<0$ can be made more justifiably, a bit of handwaving notes that the distribution is symmetric about 0 when the two half-normal distributions have the same parameter. That results in the density function of the difference being

(Sqrt[2] E^(-((d^2 θ^2)/(2 π))) θ (1 + Erf[(Abs[d] θ)/Sqrt[2 π]]))/π

A plot of this function against a huge random sample from the same distribution does not suggest the result is wrong. So borrowing from @march we have the following:

dist = TransformedDistribution[x - y, {x \[Distributed] HalfNormalDistribution[θ], 
 y \[Distributed] HalfNormalDistribution[θ]}]
sample = RandomVariate[dist /. θ -> 1, {1000000}];
f[z_] = PDF[SmoothKernelDistribution[sample, {"Adaptive", Automatic, Automatic}], z];
Plot[{(Sqrt[2] E^(-(d^2/(2 π))) Erfc[Abs[d]/Sqrt[2 π]])/π, f[d]}, {d, -5, 5},
  PlotStyle -> {{Green, Thickness[0.015]}, Blue}]

probability density function of difference of two half-normals

UPDATE

I offer the following without proof for the case where the half-normal parameters differ.

g12[a_, b_] := (D[Integrate[PDF[HalfNormalDistribution[a], u] PDF[HalfNormalDistribution[b], v],{u, 0, ∞}, {v, 0, d + u}, Assumptions -> {a > 0, b > 0, d >= 0}], d]) /. d -> Abs[d]
g[a_, b_] = Piecewise[{{g12[a, b], d <= 0}, {g12[b, a], d > 0}}]

θ1 = 4;
θ2 = 0.75;
dist = TransformedDistribution[
  x - y, {x \[Distributed] HalfNormalDistribution[θ1], 
  y \[Distributed] HalfNormalDistribution[θ2]}]
sample = RandomVariate[dist, {1000000}];
f[z_] = PDF[SmoothKernelDistribution[sample, {"Adaptive", Automatic, Automatic}], z];
Plot[{g[θ1, θ2], f[d]}, {d, -5, 5}, 
  PlotStyle -> {{Green, Thickness[0.015]}, Blue}, PlotRange -> All]

Unequal half-normal parameters

So in general the density function when the half-normal paramaters are not equal is given by

Density function for unequal parameters

$\endgroup$
  • $\begingroup$ Aren't your two expressions for positive and negative d the same? $\endgroup$ – Eckhard Jan 6 '17 at 21:03
  • 1
    $\begingroup$ @Eckhard. I don't think so. The numerators in the Erfc function differ in that one has $\theta1$ and the other has $\theta2$. $\endgroup$ – JimB Jan 6 '17 at 21:06
  • $\begingroup$ Indeed. Hasty reading on my part. $\endgroup$ – Eckhard Jan 6 '17 at 21:09
  • 1
    $\begingroup$ @Eckhard. Review comments are always welcome. I promise that I'll mess something up real soon. $\endgroup$ – JimB Jan 6 '17 at 21:11
6
$\begingroup$

To get the distribution, do

dist = TransformedDistribution[x - y, {x \[Distributed] HalfNormalDistribution[θ], y \[Distributed] HalfNormalDistribution[θ]}]

(Note: \[Distributed] can be quickly entered in Mathematica using Esc+dist+Esc.)

Ideally, we'd like to do

PDF[dist, z]

to get the pdf as a function of z, but this returns unevaluated, even in the case where you specify a value for θ.

For that reason, we take a numerical approach by generating a pseudo-random sample from the distribution and using SmoothKernelDistribution (as suggested by JimBaldwin in a comment under the OP) to generate a numerical approximation of the pdf, as follows.

sample = RandomVariate[dist /. θ -> 1, {1000000}];
f[z_] = PDF[SmoothKernelDistribution[sample], z];
Plot[f[z], {z, -5, 5}]

enter image description here

As suggested by a comment by the OP below, for different parameters, just do

dist = TransformedDistribution[x - y, {x \[Distributed] HalfNormalDistribution[θ1], y \[Distributed] HalfNormalDistribution[θ2]}];

and then

sample = RandomVariate[dist /. {\[Theta]1 -> 4, \[Theta]2 -> 0.75}, {1000000}];
f[z_] = PDF[SmoothKernelDistribution[sample], z];
Plot[f[z], {z, -5, 5}]

enter image description here

$\endgroup$
  • $\begingroup$ Thank you, March. In the case you specify, the random variables have the same theta. Can you show the case where theta can be different? I note that in "dist", theta is identical, and that in your sampling, you specify "theta arrow 1". I suspect that if I add subscripts to the thetas, and add another "theta arrow", for example "theta arrow 0.5", I get the answer I am looking for, but I have been unable to find the correct syntax. So your help is much appreciated. $\endgroup$ – user120911 Jan 6 '17 at 9:04

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.