próbkowania od odwrotny rozkład gamma w R

Question

W celu próbki z odwrotny rozkład gamma w R, jest następujący poprawny sposób to zrobić:

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

a = 4 
b = 9 

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

Answer 1

Chociaż @Dason i @Stephane już zauważyli, że Twój podejście to jest ważne, istnieje kilka pakietów w R, które to zrobić (znaleziono googlowania dla r inverse gamma:

• rigamma {pscl},
• rinvgamma {MCMCpack}

Patrz także stronę wikipedia dla gamma distribution i inverse gamma distribution dla funkcji gęstości prawdopodobieństwa obu rozkładów:

do dystrybucji gamma w stosunku:

dla osób odwrotna gamma.

Answer 2

Poniższy kod jest przykładem do porównania symulacji z odwrotnej gammy z różnych pakietów R @ user2005253 i @Stephane.

@Paul Hiemstra nie jestem pewien co do 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