伽马密度图：dgamma 在自己的函数返回错误时工作

Question

我正在尝试为参数 (shape) alfa=3457 和 (rate) 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) returns 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)

伽马密度图：dgamma 在自己的函数返回错误时工作

Gamma density plot: dgamma working while own function returning error

plot

r

gamma-function