0

I have created an iid random sample of size n=20 from a standard normal distribution in R. I have a monte-carlo simulation with m=10^5 scenarios. I am trying to find the probability that the excess kurtosis of the sample lies between -1.5 and 1.5. So far, I have this code:

set.seed(1234)
n<-10
m<-10^5

sample.normal<- rep(0,m)
for(i in (1:m)){
x<-rnorm(n,0,1)
sample.normal[i]<-kurtosis(x,na.rm=FALSE)-3
}

#finding p(-1.5<y2<1.5)
y<-sample.normal
Z<-ecdf(y)
Z(1.5)-Z(-1.5)

However, now I want to change the sample size to a range of values between n=10 and n=200 to show how as the sample size increases, more values lie between -1.5 and 1.5. So I want to calculate the Z(1.5)-Z(-1.5) at 10,20,30,40....200 (bearing in mind m=10^5) (I am ultimately showing that with a larger sample size, excess kurtosis tends to that of a standard normal which is zero) How do I do this using a loop or a different method? Any suggestions would be appreciated.

Chris
  • 6,302
  • 1
  • 27
  • 54
Chloe Stats
  • 1
  • 2

1 Answers1

3

If you want to automate your code to other values of n, the following function simplifies the simulation and then it is called in a sapply loop for n = 10, 20, ..., 200.

Note that one of the arguments is rdist, with default rnorm. If you want to pass arguments to the RNG function rdist use the dots ... argument.

library(e1071)

set.seed(1234)

simKurtosis <- function(n, m = 10^5, rdist = rnorm, ...){
  y <- replicate(m, kurtosis(rdist(n, ...)) - 3)
  Z <- ecdf(y)
  Z(1.5) - Z(-1.5)
}

ex_kurt <- sapply(seq(10, 200, by = 10), simKurtosis)
Rui Barradas
  • 70,273
  • 8
  • 34
  • 66
  • @ChloeStats It saves the results of the simulations in the vector `ex_kurt`. You can then inspect it. According to the help page, only `type = 2` guarantees unbiased kurtosis and the default is `type = 3`. maybe you would want to change that. – Rui Barradas Nov 02 '18 at 21:48