# Exit Poll (N = 3889)
X<-c() # vote outcome, with x = 1 for Jerry Brown and x = 0 for all other responses
N=0 # sample size for the Exit Poll (break at N = 3889)
set.seed(123456) # ...just for replication purpose

repeat{
N=N+1
x <- sample(c(0,1),size=1,replace=TRUE,prob=c(0.462,0.538)) # sample from population distribution
X<-c(X,x)
if(N==3889)break()
}

X

# number of votes for Brown
sum(X)
# sample proportion for Brown
sum(X)/3889
#--------------------------------------------------
#--------------------------------------------------
#--------------------------------------------------
#--------------------------------------------------
#--------------------------------------------------
#--------------------------------------------------
#--------------------------------------------------
# Repeat 1,000 Exit Polls on N = 3,889 participants

sampling.distribution<-c() # the 1,000 random sample proportion "for Brown"
R = 0 # number of replication
set.seed(123456)
repeat{
R=R+1

# Exit Poll (N = 3889)	    
		X<-c()
		N=0
set.seed=123456
		repeat{
		N=N+1
		x <- sample(c(0,1),size=1,replace=TRUE,prob=c(0.462,0.538)) # population probabilites
		X<-c(X,x)
		if(N==3889)break()
		}

		# sample proportion for Brown
sampling.distribution<-c(sampling.distribution,sum(X)/3889)
cat("R=",R,"\n")
if(R==1000)break()
	
	
}
#--------------------------------------------------
# Histogram of sampling distribution of proportion:
hist(sampling.distribution,freq=FALSE,col="blue",border="white")
points(density(sampling.distribution),type="l",lwd=2,col="gray80")
mean(sampling.distribution) # p = 0.538 
sd(sampling.distribution) # sqrt(p*(1-p)/3889) = 0.007994534

#--------------------------------------------------
#--------------------------------------------------
#--------------------------------------------------
#--------------------------------------------------
#--------------------------------------------------
#--------------------------------------------------
#--------------------------------------------------

###############################################
#     Central Limit Theorem on binomial       #
#											  #
# Repeat 1,000 Exit Polls on "N" participants #
###############################################

sampling.distribution<-c() # the 1,000 random sample proportion "for Brown"
R = 0 # number of replication
set.seed(123456)
repeat{
R=R+1
N = 2 # ...try to change N value
# Sample one value from binomial with parameters N and p = 0.538
		X<-rbinom(1,size=N,prob=0.538)

		# sample proportion for Brown
sampling.distribution<-c(sampling.distribution,X/N)
cat("R=",R,"\n")
if(R==1000)break()
	
	
}
###################################################
#### Briefly ######################################
###################################################
# R = 1000 # 1,000 replicates
# N = 30 # ...try to change N value
# set.seed(123456)

# sampling.distribution<-rbinom(R,N,prob=0.538)/N
###################################################
#--------------------------------------------------
# Histogram of sampling distribution of proportion:
plot(density(sampling.distribution),type="l",lwd=2,col="red",bty="n",xlim=c(0,1))
polygon(density(sampling.distribution), col = rgb(1, 0, 0, alpha = 0.5))
abline(v=0.538,lwd=2)
mean(sampling.distribution)
p = 0.538 ;p
sd(sampling.distribution) 
sqrt(p*(1-p)/N)
#--------------------------------------------------
#--------------------------------------------------
#--------------------------------------------------
#--------------------------------------------------
#--------------------------------------------------
#--------------------------------------------------
#--------------------------------------------------


###############################################
#     Central Limit Theorem on uniform        #
#			           			              #
# Repeat 1,000 Uniform on "N" participants    #
############################################### 

sampling.distribution<-c()
R = 0 # number of replication
set.seed(123456)
repeat{
R=R+1
N = 2 # ...try to change N value
# Random sample from uniform distribution (max = 20, min = 1) with p = 1/20    
		X<-sample(c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20),size=N,replace = TRUE,prob=c(rep(1/20,20)))

sampling.distribution<-c(sampling.distribution,mean(X))
cat("R=",R,"\n")
if(R==1000)break()
	
	
}
#--------------------------------------------------
# Histogram of sampling distribution of proportion:
plot(density(sampling.distribution),type="l",lwd=2,col="red",bty="n",xlim=c(0,20))
polygon(density(sampling.distribution), col = rgb(1, 0, 0, alpha = 0.5))
mean(sampling.distribution)
mu = (20+1)/2 ;mu
sd(sampling.distribution) 
sqrt( (20-1)^2/12 *1/N) #standard error!!!
#--------------------------------------------------



###############################################
#     Central Limit Theorem on "bimodal"      #
#							    #
# Repeat 1,000 Bimodal on "N" participants    #
###############################################

sampling.distribution<-c()
R = 0 # number of replication
set.seed(123456)
repeat{
R=R+1
N = 2 # ...try to change N value
# Random sample from bimodal distribution (max = 20, min = 1) with p = 1/20    

		X<-sample(c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20),size=N,replace = TRUE,prob=c(1/3*2/3,1/3*1/3,rep((1/3)*1/16,16),1/3*1/3,1/3*2/3))

sampling.distribution<-c(sampling.distribution,mean(X))
cat("R=",R,"\n")
if(R==1000)break()
	
	
}
#--------------------------------------------------
# Histogram of sampling distribution of proportion:
plot(density(sampling.distribution),type="l",lwd=2,col="red",bty="n",xlim=c(0,20))
polygon(density(sampling.distribution), col = rgb(1, 0, 0, alpha = 0.5))
mean(sampling.distribution)
#--------------------------------------------------
X=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20)
prob=c(1/3*2/3,1/3*1/3,rep((1/3)*1/16,16),1/3*1/3,1/3*2/3)
mu = sum(X*prob) ;mu
sd(sampling.distribution) 
sqrt( sum((X-mu)^2*prob) *1/N ) #standard error!!!
#--------------------------------------------------