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")