CODEWITHSUNDEEP
r
Roland Fiagbe
Roland Fiagbe

Reputation: 19

How to simulate data from Gumbel Bivariate Exponential Distribution

I am working on a paper and want to simulate data from Gumbel bivariate exponential distribution type I and type II.

Here's the Gumbel Bivariate distribution: Here's the pdf of Gumbel Bivariate distribution

I tried using the rgumbel function from the Gumbel package, but the package considers lambda1 = lambda2 = 1, but in my case, I want to consider different lambda values.

Also, I tried implementing the inverse CDF uniform sampling approach with the marginals but did not work because the marginals consider correlation=0.

Kindly assist if there is an R package on that or if not, what is the concept behind simulating data from the bivariate exponential distribution from scratch in R?

Upvotes: 1

Views: 530

Answers (1)

jblood94
jblood94

Reputation: 17011

The distribution of x2 conditioned on x1 is a gamma mixture with shape parameters 1 and 2 and weights and a common rate parameter that are a function of x1 and the GBVE parameters. A scheme for sampling is then

  1. Sample x1 from its marginal
  2. Sample x2 conditioned on x1

Implemented in a function (you may want to check my algebra):

rGBVE <- function(n, lambda1, lambda2, theta) {
  x1 <- rexp(n, lambda1)
  lambda12 <- lambda1*lambda2
  pprod <- lambda12*theta
  C <- exp(lambda1*x1)
  A <- (lambda12 - pprod + pprod*lambda1*x1)/C
  B <- (pprod*lambda2 + pprod^2*x1)/C
  D <- lambda2 + pprod*x1
  wExp <- A/D
  wGamma <- B/D^2
  data.frame(x1, x2 = rgamma(n, (runif(n) > wExp/(wExp + wGamma)) + 1, D))
}

This will return n samples in a data.frame.

UPDATE

To arrive at this solution, notice that when x1 is given, the density function can be written in the form (A + B*x2)/exp(logC + D*x2) = (A + B*x2)/(C*exp(D*x2)) = A/(C*D)*D*exp(-D*x2) + B/(C*D^2)*D^2*x2*exp(-D*x2), where A, B, logC, C, and D are constant for a given x1.

D*exp(-D*x2) is the PDF of an exponential distribution with a rate parameter equal to D (or a gamma distribution with a shape parameter equal to 1 and rate parameter equal to D). D^2*x2*exp(-D*x2) is the PDF of a gamma distribution with shape parameter equal to 2 and a rate parameter equal to D. A/(C*D) and B/(C*D^2) are the relative weights of the mixture distribution.

I double-checked my algebra, and I believe it is correct. The function can be streamlined a bit (updated 9/20/2024 based on a comment from JimB):

rGBVE <- function(n, lambda1, lambda2, theta) {
  x1 <- rexp(n, lambda1)
  term <- 1 + lambda1*theta*x1
  data.frame(x1, x2 = rgamma(n, (runif(n) > 1 - theta/term) + 1, lambda2*term))
}

Upvotes: 2

Related Questions