The Euler-Mascheroni constant $\gamma$ is like many other constants ($e$, $\pi$, and , $\varphi$) used in a simple pattern. That is why it shows up in many places and also in different fields that do not seem to be connected at first glance.
Below we give two 'reasons' behind the occurrence of $\gamma$ in the expression for the mean of the Gumbel distribution. (sidenote: The word 'reason' might be a bit strong. What we do is dive a bit deeper into the Gumbel distribution and see how it relates to other constructions that relate to the constant.)
We look at two expressions of the Euler-Mascheroni constant
The derivative of the gamma function (the moment generating function of the Gumbel distribution), incorporates a $\gamma$ into it. In the point $t=0$ we have $$\Gamma^{\prime}(1-\beta t) \underset{t=0}{=} -\gamma \beta$$
The difference between the n-th harmonic number, $H_n = \sum_{k=1}^n \frac{1}{k}$, and the natural logarithm of $n$ approaches the Euler-Mascheroni constant. $$\lim_{n \to \infty} H_n - \log(n) = \gamma$$
1 Derivative of the gamma function
In this section we are making a connection by computing the mean of the Gumbel distribution based on the moment generating function (MGF).
Background: the derivative of the moment generating function and it's relationship to the mean
If you know the MGF, the integral below (which is also considered as the bilateral Laplace transform with a negative argument, $\mathcal{B}(-t)$)
$$M(t) = \mathcal{B}(-t) = \int_{-\infty}^{\infty} e^{tx} f(x) \,dx$$
then you can compute
$$\mu = \int_{-\infty}^{\infty} x f(x) \,dx$$
by taking the derivative of the $M(t)$. You can see the equivalence by differentiation under the integral sign, and for $t=0$ this equals the integral expression for the mean
$$\begin{array}{rcl} M^\prime(t) &=& \int_{-\infty}^{\infty} x \underbrace{e^{tx}}_{\text{$e^{tx}= 1$ if $t=0$}} f(x) \,dx \\ \\ M^\prime(0) &=& \int_{-\infty}^{\infty} x f(x) \,dx = \mu \end{array}$$
The moment generating function of the Gumbel distribution
The moment generating function of the Gumbel distribution is
$$M(t) = \Gamma(1-\beta t) e^{\mu t}$$
so this $\Gamma(1-\beta t)$ whose derivative when $t=0$ is equal to $\gamma \beta$ is 'giving' you this constant $\gamma$ in the expression for the mean of a Gumbel distributed variable.
The connection
The connection made in this section might be a bit weak and considered as cyclical reasoning. This is effectively just like the pure mechanistic integration that you mentioned but done indirectly by converting to the Laplace domain (the MGF) and performing the integration there. But why is the MGF a gamma function?
We wrote in the first sentence of this post that there are simple patterns behind the well known constants. Behind the Gumbel distribution, and the extreme value distribution in general, there is also a simple pattern. It is described by Fisher and Tippet as following
Since the extreme member of a sample of $mn$ may be regarded
as the extreme member of a sample of $n$ of the extreme members
of samples of $m$, and since, if a limiting form exist, both of these distributions will tend to the limiting form as $m$ is increased indefinitely, it follows that the limiting distribution must be such that the extreme member of a sample of $n$ from such a distribution has itself a similar distribution.
So the extreme value distribution under taking the maximum and translating and scaling, is again an extreme value distribution of the same distribution family.
The Gumbel distribution is the special case when there is no scaling. If you take the maximum of a sample of Gumbel distributed variables, then you have again a Gumbel distribution that is translated by some value $b_n$.
(Compare this with other types of operations, like scale location family, stable distributions, or inversion. Distributions from a family that are, after some operation, still a member of the same family.)
The Gumbel satisfies the following equation:
$$ F(x)^n = F(x+b_n)$$
Let's check this by filling in the function for the Gumbel distribution
$$\begin{array}{rcl} F(x)^n &=& \left(e^{e^{-(x-\mu)/\beta}} \right)^n \\
&=& e^{ne^{-(x-\mu)/\beta}}\\
& =& e^{e^{\log(n)-(x-\mu)/\beta+}}\\
& =& e^{e^{-(x-\mu-\beta \log(n))/\beta}}\\
& =& e^{e^{-(x-\mu^\prime)/\beta}}
\end{array} $$
Thus, the operation of taking the maximum from a sample of size $n$ is effectively like shifting the distribution mean by a factor $\mu^\prime = \mu + \log(n)\beta$.
For the gamma distribution, we also have a simple pattern. It has the property
$$\Gamma(z)z = \Gamma(z+1)$$
If you multiple the function with its argument then this is like shifting the argument by 1 step.
It is in these two simple patterns that we can see the connection why the MGF of the Gumbel distribution is a gamma function. Let's repeat the two simple patterns
- The Gumbel distribution: The cumulative distribution function (CDF) follows the relation $$F(x)^n = F(x+b_n)$$
- The gamma function: The function follows the relation $$\Gamma(z)z = \Gamma(z+1)$$
My intuition tells me that we can connect these two simple patterns by means of a relationship for the moment generating function. We might do this in such a way that we do not need to explicitly compute the integral, and we could reason that the moment generating function must have the properties of the gamma function without evaluating the integral, but just by rewriting the integral. However, I do not yet get far in showing that there must be such a connection.
Moment generating function of the maximum of a sample
Work in progress
2 Difference between the n-th harmonic number and the natural logarithm of $n$
One direct connection between the Euler-Mascheroni constant and the mean of the Gumbel distribution is the following:
Let $X_k \sim Exp(1/k)$ then the limiting distribution for the sum is for $n \to \infty$
$$-\log(n) + \sum_{k=1}^n X_k \to \text{Gumbel}(\mu =0,\beta = 1)$$
For an explanation for this connection as a sum see on math.stackexchange the question "Proof that the limit of a series of exponential distributed random variables follows a Gumbel distribution"
From this point of view, as the limit for a sum of exponential distributions, the mean of the Gumbel distribution can be expressed as the limit of the sum of the means of these exponential distributions, which is the limit of $H_n - \log(n) = \gamma$
$$\lim_{n \to \infty} E\left[ -\log(n) + \sum_{k=1}^n X_k \right] = \lim_{n \to \infty} -\log(n) + \sum_{k=1}^n \frac{1}{k} = \gamma$$
The coupon collector's problem illustrates intuitively this connection between the Gumbel distribution and a sum of variables. (Intuition about the coupon collector problem approaching a Gumbel distribution)
One way to view the extreme events as a sum could be by seeing the distribution of the extreme of a period of $k$ days and the extreme of a period of $k-1$ days as the latter being the former plus an extra little bit. In the case of exponential distributions for the daily values you get that the extra little bit is also exponential distributed and has a mean that scales as $1/k$ (the extra little bits decrease in size).
