#
	for(j in 1:100){#
     for(i in 1:reps){#
	      store[i]<-sum(sample(pop,n))#
	}#
	if(k==667){k<-666}#
store2[j]<-as.numeric(summary(store==k)[3])/n#
}#
#
print(mean(store2))
white.balls <-rep(1,8000)#
red.balls <-rep(0,4000)#
#
pop<-c(white.balls,red.balls)#
#
## number of repeated samples#
reps<-100#
#
## for different n's:#
n<-100#
k<-ceiling(2/3*n)#
# store the results of repeated samples#
store<- rep(NA,reps) #
#
store2<- rep(NA,100)#
#
	for(j in 1:100){#
     for(i in 1:reps){#
	      store[i]<-sum(sample(pop,n))#
	}#
store2[j]<-as.numeric(summary(store==k)[3])/n}#
#
print(mean(store2))
dbinom(67,100,2.3)
dbinom(67,100,2/3)
white.balls <-rep(1,80000)#
red.balls <-rep(0,40000)#
#
pop<-c(white.balls,red.balls)#
#
## number of repeated samples#
reps<-100#
#
## for different n's:#
n<-100#
k<-ceiling(2/3*n)#
# store the results of repeated samples#
store<- rep(NA,reps) #
#
store2<- rep(NA,100)#
#
	for(j in 1:100){#
     for(i in 1:reps){#
	      store[i]<-sum(sample(pop,n))#
	}#
store2[j]<-as.numeric(summary(store==k)[3])/n}#
#
print(mean(store2))#
#
dbinom(67,100,2/3)
?sample
white.balls <-rep(1,80000)#
red.balls <-rep(0,40000)#
#
pop<-c(white.balls,red.balls)#
#
## number of repeated samples#
reps<-100#
#
## for different n's:#
n<-100#
k<-ceiling(2/3*n)#
# store the results of repeated samples#
store<- rep(NA,reps) #
#
store2<- rep(NA,100)#
#
	for(j in 1:100){#
     for(i in 1:reps){#
	      store[i]<-sum(sample(pop,n,replace=T))#
	}#
store2[j]<-as.numeric(summary(store==k)[3])/n}#
#
print(mean(store2))#
#
dbinom(67,100,2/3)
white.balls <-rep(1,80000)#
red.balls <-rep(0,40000)#
#
pop<-c(white.balls,red.balls)#
#
## number of repeated samples#
reps<-100#
#
## for different n's:#
n<-1000#
k<-ceiling(2/3*n)#
# store the results of repeated samples#
store<- rep(NA,reps) #
#
store2<- rep(NA,100)#
#
	for(j in 1:100){#
     for(i in 1:reps){#
	      store[i]<-sum(sample(pop,n,replace=T))#
	}#
store2[j]<-as.numeric(summary(store==k)[3])/n}#
#
print(mean(store2))#
#
dbinom(67,100,2/3)
white.balls <-rep(1,80000)#
red.balls <-rep(0,40000)#
#
pop<-c(white.balls,red.balls)#
#
## number of repeated samples#
reps<-100#
#
## for different n's:#
n<-1000#
k<-ceiling(2/3*n)#
# store the results of repeated samples#
store<- rep(NA,reps) #
#
store2<- rep(NA,100)#
#
	for(j in 1:100){#
     for(i in 1:reps){#
	      store[i]<-sum(sample(pop,n,replace=T))#
	}#
store2[j]<-as.numeric(summary(store==k)[3])/n}#
#
print(mean(store2))#
#
dbinom(667,1000,2/3)
white.balls <-rep(1,8000)#
red.balls <-rep(0,4000)#
#
pop<-c(white.balls,red.balls)#
#
## number of repeated samples#
reps<-100#
#
## for different n's:#
n<-100#
k<-ceiling(2/3*n)#
# store the results of repeated samples#
store<- rep(NA,reps) #
#
store2<- rep(NA,100)#
#
	for(j in 1:100){#
     for(i in 1:reps){#
	      store[i]<-sum(sample(pop,n))#
	}#
store2[j]<-as.numeric(summary(store==k)[3])/n}#
#
print(mean(store2))#
#
dbinom(67,100,2/3)
sample(pop)
sample(pop,12000)
sample(pop,100)
sum(sample(pop,11999))/12000
sum(sample(pop,12000))/12000
length(pop)
sum(sample(pop,12000))/12000
numballs <- 40 #
p <- 0.5      #
probs <- rep(NA,numballs)
probs
for(k in 1:numballs){#
   probs[k] <- binomialprobability(numballs,p,k)#
   }
binomialprobability <- function(n,p,k){#
 choose(n,k)*p^k*(1-p)^(n-k)#
}
for(k in 1:numballs){#
   probs[k] <- binomialprobability(numballs,p,k)#
   }
