In [1]:
# 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 III
#
# This notebook does not contain graphs of Chapter 20 that require financial data
In [2]:
%pylab inline
import scipy.linalg as npl
import scipy.optimize as opt
import scipy.special as special
import scipy.integrate as integrate
from cycler import cycler
from sklearn.isotonic import IsotonicRegression

Populating the interactive namespace from numpy and matplotlib
In [3]:
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('text.latex', preamble=r'\usepackage{mathptmx}\usepackage{amsfonts}')
plt.rc('font', family='serif', serif='New Times Roman', size=10)



#figure sizes in inches
figwidth=9./2.54
figwidth2=6/2.54
figheight=6/2.54
figdpi=200
In [4]:
#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 matsqrt(mat):
    wval, wvect = npl.eigh(mat) #diagonalize mat
    return wvect@np.diag(sqrt(wval))@wvect.transpose() #reconstruct with sqrt e.vals.#helper functions
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)
In [5]:
plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
p=0.5
lp=2*p+1+sqrt(2*(p+1)) #edges of the spectrum for inverse-Wishart
lm=2*p+1-sqrt(2*(p+1))
dd=(lp-lm)/512
x=np.arange(lm,lp+dd,dd) 
y=sqrt(abs((lp-x)*(x-lm)))/(np.pi*x*x*p*2)
plt.plot(x,y,label='inverse-Wishart $p=1/2$')
q=.5
lm=(1-sqrt(q))**2
lp=(1+sqrt(q))**2
dd=(lp-lm)/512
x=np.arange(lm,lp+dd,dd)
yW=sqrt(abs((lp-x)*(x-lm)))/(np.pi*x*q*2)
plt.plot(x,yW,label='white Wishart $q=1/2$')
plt.xlabel('$\lambda$')
plt.ylabel('$\\rho(\lambda)$')
plt.ylim(0,1.5)
plt.xlim(0,4)
plt.legend()
plt.tight_layout(pad=0.1)
#plt.savefig('fig15_1.pdf',format='pdf')
plt.show()
In [6]:
plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
x=np.arange(-3.05,0,.02)
y=2/tanh(x)-1/x
plt.plot(x,y)
plt.plot(-x,-y,color='C0')
plt.xlim(-3,3)
plt.ylim(-12,12)
plt.xlabel('$g$')
plt.ylabel('$\mathfrak{z}=R(g)+1/g$')
plt.axes([0.65,0.30,0.25,0.25])
x=np.arange(1.2,1.8,.01)
y=2/tanh(x)-1/x
plt.plot(x,y)
#plt.savefig('fig15_2.pdf',format='pdf')
plt.show()
In [7]:
def myfunc(x):#returns 2*R(x)+1/x-z it equals zero on a solution
    return 2./tanh(x)-1/x-z
def func_as_reals(x): #transform a complex function into a 2-d vector function of 2 variables
    r1, c1 = x
    #print r1,c1
    a = myfunc(complex(r1, c1))
    return [a.real, a.imag]

random.seed(myseed)
myN=1000 #size of matrix
wval, wvect = npl.eigh(wigner(myN)) #diagonalize a Wigner
data=np.arange(myN)+1
wval=1.+(1.-2.*data)/myN #wval has a flat distribution between -1 and 1
m1=wvect@diag(wval)@wvect.transpose() #reconstruct a matrix using eigenvectors
m2=m1+diag(wval) #take the free sum (m1 has random eigenvector wiz diag(wval))
wval = npl.eigvalsh(m2)

plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
xmin=-1.58 #a starting point below lambda_min
guess=np.array([1/xmin,.0001]) #below lambda_min g(z) behaves as 1/z
x=np.array([])
y=np.array([])
for xx in arange(xmin,0,0.01): #iteratively solve for 2*R(g)+1/g==z starting at z=xmin-i*eta
    z=xx-0.00000001j#z=x-i*eta where eta is a small positive number
    res=opt.root(func_as_reals,guess)
    x=np.append(x,xx)
    y=np.append(y,abs(res.x[1])/np.pi) #the density is Im(g)/pi abs is to ensure the correct root
    guess=res.x #use last value as next guess
