Gianluca's suggestion is many times what you want. In this case a little work finds a general formula for the first part of the paper.
The paper deals with asymptotic properties of some statistics but for Mathematica easily finds the exact distribution, mean, and variance. (Approximately) using the notation in the paper two random variables N and M are selected from a discrete uniform distribution of the integers 1 through x. The probability mass function (pmf) for x=6 is given by the following:
x = 6;
pmfCounts = Tally[Flatten[Table[GCD[m, n], {n, 1, x}, {m, 1, x}], 1]]
(* {{1, 23}, {2, 7}, {3, 3}, {4, 1}, {5, 1}, {6, 1}} *)
pmf = pmfCounts;
pmf[[All, 2]] = pmf[[All, 2]]/x^2;
pmf
(* {{1, 23/36}, {2, 7/36}, {3, 1/12}, {4, 1/36}, {5, 1/36}, {6, 1/36}} *)
Then we check the frequency counts at oeis.org and find sequence A242114. The general Mathematica formula for the pmf is
gcdpmf[n_, k_] := (2 Total[EulerPhi[Range[Quotient[n, k]]]] - 1)/n^2
gcdpmf[213,47]
(* 11/45369 *)
From there you can calculate means, variances, etc. If you are more interested in asymptotic properties of functions of random variables, then specific examples are helpful.
