The multivariate gamma distribution is a multivariate absolute continuous probability distribution, defined as the cumulative sum of independent gamma random variables with possibly different shape parameters \(\alpha_i > 0, i\in\{1, \dots, k\}\) and the same scale \(\beta > 0\).
Usage
Multigam(shape = 1, scale = 1)
dmultigam(x, shape, scale, log = FALSE)
rmultigam(n, shape, scale)
# S4 method for class 'Multigam,numeric'
d(distr, x, log = FALSE)
# S4 method for class 'Multigam,matrix'
d(distr, x, log = FALSE)
# S4 method for class 'Multigam,numeric'
r(distr, n)
# S4 method for class 'Multigam'
mean(x)
# S4 method for class 'Multigam'
var(x)
# S4 method for class 'Multigam'
finf(x)
llmultigam(x, shape, scale)
# S4 method for class 'Multigam,matrix'
ll(distr, x)
emultigam(x, type = "mle", ...)
# S4 method for class 'Multigam,matrix'
mle(
distr,
x,
par0 = "same",
method = "L-BFGS-B",
lower = 1e-05,
upper = Inf,
na.rm = FALSE
)
# S4 method for class 'Multigam,matrix'
me(distr, x, na.rm = FALSE)
# S4 method for class 'Multigam,matrix'
same(distr, x, na.rm = FALSE)
vmultigam(shape, scale, type = "mle")
# S4 method for class 'Multigam'
avar_mle(distr)
# S4 method for class 'Multigam'
avar_me(distr)
# S4 method for class 'Multigam'
avar_same(distr)Arguments
- shape, scale
numeric. The non-negative distribution parameters.
- x
For the density function,
xis a numeric vector of quantiles. For the moments functions,xis an object of classMultigam. For the log-likelihood and the estimation functions,xis the sample of observations.- log
logical. Should the logarithm of the probability be returned?
- n
number of observations. If
length(n) > 1, the length is taken to be the number required.- distr
an object of class
Multigam.- type
character, case ignored. The estimator type (mle, me, or same).
- ...
extra arguments.
- par0, method, lower, upper
arguments passed to optim for the mle optimization. See Details.
- na.rm
logical. Should the
NAvalues be removed?
Value
Each type of function returns a different type of object:
Distribution Functions: When supplied with one argument (
distr), thed(),p(),q(),r(),ll()functions return the density, cumulative probability, quantile, random sample generator, and log-likelihood functions, respectively. When supplied with both arguments (distrandx), they evaluate the aforementioned functions directly.Moments: Returns a numeric, either vector or matrix depending on the moment and the distribution. The
moments()function returns a list with all the available methods.Estimation: Returns a list, the estimators of the unknown parameters. Note that in distribution families like the binomial, multinomial, and negative binomial, the size is not returned, since it is considered known.
Variance: Returns a named matrix. The asymptotic covariance matrix of the estimator.
Details
The probability density function (PDF) of the multivariate gamma distribution is given by: $$ f(x; \alpha, \beta) = \frac{\beta^{-\alpha_0}}{\prod_{i=1}^k\Gamma(\alpha_i)}, e^{-x_k/\beta} x_1^{\alpha_1-1}\prod_{i=1}^k (x_i - x_{i-1})^{(\alpha_i-1)} \quad x > 0. $$
The MLE of the multigamma distribution parameters is not available in closed
form and has to be approximated numerically. This is done with optim().
Specifically, instead of solving a \((k+1)\) optimization problem w.r.t
\(\alpha, \beta\), the optimization can be performed on the shape parameter
sum \(\alpha_0:=\sum_{i=1}^k\alpha \in(0,+\infty)^k\). The default method
used is the L-BFGS-B method with lower bound 1e-5 and upper bound Inf.
The par0 argument can either be a numeric (satisfying
lower <= par0 <= upper) or a character specifying the closed-form estimator
to be used as initialization for the algorithm ("me" or "same" - the
default value).
References
Mathal, A. M., & Moschopoulos, P. G. (1992). A form of multivariate gamma distribution. Annals of the Institute of Statistical Mathematics, 44, 97-106.
Oikonomidis, I. & Trevezas, S. (2025), Moment-Type Estimators for the Dirichlet and the Multivariate Gamma Distributions, arXiv, https://arxiv.org/abs/2311.15025
Examples
# -----------------------------------------------------
# Multivariate Gamma Distribution Example
# -----------------------------------------------------
# Create the distribution
a <- c(0.5, 3, 5) ; b <- 5
D <- Multigam(a, b)
# ------------------
# dpqr Functions
# ------------------
d(D, c(0.3, 2, 10)) # density function
#> [1] 0.0001605858
# alternative way to use the function
df <- d(D) ; df(c(0.3, 2, 10)) # df is a function itself
#> [1] 0.0001605858
x <- r(D, 100) # random generator function
# ------------------
# Moments
# ------------------
mean(D) # Expectation
#> [1] 2.5 17.5 42.5
var(D) # Variance
#> [1] 12.5 87.5 212.5
finf(D) # Fisher Information Matrix
#> shape1 shape2 shape3 scale
#> shape1 4.934802 0.0000000 0.000000 0.20
#> shape2 0.000000 0.3949341 0.000000 0.20
#> shape3 0.000000 0.0000000 0.221323 0.20
#> scale 0.200000 0.2000000 0.200000 0.34
# List of all available moments
mom <- moments(D)
mom$mean # expectation
#> [1] 2.5 17.5 42.5
# ------------------
# Point Estimation
# ------------------
ll(D, x)
#> [1] -358.8854
llmultigam(x, a, b)
#> [1] -358.8854
emultigam(x, type = "mle")
#> $shape
#> [1] 0.5167311 3.0827569 5.0220514
#>
#> $scale
#> [1] 5.056676
#>
emultigam(x, type = "me")
#> $shape
#> [1] 0.5933877 3.4603289 5.6212319
#>
#> $scale
#> [1] 4.506105
#>
emultigam(x, type = "same")
#> $shape
#> [1] 0.5316584 3.1003557 5.0364629
#>
#> $scale
#> [1] 5.029296
#>
mle(D, x)
#> $shape
#> [1] 0.5167311 3.0827569 5.0220514
#>
#> $scale
#> [1] 5.056676
#>
me(D, x)
#> $shape
#> [1] 0.5933877 3.4603289 5.6212319
#>
#> $scale
#> [1] 4.506105
#>
same(D, x)
#> $shape
#> [1] 0.5316584 3.1003557 5.0364629
#>
#> $scale
#> [1] 5.029296
#>
e(D, x, type = "mle")
#> $shape
#> [1] 0.5167311 3.0827569 5.0220514
#>
#> $scale
#> [1] 5.056676
#>
mle("multigam", x) # the distr argument can be a character
#> $shape
#> [1] 0.5167311 3.0827569 5.0220514
#>
#> $scale
#> [1] 5.056676
#>
# ------------------
# Estimator Variance
# ------------------
vmultigam(a, b, type = "mle")
#> shape1 shape2 shape3 scale
#> shape1 0.2355724 0.4114692 0.7342357 -0.8125161
#> shape2 0.4114692 7.6734815 9.1744629 -10.1525963
#> shape3 0.7342357 9.1744629 20.8894190 -18.1165398
#> scale -0.8125161 -10.1525963 -18.1165398 20.0480307
vmultigam(a, b, type = "me")
#> shape1 shape2 shape3 scale
#> shape1 0.5444444 1.1 1.944444 -2.111111
#> shape2 1.1000000 14.6 20.000000 -21.000000
#> shape3 1.9444444 20.0 39.444444 -36.111111
#> scale -2.1111111 -21.0 -36.111111 37.777778
vmultigam(a, b, type = "same")
#> shape1 shape2 shape3 scale
#> shape1 0.3813248 0.1212822 0.3132515 -0.4799147
#> shape2 0.1212822 8.7276929 10.2128426 -11.2128216
#> shape3 0.3132515 10.2128426 23.1325478 -19.7991821
#> scale -0.4799147 -11.2128216 -19.7991821 21.4658137
avar_mle(D)
#> shape1 shape2 shape3 scale
#> shape1 0.2355724 0.4114692 0.7342357 -0.8125161
#> shape2 0.4114692 7.6734815 9.1744629 -10.1525963
#> shape3 0.7342357 9.1744629 20.8894190 -18.1165398
#> scale -0.8125161 -10.1525963 -18.1165398 20.0480307
avar_me(D)
#> shape1 shape2 shape3 scale
#> shape1 0.5444444 1.1 1.944444 -2.111111
#> shape2 1.1000000 14.6 20.000000 -21.000000
#> shape3 1.9444444 20.0 39.444444 -36.111111
#> scale -2.1111111 -21.0 -36.111111 37.777778
avar_same(D)
#> shape1 shape2 shape3 scale
#> shape1 0.3813248 0.1212822 0.3132515 -0.4799147
#> shape2 0.1212822 8.7276929 10.2128426 -11.2128216
#> shape3 0.3132515 10.2128426 23.1325478 -19.7991821
#> scale -0.4799147 -11.2128216 -19.7991821 21.4658137
v(D, type = "mle")
#> shape1 shape2 shape3 scale
#> shape1 0.2355724 0.4114692 0.7342357 -0.8125161
#> shape2 0.4114692 7.6734815 9.1744629 -10.1525963
#> shape3 0.7342357 9.1744629 20.8894190 -18.1165398
#> scale -0.8125161 -10.1525963 -18.1165398 20.0480307