1
$\begingroup$

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?

$\endgroup$
2
  • $\begingroup$ 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). $\endgroup$
    – flinty
    Nov 20 '20 at 13:46
  • $\begingroup$ 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. $\endgroup$
    – JimB
    Nov 20 '20 at 18:44
2
$\begingroup$
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], 
  GammaDistribution[26.8782, 0.0793163]}] *)

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

enter image description here

$\endgroup$
1
$\begingroup$

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

Histogram, numerical approximation, and skew normal approximation

$\endgroup$

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.