English: A reconstruction of Figure 1 of Moments of Wishart-Laguerre and Jacobi ensembles of random matrices (Livan 2011)

https://arxiv.org/pdf/1103.2638.pdf

```python import numpy as np import matplotlib.pyplot as plt

1. 1 for LOE, 2 for LUE, 4 for LSE

betas = 1, 2, 4 NMs = [(8, 15)]

1. Choose number of samples

Nmatr = 100000 Es = {} for n, m in NMs:

   for beta in betas:
       if beta == 1:  # Wishart Orthogonal Ensemble
           X = np.random.randn(Nmatr, n, m)
           M = np.einsum('ijk,ilk->ijl', X, X)
           E = np.linalg.eigvals(M.reshape(Nmatr, n, n)).flatten()
       elif beta == 2:  # Wishart Unitary Ensemble
           X_real = np.random.randn(Nmatr, n, m)
           X_imag = np.random.randn(Nmatr, n, m)
           X = X_real + 1j * X_imag
           M = np.einsum('ijk,ilk->ijl', X, X.conjugate())
           E = np.linalg.eigvals(M.reshape(Nmatr, n, n)).flatten()
       elif beta == 4:  # Wishart Symplectic Ensemble
           A = np.random.randn(Nmatr, n,m) + 1j * np.random.randn(Nmatr, n,m)
           B = np.random.randn(Nmatr, n,m) + 1j * np.random.randn(Nmatr, n,m)
           X = np.block([[A, B],[-np.conj(B), np.conj(A)]])
           M = np.einsum('ijk,ilk->ijl', X, X.conjugate())
           E = np.linalg.eigvals(M.reshape(Nmatr, 2 * n, 2 * n)).flatten()
       Es[(n, m, beta)] = E

for n, m in NMs:

   plt.figure(figsize=(16, 8))

   legends = {1: "LOE", 2:"LUE", 4:"LSE"}
   colors={1:"blue", 2:"red", 4:"green"}
   for beta in betas:
       color=colors[beta]
       E = Es[(n, m, beta)]
       xs = np.real(E) / np.sqrt(beta)

       bin_heights, bin_borders, _ = plt.hist(xs, bins=500, density=True, color=color, alpha=0.1)
       bin_centers = bin_borders[:-1] + np.diff(bin_borders) / 2

       # Compute sliding window average
       window_size = 5
       window = np.ones(window_size) / window_size
       smoothed_heights = np.convolve(bin_heights, window, mode='same')

       # Plot sliding window average
       plt.plot(bin_centers, smoothed_heights, label=legends[beta], color=color)

   # Add plot labels and title
   plt.xlabel('x', fontsize=14)
   plt.ylabel('ρ(x)', fontsize=14)
   plt.title(r'Eigenvalues $/\sqrtTemplate:\beta$, with (N, M) = {}'.format((n, m)), fontsize=18)
   plt.grid(True)
   plt.legend()
   plt.show()

```