SkewNormalDistribution with a mean value

Question

Since it got incredibly good answers the last time i posted a question, i try it again! ;-)

I programmed a chemical/physical simulation in mathematica, where i can use residence time distribution as an input parameter. No I'd like to compare different theoretical distributions. As a criteria, the distributions need to have the same mean value μ. It's easy to do it with NormalDistributions as follows:

dist1 = RandomVariate[NormalDistribution[0.216208304, 0.0025], 10^4];

Show[Histogram[dist1, 14, "ProbabilityDensity"], AxesLabel->{"\[Tau]/s","\!\(\*OverscriptBox[\(V\), \(.\)]\)/\!\(\*SuperscriptBox[\(m\),\\(3\)]\)\!\(\*SuperscriptBox[\(s\), \(-1\)]\)"}]

HistogramList[dist1, 14, "ProbabilityDensity"]

Therewith i obtaine a Histogram plot where i can export the container values and size to my simulation.

More difficulties i had when i tried to create a skew distribution with the same mean value as in the normal distribution, means: 0.216208304 I tried it with the following code:

CSTR = RandomVariate[SkewNormalDistribution[0.216208304, 1, 1000], 10^4];

But the location parameter μ in the SkewNormalDistribution doesn't seems to be the mean value. If i create to mean value of the CSTR-List with Mean[CSTR], i get 1.01582 - so totally wrong.

Anybody has an idea how i can create a skew distribution with a mean value of 0.2162... ?

Thanks in advance and a happy new year! Dani

Answer 1

We have

Mean[SkewNormalDistribution[loc[1, 1000, 0.216208304], 1, 1000]]

0.2162083

where

loc[σ_, α_, μ_] = μ - (Sqrt[2/π] σ α)/Sqrt[1 + α^2];

is obtained from

First[Solve[Mean[SkewNormalDistribution[l, σ, α]] == μ, l] // FullSimplify]

$\left\{l\to \mu -\frac{\sqrt{\frac{2}{\pi }} \alpha \sigma }{\sqrt{\alpha ^2+1}}\right\}$

Actually the code below gives a closed form for the mean when discarding simulations less than $0$. Solve can't determine any parameter given the other, so you need to use a numeric function like FindRoot.

expr = With[{anti = Integrate[
                      PDF[SkewNormalDistribution[a, σ, α], x] (x - a)/
                       (1 - CDF[SkewNormalDistribution[a, σ, α], 0]), x]},
  Limit[anti, x -> ∞, Assumptions -> {{(a Sqrt[1 + α^2])/(Sqrt[2] σ),
                       (a α)/(Sqrt[2] σ)} ∈ Reals, σ > 0, Sqrt[1 + α^2] > 0}] -
  Limit[anti, x -> 0] + a] // FullSimplify

$\frac{\sqrt{\frac{2}{\pi }} \sigma \left(e^{-\frac{a^2}{2 \sigma ^2}} \text{erfc}\left(\frac{a \alpha }{\sqrt{2} \sigma }\right)+\frac{\alpha \text{erf}\left(\frac{a \sqrt{\alpha ^2+1}}{\sqrt{2} \sigma }\right)}{\sqrt{\alpha ^2+1}}\right)}{4 T\left(\frac{a}{\sigma },\alpha \right)+\text{erf}\left(\frac{a}{\sqrt{2} \sigma }\right)+1}+\frac{\sqrt{\frac{2}{\pi }} \alpha \sigma }{\sqrt{\alpha ^2+1} \left(4 T\left(\frac{a}{\sigma },\alpha \right)+\text{erf}\left(\frac{a}{\sqrt{2} \sigma }\right)+1\right)}+a$

Answer 2

SeedRandom[0];

dist1 = NormalDistribution[0.216208304, 0.0025];

Random data drawn from the NormalDistribution

data1 = RandomVariate[dist1, 10^4];

To find a SkewNormalDistribution that best fits the random data

dist2 = EstimatedDistribution[data1, SkewNormalDistribution[m, s, a]]

(* SkewNormalDistribution[0.215064, 0.0027407, 0.597653] *)

Since you also indicated that you want the values of the SkewNormalDistribution to be restricted to nonnegative values, use TruncatedDistribution.

Using the parameters of the SkewNormalDistribution (dist2) as initial estimates of the parameters for the TruncatedDistribution

dist3 = EstimatedDistribution[data1, 
  TruncatedDistribution[{0, Infinity}, SkewNormalDistribution[m, s, a]], 
  Transpose[{{m, s, a}, List @@ dist2}]]

(* TruncatedDistribution[{0, ∞}, 
 SkewNormalDistribution[0.215063, 0.00274102, 0.598107]] *)

Random data drawn from the TruncatedDistribution

data2 = RandomVariate[dist3, 10^4];

Column[{
  Show[
   Histogram[data1, 25, "PDF",
    PlotLabel -> dist1,
    ImageSize -> Large],
   Plot[Tooltip[PDF[dist1, x]], {x, Min[data1], Max[data1]}]],
  Show[
   Histogram[data2, 25, "PDF",
    PlotLabel -> dist3,
    ImageSize -> Large],
   Plot[Tooltip[PDF[dist3, x]], {x, 0, Max[data2]},
    PlotRange -> All]]}]

Comparing the Mean of the TruncatedDistribution with that of the NormalDistribution

Abs[Mean[data2] - Mean[dist1]]

(* 0.000020577 *)

Answer 3

You can re-parametrize SkewNormalDistribution so that the first parameter is its mean:

snd[μ_, σ_, α_]:= SkewNormalDistribution[μ - Mean[SkewNormalDistribution[0, σ, α]], σ, α]

Mean[snd[μ, σ, α]]

μ

Mean[snd[0.216208304, 1, 1000]]

0.216208304

Variance and Skewness are not affected by this reparametrization:

Variance[snd[μ, σ, α]] == Variance[SkewNormalDistribution[μ, σ, α]]

True

Skewness[snd[μ, σ, α]] == Skewness[SkewNormalDistribution[μ, σ, α]]

True