# Plotting difference between two Half Normal Distributions

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?

• If I made the wrong changes to your question, please change it back. – JimB Jan 5 '17 at 22:46
• 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. – march Jan 5 '17 at 23:25
• @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. – JimB Jan 5 '17 at 23:30
• @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. – march Jan 5 '17 at 23:46
• @JimBaldwin : "...Maybe one needs to brute force the integration of the joint probability distribution..." - that's gonna be ugly... – ciao Jan 6 '17 at 0:34

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


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]


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

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

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


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


• 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. – user120911 Jan 6 '17 at 9:04