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

# stat 131, spring 2020 (david draper; 22 may 2020)

# R code for watching what happens as the parameters are varied
# in the binomial and poisson discrete distributions

# i'm using the histogram representation of the PMFs
# of these distributions here, because spike plots
# are more tedious to make in R than histograms

# i simulated M = 1,000,000 IID draws from each distribution
# and make a histogram of the resulting draws

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

# preamble

set.seed( 8237424 )

M <- 1000000

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

# (1) binomial, with the tay-sachs case study parameters
# ( n = 5, p = 0.25 ) defining the base distribution

# first, let's hold p constant at 0.25 and make n grow from 5 to 200

par( mfrow = c( 4, 1 ) )

p <- 0.25

n <- 5

x.star <- rbinom( M, n, p ) 

hist( x.star, breaks = ( -0.5 ):( n + 0.5 ),
  prob = T, main = 'Binomial ( n = 25, p = 0.25 )', xlab = 'x',
  xlim = c( -0.5, 70 ), ylim = c( 0, 0.4 ) )

n <- 25

x.star <- rbinom( M, n, p ) 

hist( x.star, breaks = ( -0.5 ):( n + 0.5 ),
  prob = T, main = 'Binomial ( n = 25, p = 0.25 )', xlab = 'x',
  xlim = c( -0.5, 70 ), ylim = c( 0, 0.4 ) )

n <- 50

x.star <- rbinom( M, n, p ) 

hist( x.star, breaks = ( -0.5 ):( n + 0.5 ),
  prob = T, main = 'Binomial ( n = 25, p = 0.25 )', xlab = 'x',
  xlim = c( -0.5, 70 ), ylim = c( 0, 0.4 ) )

n <- 200

x.star <- rbinom( M, n, p ) 

hist( x.star, breaks = ( -0.5 ):( n + 0.5 ),
  prob = T, main = 'Binomial ( n = 25, p = 0.25 )', xlab = 'x',
  xlim = c( -0.5, 70 ), ylim = c( 0, 0.4 ) )

# as p is held constant and n grows,

#   --> what happens to the center of the distribution?

# the expected value of a Binomial( n, p ) distribution is 

#   --> what happens to the spread of the distribution?

# the standard deviation of a Binomial( n, p ) distribution is 

#   --> what happens to the shape of the distribution?

# 

par( mfrow = c( 1, 1 ) )

# next, let's hold n constant at 5 and make p grow from 0.01 to 0.99

par( mfrow = c( 5, 1 ) )

p <- 0.01

n <- 5

x.star <- rbinom( M, n, p ) 

hist( x.star, breaks = ( -0.5 ):( n + 0.5 ),
  prob = T, main = 'Binomial ( n = 5, p = 0.01 )', xlab = 'x',
  xlim = c( -0.5, n + 0.5 ) )

p <- 0.25

x.star <- rbinom( M, n, p ) 

hist( x.star, breaks = ( -0.5 ):( n + 0.5 ),
  prob = T, main = 'Binomial ( n = 5, p = 0.25 )', xlab = 'x',
  xlim = c( -0.5, n + 0.5 ) )

p <- 0.5

x.star <- rbinom( M, n, p ) 

hist( x.star, breaks = ( -0.5 ):( n + 0.5 ),
  prob = T, main = 'Binomial ( n = 5, p = 0.5 )', xlab = 'x',
  xlim = c( -0.5, n + 0.5 ) )

p <- 0.75

x.star <- rbinom( M, n, p ) 

hist( x.star, breaks = ( -0.5 ):( n + 0.5 ),
  prob = T, main = 'Binomial ( n = 5, p = 0.75 )', xlab = 'x',
  xlim = c( -0.5, n + 0.5 ) )

p <- 0.99

x.star <- rbinom( M, n, p ) 

hist( x.star, breaks = ( -0.5 ):( n + 0.5 ),
  prob = T, main = 'Binomial ( n = 5, p = 0.99 )', xlab = 'x',
  xlim = c( -0.5, n + 0.5 ) )

# as n is held constant and p grows,

#   --> what happens to the center of the distribution?

# the expected value of a Binomial( n, p ) distribution is 

#   --> what happens to the spread of the distribution?

# the standard deviation of a Binomial( n, p ) distribution is 

#   --> what happens to the shape of the distribution?

par( mfrow = c( 1, 1 ) )

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

# (2) poisson, as lambda ranges from 0.1 to 25

par( mfrow = c( 3, 1 ) )

lambda <- 0.1

x.star <- rpois( M, lambda ) 

hist( x.star[ x.star <= 45 ], breaks = -0.5: 45.5,
  prob = T, main = 'Poisson ( lambda = 0.1 )', xlab = 'x',
  xlim = c( -0.5, 40 ), ylim = c( 0, 0.9 ) )

lambda <- 5

x.star <- rpois( M, lambda ) 

hist( x.star[ x.star <= 45 ], breaks = -0.5: 45.5,
  prob = T, main = 'Poisson ( lambda = 5 )', xlab = 'x',
  xlim = c( -0.5, 40 ), ylim = c( 0, 0.9 ) )

lambda <- 25

x.star <- rpois( M, lambda ) 

hist( x.star[ x.star <= 45 ], breaks = -0.5:45.5,
  prob = T, main = 'Poisson ( lambda = 25 )', xlab = 'x',
  xlim = c( -0.5, 40 ), ylim = c( 0, 0.9 ) )

# as lambda > 0 grows,

#   --> what happens to the center of the distribution?
#       the center also grows

# the expected value of a Poisson( lambda ) distribution is lambda

#   --> what happens to the spread of the distribution?
#       the spread also grows

# the standard deviation of a Poisson( lambda ) distribution is sqrt( lambda )

#   --> what happens to the shape of the distribution?
#       the skewness diminishes: the distribution becomes more symmetric

# the skewness of a Poisson( lambda ) distribution is 1 / sqrt( lambda )

par( mfrow = c( 1, 1 ) )

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