I know I am late to the party here, but a couple of people mentioned different Basel Problem solutions by Pace and Calabi These solutions, along with another one by Zagier and Kontsevich all have general versions that induce different probabilistic viewpoints of $\zeta(2k).$ I reference my publication for details: https://www.ams.org/journals/qam/2018-76-03/S0033-569X-2018-01499-3/ or https://arxiv.org/abs/1710.03637.
First, Pace considers i.i.d. nonnegative Cauchy Random Variables $X,Y$ and their quotient $Z=X/Y.$ It turns out $Z$ has density function:
$$f_Z(z)= \frac{4}{\pi^2} \frac{\log(z)}{z^2-1}, \quad z>0.$$ Thus,
$$\frac{\pi^2}{4} = \int_{0}^{\infty}\frac{\log(z)}{z^2-1} \ dz=\int_{0}^{1}\frac{\log(z)}{z^2-1} \ dz+\int_{1}^{\infty}\frac{\log(z)}{z^2-1} \ dz.$$ Making the substitution $t=1/z$ allows us to deduce
$$\frac{\pi^2}{8} = \int_{0}^{1}\frac{\log(z)}{z^2-1} \ dz.$$ After converting the integrand on the right hand side into a geometric series and using monotone convergence theorem to interchange sum and integral, we obtain
$$\frac{\pi^2}{8} = \sum_{n=0}^{\infty} \frac{1}{(2n+1)^2} = \frac{3}{4} \zeta(2).$$
Pace's solution generalizes. It leads to
Interpretation 1:
The sums $$S(k)=\sum_{n \geq 0} \frac{(-1)^{nk}}{(2n+1)^k}, \quad k \in \mathbb{N},$$ are obtained by considering $k$ i.i.d. nonnegative Cauchy Random Variables $X_1, \dots, X_k$ and integrating the density function of the quotient $Z=X_1/ \dots / X_k.$ Then, upon splitting the sum $\zeta(2k)$ and rearranging, we obtain $$S(2k)= \frac{2^{2k}-1} {2^{2k}}\zeta(2k).$$
With repeated use of partial fractions and Fubini's Theorem, we can obtain a density function for $f_{Z}(z).$
Interpretation 1 will yield $$f_{Z}(z)= \begin{cases} \left(\frac{2}{\pi} \right)^k \frac{k |B_k|}{4^k-2^k}\frac{\ln^{k-1}(z)}{z^2-(-1)^k} & {k \text{ even}} \\ \left(\frac{2}{\pi} \right)^k |E_{k-1}|\frac{\ln^{k-1}(z)}{z^2-(-1)^k} & {k \text{ odd}} \end{cases}, \quad z>0 $$ where $B_k$ and $E_k$ are the $k$th Bernoulli and Eulerian numbers, respectively.
Replicating Pace's steps on $\int_{0}^{\infty} f(z) \ dz$ i.e. splitting of the region of integration and using the geometric series argument, we can recover $S(k).$
Interpretation 1 also allows us to find the sums $$S(k,a)=\sum_{n \in \mathbb{Z}} \frac{(-1)^{nk}}{(an+1)^k}, \quad a>1, \quad \frac{1}{a} \notin \mathbb{N}.$$ In particular, one considers the independent random variables
$X_1,$ having density function $$f_{X_1}(x_1)= \frac{a}{\pi} \sin \left(\frac{\pi}{a} \right) \frac{1}{x_1^a + 1}, \quad x_1 \geq 0,$$ and $X_2, \dots, X_k$ all having density function $$f_{X_i}(x_i) = \frac{a}{\pi} \sin \left(\frac{\pi}{a} \right) \frac{x_i^{1-2/a}}{x_i^2 + 1}, \quad x_i \geq 0.$$
We can obtain $S(k,a)$ by integrating the density function for $Z=X_1/(X_2 \dots X_k)^{2/a}.$ Replicating Pace's arguments for $\int_{0}^{\infty} f(z) \ dz$ allows us to recover $S(k,a).$
In previous literature, Bourgade, Fujita, and Yor evaluate $S(k,a)$ in a slightly different way from here (see https://projecteuclid.org/download/pdf_1/euclid.ecp/1465224952), with the crux being differentiating under the integral sign.
Second, Calabi evaluates the double integral $$\int_{0}^{1} \int_{0}^{1} \frac{1}{1-x^2y^2} \ dx \ dy$$ in two ways. He obtains $S(2)$ with a geometric series argument and then uses the trigonometric change of variables $$x=\frac{\sin(u)}{\cos(v)}, \quad y=\frac{\sin(v)}{\cos(u)}$$ whose Jacobian determinant is implicitly the denominator of the integrand. The transformation diffeomorphically maps $(0,1)^2$ to the open isosceles triangle with base and height both $\pi/2.$ Thus $S(2)$ is geometrically the area of this triangle $\pi^2/8,$ hence solving the Basel Problem. The generalization of this is finding the integral $$\int_{(0,1)^k} \frac{1}{1-(-1)^k x_1^2 \dots x_{k}^2} \ dx_1 \dots \ dx_k, $$ which on one hand happens to be equal to $S(k)$ by a geometric series argument. On the other hand, Calabi's general change of variables $$x_i=\frac{\sin(u_i)}{\cos(u_{i+1})}, \quad \dots \quad , x_k=\frac{\sin(u_k)}{\cos(u_1)}$$ where $1 \leq i \leq k-1,$ has a Jacobian determinant which cancels with the integrand's denominator and diffeomorphically maps $(0,1)^k$ to the open convex polytope
$$\Delta^{k}=\lbrace (u_1, \dots, u_k) \in \mathbb{R}^k: u_i +u_{i+1}, u_k +u_1<\pi/2, \quad u_1, \dots, u_k >0, \quad 1 \leq i \leq k-1 \rbrace.$$ Thus $S(k)=\text{Volume}(\Delta^{k}).$ Rescaling $u_i=\frac{\pi}{2} v_i,$ and considering $k$ i.i.d. Uniform Random Variables $V_1, \dots, V_k \in (0,1),$ we get
Interpretation 2:
$$S(k)= \left(\frac{\pi}{2} \right)^k \text{Pr} \left(V_1+V_2<1, \dots V_{k-1}+V_k<1, V_k +V_1 <1 \right).$$ That is explicitly $S(k)$ is proportional to the probability that $k$ i.i.d. uniform random variables $V_1, \dots, V_k$ in $(0,1)$ have cyclically pairwise consecutive sums less than $1.$
In previous literature, Beukers, Calabi, and Kolk dissect $\Delta^k$ into congruent pyramids (see https://pdfs.semanticscholar.org/35be/01e63c0bfd32b82c97d58ccc9c35471c3617.pdf), and some users have already mentioned Elkies and Silagadze perform spectral analysis on the characteristic function of $\Delta^k$ (see (https://www.maa.org/sites/default/files/pdf/news_old/Elkies.pdf and http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.749.9234&rep=rep1&type=pdf). $\Delta^k$ turns out to be a chain polytope, which Stanley independently interprets in general (see http://dedekind.mit.edu/~rstan/pubs/pubfiles/66.pdf).
In my paper, the aforementioned probability with the uniform random variables leads to a gargantuan formula that does not involve $B_k$ or $E_k$ in any way. Rather, it is the result of imposing different conditions on $V_1, \dots, V_k.$ See my paper for details.
The result of Interpretation 2 is that
$$S(k)= \left(\frac{\pi}{4} \right)^k +\left(\frac{\pi}{4} \right)^k \sum_{n=1}^{ \left \lfloor \frac{k}{2} \right \rfloor} \sum_{\substack{(r_1, \dots, r_n) \in [k]^n: \\ |r_p-r_q| \notin \lbrace 0,1,k-1 \rbrace, \\ p,q \in [n]} } \prod_{i=1}^{n} \frac{1}{i+\sum_{j=1}^{i} \alpha_j},$$
where $[m]:= \lbrace 1, \dots, m \rbrace$ and
$$\alpha_j=2- \delta(k,2) - \sum_{m=1}^{j-1} \delta(|r_m-r_j|,2)+\delta(|r_m-r_j|,k-2)$$ and $\delta(a,b)$ is the Kronecker Delta Function. In particular, the inner sum is taken over all tuples $(r_1, \dots, r_n) \in [k]^
n$ having cyclically pairwise nonconsecutive entries.
Lastly, Zagier and Kontsevich reproduce Calabi's $\zeta(2)$ proof by evaluating $$\int_{0}^{1}\int_{0}^{1} \frac{1}{\sqrt{xy}(1-xy)}.$$ On one hand, this integral is $4S(2).$ But on the other hand, the change of variables $$x=\frac{\xi^2(\eta^2+1)}{\xi^2+1}, \quad y=\frac{\eta^2(\xi^2+1)}{\eta^2+1}$$ transforms it into
$$\int_{0}^{\infty} \int_{0}^{1/\xi} \frac{4}{(\xi^2+1)(\eta^2+1)} d \eta \ d \xi.$$
The general version of this involves evaluating
$$I=\int_{(0,1)^k} \frac{1}{\sqrt{x_1 \dots x_k}(1-x_1 \dots x_k)} \ dx_1 \dots \ dx_k,$$ which is equal to $2^k S(k).$ But the generalized transformation
$$x_i= \frac{\xi_i^2 (\xi_{i+1}^2+1)}{(\xi_{i}^2+1)}, \quad \dots \quad,x_k= \frac{\xi_k^2 (\xi_{1}^2+1)}{(\xi_{k}^2+1)},$$ and defining $\Xi_1, \dots, \Xi_k$ to be $k$ i.i.d. nonnegative Cauchy random variables, it turns out
Intepretation 3
$$S(k)= \left(\frac{\pi}{2} \right)^k \text{Pr} \left(\Xi_1\Xi_2 <1, \dots \Xi_{k-1} \Xi_k <1, \Xi_k \Xi_1 <1 \right).$$ Explicitly, $S(k)$ is proportional to the probability that $k$ i.i.d. nonnegative Cauchy random variables $\Xi_1, \dots, \Xi_k$ have cyclically pairwise consecutive products less than $1.$
The result of Interpretation 3 is the same as the result of Interpretation 2.