x=np.append(x,-x[::-1]) #law is symetric use minus negative values for the x axis
y=np.append(y,y[::-1])
plt.hist(wval,50,density=True,color=HistoColor) #plot histogram of eigenvalue of one sample
plt.plot(x,y) #plot numerically obtained theoretical curve
plt.xlabel('$\lambda$')
plt.ylabel('$\\rho(\lambda)$')
plt.tight_layout(pad=0.1)
#plt.savefig('fig15_3.pdf',format='pdf')
plt.show()
In [8]:
random.seed(myseed)
N=1000
K=40
q=0.5
prod=np.identity(N)
for i in range(K): #multiply K time left and right by sqrt of a Wishart matrix
    mats=matsqrt(wishart(N,q))
    prod=mats@prod@mats
wval = npl.eigvalsh(prod)
data=pow(wval,1./K)
plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
bin=arange(.25,1.25,.01)
plt.hist(data,bin,density=True,color=HistoColor) #plot histogram from sample matrix
x=array([0.5,0.5,1,1])
y=array([0,2,2,0])
plt.plot(x,y) #plot the asymtotic result (K -> Infinity)
tp=(q*(1-K)+sqrt(pow(q*(1-K),2)+4*K*q))/(2*K*q)
tm=(q*(1-K)-sqrt(pow(q*(1-K),2)+4*K*q))/(2*K*q)
zpok=pow((1+tp)/tp,1./K)*(1+q*tp)
zmok=pow((1+tm)/tm,1./K)*(1+q*tm)
y=array([0,10])
x=array([zmok,zmok]) #lower edge of spectrum for finite K
plt.plot(x,y,linestyle="--",color='C0')
x=array([zpok,zpok]) #upper edge of the spectrum for finite K
plt.plot(x,y,linestyle="--",color='C0')
plt.xlim(.4,1.2)
plt.ylim(0,2.5)
plt.xlabel(r'$\mu=\lambda^{1/K}$')
plt.ylabel(r'$\rho(\mu)$')
plt.tight_layout(pad=0.1)
#plt.savefig('fig16_1.pdf',format='pdf')
plt.show()
In [9]:
def cal_ldensity(alpha):
    def myfuncl(t):
        return log(t+1)-log(t)+(alpha*(t+0.5))-lz
    def funcl_as_reals(x):
        r1, c1 = x
        a = myfuncl(complex(r1, c1))
        return [a.real, a.imag]
    lxmax=sqrt(alpha)*sqrt(1+alpha/4.)+2*log(sqrt(alpha)/2+sqrt(1+alpha/4.))
    tmax=(sqrt(1+4./alpha)-1)/2.
    delta=lxmax/100.
    guess=np.array([tmax,.001])
    x=np.array([lxmax])
    y=np.array([0])
    for lx in arange(lxmax-delta,0,-delta):
        lz=lx-0.00000001j
        res=opt.root(funcl_as_reals,guess)
        x=np.append(x,lx)
        y=np.append(y,abs(res.x[1])/np.pi)
        lx-=delta
        guess=res.x
    s=lxmax/2.
    x=np.append(x,-x[::-1])
    y=np.append(y,y[::-1])
    dataz=sqrt(abs(4*s*s-x*x))/(2*np.pi*s*s)
    return x,y,dataz

plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
x,y,dataz=cal_ldensity(100)
plt.plot(x,y,label="free log-normal")
plt.plot(x,dataz,label="Wigner")
plt.legend()
plt.ylabel('$\\rho(l)$')
plt.xlabel('$l$')
ylim(0,0.012)
plt.tight_layout(pad=0.1)
#plt.savefig('fig16_2.pdf',format='pdf')
plt.show()
In [10]:
def density_exp_cor(z=1,q=0.25,b=1.25):
    p=np.zeros(5)
    p[0]=q**2
    p[1]=2*(q**2-b*z*q)
    p[2]=z**2-2*b*z*q+q**2-1
    p[3]=-2
    p[4]=-1
    sol=np.roots(p)
    return max(sol.imag)/z/np.pi
