# Finding skewness and kurtosis numerically in Mathematica

I have a function given by

f[a_, b_, g_, c_, k_] = ProbabilityDistribution[a b c k x^(c - 1)
(1 + x^c)^(k - 1) ((1 + x^c)^k - 1)^(-b - 1) (1 + g ((1 + x^c)^k -1)
^-b)^(-(a/g) - 1), {x, 0, inf},
Assumptions -> a > 0 && b > 0 && g > 0 && c > 0 && k > 0]


I want to find skewness for certain values of a,b,g,c,k.

How to find the skewness and kurtosis numerically in Mathematica?

Starting with a corrected version of your ProbabilityDistribution

f[a_, b_, g_, c_, k_] :=
ProbabilityDistribution[
a b c k x^(c - 1) (1 + x^c)^(k - 1) ((1 + x^c)^k - 1)^(-b - 1) (1 +
g ((1 + x^c)^k - 1)^-b)^(-(a/g) - 1), {x, 0, Infinity},
Assumptions -> a > 0 && b > 0 && g > 0 && c > 0 && k > 0]


One can define the numerical mean value as

mu[a_, b_, g_, c_, k_] :=
NExpectation[x, x \[Distributed] f[a, b, g, c, k]]


the numerical skewness as

skewness[a_, b_, g_, c_, k_] :=
NExpectation[((x - mu[a, b, g, c, k])/Sqrt@NExpectation[(x - mu[a, b, g, c, k])^2,
x \[Distributed] f[a, b, g, c, k]])^3, x \[Distributed] f[a, b, g, c, k]]


and the numerical kurtosis as

kurtosis[a_, b_, g_, c_, k_] :=
NExpectation[(x - mu[a, b, g, c, k])^4, x \[Distributed] f[a, b, g, c, k]]/
NExpectation[(x - mu[a, b, g, c, k])^2, x \[Distributed] f[a, b, g, c, k]]^2


Testing with 1, 2, 3, 4, 5 as the "certain values of a,b,g,c,k"

mu[1, 2, 3, 4, 5]
skewness[1, 2, 3, 4, 5]
kurtosis[1, 2, 3, 4, 5]


0.551698
-0.397584
2.94902

and comparing with a plot of the PDF

Plot[PDF[f[1, 2, 3, 4, 5]][x], {x, -0.3, 1.75}, PlotRange -> All]


Using central moments as suggested in the comment by Guess who it is.

zm[a_, b_, g_, c_, k_][m_] :=
NExpectation[(x - mu[a, b, g, c, k])^m, x \[Distributed] f[a, b, g, c, k]]


Now the skewness can be defined as

mSkewness[a_, b_, g_, c_, k_] :=
zm[a, b, g, c, k][3]/(zm[a, b, g, c, k][2])^(3/2)

mSkewness[1, 2, 3, 4, 5]


-0.397584

and the kurtosis as

mKurtosis[a_, b_, g_, c_, k_] :=
zm[a, b, g, c, k][4]/(zm[a, b, g, c, k][2])^2

mKurtosis[1, 2, 3, 4, 5]


2.94902

• Another possibility would be to just define a function for the moments, and then assemble the expressions for variance, skewness, and kurtosis in terms of those moments. Tho, this may or may not be more numerically stable, so some testing would be needed. – J. M. will be back soon Jul 1 '15 at 17:28
• One could also use the definitions based on the cumulants. – Karsten 7. Jul 1 '15 at 17:56
• @wolfies Thank you for the clarifying and explanatory comment to my earlier answer. – Karsten 7. Jul 1 '15 at 18:08
• @Karsten7.thanks a lot both solutions are working just fine. – SA-255525 Jul 2 '15 at 4:02