Simulations in R

Question

I'm trying to write a function such as to obtain a test statistic for a vector of size n over 10 simulations. I wrote the following code but I'm not getting the result I need, how can I fix this?

skw=function(n,nsims){
 t=numeric(nsims)
for (i in 1:nsims) {
     x=rnorm(n)
     t[i]=skwness(x)
     zscore=t/(6/n)
return(zscore)
}
}

where:

skwness=function(x){
  n=length(x)
  skew.stat=(1/(n))*(1/(sd(x)^3))*(sum((x-mean(x))^3))
return(skew.stat)
}

Thanks!

Answer 1

You have a couple of issues. The major one is that return should be outside the for loop. Also, you should define t and zscore as vectors, and x as a list.

I think this will work.

Side note: t seems unnecessary in this function. You could get away with using
zscore[i] <- skwness(x[[i]])/(6/n) and get rid of t all together

skwness <- function(x){
  n <- length(x)
  skew.stat <- (1/(n))*(1/(sd(x)^3))*(sum((x-mean(x))^3))
  return(skew.stat)
}

skw <- function(n, nsims){
  t <- zscore <- numeric(nsims)
  x <- vector("list", nsims)

  for (i in 1:nsims) {
    x[[i]] <- rnorm(n)
    t[i] <- skwness(x[[i]])
    zscore[i] <- t[i]/(6/n)
  }

  return(zscore)
}

Giving it a go:

> x <- rnorm(100)
> skwness(x)
[1] 0.2332121
> skw(100, 10)
 [1]  0.6643582 -1.6963196 -2.9192317 -2.7166170  4.9255001  0.0773482  3.9171435
 [8] -3.3993994 -2.4258642  0.7075989