Quiz 8 Scripts

1. Study the lecture notes and screen cast for Simulation, Part 1: set.seed() (2:50).

Set the seed for random number generation to 23. Generate three random numbers from the standard normal distribution and find their sum.

set.seed(23); nums = rnorm(3);
nums

[1]  0.1932123 -0.4346821  0.9132671

sum(nums)

[1] 0.6717973

# answer is 0.6717973, since it's the sum of 0.1932123, -0.4346821, and  0.9132671

2. Study Part 2: replicate() (3:50).

Consider estimating σ2 as 1n∑ni=1(Xi−X¯)2, the average squared deviation from the mean. To see why we usually divide by n−1 instead of n, do the following simulation.

Set the seed to 7. Do the following calculation N=100 times. Find a random sample of size n=4 from the standard normal distribution. Find the average squared deviation from the mean of the sample data (using the formula above). Then find the average of these 100 variance estimates.

set.seed(7)
#x = replicate(n=100, expr = mean(rnorm(4)))
rep = replicate(n=100, expr = {x = rnorm(4); sum((x-mean(x))^2)})
z = mean(rep)*1/4; z

[1] 0.7585977

# ask why this is correct
set.seed(7)
x = replicate(n=4, expr = mean(rnorm(1)))
rep = replicate(n=100, expr = {x = rnorm(4); sum((x-mean(x))^2)})
z = mean(rep)*1/4; z

[1] 0.7474647

#xtest = replicate(n=100, expr = rnorm(4))
#ytest = replicate(n=100, expr = sum((xtest-(mean(xtest)))^2))
#ztest = mean(ytest)*1/4
#ztest

3. Notice that your estimate is low. Correct it by multiplying by 43, which corresponds to using the usual formula for the sample variance, s2=1n−1∑ni=1(Xi−X¯)2. (Notice that the estimate is now close to 1, the true population variance.) What is the value of this corrected estimate?

set.seed(7)
x = replicate(n=4, expr = mean(rnorm(1)))
rep = replicate(n=100, expr = {x = rnorm(4); sum((x-mean(x))^2)})
z = mean(rep)*(1/3); z

[1] 0.9966197

0.7574938*4/3

[1] 1.009992

4. Study Part 3: normal distribution (9:26).

In the online lecture, I found the sample standard deviation of a sample of 1000 random numbers from the N(μ=7,σ=3) distribution. What was its value?

mean(w <- rnorm(n=1000, mean=7, sd=3))
sd(w)
# 3.122579

5. Study Part 4: t distribution (5:16). Study Part 5: p-value (6:15).

Set the seed to 3. Repeat the p-value simulation for the test for mean pouring time, this time using 10,000 replicates. What is the simulated p-value?

set.seed(3)
x = c(118, 121, 113, 116, 117, 112, 113)
out = t.test(x, mu =119.5)
mu = 119.5
sigma = sd(x)
n = length(x)
N = 10000
t = replicate(N, { x=rnorm(n, mean=mu, sd=sigma); (mean(x) - mu)/(sd(x)/sqrt(n)) })
more.extreme = (abs(t) > abs(out$statistic))
(simulated.p.value = sum(more.extreme) / N)

[1] 0.021

out$p.value

[1] 0.02164744