# 1D Example set.seed(88) n=25 x=seq(0,1,length=n) X=abs(outer(x,x,"-")) rho=0.1 R=rho^(X^2) L=t(chol(R+diag(n)*.Machine$double.eps*100)) mu=0 Z=rnorm(n,mean=0,sd=1) Y=L%*%Z+mu # plot plot(x,Y,xlab="x",ylab="Y(x)",type='b',lwd=2,pch=20) # Generate a few realizations m=10 Ymat=matrix(0,nrow=m,ncol=n) for(i in 1:m) { Z=rnorm(n,mean=0,sd=1) Ymat[i,]=L%*%Z+mu } plot(x,Ymat[1,],xlab="x",ylab="Y(x)",type='b', lwd=2,pch=20,ylim=range(Ymat)) for(i in 2:m) lines(x,Ymat[i,],col=i,lwd=2) # Change the correlation parameter set.seed(88) rho=1e-4 R=rho^(X^2) L=t(chol(R+diag(n)*.Machine$double.eps*100)) mu=0 Z=rnorm(n,mean=0,sd=1) Y=L%*%Z+mu plot(x,Y,xlab="x",ylab="Y(x)",type='b',lwd=2,pch=20) # 2D Example set.seed(88) n=25 x=as.matrix(expand.grid(seq(0,1,length=n), seq(0,1,length=n))) X=abs(outer(x[,1],x[,1],"-")) rho=0.3 R=rho^(X^2) X=abs(outer(x[,2],x[,2],"-")) rho=1e-15 R=R*rho^(X^2) L=t(chol(R+diag(n^2)*.Machine$double.eps*100)) mu=0 Z=rnorm(n^2,mean=0,sd=1) Y=L%*%Z+mu # Make a nice plot library(rgl) par3d(cex=0.5) persp3d(matrix(Y,n,n),col="grey",xlab="x1",ylab="x2",zlab="Y",box=FALSE) plot3d(x[,1],x[,2],Y,col="black",type='s',radius=0.01,add=TRUE) # Dense 2D grid X=as.matrix(expand.grid(seq(0,1,length=25), seq(0,1,length=25))) dim(X) plot(X,pch=20,xlab="x1",ylab="x2") # Kronecker Product Covariances set.seed(88) n=25 x1=seq(0,1,length=n) x2=seq(0,1,length=n) X1=abs(outer(x1,x1,"-")) rho=0.3 R1=rho^(X1^2) X2=abs(outer(x2,x2,"-")) rho=1e-15 R2=rho^(X2^2) RR=R2%x%R1 #kronecker product sum(abs(RR-R)) # no difference LL=t(chol(R2+diag(n)*.Machine$double.eps*100))%x% t(chol(R1+diag(n)*.Machine$double.eps*100)) mu=0 Z=rnorm(n^2,mean=0,sd=1) Y=LL%*%Z+mu # Plot par3d(cex=0.5) persp3d(matrix(Y,n,n),col="grey",xlab="x1",ylab="x2", zlab="Y",box=FALSE) plot3d(x[,1],x[,2],Y,col="black",type='s',radius=0.01, add=TRUE)