probs
probs2<-rep(NA,40)#
#
for(i in 1:40){#
  probs2[i]<-dbinom(i,40,0.5)#
}
(withinone <- sum(probs[19:21]))
(withintwo <- sum(probs[18:22]))
(withinthree <- sum(probs[17:23]))
(withinfour <- sum(probs[16:24]))
(withinfive <- sum(probs[15:25]))
(withinsix <- sum(probs[14:26]))
mean.index <- 20#
intervals <- rep(NA,19)#
#
for(i in 1:19){#
  indices <- seq(mean.index-i,mean.index+i,by=1)#
  range <- probs[indices]#
  intervals[i] <- sum(range)#
}#
#
conf.intervals <- data.frame(margin=rep(1:19),probability=intervals)#
#
conf.intervals40 <- conf.intervals#
#
print(head(conf.intervals))
plot(conf.intervals$margin,conf.intervals$probability,#
     type="b",xlab="Margins",ylab="Probability",main="Sample size 40")#
#
segments(0,.95,5.7,.95,col="red")#
segments(5.7,0,5.7,.95,col="red")
numballs <- 400 #
p <- 0.5      #
probs <- rep(NA,numballs)#
for(k in 1:numballs){#
   currentk <- binomialprobability(numballs,p,k)#
   probs[k] <- currentk#
   }#
#
mean.index <- 200#
intervals <- rep(NA,199)#
#
for(i in 1:199){#
  indices <- seq(mean.index-i,mean.index+i,by=1)#
  range <- probs[indices]#
  intervals[i] <- sum(range)#
}#
#
conf.intervals <- data.frame(margin=rep(1:199),probability=intervals)#
#
conf.intervals400 <- conf.intervals#
#
print(head(conf.intervals))
plot(conf.intervals$margin,conf.intervals$probability,type="b",xlab="Margins",ylab="Probability",main="Sample size 400")
segments(-6,.95,19,.95,col="red")#
segments(19,0,19,.95,col="red")
## computes the margins for any number of balls and any p:#
compute_margins <- function(numballs,p){#
  probs <- rep(NA,numballs)#
 for(k in 1:numballs){#
    currentk <- binomialprobability(numballs,p,k)#
    probs[k] <- currentk#
    }#
#
mean.index <- numballs*p#
max.margin <- numballs*p-1#
#max.margin <- numballs#
intervals <- rep(NA,max.margin)#
#
for(i in 1:max.margin){#
  indices <- seq(mean.index-i,mean.index+i,by=1)#
  range <- probs[indices]#
  intervals[i] <- sum(range)#
}#
#
conf.intervals <- data.frame(margin=rep(1:max.margin),probability=intervals)#
#
return(conf.intervals)#
}
plotcis <- function(numballs,p,color="black",margin,maintitle,interval=TRUE){#
probs <- rep(NA,numballs+1)#
for(k in 0:numballs){#
currentk <- binomialprobability(numballs,p,k)#
  probs[k+1] <- currentk#
}#
proportions <- 0:numballs/numballs#
#
plot(proportions,probs,type="l",col="black",xlab="Proportions",#
     ylab="Probability",main=maintitle)#
#
if(interval==TRUE){#
# draw 95% confidence intervals#
segments(proportions[(numballs/2+1)-margin],-0.5,#
         proportions[(numballs/2+1)-margin],.06,col=color,lty=1,lwd=2)#
segments(proportions[(numballs/2+1)+margin],-0.5,#
         proportions[(numballs/2+1)+margin],.06,col=color,lty=1,lwd=2)#
}#
}
## as smple size goes up, CIs become tighter and tighter:#
op <- par(mfrow=c(1,2),pty = "s")#
plotcis(40,.5,margin=5,maintitle="Sample size 40")#
plotcis(400,.5,margin=19,maintitle="Sample size 400")
sums<-rep(NA,21)#
#
for(i in 0:20){#
  sums[i+1]<-dbinom(i,40,0.5)#
}#
#
print(sums)#
#
sum(sums)#
#
####################################################
### chunk number 38: #
####################################################
#
sum(dbinom(0:20,40,0.5))#
#
####################################################
### chunk number 39: #
####################################################
#
pbinom(20,40,0.5)
newfunction <- function(x,mu,sigma){#
  1/(sqrt(2*pi)*sigma)*exp(1)^(-((x - mu)^2/(2*sigma^2)))}#
#
####################################################
### chunk number 41: #
####################################################
#
plotcis(40,.5,40,margin=20,#
        maintitle="Comparing the binomial and normal distributions",#
        interval=FALSE)
dbinom(400,600,2/3)
dbinom(667,1000,2/3)
rbinom(1,1,0.5)
rbinom(3,1,0.5)
rbinom(1,2,0.5)
rbinom(1,3,0.5)
rbinom(1,4,0.5)
rbinom(1,1,0.5)
for(i in 1:100){#
	print(mean(rbinom(3,1,0.5)))#
	}
for(i in 1:100){#
	print(sum(rbinom(3,1,0.5)))#
	}
st<-rep(NA,100)#
#
for(i in 1:100){#
	st[i]<-print(sum(rbinom(3,1,0.5)))#
	}
