1
$\begingroup$

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

$\endgroup$
1
$\begingroup$

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 *)
| improve this answer | |
$\endgroup$
  • $\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$ – zeraoulia rafik Jan 2 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 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$ – zeraoulia rafik Jan 2 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$ – zeraoulia rafik Jan 2 at 21:27
4
$\begingroup$
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}]

enter image description here

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 *)
| improve this answer | |
$\endgroup$
  • $\begingroup$ Great answer , Thanks so much , I up vote $\endgroup$ – zeraoulia rafik Jan 1 at 21:04
1
$\begingroup$

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?

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

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

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

| improve this answer | |
$\endgroup$
  • $\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 at 21:26
  • $\begingroup$ @Bill, it is just a wrong typo i have edited it $\endgroup$ – zeraoulia rafik Jan 1 at 21:31

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.