bivariate gamma

function Q = bivgmarnd(N, a1, b1, a2, b2, rho)

% X1 ~ Gamma(a1, b1),  X2 ~ Gamma(a2, b2),  X1 is independent of X2
% Let Y1 = X1 + p * X2, Y2 = p * X1 + X2
% E(Y1) = a1 * b1 + p * a2 * b2, E(Y2) = p*a1*b1 + a2*b2
% Var(Y1) = Var(X1) + p^2 * Var(X2) = a1*b1^2+a2*(b2*p)^2
% Var(Y2) = a1 * (p*b1)^2 + a2*b2^2
% Cov(Y1, Y2) = Cov(X1 + p*X2, p*X1+X2) = p*Var(X1) + p*Var(X2) + 2*p*Cov(X1,X2)
%   = p*a1*b1^2 + p*a2*b2^2  (Cov(X1,X2) = 0 because X1 is indep. of X2)
%
% examples:
% corr(bivgmarnd(1e6, 3, 2, 3, 1, 0.6))
% % ans =
% %     1.0000    0.6015
% %     0.6015    1.0000
% corr(bivgmarnd(1e6, 3, 2, 3, 1, 0.9))
% % ans =
% %     1.0000    0.9001
% %     0.9001    1.0000
% corr(bivgmarnd(1e6, 3, 2, 3, 1, -0.3))
% % ans =
% %     1.0000    -0.3016
% %    -0.3016     1.0000
syms p;
p_ = double(solve((p*a1*b1^2 + p*a2*b2^2) / (sqrt(a1*b1^2+a2*(b2*p)^2) * sqrt(a1 * (p*b1)^2 + a2*b2^2)) == rho, p));
rho_ = (p_*a1*b1^2 + p_*a2*b2^2) ./ (sqrt(a1*b1^2+a2*(b2*p_).^2) .* sqrt(a1 * (p_*b1).^2 + a2*b2^2));
validate = find(abs(rho_ - rho) < 1e-3);
if isempty(validate)
  error('cannot find the real solution for (p1, p2)')
else
  p_ = p_(validate(1));
end

X1 = gamrnd(a1, b1, N, 1);
X2 = gamrnd(a2, b2, N, 1);

Q = [X1, X2] * [1, p_; p_, 1];
end