In [2]:
# A first Course in Random Matrix Theory
# by Marc Potters and Jean-Philippe Bouchaud
# Cambridge University Press 2020
# Python code used to generate graphs in Part I
In [3]:
%pylab inline
import scipy.linalg as npl
from cycler import cycler
from sklearn.isotonic import IsotonicRegression

Populating the interactive namespace from numpy and matplotlib
In [4]:
ColorGraphs=True
if(ColorGraphs):
    HistoColor='powderblue'
else:
    plt.rcParams['axes.prop_cycle'] = cycler(color=['#000000', '#777777','#bbbbbb'])
    HistoColor='#777777'
plt.rc('legend', fontsize=8)
plt.rc('font', family='serif', size=10)
plt.rc('lines', linewidth=1, color='black')


#uncomment if latex is installed
plt.rc('text', usetex=True)
plt.rc('font', family='serif', serif='New Times Roman', size=10)
plt.rc('text.latex', preamble=r'\usepackage{mathptmx}\usepackage{amsfonts}')

#figure sizes in inches
figwidth=9./2.54
figwidth2=6/2.54
figheight=6/2.54
figdpi=200
In [5]:
#helper functions
myseed=180569
def wigner(myN,sigma2=1.):
    mat=np.random.randn(myN, myN) #non-symmetric Gaussian matrix
    return sqrt(sigma2/(2.*myN))*(mat+mat.transpose()) #symmetrize and normalize
def wishart(myN,q=0.5):
    myT=int(myN/q)
    mat=np.random.randn(myN, myT) #non-symmetric Gaussian matrix
    return 1./myT*mat@mat.transpose() #symmetrize and normalize
def myg(values,z):#returns the sample Stieltjes transform of data values evaluated at points z (must be 1d array)
    return np.average(1./(z[:,newaxis]-values[newaxis,:]),axis=1)
def ginibre(myN,sigma2=1.):
    return sqrt(sigma2/myN)*np.random.randn(myN, myN) #non-symmetric Gaussian matrix
def matsqrt(mat):
    wval, wvect = npl.eigh(mat) #diagonalize mat
    return wvect@np.diag(sqrt(wval))@wvect.transpose() #reconstruct with sqrt e.vals.
