#Below we see some pictures to help explain sum of independent random variables. 

# First you can compare the distribution of the average that one obtains running 

N1=10000
N2=1000
N3=100

# trials with a chance of success 

p=2/3

# You may change the values as you please.  Try to understand the role of  1/sqrt(N) and its relationship to "Diminishing returns".
 
###### Code

X=seq(0,N1,by=1)
Y=dbinom(X,N1,p)
plot(X/N1,N1*Y,cex=.5,xlab="Average",ylab="Probability Density",main=paste("The average as N increases, p = ",round(p,2),sep=""))
X=seq(0,N2,by=1)
Y=dbinom(X,N2,p)
points(X/N2,N2*Y,col="red",cex=.5)
X=seq(0,N3,by=1)
Y=dbinom(X,N3,p)
points(X/N3,N3*Y,col="blue",cex=.5)

#############


#Next we compare the standardized sums.  Once again we compare the situation for different N

N1=10000
N2=100
N3=10

# with chance of success 

p=2/3

#In this case we are exploring the Central Limit Theorem, that asserts that the distribution of these standardized sums will approach the standard normal as N increases. 

############ Codee

#We need to choose a window of dev standard deviations form the mean
 
devs=3.5

X=seq(0,N1)
X2<-c(X-1/2,X+1/2)
dim(X2)=c(N1+1,2)
X2=aperm(X2,c(2,1))
dim(X2)=c(1,2*(N1+1))
Y=dbinom(X,N1,p)
Y2<-c(Y,Y)
dim(Y2)=c(N1+1,2)
Y2=aperm(Y2,c(2,1))
dim(Y2)=c(1,2*(N1+1))
Z0=(X2-N1*p)/sqrt(N1*p*(1-p))
Z=Z0[Z0>-devs & Z0<devs]
Y=Y2[Z0>-devs & Z0<devs]
plot(Z,sqrt(N1*p*(1-p))*Y,type="l",col="blue",xlab="Standard deviations from mean",ylab="Probability Density",main=paste("Standardized sum as N increases, p = ",round(p,2),sep=""))

X=seq(0,N2)
X2<-c(X-1/2,X+1/2)
dim(X2)=c(N2+1,2)
X2=aperm(X2,c(2,1))
dim(X2)=c(1,2*(N2+1))
Y=dbinom(X,N2,p)
Y2<-c(Y,Y)
dim(Y2)=c(N2+1,2)
Y2=aperm(Y2,c(2,1))
dim(Y2)=c(1,2*(N2+1))
Z0=(X2-N2*p)/sqrt(N2*p*(1-p))
Z=Z0[Z0>-devs & Z0<devs]
Y=Y2[Z0>-devs & Z0<devs]
points(Z,sqrt(N2*p*(1-p))*Y,type="l",col="red")

X=seq(0,N3)
X2<-c(X-1/2,X+1/2)
dim(X2)=c(N3+1,2)
X2=aperm(X2,c(2,1))
dim(X2)=c(1,2*(N3+1))
Y=dbinom(X,N3,p)
Y2<-c(Y,Y)
dim(Y2)=c(N3+1,2)
Y2=aperm(Y2,c(2,1))
dim(Y2)=c(1,2*(N3+1))
Z0=(X2-N3*p)/sqrt(N3*p*(1-p))
Z=Z0[Z0>-devs & Z0<devs]
Y=Y2[Z0>-devs & Z0<devs]
points(Z,sqrt(N3*p*(1-p))*Y,type="l",col="orange")

# When your done you stick in the standard normal to compare:

X=seq(-devs,devs,by=.01)
Y=dnorm(X,0,1)
points(X,Y,type="l")

##################

#Recall,  for the standardized sum of a success failure experiment one as the rule of thumb: N is big enough to apply the CLT provided   min(p*N,(1-p)*N)>=10. We Explore this for different N and p: 

p=1/20
N=200


########## Code

X=seq(0,N)
X2<-c(X-1/2,X+1/2)
dim(X2)=c(N+1,2)
X2=aperm(X2,c(2,1))
dim(X2)=c(1,2*(N+1))
Y=dbinom(X,N,p)
Y2<-c(Y,Y)
dim(Y2)=c(N+1,2)
Y2=aperm(Y2,c(2,1))
dim(Y2)=c(1,2*(N+1))
Z0=(X2-N*p)/sqrt(N*p*(1-p))
Z=Z0[Z0>-devs & Z0<devs]
Y=Y2[Z0>-devs & Z0<devs]
plot(Z,sqrt(N*p*(1-p))*Y,type="l",col="blue",xlab="Standard deviations from mean",ylab="Probability Density",main=paste("The standardized sum, N=",N,"p=",p,sep=" "))

X=seq(0,N,by=1)
Y=dbinom(X,N,p)
Z0=(X-N*p)/sqrt(N*p*(1-p))
Z=Z0[Z0>-devs & Z0<devs]
Y=Y[Z0>-devs & Z0<devs]
points(Z,sqrt(N*p*(1-p))*Y,col="orange")

X=seq(-devs,devs,by=.01)
Y=dnorm(X,0,1)
points(X,Y,type="l")


############