r - 伽玛密度图：dgamma 工作，而自己的函数返回错误

Question

我正在尝试为参数（形状）alfa=3457 和（速率）beta=84 绘制 Gamma 密度函数。如果我这样做使用：

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 告诉我 gamma(s) 返回 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)