In [6]:
plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
mat=np.array([[1,-0.2,0.2],[-.3,2,-.2],[0,1.1,3]])
eig=npl.eigvals(mat)
x=eig.real
y=eig.imag
plot(x,y,marker='x',linestyle='none',color='C1')
x=np.diag(mat)
y=0*x
plot(x,y,marker='.',linestyle='none',color='C1')
theta=arange(0,np.pi*2.02,np.pi*.02)
xc=1
r=.4
x=xc+r*cos(theta)
y=r*sin(theta)
plot(x,y,color='C0')
xc=2
r=.5
x=xc+r*cos(theta)
y=r*sin(theta)
plot(x,y,color='C0')
xc=3
r=1.1
x=xc+r*cos(theta)
y=r*sin(theta)
plot(x,y,color='C0')
plt.xlim(0,4.2)
plt.ylim(-1.4,1.4)
plt.xlabel('Re\,($\lambda$)')
plt.ylabel('Im\,($\lambda$)')
plt.tight_layout(pad=0.1)
#plt.savefig('fig1_1.pdf',format='pdf')
plt.show()
In [7]:
plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
x=np.arange(-5,5,.01)
y=.5/(np.pi*(x**2+.25))
plt.plot(x,y)
plt.annotate(text='', xy=(.5,.30), xytext=(-.5,.30), arrowprops=dict(arrowstyle='<->'))
plt.text(0,0.33,'$2\eta$',horizontalalignment='center')
plt.xlabel('$x$')
plt.ylabel('$K_\eta(x)$')
plt.xlim(-5,5)
plt.ylim(0,0.68)
plt.tight_layout(pad=0.1)
#plt.savefig('fig2_2.pdf',format='pdf')
plt.show()
In [8]:
plt.figure(figsize=(figwidth2,figheight),dpi=figdpi)
wval= npl.eigvalsh(wigner(400)) #diagonalize a Wigner
myx=np.arange(-3,3,.002)
my100=myg(wval,myx-1j)
my10=myg(wval,myx-.05j)
my1=myg(wval,myx-.0025j)
x=np.arange(-2,2.01,.01)
y=sqrt(abs(4-x*x))/2
plt.plot(myx,my10.imag,label='$\eta=1/\sqrt{N}$')
plt.plot(myx,my100.imag,label='$\eta=1$')
plt.plot(x,y,label='analytic',color='C0',linestyle="--")
plt.xlabel('$x$')
plt.ylabel('Im\,$g(x-{\\rm i}\eta)$')
plt.xlim(-3,3)
plt.ylim(0,1.15)
plt.legend()
plt.tight_layout(pad=0.1)
#plt.savefig('fig2_3a.pdf',format='pdf')
plt.show()
In [9]:
plt.figure(figsize=(figwidth2,figheight),dpi=figdpi)
wval= npl.eigvalsh(wigner(400)) #diagonalize a Wigner
myx=np.arange(-.12,.12,.0001)
my100=myg(wval,myx-1j)
my10=myg(wval,myx-.05j)
my1=myg(wval,myx-.0025j)
x=myx
y=sqrt(abs(4-x*x))/2
plt.plot(myx,my1.imag,label='$\eta=1/N$')
plt.plot(myx,my10.imag,label='$\eta=1/\sqrt{N}$')
plt.plot(x,y,label='analytic',color='C0',linestyle="--")
plt.xlabel('$x$')
plt.ylabel('${\\rm Im}\, g(x-{\\rm i}\eta)$')
plt.xlim(-.12,.12)
plt.ylim(0,2.80)
plt.legend()
plt.tight_layout(pad=0.1)
#plt.savefig('fig2_3b.pdf',format='pdf')
plt.show()
In [10]:
plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
x=np.arange(-2,2.01,.01)
y=sqrt(abs(4-x*x))/(2*np.pi)
plt.plot(x,y)
plt.xlabel('$\lambda$')
plt.ylabel('$\\rho(\lambda)$')
plt.xlim(-3,3)
plt.ylim(0,0.345)
plt.tight_layout(pad=0.1)
#plt.savefig('fig2_4.pdf',format='pdf')
plt.show()
In [11]:
plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
random.seed(myseed)
myN=1000
gi=ginibre(myN)
ei=npl.eigvals(gi)
plt.scatter(ei.real,ei.imag,s=1)
theta=arange(0,2*np.pi+.01,.01)
x=cos(theta)
y=sin(theta)
plt.plot(x,y,linestyle='--',color='C1')
plt.xlim(-1.55,1.55)
plt.ylim(-1.1,1.1)
plt.xlabel('Re$\lambda$')
plt.ylabel('Im$\lambda$')
plt.tight_layout(pad=0.1)
#plt.savefig('fig3_1.pdf',format='pdf')
plt.show()
In [12]:
plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
q=.5
lm=(1-sqrt(q))**2
lp=(1+sqrt(q))**2
dd=(lp-lm)/512
x=np.arange(lm,lp+dd,dd)
y=sqrt(abs((lp-x)*(x-lm)))/(np.pi*x*q*2)
plt.plot(x,y,label='$q=1/2$')
q=2
lm=(1-sqrt(q))**2
lp=(1+sqrt(q))**2
dd=(lp-lm)/512
x=np.arange(lm,lp+dd,dd)
y=sqrt(abs((lp-x)*(x-lm)))/(np.pi*x*q*2)
plt.plot(x,y,'C1',label='$q=2$')
x=(0,0)
y=(0,2)
plt.plot(x,y,'C1')
plt.xlabel('$\lambda$')
plt.ylabel('$\\rho(\lambda)$')
plt.ylim(0,1)
plt.legend()
plt.tight_layout(pad=0.1)
#plt.savefig('fig4_2.pdf',format='pdf')
plt.show()
In [13]:
plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
x=np.arange(0.001,4.01,.01)
q=0.5
y=(x+(q-1)*log(x))/q
plt.plot(x,y,label="$q=1/2$")
q=2.
y=(x+(q-1)*log(x))/q
plt.plot(x,y,label="$q=2$")
plt.xlabel('$x$')
plt.ylabel('$V(x)$')
plt.xlim(-0.1,4.)
plt.ylim(-2,6.0)
plt.legend()
plt.tight_layout(pad=0.1)
#plt.savefig('fig5_1.pdf',format='pdf')
plt.show()
In [14]:
plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
x=arange(-2.2,2.4,.1)
y=x*x/2.
plt.plot(x,y,label='$V(x)$')
myN=20
wval= npl.eigvalsh(wigner(myN))
yy=wval*wval/2.
plt.plot(wval,yy,marker='.',linestyle='none',color='C1',label='$N=20$')
plt.xlim(-2.1,2.1)
plt.ylim(-0.1,2.2)
plt.xlabel('$x$')
plt.ylabel('$V(x)$')
plt.legend()
plt.tight_layout(pad=0.1)
#plt.savefig('fig5_2.pdf',format='pdf')
plt.show()
In [15]:
plt.figure(figsize=(figwidth2,figheight),dpi=figdpi)
def rho(x,gamma=0):
    a2=(sqrt(1+12*gamma)-1)/(6*gamma)
    q=gamma*x**2+1+2*gamma*a2
    return q*sqrt(np.maximum(4*a2-x**2,0))/(2*np.pi)