density_exp_corV=np.vectorize(density_exp_cor,excluded=[1,2])
plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
tvar=0.25
b=25
q=tvar/b
x=np.arange(0.01,3.3,.01)
y=density_exp_corV(x,q,b)
plt.plot(x,y,label='$q=0.01\quad b=25$',linestyle=':')
b=1.5
q=tvar/b
x=np.arange(0.01,3.3,.01)
y=density_exp_corV(x,q,b)
plt.plot(x,y,label='$q=0.17\quad b=1.5$',linestyle='--')
plt.legend()
b=1
q=tvar/b
x=np.arange(0.01,3.3,.01)
y=density_exp_corV(x,q,b)
plt.plot(x,y,label='$q=0.25\quad b=1$')
plt.xlabel('$\lambda$')
plt.ylabel('$\\rho(\lambda)$')
plt.xlim(0,2.5)
plt.ylim(0,1.2)
plt.legend()
plt.tight_layout(pad=0.1)
#plt.savefig('fig17_1.pdf',format='pdf')
plt.show()
In [11]:
plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
aa=arange(0,1,.001)
bb=(aa*aa+1)/(1-aa*aa)
q=0.5
qs=1-sqrt(q*q*(bb*bb-1)+1)+bb*q
plot(aa,qs,label="q=0.5")
qs=q*(1+2*aa*aa)
plot(aa,qs,color='C0',linestyle="--")
q=0.1
qs=1-sqrt(q*q*(bb*bb-1)+1)+bb*q
plot(aa,qs,label="q=0.1",color='C1')
qs=q*(1+2*aa*aa)
plot(aa,qs,color='C1',linestyle="--")
plt.legend()
plt.xlabel('$a$')
plt.ylabel('$q^*$')
plt.ylim(0,1.1)
plt.tight_layout(pad=0.1)
#plt.savefig('fig17_2.pdf',format='pdf')
plt.show()
In [12]:
def lim_density_exp_cor(z=1,ss=0.5):
    p=np.zeros(4)
    p[0]=-2*ss*z
    p[1]=z**2-2*ss*z+-1
    p[2]=-2
    p[3]=-1
    sol=np.roots(p)
    return max(sol.imag)/z/np.pi
lim_density_exp_corV=np.vectorize(lim_density_exp_cor,excluded=[1,2])

plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
tvar=0.25
x=np.arange(0.01,5,.01)
y=lim_density_exp_corV(x,tvar)
plt.plot(x,y,label='$\sigma^2=0.25$',linestyle=":")
tvar=0.5
x=np.arange(0.01,5,.01)
y=lim_density_exp_corV(x,tvar)
plt.plot(x,y,label='$\sigma^2=0.5$',linestyle="--")
plt.legend()
tvar=1
x=np.arange(0.01,5,.01)
y=lim_density_exp_corV(x,tvar)
plt.plot(x,y,label='$\sigma^2=1$')
#plt.plot(x,y,label='$b=2 q=0.25$')
plt.xlabel('$\lambda$')
plt.ylabel('$\\rho(\lambda)$')
plt.xlim(0,5)
plt.ylim(0,2)
plt.legend()
plt.tight_layout(pad=0.1)
#plt.savefig('fig17_3.pdf',format='pdf')
plt.show()
In [13]:
def Phip(y,s):
    return exp(s*1j*y)*special.erfc((1+s*1j*y)/sqrt(2))
def func(y):
    return real(1j*(Phip(y,-1)-Phip(y,1))/(Phip(y,-1)+Phip(y,1)))
def func2(y):
    return imag(Phip(y,1))/real(Phip(y,1))
plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
y=arange(-10,10,.01)
xb=y+func2(y)
plt.plot(y,xb)
plt.plot(y,0.475*y,linestyle='--')
plt.xlabel('$y$')
if (plt.rcParams['text.usetex']):
    plt.ylabel('$\mathbb E[x]_y$')
else:
    plt.ylabel('$E[x]_y$')
plt.xlim(-6,6)
plt.ylim(-1,1)
plt.tight_layout(pad=0.1)
#plt.savefig('fig18_1.pdf',format='pdf')
plt.show()
In [14]:
plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
def inv_gamma(x,a,b):
    return pow(b,a)/special.gamma(a)*exp(-b/x)*pow(x,-a-1)
