Question

Consider three six-sided dice, and let random variable Y = the value of the face for each. The probability mass of function of Y is given by the following table:

y	1	2	3	4	5	6	otherwise
P(Y=y)	0.35	0.30	0.25	0.05	0.03	0.02	0

1. Roll the three dice and let random variable X = sum of the three faces. Repeat this experiment 50000 times.

1. Find the simulated probability mass function (pmf) of random variable X.
2. Find the simulated mean of the random variable X.
3. Find the simulated standard deviation of random variable X.
4. Find the simulated P(2 < X < 13).
5. Draw the histogram of the pmf of X.

1. Generate 2000 samples of size 50000 values for the random X from distribution in part1. For each of these 2000 samples calculate the mean of X.

1. Find the simulated probability that the mean is between 6.49 and 6.53 exclusive.
2. Find the simulated mean of the means.
3. Find the simulated standard deviation of the means.
4. Draw the histogram of the means.

Answer 1

R code with comments (all statements starting with # are comments)

#part 1
set.seed(123)
#set the number of experiments
N<-50000
#set the probabilities of Y
py<-c(0.35,0.30,0.25,0.05,0.03,0.02)
#sample 3 Ys N times
y<-sample(1:6,size=3*N,prob=py,replace=TRUE)
#convert y to a matrix of N rows and 3 columns
#(indicating the roll of 3 dice, N times)
y<-matrix(y,nrow=N,ncol=3)
#sum each row to get N, Xs
x<-apply(y,1,sum)
#A. Find the simulated probability mass function (pmf) of random variable X.
#get unique values of X
X<-sort(unique(x))
#initialize the variable pmf
pmfx<-numeric(length(X))
#calculate the pmf
for (i in 1:length(X)){
pmfx[i]<-sum(x==X[i])/N
}
print(data.frame(X,pmfx))
#B.Find the simulated mean of the random variable X.
sprintf('The simulated mean of X is %.4f',mean(x))
#C. Find the simulated standard deviation of random variable X.
sprintf('The simulated standard deviation of X is %.4f',sd(x))
#D. Find the simulated P(2 < X < 13).
sprintf('The simulated P(2 < X < 13)is %.4f',sum(2<x&x<13)/N)
#E. Draw the histogram of the pmf of X.
barplot(pmfx,names=X,main="Histogram of pmf of X",xlab="X",ylab="Density")

# get this output

get this plot

part 2)

R code

#Generate 2000 samples of size 50000 values
#for the random X from distribution in part1.
set.seed(123)
R<-2000
xbar<-numeric(R)
for (i in 1:R) {
   #sample 3 Ys N times
   y<-sample(1:6,size=3*N,prob=py,replace=TRUE)
   #convert y to a matrix of N rows and 3 columns
   #(indicating the roll of 3 dice, N times)
   y<-matrix(y,nrow=N,ncol=3)
   #sum each row to get N, Xs
   x<-apply(y,1,sum)
   #get the mean
   xbar[i]<-mean(x)
}
#a) the simulated probability that the mean is between 6.49 and 6.53 exclusive
sprintf('The simulated probability that the mean is between 6.49 and 6.53 exclusive is %.4f',sum(6.49<xbar&xbar<6.53)/R)
#b) Find the simulated mean of the means.
sprintf('The simulated mean of the means is %.4f',mean(xbar))
#c) Find the simulated standard deviation of the means.
sprintf('The simulated standard deviation of the means is %.4f',sd(xbar))
#d. Draw the histogram of the means.
hist(xbar,main="Histogram of means",xlab="Means")

#get this output

get this plot

As predicted by the CLT, the histogram of means is normally distributed (bell shape), even though the population (the distribution of X from part 1) is not normally distributed.