############Plotting densities to random samples###############

#Binomial example:
#Set number of Monte Carlo samples:
N=1000
#Set binomial parameters:
n=20
theta=0.5
#Generate samples:
x=rbinom(N,size=n,prob=theta)
#Plot frequency chart:
plot(table(x)/N,ylab="Relative frequency",main="Binomial distribution")
#Add theoretical probabilities to above plot
lines(0:n, dbinom(0:n,n,theta), col="blue")

#Now execute the above as a function:
binom.example=function(n=20,theta=0.5,N=1000){
x=rbinom(N,size=n,prob=theta)
plot(table(x)/N,ylab="Relative frequency",main=paste("Binomial, n=",n,"theta=",theta))
lines(0:n, dbinom(0:n,n,theta), col="blue")
                                              }

#generate several cases:
par(mfrow=c(2,2),cex=0.65)
binom.example(N=5000)
binom.example(theta=0.15,N=5000)
binom.example(theta=0.85,N=5000)
binom.example(n=50,N=5000)


#Poisson example:
#Set number of Monte Carlo samples:
N=1000
#Set poisson parameter:
lambda=2
#Generate samples:
x=rpois(N,lambda=lambda)
#Plot frequency chart:
plot(table(x)/N,ylab="Relative frequency",main="Poisson distribution")
#Add theoretical probabilities to above plot
x.max=max(x)+1
lines(0:x.max, dpois(0:x.max,lambda), col="blue")

#Now execute the above as a function:
poisson.example=function(lambda=2,N=1000){
x=rpois(N,lambda=lambda)
plot(table(x)/N,ylab="Relative frequency",main=paste("Poisson, lambda=",lambda))
x.max=max(x)+1
lines(0:x.max, dpois(0:x.max,lambda), col="blue")
                                              }

#generate several cases:
par(mfrow=c(2,2),cex=0.65)
poisson.example(N=5000)
poisson.example(lambda=4,N=5000)
poisson.example(lambda=10,N=5000)
poisson.example(lambda=20,N=5000)