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

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

    return( exp( log.density ) )

}