# Generate variables from 
# standard normal distributions
Z1 <- rnorm(1000)
Z2 <- rnorm(1000)
Z3 <- rnorm(1000)
Z4 <- rnorm(1000)

# Create the sum of squares
Y <- Z1^2 + Z2^2 + Z3^2 + Z4^2

# What distribution will Y have?
plot(density(Y),
     main="Chi-sq Simulated",
     col="navy", lwd=2,
     las=1)