gamma_from_normal: Simulate A Gamma Multivariate Distribution from a...

View source: R/gamma_from_normal.R

gamma_from_normalR Documentation

Simulate A Gamma Multivariate Distribution from a Multivariate Normal Distribution

Description

To generate a positively skewed multivariate distribution from a multivariate normal distribution, this function uses the rejection sampling method, also known as the accept-reject method.

Usage

gamma_from_normal(shape, rate, mean_vec, cov_matrix, size, c, seed_num)

Arguments

shape

A shape parameter for a gamma distribution.

rate

A rate parameter for a gamma distribution.

cov_matrix

A m x m covariance matrix for an MVN distribution, where m is the number of dimensions and must be equal to the length of the mean vector.

size

is the sample size.

c

This is a constant such that f(x)/g(x) <= c, where f(x) is the density function of the standard normal distribution, and g(x) is a gamma distribution.

seed_num

This is an integer value served as the seed number.

mean_fec

A mean vector for the multivariate normal (MVN) distribution.

References

\insertRef

Braun2016AUTTT

\insertRef

Maria2007AUTTT

\insertRef

Robert2010AUTTT

See Also

dgamma(), dnorm()

Examples

# find the constant
y <- seq(-3, 3, 0.01)
# plotting pdf of a gamma dist
shape <- 5
rate <- 5
par(mfrow=c(1,1))
plot(y, dgamma(y, shape=shape, rate = rate))
plot(y, dgamma(y, shape = shape, rate = rate)/dnorm(y), 
     ylim=c(1,6), 
     xlim=c(0, 3.5), 
     ylab = expression("f(y)/g(y)"), 
     xlab = "A Random Variable Y")
abline(h=4)
# Constant c = 4 is appropriate.
c <- 4

mean_vec <- c(-0.1, 0, 0.1) # mean
sd_vec <- c(0.95, 0.98, 1.1) # sd vector for each dimension
f_cor <- c(0.49, 0.74, 0.87) # correlation vector
# convert the cor. vector to cor.matrix
my_cormat_input <- AUTTT::to_cormatrix(cor_vec = f_cor, n_dim=3)
# covert correlation matrix to covariance matrix
my_covmat_input <- AUTTT::to_covmatrix(cor_matrix = my_cormat_input, sd_vec=sd_vec)
seed_num <- 45679

theta <- gamma_from_normal(shape = shape, 
                  rate = rate, 
                  mean_vec = mean_vec, 
                  cov_matrix = my_covmat_input, 
                  size = 300, 
                  c = c, 
                  seed_num = seed_num)
# see the shape the sample
par(mfrow = c(2,2))
for (i in 1:3){
  y <- sort(theta$X[,i])
  hist(y, prob=TRUE, ylim=c(0, 1))  
  lines(y, dgamma(y, shape = 5, scale=1/4))
}

Boklauth/AUTTT documentation built on Dec. 9, 2022, 7:37 a.m.