x=arange(-sqrt(8),sqrt(8),.01)
y=rho(x,1.)
plt.plot(x,y,label='$\gamma=1$')
y=rho(x,0.001)
plt.plot(x,y,label='$\gamma=0$')
y=rho(x,-1./12)
plt.plot(x,y,label='$\gamma=-1/12$')

plt.legend()
plt.xlabel('$\lambda$')
plt.ylabel('$\\rho(\lambda)$')
plt.tight_layout(pad=0.1)
#plt.savefig('fig5_3a.pdf',format='pdf')
plt.show()
In [16]:
plt.figure(figsize=(figwidth2,figheight),dpi=figdpi)
x=arange(-4.5,4.5,.01)
y1=x**2/2.+x**4/4.
y2=x**2/2.
y3=x**2/2.-(1./12.)*x**4/4.
plt.plot(x,y1,label='$\gamma=1$')
plt.plot(x,y2,label='$\gamma=0$')
plt.plot(x,y3,label='$\gamma=-1/12$')
xx=np.array([-sqrt(8),sqrt(8)])
yy=xx**2/2.-(1./12.)*xx**4/4.
plt.plot(xx,yy,marker='.',linestyle='none',color='C2')
plt.ylim(0,5)
plt.xlim(-4.5,4.5)
plt.xlabel('$x$')
plt.ylabel('$V(x)$')
plt.legend()
plt.tight_layout(pad=0.1)
#plt.savefig('fig5_3b.pdf',format='pdf')
plt.show()
In [17]:
plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
delta=0.01
x=arange(2,10,delta)
vmg=sqrt(x*x-4)
phi_i=np.cumsum(vmg)*delta
def phiW(x):
    sq=sqrt(x*x-4)
    return x*sq/2.-2*log((sq+x)/2.)
