# TransformedDistribution of normal and skew-normal distributions

I have the following distributions:

distC =NormalDistribution[251.563, 0.317394]
distS =SkewNormalDistribution[252.973, 0.503173, 5.59741]


I want to analyse what results when distS-distC is used. I made the following approach:

gapFkt[d_] :=
PDF[TransformedDistribution[
s - c, {c \[Distributed] distC, s \[Distributed] distS}], d]


When I try to plot it, it takes veeery long. So I tried to make some points first using a table.

pkt = ParallelTable[gapFkt[d], {d, 0., 2., .1}]


This is very slow as well and all the values are 0. , although at ~1. there is a peak of ~0.8. I can use a numerical approach. However, is there a way to do it without numerical brute force?

• tdist = TransformedDistribution[s - c, {c \[Distributed] distC, s \[Distributed] distS}]; Apparently PDF[tdist,x] is given a value of zero rather than a distribution. Looks like a bug since it should have a valid PDF (see RandomVariate[tdist, 10000] // Histogram). Nov 20 '20 at 13:46
• I don't think it's exactly a fixable bug in that there is (I think) no closed-form for the pdf. If one tries smaller means (but with the same mean difference), just the expression is returned rather than a zero: distC = NormalDistribution[0, 0.317394]; distS = SkewNormalDistribution[252.973 - 251.563, 0.503173, 5.59741]; gapFkt[d_] := PDF[TransformedDistribution[s - c, {c \[Distributed] distC, s \[Distributed] distS}], d] gapFkt[2]. One might need to compute the cdf numerically followed by the pdf.
– JimB
Nov 20 '20 at 18:44

Clear["Global*"]

distC = NormalDistribution[251.563, 0.317394] // Rationalize;

distS = SkewNormalDistribution[252.973, 0.503173, 5.59741] //
Rationalize[#, 0] &;

SeedRandom[1234]

With[{n = 25000}, data =
RandomVariate[distS, n] -
RandomVariate[distC, n]];

distD = FindDistribution[data, RandomSeeding -> 1234]

(* MixtureDistribution[{0.574031,
0.425969}, {NormalDistribution[1.60781, 0.334762],

Show[
Histogram[data, Automatic, "PDF"],
Plot[PDF[distD, d], {d, Min@data, Max@data}]]


Here are two ways to get a good approximation of the density with one of the ways being about as parsimonious as you can get: numerical approximation and fit a skew normal with the same mean, variance, and skewness (method of moments).

A brute-force definition of the pdf of the difference at a value $$z$$ is defined to be

Integrate[PDF[distC, x] PDF[distS, z + x], {x, -∞,∞}]


But there is no closed-form solution. Next one attempts to use numerical integration for particular values of $$z$$:

NIntegrate[PDF[distC, x] PDF[distS, 2 + x], {x, -∞,∞}]


But here one gets 0. for probably all values of $$z$$. The fix is just to numerically integrate over the range of values of $$x$$ that cover most of the positive density. Here is a table of values using that technique:

xLow = 251.563 - 4*0.317394
xHigh = 251.563 + 4*0.317394
pdf = Table[{z, NIntegrate[PDF[distC, x] PDF[distS, z + x], {x, xLow, xHigh}]},
{z, 0, 4, 0.01}]


It is not unreasonable to think that maybe that a linear combination of a normal and a skew normal might have approximately a skew normal distribution. So we find the parameters of a skew normal distribution that has the same mean, variance, and skewness as s - c:

dist = TransformedDistribution[s - c, {s \[Distributed] distS, c \[Distributed] distC}]
{μ0, σ0, α0} = Select[{μ, σ, α} /.
Solve[{Mean[dist] == Mean[SkewNormalDistribution[μ, σ, α]],
Variance[dist] == Variance[SkewNormalDistribution[μ, σ, α]],
Skewness[dist] == Skewness[SkewNormalDistribution[μ, σ, α]]},
{μ, σ, α}, Reals], #[[2]] > 0 &][[1]]


Now check this with a large random sample:

SeedRandom[12345];
zz = RandomVariate[distS, 100000] - RandomVariate[distC, 100000];

(* Plot them all *)
Show[Histogram[zz, "FreedmanDiaconis", "PDF"],
Plot[{pdf[z], PDF[SkewNormalDistribution[μ0, σ0, α0], z]}, {z, 0, 4},
PlotStyle -> {{Thickness[0.02], LightGray}, Red},
PlotLegends -> {"Numerical integration", "Skew normal approximation"}]]
`