n = 10
B = 10000
popn = 1:6
set.seed(36784)
draws=sample(x=popn, size=n*B, replace=TRUE)
samples=matrix(draws, nrow=B, ncol=n, byrow=TRUE)
head(samples)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 5 6 2 6 2 5 2 1 2 5
## [2,] 6 6 4 4 5 6 4 1 2 5
## [3,] 6 3 1 2 6 3 1 5 5 1
## [4,] 1 6 4 5 1 3 5 6 1 3
## [5,] 1 5 6 6 4 4 6 4 3 6
## [6,] 2 6 1 6 4 4 6 6 6 6
sim.sums = rowSums(samples)
head(sim.sums)
## [1] 36 43 33 35 45 47
EY=mean(sim.sums)
EY
## [1] 35.0819
s=seq(1:6)
mu=mean(s)
10*mu
## [1] 35
###The approximation value is 35.0819. The theoretical value is 35.//Thus there is a very small difference between the approximation and the theoretical value.
VarY=mean((sim.sums-EY)^2)
VarY
## [1] 29.32419
sig.sq=sum(s^2)/6-mu^2
10*sig.sq
## [1] 29.16667
###The approximation value is 29.32419. The theoretical value is 29.16667. //Thus there is a very small difference between the approximation and the theoretical value.
Prob=mean(sim.sums>=45)
Prob
## [1] 0.0411
###It is concluded that there is approximately a 4% probability //that the sample sum is larger than or equal to 45.
propns=table(sim.sums)/length(sim.sums)
plot(propns, type="h", main="Simulated Probability Distribution of Y",
xlab="Sum of Sampled Values (Y)", ylab="Probability")
