inverse.gamma.density <- function( lambda, alpha, beta ) {

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

  return( exp( log.density ) )

}