Ch.6　Gauss Process and Functional Data Analysis
6.1 Regression
Example 83
## Definition of (m, k)
m=function(x) 0; k=function(x,y) exp(-(x-y)^2/2)
## Definition of function gp_sample
gp.sample=function(x,m,k){
n=length(x)
m.x=m(x)
k.xx=matrix(0,n,n); for(i in 1:n)for(j in 1:n)k.xx[i,j]=k(x[i],x[j])
R=t(chol(k.xx))
u=rnorm(n)
return(as.vector(R%*%u+m.x))
}
## Generate random numbers, generate covariance matrix and compare with k.xx
x=seq(-2,2,1); n=length(x)
r=100; z=matrix(0,r,n); for(i in 1:r)z[i,]=gp.sample(x,m,k)
k.xx=matrix(0,n,n); for(i in 1:n)for(j in 1:n)k.xx[i,j]=k(x[i],x[j])
cov(z)## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.02638283 0.5235439 0.1079838 0.00681492 -0.02331743
## [2,] 0.52354392 0.8813942 0.5689148 0.12409130 0.10024847
## [3,] 0.10798383 0.5689148 1.0639257 0.62118779 0.17264526
## [4,] 0.00681492 0.1240913 0.6211878 0.91195506 0.43095029
## [5,] -0.02331743 0.1002485 0.1726453 0.43095029 0.71963669k.xx## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.0000000000 0.6065307 0.1353353 0.0111090 0.0003354626
## [2,] 0.6065306597 1.0000000 0.6065307 0.1353353 0.0111089965
## [3,] 0.1353352832 0.6065307 1.0000000 0.6065307 0.1353352832
## [4,] 0.0111089965 0.1353353 0.6065307 1.0000000 0.6065306597
## [5,] 0.0003354626 0.0111090 0.1353353 0.6065307 1.0000000000Example 84
## Definition of (m, k)
m=function(x) x[,1]-x[,2]
k=function(x,y) exp(-sum((x-y)^2)/2)
## Definition of function gp_sample
gp.sample=function(x,m,k){
n=nrow(x)
m.x=m(x)
k.xx=matrix(0,n,n); for(i in 1:n)for(j in 1:n)k.xx[i,j]=k(x[i,],x[j,])
R=t(chol(k.xx))
u=rnorm(n)
return(R%*%u+m.x)
}
## Generate random numbers, generate covariance matrix and compare with k.xx
n=5; x=matrix(rnorm(n*2),n,n)
r=100; z=matrix(0,r,n); for(i in 1:r)z[i,]=gp.sample(x,m,k)
k.xx=matrix(0,n,n); for(i in 1:n)for(j in 1:n)k.xx[i,j]=k(x[i],x[j])
cov(z)## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.0393134 0.12070296 0.26763331 0.34890842 0.47179829
## [2,] 0.1207030 0.94633809 -0.07968400 -0.01169706 0.06137032
## [3,] 0.2676333 -0.07968400 1.19382568 0.06649081 0.18494846
## [4,] 0.3489084 -0.01169706 0.06649081 0.75371975 0.17094662
## [5,] 0.4717983 0.06137032 0.18494846 0.17094662 0.93273057k.xx## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.0000000 0.6343544 0.9094558 0.8484902 0.8959100
## [2,] 0.6343544 1.0000000 0.3807027 0.3114987 0.8889401
## [3,] 0.9094558 0.3807027 1.0000000 0.9905846 0.6642492
## [4,] 0.8484902 0.3114987 0.9905846 1.0000000 0.5810156
## [5,] 0.8959100 0.8889401 0.6642492 0.5810156 1.0000000gp.1=function(x.pred){
h=array(dim=n); for(i in 1:n)h[i]=k(x.pred,x[i])
R=solve(K+sigma.2*diag(n)) ## Calculaion of O(n^3)
mm=mu(x.pred)+t(h)%*%R%*%(y-mu(x))
ss=k(x.pred,x.pred)-t(h)%*%R%*%h
return(list(mm=mm,ss=ss))
}
gp.2=function(x.pred){
h=array(dim=n); for(i in 1:n)h[i]=k(x.pred,x[i])
L=chol(K+sigma.2*diag(n)) ## Calculaion of O(n^3/3)
alpha=solve(L,solve(t(L),y-mu(x))) ## Calculaion of O(n^2)
mm=mu(x.pred)+sum(t(h)*alpha)
gamma=solve(t(L),h) ## Calculaion of O(n^2)
ss=k(x.pred,x.pred)-sum(gamma^2)
return(list(mm=mm,ss=ss))
}Example 85
sigma.2=0.2
k=function(x,y)exp(-(x-y)^2/2/sigma.2) # covariance function
mu=function(x) x # average function
n=1000; x=runif(n)*6-3; y=sin(x/2)+rnorm(n) # Data Generation
K=array(dim=c(n,n)); for(i in 1:n)for(j in 1:n)K[i,j]=k(x[i],x[j])
## Measure execution time
library(tictoc)## Warning: package 'tictoc' was built under R version 4.0.5tic(); gp.1(0); toc()## $mm
## [,1]
## [1,] 0.06461642
##
## $ss
## [,1]
## [1,] 0.002643556## 0.85 sec elapsedtic(); gp.2(0); toc()## $mm
## [1] 0.06461642
##
## $ss
## [1] 0.002643556## 0.98 sec elapsed# 3 sigma width before and after the average is noted
u.seq=seq(-3,3,0.1); v.seq=NULL; w.seq=NULL;
for(u in u.seq){res=gp.1(u); v.seq=c(v.seq,res$mm); w.seq=c(w.seq,sqrt(res$ss))}
plot(u.seq,v.seq,xlim=c(-3,3),ylim=c(-3,3),type="l")
lines(u.seq,v.seq+3*w.seq,col="blue"); lines(u.seq,v.seq-3*w.seq,col="blue")
points(x,y)
## 5 times with different samples
plot(0,xlim=c(-3,3),ylim=c(-3,3),type="n")
n=100
for(h in 1:5){
x=runif(n)*6-3; y=sin(pi*x/2)+rnorm(n)
sigma2=0.2
K=array(dim=c(n,n)); for(i in 1:n)for(j in 1:n)K[i,j]=k(x[i],x[j])
u.seq=seq(-3,3,0.1); v.seq=NULL
for(u in u.seq){res=gp.1(u); v.seq=c(v.seq,res$mm)}
lines(u.seq,v.seq,col=h+1)
}
## 6.2 Classification
Example 86
## Iris Data
df=iris
x=df[1:100,1:4]
y=c(rep(1,50),rep(-1,50))
n=length(y)
## Compute kernels with 4 covariates
k=function(x,y) exp(sum(-(x-y)^2)/2)
K=matrix(0,n,n)
for(i in 1:n)for(j in 1:n)K[i,j]=k(x[i,],x[j,])
eps=0.00001
f=rep(0,n)
g=rep(0.1,n)
while(sum((f-g)^2)>eps){
g=f ## Store the value before updating for comparison
v=exp(-y*f)
u=y*v/(1+v)
w=as.vector(v/(1+v)^2)
W=diag(w); W.p=diag(w^0.5); W.m=diag(w^(-0.5))
L=chol(diag(n)+W.p%*%K%*%W.p)
L=t(L) ## R language chol function outputs a transposed matrix
gamma=W%*%f+u
beta=solve(L,W.p%*%K%*%gamma)
alpha=solve(t(L)%*%W.m,beta)
f=K%*%(gamma-alpha)
}
as.vector(f)## [1] 2.901760 2.666188 2.736000 2.596215 2.888826 2.422990 2.712837
## [8] 2.896583 2.263840 2.722794 2.675787 2.804277 2.629917 2.129059
## [15] 1.994737 1.725577 2.502404 2.894768 2.211715 2.757889 2.580703
## [22] 2.788434 2.450147 2.598253 2.493633 2.589272 2.799560 2.854338
## [29] 2.858034 2.682198 2.663180 2.652952 2.409809 2.157029 2.738196
## [36] 2.777507 2.605493 2.848624 2.342636 2.882683 2.887406 1.561917
## [43] 2.454169 2.649399 2.407117 2.633906 2.727124 2.673216 2.749571
## [50] 2.884288 -1.870442 -2.537382 -1.932737 -2.531886 -2.579367 -2.785015
## [57] -2.381783 -1.467521 -2.486519 -2.356974 -1.600772 -2.811324 -2.381371
## [64] -2.734406 -2.330770 -2.341664 -2.614817 -2.748088 -2.372972 -2.628800
## [71] -2.251908 -2.789266 -2.345693 -2.704238 -2.712489 -2.502956 -2.162480
## [78] -2.015710 -2.840261 -2.194749 -2.474090 -2.342643 -2.755051 -2.189084
## [85] -2.415174 -2.382286 -2.275106 -2.496755 -2.695880 -2.664357 -2.648632
## [92] -2.771846 -2.807510 -1.546826 -2.814728 -2.756341 -2.834644 -2.849801
## [99] -1.182284 -2.845863pred=function(z){
kk=array(0,dim=n); for (i in 1:n)kk[i]=k(z,x[i,])
mu=sum(kk*as.vector(u)) ## Average
alpha=solve(L,W.p%*%kk); sigma2=k(z,z)-sum(alpha^2) ## Variance
m=1000; b=rnorm(m,mu,sigma2); pi=sum((1+exp(-b))^(-1))/m ## Predicted value
return(pi)
}Example 87
z=array(0,dim=4)
for(j in 1:4)z[j]=mean(x[1:50,j])
pred(z)## [1] 0.9467988for(j in 1:4)z[j]=mean(x[51:100,j])
pred(z)## [1] 0.053256086.3 Inducing variable method
Eample 88
sigma.2=0.05 #Originally, it should be presumed
k=function(x,y)exp(-(x-y)^2/2/sigma.2) # Covariance function
mu=function(x) x # Average function
n=200; x=runif(n)*6-3; y=sin(x/2)+rnorm(n) # Data Generation
eps=10^(-6)
m=100
index=sample(1:n, m, replace=FALSE)
z=x[index]
m.x=0
m.z=0
K.zz=array(dim=c(m,m)); for(i in 1:m)for(j in 1:m)K.zz[i,j]=k(z[i],z[j])
K.xz=array(dim=c(n,m)); for(i in 1:n)for(j in 1:m)K.xz[i,j]=k(x[i],z[j])
K.zz.inv=solve(K.zz+diag(rep(10^eps,m)))
lambda=array(dim=n)
for(i in 1:n)lambda[i]=k(x[i],x[i])-K.xz[i,1:m]%*%K.zz.inv%*%K.xz[i,1:m]
Lambda.0.inv=diag(1/(lambda+sigma.2))
Q=K.zz+t(K.xz)%*%Lambda.0.inv%*%K.xz ## The computation of Q does not require O(n^3)
Q.inv=solve(Q+diag(rep(eps,m)))
muu=Q.inv%*%t(K.xz)%*%Lambda.0.inv%*%(y-m.x)
dif=K.zz.inv-Q.inv
K=array(dim=c(n,n)); for(i in 1:n)for(j in 1:n)K[i,j]=k(x[i],x[j])
R=solve(K+sigma.2*diag(n)) ## O(n^3) calculations are required
gp.ind=function(x.pred){
h=array(dim=m); for(i in 1:m)h[i]=k(x.pred,z[i])
mm=mu(x.pred)+h%*%muu
ss=k(x.pred,x.pred)-h%*%dif%*%h
return(list(mm=mm,ss=ss))
} ## Using the inducing variable method
gp.1=function(x.pred){
h=array(dim=n); for(i in 1:n)h[i]=k(x.pred,x[i])
mm=mu(x.pred)+t(h)%*%R%*%(y-mu(x))
ss=k(x.pred,x.pred)-t(h)%*%R%*%h
return(list(mm=mm,ss=ss))
} ## Not using the inducing variable method
x.seq=seq(-2,2,0.1)
mmv=NULL; ssv=NULL
for(u in x.seq){
mmv=c(mmv,gp.ind(u)$mm)
ssv=c(ssv,gp.ind(u)$ss)
}
plot(0, xlim=c(-2,2),ylim=c(min(mmv),max(mmv)),type="n")
lines(x.seq,mmv,col="red")
lines(x.seq,mmv+3*sqrt(ssv),lty=3,col="red")
lines(x.seq,mmv-3*sqrt(ssv),lty=3,col="red")
x.seq=seq(-2,2,0.1)
mmv=NULL; ssv=NULL
for(u in x.seq){
mmv=c(mmv,gp.1(u)$mm)
ssv=c(ssv,gp.1(u)$ss)
}
lines(x.seq,mmv,col="blue")
lines(x.seq,mmv+3*sqrt(ssv),lty=3,col="blue")
lines(x.seq,mmv-3*sqrt(ssv),lty=3,col="blue")
points(x,y)
6.4 Karhunen-Loeve Expansion
Example 89
lambda=function(j) 4/((2*j-1)*pi)^2 ## eigenvalue
ee=function(j,x) sqrt(2)*sin((2*j-1)*pi/2*x) ## Definition of Eigenfunction
n=10; m=7
f=function(z,x){ ## Definition of Gauss Process
n=length(z)
S=0; for(i in 1:n)S=S+z[i]*ee(i,x)*sqrt(lambda(i))
return(S)
}
plot(0,xlim=c(-3,3),ylim=c(-2,2),type="n",xlab="x",ylab="f(omega,x)")
for(j in 1:m){
z=rnorm(n)
x.seq=seq(-3,3,0.001)
y.seq=NULL; for(x in x.seq)y.seq=c(y.seq,f(z,x))
lines(x.seq,y.seq,col=j)
}
title("Brown Motion")
matern=function(nu,l,r){
p=nu-1/2
S=0
for(i in 0:p)S=S+gamma(p+i+1)/gamma(i+1)/gamma(p-i+1)*(sqrt(8*nu)*r/l)^(p-i)
S=S*gamma(p+2)/gamma(2*p+1)*exp(-sqrt(2*nu)*r/l)
return(S)
}Example 90
m=10
l=0.1
for(i in 1:1)curve(matern(i-1/2,l,x),0,0.5,ylim=c(0,10),col=i+1)
for(i in 2:m)curve(matern(i-1/2,l,x),0,0.5,ylim=c(0,10),ann=FALSE,add=TRUE,col=i+1)
legend("topright",legend=paste("nu=",(1:m)+0.5),lwd=2,col=1:m,)
title("Matern Kernel (l=1)")
Example 91
rand.100=function(Sigma){
L=t(chol(Sigma)) ## Cholesky decomposition of the covariance matrix
u=rnorm(100)
y=as.vector(L%*%u) ## Generates a pair of random numbers of covariance matrix Sigma with mean 0
}
x = seq(0,1,length=100)
z = abs(outer(x,x,"-")) # compute distance matrix, d_{ij} = |x_i - x_j|
l=0.1
Sigma_OU = exp(-z/l) ## OU: Slow in matern(1/2, l, z)
y=rand.100(Sigma_OU)
plot(x,y,type="l",ylim=c(-3,3))
for(i in 1:5){
y = rand.100(Sigma_OU)
lines(x,y,col=i+1)
}
title("OU process (nu=1/2,l=0.1)")
Sigma_M=matern(3/2,l,z) ## Matern
y = rand.100(Sigma_M)
plot(x,y,type="l",ylim=c(-3,3))
for(i in 1:5){
y = rand.100(Sigma_M)
lines(x,y,col=i+1)
}
title("Matern process (nu=3/2,l=0.1)")
6.5 Functional Data Analysis
Exampele 92
library(fda)## Loading required package: splines## Loading required package: fds## Warning: package 'fds' was built under R version 4.0.5## Loading required package: rainbow## Warning: package 'rainbow' was built under R version 4.0.5## Loading required package: MASS## Warning: package 'MASS' was built under R version 4.0.5## Loading required package: pcaPP## Warning: package 'pcaPP' was built under R version 4.0.5## Loading required package: RCurl## Warning: package 'RCurl' was built under R version 4.0.5## Loading required package: deSolve## Warning: package 'deSolve' was built under R version 4.0.5##
## Attaching package: 'fda'## The following object is masked from 'package:graphics':
##
## matplotg=function(j,x){ ## Prepare p bases
if(j==1) return(1/sqrt(2*pi))
if(j%%2==0) return(cos((j%/%2)*x)/sqrt(pi))
else return(sin((j%/%2)*x)/sqrt(pi))
}
beta=function(x,y){ ## Compute the coefficients before the p bases of the function
X=matrix(0,N,p)
for(i in 1:N)for(j in 1:p)X[i,j]=g(j,x[i])
beta=solve(t(X)%*%X+0.0001*diag(p))%*%t(X)%*%y
return(drop(beta))
}
N=365; n=35; m=5; p=100; df=daily
C=matrix(0,n,p)
for(i in 1:n){x=(1:N)*(2*pi/N)-pi; y=as.vector(df[[2]][,i]); C[i,]=beta(x,y)}
res=prcomp(C)
B=res$rotation
xx=res$x
z=function(i,m,x){ ## The original function approximated by m principal components out of p bases
S=0
for(j in 1:p)for(k in 1:m)for(r in 1:p)S=S+C[i,j]*B[j,k]*B[r,k]*g(r,x)
return(S)
}
x.seq=seq(-pi,pi,2*pi/100)
plot(0,xlim=c(-pi,pi),ylim=c(-15,25),type="n",xlab="Days",ylab="Temp(C)",
main="Reconstruction for each m")
lines(x,df[[2]][,14],lwd=2)
for(m in 2:6){
lines(x.seq,z(14,m,x.seq),col=m,lty=1)
}
legend("bottom",legend=c("Original",paste("m=",2:6)), lwd=c(2,rep(1,5)), col=1:6,ncol=2)
lambda=res$sdev^2
ratio=lambda/sum(lambda)
plot(1:5,ratio[1:5],xlab="PC1 through PC5",ylab="Ratio",type="l",main="Ratio ")
h=function(coef,x){ ## Define functions with coefficients
S=0
for(j in 1:p)S=S+coef[j]*g(j,x)
return(S)
}
plot(0,xlim=c(-pi,pi),ylim=c(-1,1),type="n")
for(j in 1:3)lines(x.seq,h(B[,j],x.seq),col=j)
index=c(10,12,13,14,17,24,26,27)
others=setdiff(1:35,index)
first=substring(df[[1]][index],1,1)
plot(0,xlim=c(-25,35),ylim=c(-15,10),type="n",xlab="PC1",ylab="PC2",main="Canadian Weather")
points(xx[others,1],xx[others,2],pch=4)
points(xx[index,1], xx[index,2], pch = first, cex = 1, col =rainbow(8))
legend("bottom",legend=df[[1]][index], pch=first, col=rainbow(8),ncol=2)