# Expression for variation in the fourth moment of standard normal vs. sample size

Question:
What is the expression for the variation in the fourth (not-centralized) moment as a function of sample size for the standard normal distribution and is there a proper analytic/symbolic way to derive it?

Details:
When I use random number generator to create samples from a standard normal distribution, then I compute the mean, and I do this many times per sample-level, I can create a plot like the following: My textbooks told me that the relationship between the variation and the sample size was:

$$Err \propto \frac{1}{\sqrt{n}}$$

where $Err$ is the standard deviation of the estimate and $n$ is the sample size.

When I plug it into this, I get a visual confirmation. When I plug into a linear fit I get the following summary:

> summary(est)

Call:
lm(formula = s2 ~ I(sqrt(1/n)))

Residuals:
Min         1Q     Median         3Q        Max
-0.0097079 -0.0004087  0.0000041  0.0004380  0.0093119

Coefficients:
Estimate Std. Error  t value Pr(>|t|)
(Intercept)  -2.633e-04  5.312e-05   -4.957  8.4e-07 ***
I(sqrt(1/n))  1.005e+00  7.212e-04 1393.680  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.0009955 on 994 degrees of freedom
Multiple R-squared:  0.9995,    Adjusted R-squared:  0.9995
F-statistic: 1.942e+06 on 1 and 994 DF,  p-value: < 2.2e-16


The R-squared of 99.95% is an indicator of a pretty good fit.

What is the expression for the variation in the fourth (not-centralized) moment as a function of sample size for the standard normal distribution?

The plot of the variation in fourth moment versus sample size looks like this: When I plot Error versus square root of inverse sample count, I do not get a straight line. I could wire a stack of transformations on "n" into "glmulti" (link) but that can be a lot of work for little return.

Some "hacking" in algebra gives me a fair error expression, but I have no way of telling if it is local to my data, or more globally valid.

Here is my source code:

set.seed(5) #for reproducibility

n <- seq(from=5,to=1000,by = 1)

#number of iterations at each sample size
N <- 3000

#predeclare
store1 <- matrix(0, nrow=length(n),ncol=N)
store2 <- matrix(0, nrow=length(n),ncol=N)
mymean <- numeric()
mymom<- numeric()

nn <- matrix(0, nrow=length(n),ncol=N)

#for each sample size
for(i in 1:length(n)){

#repeat the measure "N" times
for (j in 1:N){

#store so we can separate it out later
nn[i,j] <- n[i]

#take sample
y <- rnorm(n = nn[i,j],mean = 0,sd = 1)

#compute moment
mymean <- mean(y)
mymom <- mean( (y-mean(y))^4)

#store
store1[i,j] <- mymean
store2[i,j] <- mymom
}
}

##compute variation

#predeclare
s1 <- numeric()
s2 <- numeric()

s3 <- numeric()
s4 <- numeric()

#loop
for (i in 1:length(n)){

#find elements which have a particular sample size
idx <- which(nn==n[i],arr.ind=T)

#compute statistics of first moment (aka mean)
s1[i] <- mean(store1[idx])
s2[i] <- sd(store1[idx])

#compute statistics of fourth moment (not centered
s3[i] <- mean(store2[idx])
s4[i] <- sd(store2[idx])

}

#make figures

plot(n, s2,xlab="Sample count",ylab="Error in mean")
lines(lowess(x=n,y=s2,f=0.01),col="Red",lwd=2)
grid()

plot(sqrt(1/n), s2,xlab="SQRT inverse Sample count",ylab="Error in mean")
lines(lowess(x=sqrt(1/n),y=s2,f=0.01),col="Red",lwd=2)
grid()

plot(n, s4,xlab="Sample count",ylab="Error in fourth moment")
lines(lowess(x=n,y=s4,f=0.005),col="Red",lwd=2)
grid()

plot(sqrt(1/n), s4,xlab="SQRT inverse Sample count",ylab="Error in fourth moment",
xlim=c(0,0.2),ylim=c(0,2))
lines(lowess(x=sqrt(1/n),y=s4,f=0.02),col="Red",lwd=2)
grid()

#fit to models
est <- lm(s2~I(sqrt(1/n)))
summary(est)

est2 <- lm(s2~I(sqrt(1/n))+I((1/n)^(1/4)) )
summary(est2)


Here is the algebra for the standard normal case. I am using formulas for moments from https://en.wikipedia.org/wiki/Normal_distribution#Moments Let $X_1, \dotsc, X_n$ be iid observations from a standard normal distribution. The fourth empirical moment (about zero) is given by $$\DeclareMathOperator{\E}{\mathbb{E}} \frac{1}{n} \sum_i X_i^4.$$Its expected value can be calculated: $$\E \left\{ \frac1{n} \sum_i X_i^4 \right\}= 3$$ and its variance: $$\DeclareMathOperator{\Var}{\mathbb{Var}} \Var\left\{ \frac1{n} \sum_i X_i^4 \right\} = (\frac1{n})^2\sum_i \Var X_i^4= \frac{96}{n}$$ and some details of the last calculation here: $$\Var X_i^4 = \E \left( X_i^4 -\E X_i^4 \right)^2 = \E \left( X_i^4-3\right)^2 = \E X_i^8 - 2\cdot 3 \E X_i^4 + 9 = \sigma^8\cdot 105 -6 \cdot 3 + 9= 96$$