# 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) – zeraoulia rafik Jan 2 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 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 – zeraoulia rafik Jan 2 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 – zeraoulia rafik Jan 2 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 – zeraoulia rafik Jan 1 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?

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 at 21:26