x=arange(.001,3,.01)
y=inv_gamma(x,2,1)
plt.plot(x,y,label="$a=2,b=1$",linestyle=':')
y=inv_gamma(x,2,2)
plt.plot(x,y,linestyle='--',label="$a=2, b=2$")
y=inv_gamma(x,30,30)
plt.plot(x,y,label="$a=30,b=30$")
plt.xlabel('$x$')
plt.ylabel('$P_0(x)$')
plt.xlim(0,3)
plt.legend()
plt.tight_layout(pad=0.1)
#plt.savefig('fig18_2.pdf',format='pdf')
plt.show()
In [15]:
random.seed(myseed)
N=500
T=1000
sigma_n=sqrt(1000.)
error_in=np.array([])
error_out=np.array([])
data_nonz=np.array([])
data_alpha=arange(0.05,2,.05)
true_coef=np.random.laplace(0,1,N)
data_in=np.random.randn(N, T)
data_out=np.random.randn(N, T)
noise=sigma_n*np.random.randn(T)
noise_out=sigma_n*np.random.randn(T)
yyy=true_coef@data_in+noise
yyy_out=true_coef@data_out+noise_out
from sklearn import linear_model
for alpha in data_alpha:
    if(alpha>0):
        clf = linear_model.Lasso(alpha,fit_intercept=False)
    else:
        clf = linear_model.LinearRegression(fit_intercept=False)
    clf.fit(data_in.transpose(),yyy)
    yyy_in_fit=clf.coef_@data_in
    yyy_out_fit=clf.coef_@data_out
    error_in=np.append(error_in,1./T*np.sum((yyy-yyy_in_fit)**2))
    error_out=np.append(error_out,1./T*np.sum((yyy_out-yyy_out_fit)**2))
    data_nonz=np.append(data_nonz,1.-np.count_nonzero(clf.coef_)/float(N))
fig, ax1 = plt.subplots(figsize=(figwidth,figheight),dpi=figdpi)
ax2 = ax1.twinx()
ax1.plot(data_alpha,error_out,label=r'$\mathcal{R}_{\mathrm{out}}^2$')
ax1.plot(data_alpha,error_in,label=r'$\mathcal{R}_{\mathrm{in}}^2$',linestyle='--')
ax1.set_xlabel('$b$')
ax1.set_ylabel(r'$\mathcal{R}^2$')
if (plt.rcParams['text.usetex']):
    ax2.set_ylabel(r'$f_\text{\rm zero}$')
    ax2.plot(data_alpha,data_nonz,label=r'$f_\text{\rm zero}$',linestyle=':')
else:
    ax2.set_ylabel(r'$f_{\rm zero}$')        
    ax2.plot(data_alpha,data_nonz,label=r'$f_{\rm zero}$',linestyle=':')
ax1.set_xlim(0,2)
ax2.set_ylim(0,1)
fig.legend(loc="lower right", bbox_to_anchor=(1,0), bbox_transform=ax1.transAxes)
plt.tight_layout(pad=0.1)
#plt.savefig('fig18_3.pdf',format='pdf')
plt.show()
In [16]:
random.seed(myseed)
plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
myN=2000
M1=wigner(myN)
M2=wigner(myN)
E=M1+M2
Eval,Evect=npl.eigh(E)
Mval,Mvect=npl.eigh(M1)
EvectT=Evect.transpose()
MvectT=Mvect.transpose()
eps=2/sqrt(myN)
target=-1.5
indices=np.where(abs(Mval-target) <eps)
res=pow(EvectT@MvectT[indices].transpose(),2)
fullres=sum(res,axis=1)
fullres*=myN/size(indices)
plt.plot(Eval,fullres,linestyle='none',marker='.',markersize=2,color='C1',label=r'simulation $\mu=-1.5$')
target=0.5
indices=np.where(abs(Mval-target) <eps)
res=pow(EvectT@MvectT[indices].transpose(),2)
fullres=sum(res,axis=1)
fullres*=myN/size(indices)
plt.plot(Eval,fullres,linestyle='none',marker='s',markersize=1,color='C2',label=r'simulation $\mu=0.5$')
target=-1.5
y=1/(pow(target-3*Eval/4,2)+1/2-pow(Eval,2)/16)
plt.plot(Eval,y,color='C0',label=r'theory $\mu=-1.5$')
target=0.5
y=1/(pow(target-3*Eval/4,2)+1/2-pow(Eval,2)/16)
plt.plot(Eval,y,linestyle='--',color='C0',label=r'theory $\mu=0.5$')
plt.legend()
plt.xlabel(r'$\lambda$')
plt.ylabel(r'$\Phi(\lambda,\mu)$')
plt.xlim(-3,3)
plt.ylim(0,6)
plt.tight_layout(pad=0.1)
#plt.savefig('fig19_1.pdf',format='pdf')
plt.show()
In [17]:
def TwishIwish(z,p,q):
    arg=pow(q+1-z,2)-4*(q+z*p)
    return (z-q-1-sqrt(arg)*sign(real(z-1-q-2*p)))/2/(q+z*p)
def RIEwishIwish(z,p,q):
    return real(z/pow(np.abs(1+q*TwishIwish(z,p,q)),2))

plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
z=arange(0.005,4,.01)-1e-10j
plt.plot(real(z),RIEwishIwish(z,.25,.25),label='RIE')
plt.plot(real(z),0.5*real(z)+0.5,linestyle=':',label='linear')
plt.ylim(0,3)
plt.xlim(0.,4)
plt.legend()
plt.xlabel(r'$\lambda$')
plt.ylabel(r'$\xi(\lambda)$')
plt.tight_layout(pad=0.1)
#plt.savefig('fig19_2.pdf',format='pdf')
plt.show()
In [18]:
random.seed(myseed)
sig2=.2
a=sqrt(3*sig2)
q=.40
myN=1000
values=1-a+2*a*arange(0.5/myN,1,1/myN)
C=diag(values)
C12=diag(sqrt(values))
E=C12@wishart(myN,q)@C12
Eval,Evect=npl.eigh(E)
oracle=diag(Evect.transpose()@C@Evect)
lm=Eval[0]
lp=Eval[myN-1]

def tflat(z):
    return (z/(2*a))*log((z-1+a)/(z-1-a))-1
def iterf(t,z):
    return tflat(z/(1+q*t))
zreal=arange(lp,lm-0.005,-.005)
myt=zreal*(0+0j)
t=1/(3.-1e-7j)
for i in range(size(zreal)):
    z=zreal[i]-1e-7j
    t=opt.fixed_point(iterf,t,(z,))
    myt[i]=t

def myrho(x,params,lm,lp):
    res=sqrt(np.maximum((x-lm)*(lp-x),0))*(1+params[0]*x+params[1]*x*x)/(1+params[2]*x+params[3]*x*x)
    return np.maximum(res,0)
def optfunc2(params):
    rho=myrho(Eval,params,Eval[0],Eval[myN-1])
    pcum=integrate.cumtrapz(rho,Eval,initial=0)
    pcum=pcum/pcum[myN-1]
    sampsum=(arange(0.5,myN,1))/myN
    diff=pcum-sampsum
    return sum(diff**2)
ps=np.array((.1,.1,.05,.05))
res=opt.minimize(optfunc2,ps)
rho=myrho(Eval,res.x,Eval[0],Eval[myN-1])
pcum=integrate.cumtrapz(rho,Eval,initial=0)
norm=1/pcum[myN-1]


x=arange(lm,lp+0.01,.01)
rhor=norm*myrho(x,res.x,lm,lp)
plt.figure(figsize=(figwidth2,figheight),dpi=figdpi)
bin=arange(0,3.5,.05)
plt.hist(Eval,bin,density=True,color=HistoColor,label='sample $N=1000$') 
plt.plot(zreal,imag(myt)/zreal/np.pi,linestyle='--',label='theoretical density')
plt.plot(x,rhor,color='C2',label='fitted density')
plt.legend()
plt.xlabel(r'$\lambda$')
plt.ylabel(r'$\rho(\lambda)$')
plt.xlim(0,3.5)
plt.tight_layout(pad=0.1)
#plt.savefig('fig19_3a.pdf',format='pdf')
plt.show()
In [19]:
def integrant(x):
    return norm*myrho(x,res.x,lm,lp)
def hhh(x):
    return integrate.quad(integrant, lm,lp,weight='cauchy',wvar=x)[0]
hhhV=np.vectorize(hhh)
x=arange(lm+0.01,lp,.01)
h=hhhV(x)
rhor2=norm*myrho(x,res.x,lm,lp)

plt.figure(figsize=(figwidth2,figheight),dpi=figdpi)
shrink_th=zreal/pow(np.absolute(1+q*myt),2)
plt.plot(Eval,oracle,linestyle='none',marker='.',markersize=2,color='C2',label=r'simulation oracle')
plt.plot(zreal,shrink_th,linestyle='--',label='theoretical shrinkage')
ggg=-h+1j*np.pi*rhor2
ttt=x*ggg-1
shrink=x/pow(np.absolute(1+q*ttt),2)
plt.plot(x,shrink,label='shrinkage from fit')
plt.legend()
plt.xlabel(r'$\lambda$')
plt.ylabel(r'$\xi(\lambda)$')
plt.xlim(0,3.5)
plt.ylim(0,1.6)
plt.tight_layout(pad=0.1)
#plt.savefig('fig19_3b.pdf',format='pdf')
plt.show()
In [20]:
def g_wig(lam,x,eta):
    z=x-lam-1e-7j
    return 0.5*z*(1-sqrt(1-4*pow(eta/z,2)))/pow(eta,2)