phi=phiW(x)
q=pow(sqrt(2.)-1,2)
lp=pow(1+sqrt(q),2)
lm=pow(1-sqrt(q),2)
vmg_wis=sqrt(abs((x-lp)*(x-lm)))/(x*q)
phi_wis=np.cumsum(vmg_wis)*delta
plt.plot(x,phi,linestyle='--',label='Wigner')
plt.plot(x,phi_wis,label='Wishart')
plt.legend()
plt.xlabel('$x$')
plt.ylabel('$\Phi(x)$')
plt.xlim(2,8)
plt.ylim(0,25)
plt.tight_layout(pad=0.1)
#plt.savefig('fig5_4.pdf',format='pdf')
plt.show()
In [18]:
plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
step=0.01
lp=4./sqrt(3)
x=np.arange(step*0.5,2.5,step)
y=sqrt(np.maximum((lp-x)/x,0))*(lp+2*x)/(4*np.pi)
plt.plot(x,y,label='Conditioned Wigner')
y=sqrt(np.maximum(4-x*x,0))/(2*np.pi)
plt.plot(x,y,linestyle='--',label='Wigner')
plt.xlabel('$\lambda$')
plt.ylabel('$\\rho(\lambda)$')
plt.legend()
plt.xlim(0,2.5)
plt.ylim(0,1.5)
plt.tight_layout(pad=0.1)
#plt.savefig('fig5_5.pdf',format='pdf')
plt.show()
In [19]:
#The 64 zeros of H_64(x) computed by Mathematica
data64=[-14.8862, -13.994, -13.2556, -12.5966, -11.9894, -11.4167, -10.8764, \
-10.3435, -9.88055, -9.35959, -8.89018, -8.38796, -7.96656, -7.51555, \
-7.07057, -6.65082, -6.22236, -5.80063, -5.3833, -4.97007, -4.56046, \
-4.15404, -3.75046, -3.34936, -2.95044, -2.55339, -2.15793, -1.76381, \
-1.37075, -0.978526, -0.586883, -0.195589, 0.195589, 0.586883, \
0.978526, 1.37075, 1.76381, 2.15793, 2.55339, 2.95044, 3.34936, \
3.75046, 4.15404, 4.56045, 4.97009, 5.38331, 5.80054, 6.22305, \
6.6492, 7.06917, 7.52831, 7.9621, 8.41828, 8.89073, 9.36924, 9.81121, \
10.3538, 10.876, 11.4199, 11.9896, 12.5965, 13.2557, 13.9941, 14.8862]
In [20]:
plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
bin=range(-16,18,2)
plt.hist(data64,bin,density=True,color=HistoColor) 
x=np.arange(-16,16.05,.05)
y=sqrt(abs(4*64-x*x))/(2*np.pi*64)
plt.plot(x,y)
plt.xlabel('$\lambda$')
plt.ylabel('$\\rho(\lambda)$')
plt.tight_layout(pad=0.1)
#plt.savefig('fig6_1.pdf',format='pdf')
plt.show()
In [21]:
#The 100 zeros of the Laguerre polynomial L_100^(100) computed by Mathematica
dataL100_100=[20.289, 22.9801, 25.3914, 27.6883, 29.9315, 32.1509, 34.3641, \
36.5823, 38.8132, 41.0623, 43.3339, 45.6313, 47.9571, 50.3134, \
52.7021, 55.1248, 57.5828, 60.0774, 62.6097, 65.1806, 67.7912, \
70.4424, 73.1348, 75.8695, 78.6471, 81.4684, 84.3342, 87.2453, \
90.2022, 93.2059, 96.257, 99.3564, 102.505, 105.703, 108.951, \
112.251, 115.603, 119.008, 122.466, 125.98, 129.549, 133.174, \
136.857, 140.599, 144.4, 148.262, 152.185, 156.172, 160.223, 164.339, \
168.522, 172.773, 177.094, 181.486, 185.951, 190.489, 195.105, \
199.798, 204.571, 209.426, 214.365, 219.39, 224.504, 229.709, \
235.009, 240.405, 245.902, 251.502, 257.208, 263.026, 268.958, \
275.008, 281.183, 287.487, 293.924, 300.502, 307.228, 314.107, \
321.148, 328.36, 335.753, 343.338, 351.126, 359.133, 367.372, \
375.863, 384.626, 393.686, 403.071, 412.817, 422.965, 433.567, \
444.69, 456.42, 468.873, 482.215, 496.694, 512.724, 531.099, 553.883]
In [22]:
plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
bin=range(0,600,20)
plt.hist(dataL100_100,bin,density=True,color=HistoColor) 
q=.5
T=200
lm=T*(1-sqrt(q))**2
lp=T*(1+sqrt(q))**2
dd=(lp-lm)/512.
x=np.arange(lm,lp+dd,dd)
y=sqrt(abs((lp-x)*(x-lm)))/(np.pi*x*q*2*T)
plt.plot(x,y)
plt.xlabel('$\lambda$')
plt.ylabel('$\\rho(\lambda)$')
plt.tight_layout(pad=0.1)
#plt.savefig('fig6_2.pdf',format='pdf')
plt.show()
In [23]:
def jacobi(z,q1,q2): #returns the density of a Jacobi random matrix(q1,q2) at point z
    alp=1-1./q1
    gam=1./q1+1./q2-2.
    arg=(4*(1+gam)-2*alp*gam)*z-alp*alp-pow((2+gam)*z,2)
    return sqrt(np.maximum(arg,0))/(2*np.pi*z*(1-z))

