my.dgamma.log()generates the log of the probability density of the Gamma distribution. Use?lgammato understand what it does. Note that this is exactly equal to the functiondgamma(..., log=TRUE)so use this to verify your answers.

The function to be debugged:

my.dgamma.log <- function(xx, shape, rate) { shape*log(rate) - lgamma(shape) + (shape-1)*xx - rate*xx } 

The testing code:

my.dgamma.log(seq(0.2, 5, by=0.2), shape=0.2, rate=5) dgamma(seq(0.2, 5, by=0.2), shape=0.2, rate=5, log=TRUE) my.dgamma.log (seq(0.4, 10, by=0.4), shape=0.2, rate=5) my.dgamma.log (seq(1, 10, by=0.2), shape=2, rate=6)