mp_eigval_plot.py

def proba_mp(alpha, x, lambda_m, lambda_p):
    if x == 0:
        return 1-alpha if alpha<1 else 0

    if x<lambda_m or lambda_p<x:
        return 0

    return math.sqrt((lambda_p-x) * (x-lambda_m)) / (2*x*math.pi*alpha) 
    
xmax = 5
d = 1000
X = [i*xmax/d for i in range(d)]

for alpha in [0.1, 0.2, 0.3, 0.4, 0.5, 1.0, 1.5, 2.0]:
    lambda_m = (1 - math.sqrt(alpha)) **2
    lambda_p = (1 + math.sqrt(alpha)) **2

    Y = [proba_mp(alpha, x, lambda_m, lambda_p) for x in X]
    plt.plot(X, Y, label=r'$\alpha$ = %.2f' % alpha)

plt.ylim((0,2))
plt.xlim((0,xmax))
plt.legend(loc="upper right")
plt.show()