Lấy mẫu từ một phân phối gamma nghịch đảo trong R

Question

Để lấy mẫu từ một phân phối gamma nghịch đảo trong R, là sau cách chính xác để làm điều đó:

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

a = 4 
b = 9 

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

Answer 1

Mặc dù @Dason và @Stephane đã nhận xét rằng bạn cách tiếp cận có giá trị, có một số gói trong R mà làm điều này (tìm thấy googling cho r inverse gamma:

• rigamma {pscl},
• rinvgamma {MCMCpack}

cũng Xem trang wikipedia cho gamma distribution và inverse gamma distribution cho hàm mật độ xác suất của cả các bản phân phối:

cho việc phân phối gamma so:

cho gamma nghịch đảo.

Answer 2

Mã dưới đây là ví dụ để so sánh các mô phỏng từ gamma nghịch đảo từ các gói R @ user2005253 và @Stephane.

@ Paul Hiemstra Tôi không chắc chắn về 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