График гамма-плотности: dgamma работает, а собственная функция возвращает ошибку

Question

Я пытаюсь построить функцию плотности гамма для параметров (форма) альфа =3457 и (скорость) бета =84. Если я сделаю это с помощью:

curve(dgamma(x,shape=3457,rate=84),from=35,to=50,xlab="posterior theta",ylab="density")

все работает просто отлично, и я получаю плотность по центру ~41. Хотя, если я сначала построю плотность как:

para_density = function(x,s=3457,r=84,N) ((r^s)/(gamma(s)))*x^(s-1)*exp(-r*x)

и затем нанесите на карту это с:

curve(para_density(x,s=3457,r=84,N=N),from=0,to=100)

R говорит мне, что гамма (ы) возвращает Inf, что разумно, учитывая его значение.

Почему же тогда я могу построить его с одинаковой параметризацией, используя функции plot и dgamma?

Спасибо всем заранее.

Answer 1

Ну, тогда просто используйте логарифм

Код (не проверено!)

mygamma <- function(x, s, r) {
    l <- s*log(r) - lgamma(s) + (s-1.0)*log(x) - x*r
    exp(l)
}

ОБНОВИТЬ

Да, это работает

library(ggplot2)

p <- ggplot(data = data.frame(x = 0), mapping = aes(x = x))
p <- p + stat_function(fun = function(x) dgamma(x,shape=3457,rate=84))
p <- p + stat_function(fun = function(x) mygamma(x,s=3457,r=84))
p <- p + xlim(35.0, 50.0)
print(p)