r - 从 R 中的逆伽马分布采样

Question

为了从 R 中的逆伽马分布中采样，以下是正确的方法：

#I want to sample an inverse-gamma(a,b)

a = 4
b = 9

x = 1/rgamma(1,a,b)

Answer 1

尽管@Dason 和@Stephane 已经评论说您的方法是有效的，但是 R 中有几个包可以做到这一点（发现谷歌搜索r inverse gamma：

• rigamma {pscl} ,
• rinvgamma {MCMCpack}

另请参阅 wikipedia 页面以了解gamma 分布和两种分布的概率密度函数的逆 gamma分布：

对于伽马分布与：

为逆伽玛。

Answer 2

下面的代码是一个示例，用于比较来自各种 R 包 @user2005253 和 @Stephane 的反伽马的模拟。

@Paul Hiemstra 我不确定 ringvamma{MCMCpack}

# double check implementations from various packages
library(ggplot2)
alpha = 1
rate = 0.5

# stats library ----------------------------------
library(stats) 
x.base<- 1/rgamma(10000, shape = alpha, rate = rate)
x11()
p.try0<- ggplot(data.frame(x = x.base), aes(x=x)) + geom_density() +
ggtitle(paste("Stats package: shape", alpha, "rate ", rate)) + xlim(c(0, 3))
p.try0

# invgamma library -------------------------------
library(invgamma)
sims.1<- rinvgamma(10000, shape = alpha, rate = rate)
p.try1<- ggplot(data.frame(x = sims.1), aes(x=x)) + geom_density() +
  ggtitle(paste("Package (invgamma) shape", alpha, " rate ", rate, sep = ""))+
  xlim(c(0, 3))
x11()
p.try1

detach("package:invgamma", unload = TRUE) 

# MCMCpack library -------------------------------
library(MCMCpack) # no rate argument - this works only on shape and scale.
                  #That's ok since scale = 1/rate
sims.2<- rinvgamma(10000, shape = alpha, scale = 1/rate)

p.try2<- ggplot(data.frame(x = sims.2), aes(x=x)) + geom_density() + 
  ggtitle(paste("Package MCMCpack: shape", alpha, " scale", 1/rate, sep = "")) +   
  xlim(c(0, 3))

x11()
p.try2

# Implementation of rinvgamma incorrect for MCMC pack? Because this works with
sims.3<- rinvgamma(10000, shape = alpha, scale = rate)

p.try3<- ggplot(data.frame(x = sims.2), aes(x=x)) + geom_density() +
   ggtitle(paste("again MCMCpack: here scale = rate ???")) + xlim(c(0, 3))

x11()
p.try3