# Calculating Kurtosis with Beta distribution in R

I am working on an assignment question in R. The problem I am having is to calculate the Kurtosis using Beta distribution.

So far I am able to get the kurtosis value of alpha = 3 beta = 2 using this code:

> kurtosis=(3*alpha^3*beta+3*alpha*beta^3+6*alpha^2*beta^2+alpha^3+beta^3+13*alpha^2*beta+13*alpha*beta^2+alpha^2+beta^2+14*alpha*beta)/((alpha*beta*(alpha+beta+2)*(alpha+beta+3)))


kurtosis [1] 2.892857

The question asks me to generate 1 000 000 independent identically distributed random variables following a Beta distribution. That is easy:

n=100 (I am using 100 for now as 1 000 000 pretty much just hangs my machine up for hours)
p<-rbeta(n, shape1=alpha, shape2=beta)
p

[1] 0.6264632 0.8556988 0.8947960 0.2508657 0.2478296 0.3652293 0.4058185 0.9133942 0.2521895 0.3846150 0.7311879 0.6123567 0.5687362 0.5272337 0.7321268 0.5167320 0.2628915
[18] 0.6827633 0.4034805 0.7767011 0.9246785 0.5946597 0.7077071 0.4014387 0.5738073 0.7225824 0.2736515 0.2674327 0.3449036 0.8527233 0.9142692 0.8225190 0.8081674 0.7760414
[35] 0.5947241 0.8211220 0.6825275 0.1994827 0.3839823 0.8674370 0.4383073 0.6895321 0.6948742 0.4175201 0.4870160 0.4839224 0.5868314 0.6400326 0.5956450 0.7255520 0.4290051
[52] 0.8997959 0.8925137 0.7729661 0.6168284 0.5989350 0.6112734 0.1849604 0.5710712 0.5314022 0.8364838 0.7888875 0.8285151 0.5491049 0.9450382 0.6302862 0.6035714 0.4513216
[69] 0.1588688 0.6744496 0.6010260 0.6908575 0.5219169 0.6897020 0.1486987 0.6108051 0.5032918 0.8729911 0.8216785 0.4243647 0.6615128 0.5479400 0.6326745 0.7449108 0.7917081
[86] 0.3582359 0.7631012 0.6198699 0.4402614 0.5099704 0.3586977 0.4702083 0.5465371 0.7308831 0.3773918 0.5209375 0.2836313 0.6724058 0.7975775 0.6507358


Now I am using this formula to calculate Kurtosis:

(sum(Xi-Xbar)^4/n)/(sum(Xi-Xbar)^2/n)^2


Putting this into R

for (i in 1:n) {
value[i] <- ((p[i]-mean(p))^4/n)/(((p[i]-mean(p))^2/n)^2)}


When I call value I get this return:

 value

[1] 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100
[44] 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100
[87] 100 100 100 100 100 100 100 100 100 100 100 100 100 100


I tried to do the first value of (p[1]-mean(p))^4/n)/(((p[1]-mean(p))^2/n)^2) in mathematica and i get value of 7.389056099

I found that the problem is when I try to divide numerator by denominator, is there another way to do division other than / ? Or am I approaching the problem wrong?

• Did I understand you correctly: You want to calculate the empirical kurtosis of the random variables generated by rbeta, right? So there will be finally one calculated number, which is the kurtosis? Why should there result an array of values? I am not getting this? – Stat Tistician May 5 '13 at 7:33
• see my answer, maybe this helps you. – Stat Tistician May 5 '13 at 7:43

I am not sure, if I got you right, but you want to calculate the kurtosis of the numbers generated by rbeta which you saved into the variable p?

You can use the e1071 library, and NOTICE: THE COMMAND KURTOSIS GIVES YOU THE EXCESS KURTOSIS, so you have to add 3:

library(e1071)
kurtosis(p)+3


Or you calculate it manually: First of all, you do not need a loop. Second, your loop is wrong, because it is not considering the sums correctly.

You can calculate the kurtosis manually with:

m4 <- mean((p-mean(p))^4)
kurt <- m4/(sd(p)^4)
kurt


in this simple case you do not need a loop.

I give you a more simple example (not from a beta distribution):

testvector<-c(1,2,1,2,1,2)
kurtosis(testvector)+3
m4 <- mean((testvector-mean(testvector))^4)
kurt <- m4/(sd(testvector)^4)
kurt


In both cases, the result is

[1] 0.6944444

• Thank you for your reply, that makes much more sense. The actual question is to use a monte-carlo technique to approach skewness and kurtosis using the sample p created. I was under the impression that I should have used the manual calculation for kurtosis from p and then compare it with the kurtosis value from the alpha, beta values. Your explanation is much better! I think that it may be on the right track for my original question, time for me to think more :) Thank you – netuser01 May 5 '13 at 8:47
• @netuser01 you can accept my anwer, by clicking on the button to the left of my post (the hook). – Stat Tistician May 5 '13 at 9:01