st
as.factor(st)
summary(as.factor(st))
summary(as.factor(st))/100
0.5^3
rbinom(1,1,0.5)
rbinom(3,1,0.5)
sum(rbinom(3,1,0.5))
rbinom(3,1,0.5)
st
12/100
0.5^3
summary(st)
summary(as.factor(st))
summary(as.factor(st))/100
0:1
mean(0:1)
as.factor(0:1)
means(as.factor(0:1))
mean(as.factor(0:1))
white.balls <-rep(1,8000)#
red.balls <-rep(0,4000)#
#
pop<-c(white.balls,red.balls)
pop
sample(pop,10)
white.balls <-rep(1,8000)#
red.balls <-rep(0,4000)#
#
pop<-c(white.balls,red.balls)
sample(pop,100)
sum(sample(pop,100))
green.balls<-rep(2,4000)#
white.balls <-rep(1,8000)#
red.balls <-rep(0,4000)#
#
pop<-c(white.balls,red.balls)#
pop2<-c(green.balls,white.balls,red.balls)
pop2
sum(sample(pop2,100))
sample(pop2,100)
sum(sample(pop,100))
sample(pop,100)
sum(sample(pop,100))
sum(sample(pop,100))/100
sample(pop2,100)
summary(sample(pop2,100))
summary(as.factor(sample(pop2,100)))
summary(as.factor(sample(pop,100)))
op <- par(mfrow=c(2,3),pty = "s")#
#
#
## Observe 40 drops 15,25,... times#
####################################################
### chunk number 14: #
####################################################
#
size <- 1#
p <- 0.5#
k <- 40#
# create a vector "observations"#
observations <- c(15,25,50,100,500,1000)#
n <- length(observations)#
#
####################################################
### chunk number 15: #
####################################################
#
for(j in 1:n){#
  results <- rep(NA,observations[j])#
 #
 for(i in 1:observations[j]){#
  results[i] <- sum(rbinom(k,size,p))#
 } ## end-for#
 #
  # create title on the fly:#
  title <- paste(c("Num Obs.",observations[j],sep=" ")) #
  #
  # plot histogram:#
  hist(results,#
       ylab="Frequency", #
       xlab="No. of R-stone hits",      #
       main=title)#
  Sys.sleep(3)#
#
}#end-for
size <- 1   #
p <- 0.5#
drops <- c(4,40,400,4000)#
#
op <- par(mfrow=c(2,2),pty = "s")#
#
for(num.drops in drops) {#
results <- rep(NA,100)#
for (i in 1:100){#
  results[i] <- sum(rbinom(num.drops,size,p)) #
  }#end-for#
#
maintitle <- paste(num.drops,"drops",sep=" ")#
#
hist(results,#
     xlim=range(c(0:num.drops)),#
     xlab="Number of R-stone hits",main=maintitle)#
}#end-for
drops<-rep(1:400,1)#
standard.dev <- rep(NA,400)#
#
for(j in 1:400) {#
   results <- rep(NA,100)#
   for (i in 1:100){#
      results[i] <- sum(rbinom(j,size,p)) #
      }#
   standard.dev[j] <- sd(results)#
}#
#
plot(drops,standard.dev,#
     xlim=c(1,400),#
     xlab="Number of drops",#
     ylab="Standard deviation",#
     main=expression("SD "*italic("increases")*" as we increase sample size."))
#
## for any number of drops ranging from 1 to 400:#
drops<-rep(1:400,1)#
standard.dev <- rep(NA,400)#
#
for(j in 1:400) {#
   results <- rep(NA,100)#
   for (i in 1:100){#
      results[i] <- sum(rbinom(j,size,p)) #
      }#
   standard.dev[j] <- sd(results)#
}#
#
plot(drops,standard.dev,#
     xlim=c(1,400),#
     xlab="Number of drops",#
     ylab="Standard deviation",#
     main=expression("SD "*italic("increases")*" as we increase sample size."))
size <- 1   #
p <- 0.5#
drops <- c(4,40,400,4000)#
#
op <- par(mfrow=c(2,2),pty = "s")#
#
for(num.drops in drops) {#
results <- rep(NA,100)#
for (i in 1:100){#
  results[i] <- mean(rbinom(num.drops,size,p)) #
  }#end-for#
#
maintitle <- paste(num.drops,"drops",sep=" ")#
#
hist(results,#
     xlim=range(c(0:1)),#
     xlab="Proportion of R-stone hits",main=maintitle)#
}#end-for
drops<-rep(1:400,1)#
standard.dev <- rep(NA,400)#
#
for(j in 1:400) {#
   results <- rep(NA,100)#
   for (i in 1:100){#
      results[i] <- mean(rbinom(drops[j],size,p)) #
      }#
   standard.dev[j] <- sd(results)#
}#
#
plot(drops,standard.dev,#
     xlim=c(1,400),#
     xlab="Number of drops",#
     ylab="Standard deviation",#
     main=expression("SD decreases as we increase sample size."))
st<-rep(NA,100)#
#
for(i in 1:100){#
	st[i]<-sum(rbinom(3,1,0.5))#
	}#
#
summary(as.factor(st))/100
