Movatterモバイル変換

[0]ホーム
URL:

NSGPR example

library(GPFDA)require(MASS)

Simulating data from a GP with nonstationary and nonseparablecovariance structure

In this example, we simulate data from a GP whose covariancestructure has a spatially-varying anisotropy matrix.

We simulate\(10\) independentrealizations of\(X(\boldsymbol{{t}}), \\boldsymbol{{t}} = (s_1, s_2)^\top\), from the following model:\[\begin{equation}X(\boldsymbol{{t}}) = f(\boldsymbol{{t}}) +\varepsilon(\boldsymbol{{t}}), \quad \varepsilon(\boldsymbol{{t}}) \simN(0, \sigma_\varepsilon^2),\quad \text{ where } \boldsymbol{{t}} \in \mathcal{T} \subset\left[0,1\right]^2.\end{equation}\] We assume that\(f\) is a zero-mean GP with a nonseparableand nonstationary covariance function given by\[\begin{align}\label{NScovfun}\text{Cov}\big[ f(\boldsymbol{{t}}),f(\boldsymbol{{t}}') \big] =\sigma(\boldsymbol{{t}}) \sigma(\boldsymbol{{t}}') |\boldsymbol{{A}}(\boldsymbol{{t}}) | ^{-1/4} |\boldsymbol{{A}}(\boldsymbol{{t}}') | ^{-1/4} \times\bigg| \frac{\boldsymbol{{A}}^{-1}(\boldsymbol{{t}}) +\boldsymbol{{A}}^{-1}(\boldsymbol{{t}}')}{2} \bigg| ^{-1/2} g\Big(\sqrt{Q_{\boldsymbol{{t}}\boldsymbol{{t}}'}}\Big),\end{align}\] where\(g(\cdot)\)is the exponential correlation function,\(\sigma(\boldsymbol{{t}})=1\) and\(\sigma_\varepsilon^2\) is negligible.

The elements of the varying anisotropy matrix\(\boldsymbol{{A}}(s_1,s_2)\) are:\[\begin{align}& \boldsymbol{{A}}_{11}(s_1,s_2) = \exp\big(6 \cos(10 s_1 - 5s_2)\big), \\& \boldsymbol{{A}}_{22}(s_1,s_2) = \exp\big(\sin(6 s_1^3) + \cos(6s_2^4)\big),\\& \boldsymbol{{A}}_{12}(s_1,s_2) =\{\boldsymbol{{A}}_{11}(s_1,s_2)\boldsymbol{{A}}_{22}(s_1,s_2)\}^{1/2}\rho_{12}(s_1,s_2),\\& \rho_{12}(s_1,s_2) = \tanh\big((s_1^2+s_2^2)/2\big).\end{align}\]

set.seed(123)n1<-30# sample size for input coordinate 1n2<-30# sample size for input coordinate 2n<- n1*n2# total sample size# Creating evenly spaced spatial coordinatesinput<-list()input[[1]]<-seq(0,1,length.out = n1)input[[2]]<-seq(0,1,length.out = n2)inputMat<-as.matrix(expand.grid(input))# inputs in matrix form# Creating the varying anisotropy matrix with spatially-varying parametersA11<-function(s1,s2){exp(6*cos(10*s1-5*s2))}A22<-function(s1,s2){exp(sin(6*s1^3)+cos(6*s2^4))}A12<-function(s1,s2){sqrt(A11(s1,s2)*A22(s1,s2))*tanh((s1^2+s2^2)/2)}A_List<-list()R12_vec<-rep(NA, n)for(iin1:n){  s1<- inputMat[i,1]  s2<- inputMat[i,2]  A_i11<-A11(s1=s1,s2=s2)  A_i22<-A22(s1=s1,s2=s2)  A_i12<-A12(s1=s1,s2=s2)  A_i<-matrix(NA,2,2)  A_i[1,1]<- A_i11  A_i[2,2]<- A_i22  A_i[1,2]<- A_i12  A_i[2,1]<- A_i12  A_List[[i]]<- A_i  R12_vec[i]<- A_i12/sqrt(A_i11*A_i22)}# Constructing the (n x n) covariance matrix KScaleDistMats<-calcScaleDistMats(A_List=A_List,coords=inputMat)Scale.mat<- ScaleDistMats$Scale.matDist.mat<- ScaleDistMats$Dist.matcorrModel<-"pow.ex"gamma<-1K<- Scale.mat*unscaledCorr(Dist.mat=Dist.mat,corrModel=corrModel,gamma=gamma)diag(K)<-diag(K)+1e-8# Generate response surfacesmeanFunction<-rep(0, n)nrep<-10response<-t(mvtnorm::rmvnorm(nrep, meanFunction, K))

Estimating NSGP covariance function

The estimation of the NSGP covariance function is done by using thensgpr function.

We want a nonseparable covariance structure (sepCov=F),assuming the process has unit variance(unitSignalVariance=T) and a zero noise variance(zeroNoiseVariance=T).

The unconstrained parameters associated to the varying anisotropymatrix will be modeled via\(L=M=6\)B-spline basis functions with evenly spaced knots. To allow theparameters of the anisotropy matrix to change along both directions\(s_1\) and\(s_2\), we setwhichTau = c(T,T).

The estimated hyperparameters (B-spline coefficients) from thensgpr function are saved in the objecthp.

### NOT RUN# fit <- nsgpr(response = response,#               input = input,#               corrModel = corrModel,#               gamma = gamma,#               whichTau = c(T,T),#               absBounds = 8,#               nBasis = 6,#               nInitCandidates = 1000,#               cyclic = c(F,F),#               unitSignalVariance = T,#               zeroNoiseVariance = T,#               sepCov = F)## Taking ML estimates of B-spline coefficients# hp <- fit$MLEsts###  end NOT RUNhp<-c(5.60699,1.14865,-6.61148,8,-4.44439,-2.14159,3.98823,5.42213,-8,6.91389,-0.17032,-6.29215,0.48661,8,-2.84537,-1.59173,8,-2.55939,-8,-0.95508,7.58644,-8,4.86161,8,-4.08238,-8,8,-2.74033,-3.5963,8,-1.30662,-8,2.05985,4.27012,-7.26238,4.76814,0.67189,0.62838,0.42846,1.33653,-0.40139,0.43309,-0.74954,0.79752,-1.159,3.31145,-0.87894,-0.99702,1.72585,-0.29033,2.41099,0.46065,-0.95101,2.0856,-1.39188,0.57495,-0.39334,3.24682,1.09657,-1.90269,1.5223,-2.87977,1.54015,-3.89324,-1.67514,-0.18049,3.93999,-1.50017,1.18639,2.16521,-6.93837,-1.88978,-2.06261,-0.80103,3.17891,-3.68428,1.64603,0.58847,0.86276,-2.38605,3.99548,-0.52069,1.33181,-0.10872,-0.43153,-4.91787,2.56248,1.90786,-5.382,0.42587,1.43742,0.54047,-0.23679,2.63721,-0.11159,0.57184,-2.06124,1.82977,-6.13951,4.22915,-0.29822,-2.69144,-7.61238,5.1746,-5.04487,7.5953,0.16564,-1.16696,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-23.02585)

Prediction

Based on the estimated NSGP model, we want to make predictions attest input points:

# Creating test input pointsn1test<- n1*2n2test<- n2*2inputNew<-list()inputNew[[1]]<-seq(0,1,length.out = n1test)inputNew[[2]]<-seq(0,1,length.out = n2test)inputMatTest<-as.matrix(expand.grid(inputNew[[1]], inputNew[[2]]))pred<-nsgprPredict(hp=hp,response=response,input=input,inputNew=inputNew,noiseFreePred=F,nBasis=6,corrModel=corrModel,gamma=gamma,cyclic=c(F,F),whichTau=c(T,T))

Below we plot the first realisation and then the correspondingpredictions for the new test input locations:

zlim<-range(c(pred$pred.mean, response))plotImage(response = response,input = inputMat,realisation =1,n1 = n1,n2 = n2,zlim = zlim,main ="observed")

plotImage(response = pred$pred.mean,input = inputMatTest,realisation =1,n1 = n1test,n2 = n2test,zlim = zlim,main ="prediction")

Plots of covariance functions

To calculate the covariance matrix based on thehpestimates found above:

FittedCovMat<-nsgpCovMat(hp=hp,input=input,corrModel=corrModel,gamma=gamma,nBasis=6,cyclic=c(F,F),whichTau=c(T,T),calcCov=T)$Cov

We will plot slices of the true and fitted covariance functionscentred at specific spatial locations (the\(8\)th observed point of\(s_1\) and the\(10\)th of\(s_2\)). The idea is to obtain image plotsof\(k(s_1-0.24, s_2-0.31; 0.24,0.31)\).

# centre points for the covariance functionsinput1cent<- input[[1]][8]input2cent<- input[[2]][10]# half-widthmaxdist<-0.2zlim<-range(c(K, FittedCovMat))+0.2*c(-1,1)centre_idx<-which(inputMat[,1]==input1cent& inputMat[,2]==input2cent)other_idx<-which(inputMat[,1]>input1cent-maxdist&                     inputMat[,1]<input1cent+maxdist&                     inputMat[,2]>input2cent-maxdist&                     inputMat[,2]<input2cent+maxdist)other_idx_le<-length(other_idx)# Locations of the slice:sliceInputMat<- inputMat[other_idx,]-cbind(rep(input1cent, other_idx_le),rep(input2cent, other_idx_le))# Slice of the true covariance matrixsliceCov<- K[centre_idx, other_idx]nEach<-sqrt(length(sliceCov))sliceTrueCov<-c(t(matrix(sliceCov,byrow=T,ncol=nEach)))plotImage(response = sliceTrueCov,input = sliceInputMat,n1 = nEach,n2 = nEach,zlim = zlim,main ="true covariance function")

# Slice of the estimated covariance matrixsliceCov<- FittedCovMat[centre_idx, other_idx]nEach<-sqrt(length(sliceCov))sliceFittedCov<-c(t(matrix(sliceCov,byrow=T,ncol=nEach)))plotImage(response = sliceFittedCov,input = sliceInputMat,n1 = nEach,n2 = nEach,zlim = zlim,main ="estimated covariance function")

[8]ページ先頭
©2009-2025 Movatter.jp