#compute numerically a random Jacobi matrix of size myN
myN=1000
q1=1./5
q2=1./2
w1=(1./q1)*wishart(myN,q1)
w2=(1./q2)*wishart(myN,q2)
ep=w1+w2
sw1=matsqrt(inv(ep))
jac=sw1@w1@sw1
wval= npl.eigvalsh(jac)

plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
plt.hist(wval,40,density=True,color=HistoColor)
x=arange(0.0025,1,.005)
y=jacobi(x,q1,q2)
plt.plot(x,y)
plt.xlim(0,1)
plt.xlabel('$\lambda$')
plt.ylabel('$\\rho(\lambda)$')
plt.tight_layout(pad=0.1)
#plt.savefig('fig7_1.pdf',format='pdf')
plt.show()
In [24]:
plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
x=arange(.0025,1,.005)
q1=1./10.
q2=q1
y=jacobi(x,q1,q2)
plt.plot(x,y,linestyle=':',label="$c=10$")
q1=1./2.
q2=q1
y=jacobi(x,q1,q2)
plt.plot(x,y,linestyle='--',label="$c=2$")
q1=.999
q2=q1
y=jacobi(x,q1,q2)
plt.plot(x,y,label="arcsine law")
plt.xlim(0,1)
plt.ylim(0,3.4)
plt.xlabel('$\lambda$')
plt.ylabel('$\\rho(\lambda)$')
plt.legend()
plt.tight_layout(pad=0.1)
#plt.savefig('fig7_2.pdf',format='pdf')
plt.show()
In [25]:
#Here are the 200 zeros of the polynomial Jacobi_200^(200,600) computed by Mathematica
data_jacobi_200_200_600=[-0.4564354946708801 , -0.4372304554321497 , -0.42118358492456714 ,  -0.4067595135299716 ,\
-0.3933762353651123 , -0.3807359932005732 ,  -0.3686614309336205 , -0.3570365061327438 , -0.345780257620538 , \
-0.33483349216609265 , -0.3241513587454109 ,  -0.3136989037091532 , -0.3034482599163859 , -0.2933767890571437 ,\
-0.28346580868459254 , -0.2736996931120483 ,  -0.2640652218379409 , -0.2545510967967356 ,  -0.24514757774924237 ,\
-0.2358462022141621 ,  -0.22663956710136124 , -0.21752115617160717 ,  -0.2084852020676825 , -0.19952657479474334 ,\
-0.1906406906943753 , -0.1818234374818934 ,  -0.17307111200735964 , -0.16438036819267035 ,  -0.15574817317962417 ,\
-0.14717177015776575 ,  -0.13864864666763954 , -0.13017650742389258 ,  -0.12175325089392135 , -0.1133769490161133 ,\
-0.10504582955779079 , -0.09675826070447353 ,  -0.08851273754476698 , -0.08030787017332866 ,  -0.0721423731811792 ,\
-0.06401505634055128 ,  -0.05592481632237371 , -0.04787062930980953 ,  -0.03985154439212301 , -0.03186667764041304 ,\
-0.023915206781106833 , -0.015996366395100103 ,  -0.00810944358049031 , -0.00025377402532386574 ,  0.007571261556056406 ,\
0.015366240663395668 ,  0.023131701962795594 , 0.030868147981592604 ,  0.038576047572775926 , 0.046255838172611315 ,  \
0.05390792787232026 , 0.06153269732222852 , 0.06913050148468444 ,  0.07670167125020932 , 0.08424651492973909 , \
0.0917653196344154 ,  0.09925835255315635 , 0.1067258621371583 , 0.11416807919953174 ,  0.12158521793743596 , \
0.12897747688333808 , 0.13634503979136595 ,  0.14368807646414314 , 0.1510067435249767 , 0.1583011851398088 ,  \
0.16557153369293115 , 0.17281791042009562 , 0.1800404260023239 ,  0.18723918112342652 , 0.1944142669939751 , \
0.20156576584423508 ,  0.2086937513883516 , 0.21579828926188804 , 0.2228794374346421 ,  0.2299372466005075 , \
0.23697176054600544 , 0.24398301649898235 ,  0.25097104545885135 , 0.25793587250964967 , 0.26487751711708707 , \
0.27179599341067223 , 0.27869131045192386 , 0.2855634724896002 ,  0.2924124792028137 , 0.29923832593283856 , \
0.3060410039043597 ,  0.31282050043686505 , 0.31957679914683484 , 0.32630988014133827 ,  0.33301972020361165 ,\
0.3397062929711555 , 0.3463695691068554 ,  0.3530095164636035 , 0.35962610024286973 , 0.3662192831476472 ,  \
0.3727890255301763 , 0.3793352855348286 , 0.38585801923651697 ,  0.39235718077497994 , 0.3988327224852752 , \
0.405284595024803 ,  0.4117127474971707 , 0.418117127573198 , 0.424497681609356 ,  0.43085435476392464 , \
0.43718709111114895 , 0.44349583375366924 ,  0.4497805249334983 , 0.4560411061418172 , 0.4622775182278601 , \
0.4684897015071616 , 0.47467759586944086 , 0.48084114088640245 ,  0.48698027591973897 , 0.4930949402296273 ,\
0.49918507308402055 ,  0.5052506138690473 , 0.5112915022008443 , 0.5173076780391606 ,  0.5232990818030925 , \
0.529265654489325 , 0.5352073377932791 ,  0.5411240742335883 , 0.5470158072803609 , 0.5528824814877097 ,  \
0.5587240426310727 , 0.5645404378498883 , 0.5703316157962315 ,  0.5760975267900723 , 0.5818381229818733 , \
0.5875533585233088 ,  0.5932431897469611 , 0.5989075753559306 , 0.6045464766243881 ,  0.6101598576102025 ,\
0.6157476853808938 , 0.6213099302542902 ,  0.6268465660554234 , 0.6323575703913572 , 0.6378429249458425 ,  \
0.6433026157959053 , 0.6487366337527244 , 0.654144974729434 ,  0.6595276401388149 , 0.6648846373242017 , \
0.6702159800273607 ,  0.6755216888975816 , 0.6808017920467871 , 0.6860563256561184 ,  0.6912853346402075 , \
0.6964888733762278 , 0.7016670065058384 ,  0.7068198098193358 , 0.7119473712327365 , 0.7170497918701697 ,  \
0.72212718726592 , 0.7271796887027864 , 0.7322074447061907 ,  0.7372106227167815 , 0.7421894109682571 , \
0.7471440206019263 ,  0.7520746880553447 , 0.7569816777694403 , 0.7618652852672222 ,  0.7667258406678364 , \
0.7715637127129586 , 0.7763793133989894 ,  0.7811731033291819 , 0.7859455979259375 , 0.7906973746767121 ,  \
0.7954290816295813 , 0.8001414474096127 , 0.8048352930991168 ,  0.8095115464196567 , 0.8141712587799824 , \
0.8188156259242058 ,  0.8234460131466187 , 0.8280639863604254 , 0.8326713507578033 ,  0.8372701994405017 ,\
0.8418629753315873 , 0.8464525510575275 ,  0.8510423335755142 , 0.855636403555115 , 0.8602397046799187 ,  \
0.864858306521436 , 0.8694997791241686 , 0.874173743244184 ,  0.8788927085123823 , 0.8836734079418276 , \
0.8885390431756565 ,0.8935233392768819 , 0.8986785948325702 , 0.9040939587215077 ,  0.9099465403733097 ,\
0.9167120707867508 ]
In [26]:
def jacobiC(z,q1,q2): #compute the density of a centered-range Jacobi ramdom matrix
    ap=1./q1+1./q2-2.
    am=1./q2-1./q1
    arg=pow(ap+2,2)*z*z+2*ap*am*z+(pow(am,2)-4*(ap+1))
    return sqrt(np.maximum(-arg,0))/(2*np.pi*(1-z*z))
plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
bin=arange(-0.5,1,.025)
plt.hist(data_jacobi_200_200_600,bin,density=True,color=HistoColor) 
q1=1./4.
q2=1./2.
x=arange(-.9925,1,.005)
y=jacobiC(x,q1,q2)
plt.plot(x,y)
plt.xlim(-1,1)
plt.xlabel('$\lambda$')
plt.ylabel('$\\rho(\lambda)$')
plt.tight_layout(pad=0.1)
#plt.savefig('fig7_3.pdf',format='pdf')
plt.show()