The Gamma-Gamma Distribution

Description

Density and random generation for the Gamma-Gamma distribution with parameters shape1, rate1, and shape2.

Usage

1
2
  dggamma(x, shape1, rate1, shape2)
  rggamma(n, shape1, rate1, shape2)

Arguments

x

Vector. Quantiles.

n

Scalar. Number of random variates to generate (sample size).

shape1, rate1

Vector. Shape and rate parameters for y-distribution. Must be strictly positive.

shape2

Vector. Shape parameter for conditional x-distribution. Must be a positive integer.

Details

A Gamma-Gamma distribution with parameters shape1 = a, rate1 = r and shape2 = b has density

f(x) = [(r^a)/(Gamma(a))][Gamma(a+b)/Gamma(b)] [x^(b-1)/(r+x)^(a+b)]

for x > 0 where a,r > 0 and b = 1,2,….

The distribution is generated using the following scheme:

  1. Generate Y ~ Gamma(shape=shape1,rate=rate1).

  2. Generate X ~ Gamma(shape=shape2,rate=Y).

Then, X follows a Gamma-Gamma distribution.

Value

dggamma gives the density, and rggamma gives random variates.

References

Bernardo JM, Smith AFM. (1994) Bayesian Theory. Wiley, New York.

See Also

dgamma

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
############################################################
# Construct a plot of the density function with median and
# quantiles marked.

# define parameters
shape1 <- 4
rate1 <- 4
shape2 <- 20

# construct density plot
x <- seq(0.1,150,0.1)
plot(dggamma(x,shape1,rate1,shape2)~x,
     type="l",lwd=2,main="",xlab="x",ylab="Density f(x)")
     
# determine median and quantiles
set.seed(123)
X <- rggamma(5000,shape1,rate1,shape2)
quants <- quantile(X,prob=c(0.25,0.5,0.75))

# add quantities to plot
abline(v=quants,lty=c(3,2,3),lwd=2)


############################################################
# Consider the following set-up:
#   Let x ~ N(theta,sigma2), sigma2 is unknown variance.
#   Consider a prior on theta and sigma2 defined by
#      theta|sigma2 ~ N(mu,(r*sigma)^2)
#      sigma2 ~ InverseGamma(a/2,b/2), (b/2) = rate.
#
#   We want to generate random variables from the marginal
#   (prior predictive) distribution of the sufficient
#   statistic T = (xbar,s2) where the sample size n = 25.

# define parameters
a <- 4
b <- 4
mu <- 1 
r <- 3
n <- 25


# generate random variables from Gamma-Gamma
set.seed(123)
shape1 <- a/2
rate1 <- b
shape2 <- 0.5*(n-1)

Y <- rggamma(5000,shape1,rate1,shape2)

# generate variables from a non-central t given Y
df <- n+a-1
scale <- (Y+b)*(1/n + r^2)/(n+a-1)

X <- rt(5000,df=df)*sqrt(scale) + mu

# the pair (X,Y) comes from the correct marginal density

# mean of xbar and s2, and xbar*s2
mean(X)
mean(Y)
mean(X*Y)