This distribution has the
bivariate gamma distribution
bigamma.mckay
as a special case.
Let \(Q > 1\) be the number of columns of the
response matrix y.
Then the
joint probability density function is given by
$$f(y_1,\ldots,y_Q; b, s_1, \ldots, s_Q) =
y_1^{s_1} (y_2 - y_1)^{s_2}
\cdots (y_Q - y_{Q-1})^{s_Q}
\exp(-y_Q / b) / [b^{s_Q^*}
\Gamma(s_1) \cdots \Gamma(s_Q)]$$
for \(b > 0\),
\(s_1 > 0\), ...,
\(s_Q > 0\) and
\(0<y_1< y_2<\cdots<y_Q<\infty\).
Also,
\(s_Q^* = s_1+\cdots+s_Q\).
Here, \(\Gamma\) is
the gamma function,
By default, the linear/additive predictors are
\(\eta_1=\log(b)\),
\(\eta_2=\log(s_1)\),
...,
\(\eta_M=\log(s_Q)\).
Hence \(Q = M - 1\).
The marginal distributions are gamma,
with shape parameters
\(s_1\) up to \(s_Q\), but they have a
common scale parameter \(b\).
The fitted value returned
is a matrix with columns equalling
their respective means;
for column \(j\) it is
sum(shape[1:j]) * scale.
The correlations are always positive;
for columns \(j\) and \(k\)
with \(j < k\),
the correlation is
sqrt(sum(shape[1:j]) /sum(shape[1:k])).
Hence the variance of column \(j\)
is sum(shape[1:j]) * scale^2.