y <- c( 1, 2, 1, 1, 4, 1, 2, 2, 0, 3, 6, 2, 1, 3 )

print( n <- length( y ) )

14

stem( y, scale = 2 )

  The decimal point is at the |

  0 | 0
  1 | 00000
  2 | 0000
  3 | 00
  4 | 0
  5 | 
  6 | 0

setwd( "C:/DD/Teaching/AMS-132/Winter-2017" )

pdf( "homework-1-problem-2-a.pdf" )

hist( y, breaks = ( 0:7 ) - 0.5, prob = T, 
  xlab = 'LoS', main = '' )

dev.off( )

mean( y )

2.071429

sort( y )

# 0 1 1 1 1 1 2 2 2 2 3 3 4 6

print( s <- sum( y ) )

29

print( theta.hat <- s / n )
 
2.071429

alpha <- 0.05

qnorm( 1 - alpha / 2 )

1.959964

print( se.theta.hat <- sqrt( theta.hat / n ) )

0.3846546

print( CI.95 <- c( theta.hat - qnorm( 1 - alpha / 2 ) * se.theta.hat,
  theta.hat + qnorm( 1 - alpha / 2 ) * se.theta.hat ) )

# 1.317519 2.825338

theta.grid <- seq( 1, 4, length = 500 )

likelihood <- theta.grid^s * exp( - n * theta.grid )

par( mfrow = c( 1, 2 ) )

plot( theta.grid, likelihood, type = 'l', xlab = 'theta',
  ylab = 'Likelihood', lwd = 2 )

plot( theta.grid, log( likelihood ), type = 'l', xlab = 'theta',
  ylab = 'Log Likelihood', lwd = 2 )

pdf( "homework-1-problem-2-c.pdf" )

par( mfrow = c( 2, 1 ) )

plot( theta.grid, likelihood, type = 'l', xlab = 'theta',
  ylab = 'Likelihood', lwd = 2 )

plot( theta.grid, log( likelihood ), type = 'l', xlab = 'theta',
  ylab = 'Log Likelihood', lwd = 2 )

dev.off( )

table( y )

y
0 1 2 3 4 6 
1 5 4 2 1 1

table( y ) / n

y
         0          1          2          3          4          6 
0.07142857 0.35714286 0.28571429 0.14285714 0.07142857 0.07142857

dpois( 0:7, theta.hat )

0.126005645 0.261011693 0.270333539 0.186658872 0.096662630 0.040045947
0.013825386 0.004091186

1 - sum( dpois( 0:6, theta.hat ) )

0.005456286

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

alpha <- 0.01

beta <- 0.01

s <- 29

n <- 14

# posterior is gamma( alpha + s, beta + n )
# here this is gamma( 29.01, 14.01 )

theta.grid <- seq( 0, 4, length = 500 )

# we want shape, rate not shape, scale

pdf( "homework-1-problem-2-j.pdf" )

par( mfrow = c( 1, 1 ) )

plot( theta.grid, dgamma( theta.grid, shape = alpha + s, 
  rate = beta + n ), type = 'l', lwd = 2, xlab = 'theta',
  ylab = 'Density' )

# abline( v = ( alpha + s ) / ( beta + n ), col = 'red', lwd = 2 )

lines( theta.grid, dgamma( theta.grid, shape = alpha,
  rate = beta ), lty = 2, lwd = 2, col = 'blue')

lines( theta.grid, dgamma( theta.grid, shape = s + 1,
  rate = n ), lty = 3, lwd = 2, col = 'red')

text( 0.5, 0.8, 'prior', lwd = 2, col = 'blue' )

text( 1.25, 0.9, 'posterior', lwd = 2, col = 'black' )

text( 3, 0.7, 'likelihood', lwd = 2, col = 'red' )

dev.off( )

alpha.posterior <- alpha + s

beta.posterior <- beta + n

print( posterior.mean <- alpha.posterior / beta.posterior )

2.070664

print( theta.hat <- s / n )

2.071429

print( posterior.sd <- sqrt( alpha.posterior ) / beta.posterior )

0.3844463

print( se.theta.hat <- sqrt( theta.hat / n ) )

0.3846546

c( theta.hat - 1.96 * se.theta.hat, theta.hat + 1.96 * se.theta.hat )

c( qgamma( 0.025, alpha + s, beta + n ), 
  qgamma( 0.975, alpha + s, beta + n ) )


          method             estimate    SE/SD     95% interval
 
   large-sample ML            2.0714    0.3847  ( 1.3175, 2.8254 )
Bayes with diffuse prior      2.0707    0.3844  ( 1.3869, 2.8894 )