I have tried to compute Skewness and variance and moment using the code below, I have tried first to create a distribution from my PDF using the command distr = [PDF, {x,x_min , x_max}], but when i tried to plug my PDF there to get Variance I received the error message "Aborted" computation. Any help?

This is my attempt to get Variance:

distr = ProbabilityDistribution[(Pi! + 6 - Pi)(z - σ)^2/(Sqrt[(1 +  0.0025μ^2)*2*Pi]) Exp[-(z - σ)^2/((1 +0.0025 μ^2))*
  Sqrt[(1 +0.0025  μ^2)*2*Pi]*Erf[(z - σ)^2/((1 + 0.0025μ^2)) Sqrt[(1 +0.0025 μ^2)*2*Pi]]],{z,-Infinity,Infinity}/.{σ->0,μ->0.05}]

Note: This question is related to my question here


4 Answers 4


First we start with the function you supplied that needs an appropriate multiplier to turn it into a probability density function:

$$ F(z,\mu,\sigma)=\frac{2 (z-\sigma )^2 \exp \left(-\frac{(z-\sigma )^2 \sqrt{\left(1+0.25 \mu ^2\right) 2 \pi } \text{erf}\left(\frac{(z-\sigma )^2 \sqrt{\left(1+0.25 \mu ^2\right) 2 \pi }}{1+0.25 \mu ^2}\right)}{1+0.25 \mu ^2}\right)}{\pi ^2 \sqrt{\left(1+0.25 \mu ^2\right) 2 \pi }} $$

As suggested by @Xi'an (https://stats.stackexchange.com/questions/442784/mean-and-variance-of-a-non-standard-pdf) you should replace

$$\frac{(z-\sigma)^2 \sqrt{\left(1+0.25 \mu ^2\right) 2 \pi}}{1+0.25 \mu ^2}$$


$$\frac{(z-\sigma)^2 \sqrt{2 \pi}}{\sqrt{1+\mu^2/4}}$$

This gives us

$$ F(z,\mu,\sigma)=\frac{1}{\pi^3} \frac{(z-\sigma )^2 \sqrt{2 \pi }}{\sqrt{1+ \mu^2/4}} \exp\left(-\frac{(z-\sigma )^2 \sqrt{2 \pi }}{\sqrt{1+ \mu^2/4}} \text{erf}\left(\frac{(z-\sigma )^2 \sqrt{2 \pi }}{\sqrt{1+ \mu^2/4}}\right) \right)$$

Note that $F(z,\mu,\sigma)$ is symmetric about $\sigma$ so that the constant of integration will not depend on $\sigma$ but will depend on $\mu$. This also means that all odd central moments (assuming they exist) will be zero.

We could just use numerical integration to find the integration constant for an assorted set of values for $\mu$ but that turns out to be unnecessary. We perform a change of variables where

$$h = (z-\sigma)\sqrt[4]{2\pi}/\sqrt[4]{1+\mu^2/4} $$

So $dz =\frac{\sqrt[4]{\frac{\mu ^2}{4}+1}}{\sqrt[4]{2 \pi }} dh$.

The constant of integration becomes

$$c=\pi^3\frac{\sqrt[4]{2 \pi }}{\sqrt[4]{1+ \mu^2/4}}\frac{1}{\int_{-\infty}^{\infty} h^2 \exp\left(-h^2 \text{erf}\left(h^2\right)\right)dh}$$

The reciprocal of the integral in the constant of integration is

1/NIntegrate[h^2 Exp[-h^2 Erf[h^2]], {h, -\[Infinity], \[Infinity]}]
(* 1.0105750026505362 *)

resulting in the constant of integration being approximately

$$c\simeq 1.0105750026505362 \times \pi^3\frac{\sqrt[4]{2 \pi }}{\sqrt[4]{1+ \mu^2/4}}$$

The pdf can be written as

$$1.0105750026505362 \times \frac{\sqrt[4]{2 \pi }}{\sqrt[4]{1+ \mu^2/4}} \frac{(z-\sigma )^2 \sqrt{2 \pi }}{\sqrt{1+ \mu^2/4}} \exp\left(-\frac{(z-\sigma )^2 \sqrt{2 \pi }}{\sqrt{1+ \mu^2/4}} \text{erf}\left(\frac{(z-\sigma )^2 \sqrt{2 \pi }}{\sqrt{1+ \mu^2/4}}\right) \right)$$

pdf[z_, σ_, μ_] := Module[{h, k},
  k = (2 π/(1 + μ^2/4))^(1/4);
  h = k (z - σ);
  1.0105750026505362*k*h^2 Exp[-h^2 Erf[h^2]]]

We could then determine the central moments for a specific value of μ (say $\mu=1/20$):

NIntegrate[(z - 1)^2 pdf[z, 1, 1/20], {z, -∞, ∞}]
(* 0.565411 *)
NIntegrate[(z - 1)^4 pdf[z, 1, 1/20], {z, -∞, ∞}]
(* 0.545414 *)
NIntegrate[(z - 1)^6 pdf[z, 1, 1/20], {z, -∞, ∞}]
(* 0.751608 *)

But the $i$-th central moment for any particular value of μ will be

centralMoment[i_, μ_] := 
 If[OddQ[i], 0, 
  1.0105750026505362 (2 π/(1 + μ^2/4))^(-i/4) 
  NIntegrate[h^(2 + i) Exp[-h^2 Erf[h^2]], {h, -∞, ∞}]]

Here is the variance:

FullSimplify[centralMoment[2, μ], μ > 0]
0.282617 Sqrt[4 + μ^2]

As a check:

centralMoment[2, 1/20]
(* 0.565411 *)
centralMoment[4, 1/20]
(* 0.545414 *)
centralMoment[6, 1/20]
(* 0.751608 *)
  • $\begingroup$ Nice answer Jim ,according to the even moments ( Moment generating function) it is good to seek about its closed form since \mu lie in narrow region(0,1) $\endgroup$ Jan 2, 2020 at 20:32
  • $\begingroup$ It would be good to have a closed if indeed a closed form exists for the moment generating function. What makes you think it exists? It would seem to me that finding a moment generating function is beyond your original question. $\endgroup$
    – JimB
    Jan 2, 2020 at 21:21
  • $\begingroup$ I have asked that question in cross validated , I have asked about Laplace transform of the PDF which it does mean the moment generating function $\endgroup$ Jan 2, 2020 at 21:23
  • $\begingroup$ And i don't think so about existence we may use approximation such us PAde and Fadi BRUNO FORMULA for the representation $\endgroup$ Jan 2, 2020 at 21:27
distr = With[{σ = 0, μ = 1/20},
  ProbabilityDistribution[(Pi! + 6 - 
      Pi) (z - σ)^2/(Sqrt[(1 + μ^2/400)*2*
        Pi]) Exp[-(z - σ)^2/((1 + μ^2/400))*
      Sqrt[(1 + μ^2/400)*2*Pi]*
      Erf[(z - σ)^2/((1 + μ^2/400)) Sqrt[(1 + μ^2/400)*2*
          Pi]]], {z, -Infinity, Infinity}]]

Checking whether the distribution is properly normalized

NIntegrate[PDF[distr, z], {z, -Infinity, Infinity},
 WorkingPrecision -> 12]

(* 0.999359291501 *)

Since it is not, redefine as

a = With[{σ = 0, μ = 1/20},
  NIntegrate[(z - σ)^2/(Sqrt[(1 + μ^2/400)*2*
        Pi]) Exp[-(z - σ)^2/((1 + μ^2/400))*
      Sqrt[(1 + μ^2/400)*2*Pi]*
      Erf[(z - σ)^2/((1 + μ^2/400)) Sqrt[(1 + μ^2/400)*2*
          Pi]]], {z, -Infinity, Infinity}, WorkingPrecision -> 12]]

(* 0.0994734762093 *)

distr2 = With[{σ = 0, μ = 1/20},
   ProbabilityDistribution[(z - σ)^2/(Sqrt[(1 + μ^2/400)*2*
         Pi]) Exp[-(z - σ)^2/((1 + μ^2/400))*
        Sqrt[(1 + μ^2/400)*2*Pi]*
        Erf[(z - σ)^2/((1 + μ^2/400)) Sqrt[(1 + μ^2/400)*2*
            Pi]]]/a, {z, -Infinity, Infinity}]];


NIntegrate[PDF[distr2, z], {z, -Infinity, Infinity},
 WorkingPrecision -> 12]

(* 1.00000000000 *)

Plot[PDF[distr2, z], {z, -2.5, 2.5}]

enter image description here


(* 0 *)

This is to as expected from the symmetry of the PDF

PDF[distr2, -x] == PDF[distr2, x]

(* True *)

The modes are at

NSolve[{D[PDF[distr2, z], z] == 0, 0 < Sqrt[z^2] < 2}, z]

(* {{z -> -0.570721333298}, {z -> 0.570721333298}} *)

Since the mean is zero, the Variance is just the second Moment

var = NExpectation[z^2, z \[Distributed] distr2, WorkingPrecision -> 12]

(* 0.565236192459 *)

From the symmetry, the odd moments are zero; consequently, the Skewness is zero

NExpectation[z^3, z \[Distributed] distr2, WorkingPrecision -> 12]/var^(3/2)

(* 0.*10^-13 *)
  • $\begingroup$ Great answer , Thanks so much , I up vote $\endgroup$ Jan 1, 2020 at 21:04

I think your replacement pattern might be incorrectly positioned.

Try this

  (1+0.0025 μ^2)*Sqrt[(1+0.0025μ^2)*2*Pi]*Erf[(z-σ)^2/

Can you now use the standard formula from statistics to calculate skewness from these?

Because of the history of your attempts at this question, now you need to think very very carefully about this. Are the calculations I did correct? Did I make any mistakes? Is the result correct? How can you check this result several different ways? From the plot it appears your function is symmetric about zero. If so then a mean of zero makes sense. But what can you do to test this and verify the function is symmetric and thus the mean should be exactly zero. What about that value for the variance? Is that result correct? How can you tell if it is correct. If your function were a gaussian bell curve then we might have the intuition and background to look at the size and shape of the bell and have more confidence whether the calculated variance was correct or not. That value seems awfully small given the size and shape of your plot. And the std is similarly small. Those values worry me. How can you have any confidence they are correct? Just jumping with glee because you finally got a number can be a very very bad mistake. How are you going to be able to convince yourself that your number is correct? Because the plot looks symmetric I have an intuition about what the skewness might be. Will your calculations be correct, whether they give the result that I expect or not?

Someone much much wiser than I once said to me "Delay as long as you possibly can thinking you know what the answer is, because the moment you do that everything I have taught you is thrown out the window."

  • $\begingroup$ In my attempt now i have got the variance is 30.886 $\endgroup$ Jan 1, 2020 at 20:27
  • $\begingroup$ yes coorect , I have fot also 0.56.., thank us so much now am very happy about my PDF $\endgroup$ Jan 1, 2020 at 20:29
  • $\begingroup$ I think we can get Variance for Mu in (0,1) $\endgroup$ Jan 1, 2020 at 20:30
  • $\begingroup$ the variance is equal :1/ sqrt(pi) $\endgroup$ Jan 1, 2020 at 20:50
  • $\begingroup$ Yes that function is symetric , I have published that function before in journal , just i want this PDF to check its application in plasma physics this is My Goal to build that PDF , and am sure that all what you have got is correct $\endgroup$ Jan 1, 2020 at 20:53

This just a hint , The nice answer given by @Bill is that Variance equal $\frac{1}{\sqrt\pi}=\dfrac{\operatorname{erf}(x)}{\int_{-\infty}^{\infty}\exp {(-t^2)}dt}$ and that PDF is symetric

  • $\begingroup$ I am confused. Where did x come from? Mathematica says Erf[x]*Integrate[E^-t^2,{t,-Infinity,Infinity}] == Sqrt[Pi]*Erf[x] $\endgroup$
    – Bill
    Jan 1, 2020 at 21:26
  • $\begingroup$ @Bill, it is just a wrong typo i have edited it $\endgroup$ Jan 1, 2020 at 21:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.