#### Theorem 1 #### for ising models on a lattice in a high temp regime, sample covariance and correlation have the same asymptotic distribution N(0,\nu^2) with \nu^2 increasing with the inverse temp parameter beta install.packages("potts",repos="http://cran.us.r-project.org") install.packages("latex2exp",repos="http://cran.us.r-project.org") library(potts) library(latex2exp) ######### NOTES ############ # see https://cran.r-project.org/web/packages/potts/potts.pdf # # q=2 # critical point at log(1+\sqrt(q))=0.8813736 ######### SIMS ######### ################ plots for paper ######################## beta_vec<-c(0.1,0.3,0.5,0.7,0.8,0.83,0.85,0.87,0.8813736) nsim=1000 nrow<-100 ncol<-100 n=nrow*ncol covariance<-matrix(999,nrow=length(beta_vec),ncol=nsim) correlation<-matrix(999,nrow=length(beta_vec),ncol=nsim) for(b in 1:length(beta_vec)){ empirical_corr<-rep(9,nsim) empirical_cov<-rep(9,nsim) ncolor<-2 #beta<-log(1+sqrt(ncolor)) #CUTOFF = 0.8813736 beta<-beta_vec[b] theta<-c(rep(0,ncolor),beta) for(l in 1:nsim){ x<-matrix(1,nrow=nrow,ncol=ncol) foo<-packPotts(x,ncolor) out1<- (-1)^unpackPotts(potts(foo,theta,boundary="free",nbatch=50,)$final) out2<- (-1)^unpackPotts(potts(foo,theta,boundary="free",nbatch=50,)$final) empirical_corr[l]<-cor(as.vector(out1),as.vector(out2)) empirical_cov[l]<-cov(as.vector(out1),as.vector(out2)) #print(l) } correlation[b,]<-sqrt(n)*empirical_corr covariance[b,]<-sqrt(n)*empirical_cov print(beta) } typeIerror<-rep(999, length(b)) for(b in 1:length(beta_vec)){ typeIerror[b]<-1-sum(abs(correlation[b,])<1.96)/length(correlation[b,]) } pdf("ising_lattice_corr.pdf") #note: headers are cut off when saving directly, instead just run par(mfrow=c(3,3)) and save as pdf par(mfrow=c(3,3)) hist(correlation[1,],main=paste("\U03B2 =",beta_vec[1]*2,", sd =",round(sqrt(var(correlation[1,])),digits=2), ", \U03B1 =",typeIerror[1]), xlab=expression(hat(rho))) hist(correlation[2,],main=paste("\U03B2 =",beta_vec[2]*2,", sd =",round(sqrt(var(correlation[2,])),digits=2), ", \U03B1 =",typeIerror[2]), xlab=expression(hat(rho))) hist(correlation[3,],main=paste("\U03B2 =",beta_vec[3]*2,", sd =",round(sqrt(var(correlation[3,])),digits=2), ", \U03B1 =",typeIerror[3]), xlab=expression(hat(rho))) hist(correlation[4,],main=paste("\U03B2 =",beta_vec[4]*2,", sd =",round(sqrt(var(correlation[4,])),digits=2), ", \U03B1 =",typeIerror[4]), xlab=expression(hat(rho))) hist(correlation[5,],main=paste("\U03B2 =",beta_vec[5]*2,", sd =",round(sqrt(var(correlation[5,])),digits=2), ", \U03B1 =",typeIerror[5]), xlab=expression(hat(rho))) hist(correlation[6,],main=paste("\U03B2 =",beta_vec[6]*2,", sd =",round(sqrt(var(correlation[6,])),digits=2), ", \U03B1 =",typeIerror[6]), xlab=expression(hat(rho))) hist(correlation[7,],main=paste("\U03B2 =",beta_vec[7]*2,", sd =",round(sqrt(var(correlation[7,])),digits=2), ", \U03B1 =",typeIerror[7]), xlab=expression(hat(rho))) hist(correlation[8,],main=paste("\U03B2 =",beta_vec[8]*2,", sd =",round(sqrt(var(correlation[8,])),digits=2), ", \U03B1 =",typeIerror[8]), xlab=expression(hat(rho))) hist(correlation[9,],main=paste("\U03B2=\U03B2_c, sd =",round(sqrt(var(correlation[9,])),digits=2), ", \U03B1 =",typeIerror[9]), xlab=expression(hat(rho))) dev(off) typeIcov<-rep(999, length(b)) for(b in 1:length(beta_vec)){ typeIcov[b]<-1-sum(abs(correlation[b,])<1.96)/length(correlation[b,]) } pdf("ising_lattice_covariance.pdf") par(mfrow=c(3,3)) hist(covariance[1,],main=paste("\U03B2 =",beta_vec[1]*2,", sd=",round(sqrt(var(covariance[1,])),digits=2), ", \U03B1 =",typeIcov[1]), xlab=expression(hat(T))) hist(covariance[2,],main=paste("\U03B2 =",beta_vec[2]*2,", sd=",round(sqrt(var(covariance[2,])),digits=2), ", \U03B1 =",typeIcov[2]), xlab=expression(hat(T))) hist(covariance[3,],main=paste("\U03B2 =",beta_vec[3]*2,", sd=",round(sqrt(var(covariance[3,])),digits=2), ", \U03B1 =",typeIcov[3]), xlab=expression(hat(T))) hist(covariance[4,],main=paste("\U03B2 =",beta_vec[4]*2,", sd=",round(sqrt(var(covariance[4,])),digits=2), ", \U03B1 =",typeIcov[4]), xlab=expression(hat(T))) hist(covariance[5,],main=paste("\U03B2 =",beta_vec[5]*2,", sd=",round(sqrt(var(covariance[5,])),digits=2), ", \U03B1 =",typeIcov[5]), xlab=expression(hat(T))) hist(covariance[6,],main=paste("\U03B2 =",beta_vec[6]*2,", sd=",round(sqrt(var(covariance[6,])),digits=2), ", \U03B1 =",typeIcov[6]), xlab=expression(hat(T))) hist(covariance[7,],main=paste("\U03B2 =",beta_vec[7]*2,", sd=",round(sqrt(var(covariance[7,])),digits=2), ", \U03B1 =",typeIcov[7]), xlab=expression(hat(T))) hist(covariance[8,],main=paste("\U03B2 =",beta_vec[8]*2,", sd=",round(sqrt(var(covariance[8,])),digits=2), ", \U03B1 =",typeIcov[8]), xlab=expression(hat(T))) hist(covariance[9,],main=paste("\U03B2=\U03B2_c, sd=",round(sqrt(var(covariance[9,])),digits=2), ", \U03B1 =",typeIcov[9]), xlab=expression(hat(T))) dev.off() par(mfrow=c(3,3)) hist(covariance[1,],main=paste0("beta=",beta_vec[1],", sd=",round(sqrt(var(covariance[1,])),digits=2))) hist(covariance[2,] ,main=paste0("beta=",beta_vec[2],", sd=",round(sqrt(var(covariance[2,])),digits=2))) hist(covariance[3,] ,main=paste0("beta=",beta_vec[3],", sd=",round(sqrt(var(covariance[3,])),digits=2))) hist(covariance[4,] ,main=paste0("beta=",beta_vec[4],", sd=",round(sqrt(var(covariance[4,])),digits=2))) hist(covariance[5,] ,main=paste0("beta=",beta_vec[5],", sd=",round(sqrt(var(covariance[5,])),digits=2))) hist(covariance[6,] ,main=paste0("beta=",beta_vec[6],", sd=",round(sqrt(var(covariance[6,])),digits=2))) hist(covariance[7,] ,main=paste0("beta=",beta_vec[7],", sd=",round(sqrt(var(covariance[7,])),digits=2))) hist(covariance[8,] ,main=paste0("beta=",beta_vec[8],", sd=",round(sqrt(var(covariance[8,])),digits=2))) hist(covariance[9,] ,main=paste0("beta=",beta_vec[9],", sd=",round(sqrt(var(covariance[9,])),digits=2))) dev.off() ############ test plots ############################# beta_vec<-c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.81,0.83,0.85,0.86,0.87,0.88,0.8813736,0.9,0.95) nsim=1000 nrow<-100 ncol<-100 n=nrow*ncol covariance<-matrix(999,nrow=length(beta_vec),ncol=nsim) correlation<-matrix(999,nrow=length(beta_vec),ncol=nsim) for(b in 1:length(beta_vec)){ empirical_corr<-rep(9,nsim) empirical_cov<-rep(9,nsim) ncolor<-2 #beta<-log(1+sqrt(ncolor)) #CUTOFF = 0.8813736 beta<-beta_vec[b] theta<-c(rep(0,ncolor),beta) for(l in 1:nsim){ x<-matrix(1,nrow=nrow,ncol=ncol) foo<-packPotts(x,ncolor) out1<- (-1)^unpackPotts(potts(foo,theta,boundary="free",nbatch=50,)$final) out2<- (-1)^unpackPotts(potts(foo,theta,boundary="free",nbatch=50,)$final) empirical_corr[l]<-cor(as.vector(out1),as.vector(out2)) empirical_cov[l]<-cov(as.vector(out1),as.vector(out2)) #print(l) } correlation[b,]<-sqrt(n)*empirical_corr covariance[b,]<-sqrt(n)*empirical_cov print(beta) #pdf(paste0("beta=",beta,"_lattice_corrXY.pdf")) #pdf("beta=0.6_lattice_corrXY.pdf") # Creating a plot #hist(sqrt(n)*empirical_corr) # Closing the graphical device #dev.off() } pdf(paste0("Thm 3.3 (correlation) n=",n,".pdf")) par(mfrow=c(3,3)) hist(correlation[1,],main=paste0("beta=",beta_vec[1],", var=",round(\var(correlation[1,]),digits=2))) hist(correlation[2,] ,main=paste0("beta=",beta_vec[2],", var=",round(var(correlation[2,]),digits=2))) hist(correlation[3,] ,main=paste0("beta=",beta_vec[3],", var=",round(var(correlation[3,]),digits=2))) hist(correlation[4,] ,main=paste0("beta=",beta_vec[4],", var=",round(var(correlation[4,]),digits=2))) hist(correlation[5,] ,main=paste0("beta=",beta_vec[5],", var=",round(var(correlation[5,]),digits=2))) hist(correlation[6,] ,main=paste0("beta=",beta_vec[6],", var=",round(var(correlation[6,]),digits=2))) hist(correlation[7,] ,main=paste0("beta=",beta_vec[7],", var=",round(var(correlation[7,]),digits=2))) hist(correlation[8,] ,main=paste0("beta=",beta_vec[8],", var=",round(var(correlation[8,]),digits=2))) hist(correlation[9,] ,main=paste0("beta=",beta_vec[9],", var=",round(var(correlation[9,]),digits=2))) par(mfrow=c(3,3)) hist(correlation[10,] ,main=paste0("beta=",beta_vec[10],", var=",round(var(correlation[10,]),digits=2))) hist(correlation[11,] ,main=paste0("beta=",beta_vec[11],", var=",round(var(correlation[11,]),digits=2))) hist(correlation[12,] ,main=paste0("beta=",beta_vec[12],", var=",round(var(correlation[12,]),digits=2))) hist(correlation[13,] ,main=paste0("beta=",beta_vec[13],", var=",round(var(correlation[13,]),digits=2))) hist(correlation[14,] ,main=paste0("beta=",beta_vec[14],", var=",round(var(correlation[14,]),digits=2))) hist(correlation[15,] ,main=paste0("beta=",beta_vec[15],", var=",round(var(correlation[15,]),digits=2))) hist(correlation[16,] ,main=paste0("beta=",beta_vec[16],", var=",round(var(correlation[16,]),digits=2))) hist(correlation[17,] ,main=paste0("beta=",beta_vec[17],", var=",round(var(correlation[17,]),digits=2))) hist(correlation[18,] ,main=paste0("beta=",beta_vec[18],", var=",round(var(correlation[18,]),digits=2))) dev.off() pdf(paste0("Thm 3.3 (covariance) n=",n,".pdf")) par(mfrow=c(3,3)) hist(covariance[1,],main=paste0("beta=",beta_vec[1],", var=",round(var(covariance[1,]),digits=2))) hist(covariance[2,] ,main=paste0("beta=",beta_vec[2],", var=",round(var(covariance[2,]),digits=2))) hist(covariance[3,] ,main=paste0("beta=",beta_vec[3],", var=",round(var(covariance[3,]),digits=2))) hist(covariance[4,] ,main=paste0("beta=",beta_vec[4],", var=",round(var(covariance[4,]),digits=2))) hist(covariance[5,] ,main=paste0("beta=",beta_vec[5],", var=",round(var(covariance[5,]),digits=2))) hist(covariance[6,] ,main=paste0("beta=",beta_vec[6],", var=",round(var(covariance[6,]),digits=2))) hist(covariance[7,] ,main=paste0("beta=",beta_vec[7],", var=",round(var(covariance[7,]),digits=2))) hist(covariance[8,] ,main=paste0("beta=",beta_vec[8],", var=",round(var(covariance[8,]),digits=2))) hist(covariance[9,] ,main=paste0("beta=",beta_vec[9],", var=",round(var(covariance[9,]),digits=2))) par(mfrow=c(3,3)) hist(covariance[10,] ,main=paste0("beta=",beta_vec[10],", var=",round(var(covariance[10,]),digits=2))) hist(covariance[11,] ,main=paste0("beta=",beta_vec[11],", var=",round(var(covariance[11,]),digits=2))) hist(covariance[12,] ,main=paste0("beta=",beta_vec[12],", var=",round(var(covariance[12,]),digits=2))) hist(covariance[13,] ,main=paste0("beta=",beta_vec[13],", var=",round(var(covariance[13,]),digits=2))) hist(covariance[14,] ,main=paste0("beta=",beta_vec[14],", var=",round(var(covariance[14,]),digits=2))) hist(covariance[15,] ,main=paste0("beta=",beta_vec[15],", var=",round(var(covariance[15,]),digits=2))) hist(covariance[16,] ,main=paste0("beta=",beta_vec[16],", var=",round(var(covariance[16,]),digits=2))) hist(covariance[17,] ,main=paste0("beta=",beta_vec[17],", var=",round(var(covariance[17,]),digits=2))) hist(covariance[18,] ,main=paste0("beta=",beta_vec[18],", var=",round(var(covariance[18,]),digits=2))) dev.off() ######### exploratory sims 9/19/23 ######### beta_vec<-c(0.62,0.63,0.64,0.65,0.66,0.67,0.68,0.69,0.7,0.71,0.72,0.73,0.74,0.75,0.76,0.77,0.78,0.79) nsim=1000 nrow<-1000 ncol<-1000 n=nrow*ncol covariance<-matrix(999,nrow=length(beta_vec),ncol=nsim) correlation<-matrix(999,nrow=length(beta_vec),ncol=nsim) for(b in 1:length(beta_vec)){ empirical_corr<-rep(9,nsim) empirical_cov<-rep(9,nsim) ncolor<-2 #beta<-log(1+sqrt(ncolor)) #CUTOFF = 0.8813736 beta<-beta_vec[b] theta<-c(rep(0,ncolor),beta) for(l in 1:nsim){ x<-matrix(1,nrow=nrow,ncol=ncol) foo<-packPotts(x,ncolor) out1<- (-1)^unpackPotts(potts(foo,theta,boundary="free",nbatch=50,)$final) out2<- (-1)^unpackPotts(potts(foo,theta,boundary="free",nbatch=50,)$final) empirical_corr[l]<-cor(as.vector(out1),as.vector(out2)) empirical_cov[l]<-cov(as.vector(out1),as.vector(out2)) #print(l) } correlation[b,]<-sqrt(n)*empirical_corr covariance[b,]<-sqrt(n)*empirical_cov print(beta) #pdf(paste0("beta=",beta,"_lattice_corrXY.pdf")) #pdf("beta=0.6_lattice_corrXY.pdf") # Creating a plot #hist(sqrt(n)*empirical_corr) # Closing the graphical device #dev.off() } pdf(paste0("Thm 3.3, correlation, beta near .7, n=",n,".pdf")) par(mfrow=c(3,3)) hist(correlation[1,],main=paste0("beta=",beta_vec[1],", var=",round(var(correlation[1,]),digits=2))) hist(correlation[2,] ,main=paste0("beta=",beta_vec[2],", var=",round(var(correlation[2,]),digits=2))) hist(correlation[3,] ,main=paste0("beta=",beta_vec[3],", var=",round(var(correlation[3,]),digits=2))) hist(correlation[4,] ,main=paste0("beta=",beta_vec[4],", var=",round(var(correlation[4,]),digits=2))) hist(correlation[5,] ,main=paste0("beta=",beta_vec[5],", var=",round(var(correlation[5,]),digits=2))) hist(correlation[6,] ,main=paste0("beta=",beta_vec[6],", var=",round(var(correlation[6,]),digits=2))) hist(correlation[7,] ,main=paste0("beta=",beta_vec[7],", var=",round(var(correlation[7,]),digits=2))) hist(correlation[8,] ,main=paste0("beta=",beta_vec[8],", var=",round(var(correlation[8,]),digits=2))) hist(correlation[9,] ,main=paste0("beta=",beta_vec[9],", var=",round(var(correlation[9,]),digits=2))) par(mfrow=c(3,3)) hist(correlation[10,] ,main=paste0("beta=",beta_vec[10],", var=",round(var(correlation[10,]),digits=2))) hist(correlation[11,] ,main=paste0("beta=",beta_vec[11],", var=",round(var(correlation[11,]),digits=2))) hist(correlation[12,] ,main=paste0("beta=",beta_vec[12],", var=",round(var(correlation[12,]),digits=2))) hist(correlation[13,] ,main=paste0("beta=",beta_vec[13],", var=",round(var(correlation[13,]),digits=2))) hist(correlation[14,] ,main=paste0("beta=",beta_vec[14],", var=",round(var(correlation[14,]),digits=2))) hist(correlation[15,] ,main=paste0("beta=",beta_vec[15],", var=",round(var(correlation[15,]),digits=2))) hist(correlation[16,] ,main=paste0("beta=",beta_vec[16],", var=",round(var(correlation[16,]),digits=2))) hist(correlation[17,] ,main=paste0("beta=",beta_vec[17],", var=",round(var(correlation[17,]),digits=2))) hist(correlation[18,] ,main=paste0("beta=",beta_vec[18],", var=",round(var(correlation[18,]),digits=2))) dev.off() pdf(paste0("Thm 3.3, covariance, beta near .7, n=",n,".pdf")) par(mfrow=c(3,3)) hist(covariance[1,],main=paste0("beta=",beta_vec[1],", var=",round(var(covariance[1,]),digits=2))) hist(covariance[2,] ,main=paste0("beta=",beta_vec[2],", var=",round(var(covariance[2,]),digits=2))) hist(covariance[3,] ,main=paste0("beta=",beta_vec[3],", var=",round(var(covariance[3,]),digits=2))) hist(covariance[4,] ,main=paste0("beta=",beta_vec[4],", var=",round(var(covariance[4,]),digits=2))) hist(covariance[5,] ,main=paste0("beta=",beta_vec[5],", var=",round(var(covariance[5,]),digits=2))) hist(covariance[6,] ,main=paste0("beta=",beta_vec[6],", var=",round(var(covariance[6,]),digits=2))) hist(covariance[7,] ,main=paste0("beta=",beta_vec[7],", var=",round(var(covariance[7,]),digits=2))) hist(covariance[8,] ,main=paste0("beta=",beta_vec[8],", var=",round(var(covariance[8,]),digits=2))) hist(covariance[9,] ,main=paste0("beta=",beta_vec[9],", var=",round(var(covariance[9,]),digits=2))) par(mfrow=c(3,3)) hist(covariance[10,] ,main=paste0("beta=",beta_vec[10],", var=",round(var(covariance[10,]),digits=2))) hist(covariance[11,] ,main=paste0("beta=",beta_vec[11],", var=",round(var(covariance[11,]),digits=2))) hist(covariance[12,] ,main=paste0("beta=",beta_vec[12],", var=",round(var(covariance[12,]),digits=2))) hist(covariance[13,] ,main=paste0("beta=",beta_vec[13],", var=",round(var(covariance[13,]),digits=2))) hist(covariance[14,] ,main=paste0("beta=",beta_vec[14],", var=",round(var(covariance[14,]),digits=2))) hist(covariance[15,] ,main=paste0("beta=",beta_vec[15],", var=",round(var(covariance[15,]),digits=2))) hist(covariance[16,] ,main=paste0("beta=",beta_vec[16],", var=",round(var(covariance[16,]),digits=2))) hist(covariance[17,] ,main=paste0("beta=",beta_vec[17],", var=",round(var(covariance[17,]),digits=2))) hist(covariance[18,] ,main=paste0("beta=",beta_vec[18],", var=",round(var(covariance[18,]),digits=2))) dev.off() ########## MRF sims ####### library(MASS) library(latex2exp) n=200 v=c(1/sqrt(2), -1/sqrt(2), rep(0,n-2)) lambda=diag(rep(1,n))+n^3*v%*%t(v) nsim=1000 correlation=rep(99,nsim) covariance=rep(99,nsim) for(i in 1:nsim){ x=mvrnorm(1,rep(0,n),lambda) y=mvrnorm(1,rep(0,n),lambda) correlation[i]=cor(x,y) covariance[i]=cov(x,y) } hist(correlation, main=paste("histogram of","\U03C1", "estimated with n=200")) A=rchisq(1000,1) B=rnorm(1000,0,1) hist(A*B) var(A*B) ##### COROLLARY 1 to Thm 4 n=200 nsim=1000 rho_vec<-c(.1,.2,.3,.4,.5,.6) covariance<-matrix(999,nrow=length(rho_vec),ncol=nsim) correlation<-matrix(999,nrow=length(rho_vec),ncol=nsim) for(i in 1:length(rho_vec)){ rho<-rho_vec[i] sigma=matrix(rho,nrow=n,ncol=n)+diag(1-rho,nrow=n,ncol=n) empirical_corr=rep(99,nsim) empirical_cov=rep(99,nsim) for(i in 1:nsim){ x=mvrnorm(1,rep(0,n),sigma) y=mvrnorm(1,rep(0,n),sigma) empirical_corr[i]=cor(x,y) empirical_cov[i]=cov(x,y) } correlation[i,]<-sqrt(n)*empirical_corr covariance[i,]<-sqrt(n)*empirical_cov print(rho) } par(mfrow=c(2,3)) hist(correlation[1,],main=paste("histogram of sample correlation,","\U03C1=",rho_vec[1]*2), xlab=expression(hat(rho))) hist(correlation[2,],main=paste("\U03B2 =",beta_vec[2]*2,", sd=",round(sqrt(var(correlation[2,])),digits=2)), xlab=expression(hat(rho))) hist(correlation[3,],main=paste("\U03B2 =",beta_vec[3]*2,", sd=",round(sqrt(var(correlation[3,])),digits=2)), xlab=expression(hat(rho))) hist(correlation[4,],main=paste("\U03B2 =",beta_vec[4]*2,", sd=",round(sqrt(var(correlation[4,])),digits=2)), xlab=expression(hat(rho))) hist(correlation[5,],main=paste("\U03B2 =",beta_vec[5]*2,", sd=",round(sqrt(var(correlation[5,])),digits=2)), xlab=expression(hat(rho))) hist(correlation[6,],main=paste("\U03B2 =",beta_vec[6]*2,", sd=",round(sqrt(var(correlation[6,])),digits=2)), xlab=expression(hat(rho))) hist(sqrt(n)*correlation, main="histogram of sample correlation estimated with n=200") sd(sqrt(n)*correlation) hist(sqrt(n)*covariance, main="histogram of sample covariance estimated with n=200") sd(sqrt(n)*covariance) hist(rnorm(1000,mean=0,sd=.5))