def g_wigFull(x):
    eta=1/sqrt(myN)
    return average(g_wig(Eval[newaxis,:],x[:,newaxis],eta),axis=1)

x=arange(0,4,.02)
ggg1=myg(Eval,x-1j/sqrt(myN)) #Sample Stieltjes transform evaluated with imaginary part =-1/sqrt(N)
gggw=g_wigFull(x)

plt.figure(figsize=(figwidth2,figheight),dpi=figdpi)
plt.plot(zreal,imag(myt)/zreal/np.pi,linestyle='--',label='theoretical density')
plt.plot(x,imag(gggw)/np.pi,color='C1',linestyle=':',label='Wigner kernel')
plt.plot(x,imag(ggg1)/np.pi,color='C2',label='Cauchy kernel')
plt.legend()
plt.xlabel(r'$\lambda$')
plt.ylabel(r'$\rho(\lambda)$')
plt.xlim(0,3.5)
plt.ylim(0,1.2)
plt.tight_layout(pad=0.1)
#plt.savefig('fig19_4a.pdf',format='pdf')
plt.show()
In [21]:
plt.figure(figsize=(figwidth2,figheight),dpi=figdpi)
ttt1=x*ggg1-1
shrink1=x/pow(np.absolute(1+q*ttt1),2)
tttw=x*gggw-1
shrinkw=x/pow(np.absolute(1+q*tttw),2)
ir = IsotonicRegression()
shrinkw_iso=ir.fit_transform(x,shrinkw)
plt.plot(zreal,shrink_th,linestyle='--',label='theoretical shrinkage')
plt.plot(x[3:160],shrinkw[3:160],color='C1',linestyle=':',label='Wigner kernel')
plt.plot(x[3:160],shrinkw_iso[3:160],color='C2',label='isotonic Wigner kernel')
plt.legend(loc='lower right')
plt.xlabel(r'$\lambda$')
plt.ylabel(r'$\xi(\lambda)$')
plt.xlim(0,3.5)
plt.ylim(0,1.6)
plt.tight_layout(pad=0.1)
#plt.savefig('fig19_4b.pdf',format='pdf')
plt.show()
In [22]:
random.seed(myseed)
q=0.4
nblock=10
Tblock=int(myN/q/nblock)
hhh0=np.random.randn(nblock,Tblock,myN)
hhhC=matmul(hhh0,C12)
allin=np.array([])
allout=np.array([])
for j in range(nblock):
        Emat=np.zeros((myN,myN))
        for i in range(nblock):
            if (i==j):
                continue
            Emat=Emat+matmul(hhhC[i].transpose(),hhhC[i])
        Emat=Emat/((nblock-1)*Tblock)
        Eval1,Evect1=eig(Emat)
        Emat2=hhhC[j].transpose()@hhhC[j]/Tblock
        xval=diag(Evect1.transpose()@Emat2@Evect1)
        allin=np.append(allin,Eval1)
        allout=np.append(allout,xval)

ind=np.argsort(allin)
allin_s=allin[ind]
allout_s=allout[ind]
ir = IsotonicRegression()
allout_iso=ir.fit_transform(allin_s,allout_s)


plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
plt.plot(allin,allout,linestyle='none',marker='.',markersize=1,color='C2',label='cross-validation')
plt.plot(zreal,shrink_th,linestyle='--',label='theoretical shrinkage',color='C0')
plt.plot(allin_s,allout_iso,color='C1',label='isotonic cross-validation')
plt.legend()
plt.xlabel(r'$\lambda$')
plt.ylabel(r'$\xi(\lambda)$')
plt.xlim(0,3.4)
plt.ylim(0.3,1.8)
plt.tight_layout(pad=0.1)
#plt.savefig('fig19_5.pdf',format='pdf')
plt.show()
In [23]:
plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
q=0.5
p=0.5
N=100
truef=1./(N*(1+p))
truef=1
Einf=truef*(1-q)
Eoutf=truef/(1-q)
Xiinf=truef*(1-p*q/(p+q*1+p))
Xioutf=truef*(1+p*q*q/((p+q)*(1+p)))
Xioutf=truef*(1+p*q/(p+q))
truef,Einf,Eoutf,Xiinf,Xioutf
datag=arange(0,2,.01)
data=Einf*datag*datag
plot(data,datag,label="$\mathcal{R}_{\mathrm{in}}^2(\mathbf{E})$",color='C1',linestyle=":")
data=Xiinf*datag*datag
plot(data,datag,label="$\mathcal{R}_{\mathrm{in}}^2(\Xi)$",color='C2',linestyle=":")
data=truef*datag*datag
plot(data,datag,label="$\mathcal{R}_{\mathrm{true}}^2$",color='C0')
data=Xioutf*datag*datag
plot(data,datag,label="$\mathcal{R}_{\mathrm{out}}^2(\Xi)$",color='C2')
data=Eoutf*datag*datag
plot(data,datag,label="$\mathcal{R}_{\mathrm{out}}^2(\mathbf{E})$",color='C1')
plt.xlabel('$\mathcal{R}^2$')
plt.ylabel('$\mathcal{G}$')
plt.xlim(0,2)
plt.ylim(0,1.5)
plt.legend()
plt.tight_layout(pad=0.1)
#plt.savefig('fig20_1.pdf',format='pdf')
plt.show()
In [24]:
z=16
lim=6
def airy_int(t):
    return 1j*(t+pow(t,3)/(3*z))
x = np.linspace(-6, 6, 60)
y = np.linspace(-6, 6, 60)
X, Y = np.meshgrid(x, y)
Z = real(airy_int(X+1j* Y))
Z2 = imag(airy_int(X+1j* Y))
tst=sqrt(-z+0j)

plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
C=plt.contour(X, Y, Z,(-12,-9,-6,-4,-3,-2,-1,0,1,2,3,4,6,9,12), colors='C0');
C1=plt.contour(X, Y, Z2,(0,), colors='C1');
plt.plot((real(tst)),(imag(tst)),marker='o',color='C0',markersize=3)
plt.plot((real(tst)),-(imag(tst)),marker='s',color='C0',markersize=3)
plt.clabel(C1, inline=True, fontsize=8,fmt='%2.0fi')
plt.clabel(C, inline=True, fontsize=8,fmt='%2.0f')
plt.xlabel(r'$x$')
plt.ylabel(r'$y$')
plt.tight_layout(pad=0.1)
#plt.savefig('figA_1.pdf',format='pdf')
plt.show()
In [25]:
plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
ymax=6
a=.25
npt=2000
lp=(1+a)/(1-a)
lm=1/lp
delta=(lp-lm)/npt
x=np.arange(lm+delta,lp,delta)
y=1/sqrt(abs((lp-x)*(x-lm)))/(np.pi*x)
plt.plot(x,y,label='$a=0.25$')
x=[lm,lm]
y=[0,ymax]
plt.plot(x,y,linestyle='dashed',color='C0')
x=[lp,lp]
y=[0,ymax]
plt.plot(x,y,linestyle='dashed',color='C0')
a=.5
lp=(1+a)/(1-a)
lm=1/lp
delta=(lp-lm)/npt
x=np.arange(lm+delta,lp,delta)
y=1/sqrt(abs((lp-x)*(x-lm)))/(np.pi*x)
y[0]=ymax
x[npt-2]=lp
y[npt-2]=ymax
plt.plot(x,y,label='$a=0.5$')
x=[lm,lm]
y=[0,ymax]
plt.plot(x,y,linestyle='dashed',color='C1')
x=[lp,lp]
y=[0,ymax]
plt.plot(x,y,linestyle='dashed',color='C1')
a=.75
lp=(1+a)/(1-a)
lm=1/lp
delta=(lp-lm)/npt
x=np.arange(lm+delta,lp,delta)
y=1/sqrt(abs((lp-x)*(x-lm)))/(np.pi*x)
y[0]=ymax
x[npt-1]=lp
y[npt-1]=ymax
plt.plot(x,y,label='$a=0.75$')
x=[lm,lm]
y=[0,ymax]
plt.plot(x,y,linestyle='dashed',color='C2')
x=[lp,lp]
y=[0,ymax]
plt.plot(x,y,linestyle='dashed',color='C2')
plt.xlabel('$\lambda$')
plt.ylabel('$\\rho(\lambda)$')
#plt.xlim(-3,3)
plt.ylim(0,ymax)
plt.legend(loc=(.6,.65))
plt.tight_layout(pad=0.1)
#plt.savefig('figA_2.pdf',format='pdf')
plt.show()
In [ ]: