# Compute moments, variance, skewness, mean from the PDF of a non standard distribution

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

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

with

$$\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 *)

• 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) Jan 2, 2020 at 20:32
• 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.
– JimB
Jan 2, 2020 at 21:21
• 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 Jan 2, 2020 at 21:23
• And i don't think so about existence we may use approximation such us PAde and Fadi BRUNO FORMULA for the representation 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}]];


Checking

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

(* 1.00000000000 *)

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


Mean[distr2]

(* 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 *)

• Great answer , Thanks so much , I up vote Jan 1, 2020 at 21:04

I think your replacement pattern might be incorrectly positioned.

Try this

d=(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]]]/.{σ->0,μ->0.05}
Plot[d,{z,-2,2},PlotRange->All]
mean=NIntegrate[z*d,{z,-Infinity,Infinity}]
variance=NIntegrate[(z-mean)^2*d,{z,-Infinity,Infinity}]
std=Sqrt[variance]


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

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."

• In my attempt now i have got the variance is 30.886 Jan 1, 2020 at 20:27
• yes coorect , I have fot also 0.56.., thank us so much now am very happy about my PDF Jan 1, 2020 at 20:29
• I think we can get Variance for Mu in (0,1) Jan 1, 2020 at 20:30
• the variance is equal :1/ sqrt(pi) Jan 1, 2020 at 20:50
• 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 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

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