# unbuffer the output

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

print( epsilon <- 0.001 )

ln.poisson.gamma <- function( y, alpha, beta ) {

  lgamma( alpha + y ) + alpha * log( beta / ( beta + 1 ) ) + 
    y * log( 1 / ( beta + 1 ) ) - lgamma( alpha ) - lgamma( y + 1 )

}

step1 <- function( y, epsilon ) {

  n <- length( y )

  s <- sum( y )

  als <- mean( ln.poisson.gamma( y, epsilon + s, epsilon + n ) )

  return( c( n, s, als ) )

}

print( step1.result <- step1( y, epsilon ) )

step2 <- function( n, s, epsilon, als, m1 ) {

  lambda <- rgamma( m1, epsilon + s, epsilon + n )

  ls <- rep( 0, m1 )

  for ( i in 1:m1 ) {

    y.star <- rpois( n, lambda[ i ] )

    s.star <- sum( y.star )

    ls[ i ] <- mean( ln.poisson.gamma( y.star, epsilon + s.star, 
      epsilon + n ) )

  }

  uata <- sum( ls <= als ) / m1

  write( ls, "ls.out" )

  return( uata )

}

m1 <- 1000

print( step2.result <- step2( step1.result[ 1 ], 
  step1.result[ 2 ], epsilon, step1.result[ 3 ], m1 ) )

v.ls <- scan( "ls.out" )

hist( v.ls, nclass = 20, probability = T,
  main = '', xlab = 'uncalibrated log score' )

abline( v = step1.result[ 3 ] )

step3 <- function( y, epsilon, m1, m2 ) {

  step1.result <- step1( y, epsilon )

  n <- step1.result[ 1 ] 

  s.actual <- step1.result[ 2 ]

  uata <- step2( step1.result[ 1 ], step1.result[ 2 ], epsilon, 
    step1.result[ 3 ], m1 ) 

  v.p <- rep( 0, m2 )

  for ( j in 1:m2 ) {

    lambda.star <- rgamma( 1, epsilon + s.actual, epsilon + n )

    y.sim <- rpois( n, lambda.star )

    step1.result <- step1( y.sim, epsilon )

    v.p[ j ] <- step2( step1.result[ 1 ], step1.result[ 2 ], 
      epsilon, step1.result[ 3 ], m1 )

    if ( ( j %% 10 ) == 0 ) print( j )

  }

  aata <- sum( v.p <= uata ) / m2

  write( v.p, "v.p.out" )

  return( aata )

}

m2 <- 100

print( step3.result <- step3( y, epsilon, m1, m2 ) )

v.p <- scan( "v.p.out" )

hist( v.p, nclass = 20, probability = T, xlim = c( 0, 1 ),
  main = '', xlab = 'calibrated tail areas' )

abline( v = step2.result )

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

n <- 10

e <- rnorm( n, 0.0, 0.5 )

mu <- 0

lambda <- rep( 0, n )

y <- rep( 0, n )

for ( i in 1:n ) {

  lambda[ i ] <- exp( mu + e[ i ] )

  y[ i ] <- rpois( 1, lambda[ i ] )

}

print( y <- sort( y ) )

[1] 0 0 0 1 1 1 2 3 4 4

var( y ) / mean( y )

hist( y, breaks = -0.5 + ( 0:5 ), probability = T, main = '' )

print( step1.result <- step1( y, epsilon ) )

print( step2.result <- step2( step1.result[ 1 ], step1.result[ 2 ], 
  epsilon, step1.result[ 3 ], m1 ) )

v.ls <- scan( "ls.out" )

hist( v.ls, nclass = 20, probability = T,
  main = '', xlab = 'uncalibrated log score' )

abline( v = step1.result[ 3 ] )

m2 <- 1000

print( step3.result <- step3( y, epsilon, m1, m2 ) )

v.p <- scan( "v.p.out" )

hist( v.p, nclass = 20, probability = T, xlim = c( 0, 1 ),
  main = '', xlab = 'calibrated tail areas' )

abline( v = step2.result )