General beta simulation of Circular ensembles

[adrienhardy]

One alternative way to sample from the C\beta E for general $beta\geq 0$ goes as follows: Sample $y_1,\ldots,yn$ from a Jacobi ensemble on $[0,1]$ with parameters $a=b=1$, thus from the density $$ \prod{i<j}|y_i-yj|^\beta\prod{i=1}^n 1_0,1, $$ and then return $x_j:=e^{2i\pi x_j}$. This has the same laws as the C\beta E for arbitrary $\beta\geq 0$.

[guilgautier]

Thanks @adrienhardy, that's a good point! It is a just thought or do you have reference for it ?

Not only that would mean there's a quindiagonal model for circular ensemble (see KiNe04 and/or the docs) which is easy to sample only for $\beta\in \mathbb{N}^{*}$, but also a tridiagonal model through the Jacobi ensemble with reference measure $Beta(1,1)=Unif_[0,1]$ which is available for all $\beta\geq 0$ :)

[guilgautier]

After a quick implementation, one gets an accumulation of points around (1,0).

At first glance, the mapping Jacobi with Beta(1,1)=Unif[0,1] to the unit circle using x |->exp(2 i pi x) doesn't seem to the right thing to do when expecting expecting a Circular ensemble. With Jacobi the points are confined in [0,1] but the repulsive interaction send the could of points close the boundary. In the limit N->oo Jacobi "behaves like arcsine law" and the above mapping would send may points around (0,1).

Here, we wish to send points onto the circle from [0,1], whereas (if I remember well) KiNe04 "do" the opposite: they project points in interaction on the unit circle onto the real line.

Here is the code if you wish to reproduce the figure above

from dppy.beta_ensembles import *

N, beta = 30, 2

# Circular
circular = BetaEnsemble('circular', beta=beta)
circular.sample(**{'size':40})
circular.plot()

# From Jacobi to Circular ?
# Map [0,1] onto the unit circle
# x |-> exp(2 i pi x)
jacobi = BetaEnsemble('jacobi', beta=beta)

## Sample from Jacobi ensemble with mu=Beta(1,1)=Unif_[0,1]
for sampl_mode in ('banded', 'full'):

    if sampl_mode == 'banded': # With tridiagonal matrix
        sampl_params = {'a':1, 'b':1, 'size':N}

    elif sampl_mode == 'full': # With full matrix (draw with Haar measure on U_N)
        M_1 = int(N-1+2/beta)
        M_2 = M_1
        sampl_params = {'M_1':M_1, 'M_2':M_2, 'N':N}

    jacobi.sample(sampl_mode, **sampl_params)
    jacobi_pts = jacobi.list_of_samples[-1]

    # Mapping onto the circle
    circle_pts = np.exp(2*1j*np.pi*jacobi_pts)

    # Display the points
    plt.scatter(circle_pts.real, circle_pts.imag)
    plt.title(r'Circular from Jacobi with $\beta$={}, Jacobi drawn using {} matrix model'.format(beta, sampl_mode))
    plt.axis('equal')
    plt.show()

[adrienhardy]

You're right, I just did my computations again and I realize that I got everything backward... Please don't tell too many people that I'm really bad at making changes of variables ;-)

[guilgautier]

All right, I'm closing this issue. Feel free to reopen it if you find a way to go from Jacobi to Circular!

PS: I won't tell anyone :no_mouth:, however this thread is public ...