This page is designed to explain `?outermod`

and
`?outerbase`

.

`library(outerbase)`

A three dimensional case study is sufficient for understanding what
`outermod`

and `outerbase`

are and how to
manipulate them.

```
= 30
sampsize = 3
d = seq(1/(2*sampsize),1-1/(2*sampsize),1/sampsize)
design1d = cbind(design1d,sample(design1d),sample(design1d))
x = obtest_borehole3d(x) y
```

Covariance functions are an important building block of Gaussian
process inference. This package uses a custom class to represent
covariance function. See `?covf`

for more information on the
base class.

Creating instances is done through the `?methods::new`

call with the class name listed inside.

`= new(covf_mat25) corf `

The `cov`

method builds covariances matrices. Below is an
example of calling this method. Note that these are designed to be
single dimensional covariance functions.

```
= x[1:5,1]
xred print(corf$cov(xred,xred),3)
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 1.000 1.000 0.999 0.998 0.997
#> [2,] 1.000 1.000 1.000 0.999 0.998
#> [3,] 0.999 1.000 1.000 1.000 0.999
#> [4,] 0.998 0.999 1.000 1.000 1.000
#> [5,] 0.997 0.998 0.999 1.000 1.000
```

Hyperparameters are important for almost all covariance functions.
They control the general shape and behavior of the covariance function.
They are stored in the `covf`

class in the (editable) field
`hyp`

.

```
$hyp
corf#> [,1]
#> [1,] 0
```

You can see the effect of alternating `hyp`

below on this
correlation function.

```
$hyp = c(-0.5)
corfplot(x[,1],corf$cov(x[,1],0.5), type='l',
ylab='correlation with 0.5', xlab='input')
$hyp = c(-0.25)
corflines(x[,1],corf$cov(x[,1],0.5), type='l')
$hyp = c(0)
corflines(x[,1],corf$cov(x[,1],0.5), type='l')
```

Gaussian processes have long been shown to be top performers for near interpolation. For more information on general Gaussian processes, see the textbooks Gaussian Processes for Machine Learning or Surrogates, among others.

The idea is to represent a surface as a realization of a Gaussian process controlled by a covariance function. An outer product of covariance functions can do this job in three dimensions. This means we need to first build covariance functions.

```
= new(covf_mat25)
corf1 = new(covf_mat25)
corf2 = new(covf_mat25)
corf3 $hyp = c(-0.5) # just setting them all to the same
corf1$hyp = c(-0.5) # hyperparameter for now
corf2$hyp = c(-0.5) corf3
```

And then multiply them to calculate the covariance between two sets of points.

```
= function(x1,x2){
covftot $cov(x1[,1],x2[,1])*
corf1$cov(x1[,2],x2[,2])*
corf2$cov(x1[,3],x2[,3])
corf3
}= covftot(x,x) #total correlation matrix cormattot
```

The goal of Gaussian process inference is to take our data and
predict at some number of points, say `1000`

points.

```
= 1000
testsampsize = matrix(runif(testsampsize*d),ncol=d) xtest
```

The predictor follows from typical formulas (assuming y is zero mean, see textbooks).

`= covftot(xtest,x) %*% solve(cormattot,y) yhat `

This gives prediction accuracy that can be summarized below.

```
= obtest_borehole3d(xtest)
ytest plot(yhat, ytest, ylab="actual", xlab="prediction")
hist(ytest-yhat, main="test residuals",
xlab = "test residuals")
```

We can also use this framework to get predictive variances. These equations will not be explained in this documentation for brevity. One point here is that is does pretty well! The plot below looks standard Normal enough.

```
= as.double(t(y)%*% solve(cormattot,y)/length(y))
sigma2hat
= sigma2hat*(covftot(xtest,xtest)-t(covftot(x,xtest))%*%
varpred solve(cormattot,covftot(x,xtest)))
hist((ytest-yhat)/sqrt(diag(varpred)),
main="standarized test residuals",
xlab = "standarized test residuals")
```

The main complaints about Gaussian process inference are stability and computation speed. This package is designed to reduce those concerns.

The core classes in this package are `?outermod`

and
`?outerbase`

. An `outermod`

instance contains all
the information to build an `outerbase`

instance, but does
not build the objects corresponding to a specific `x`

. An
`outerbase`

instance is used build inference at a specific
`x`

.

An instance of the class `outermod`

is designed to hold
the information needed to create a basis matrix. An
`outermod`

instance is created using `new`

command.

`= new(outermod) om `

The first step is to set the covariance functions and the knots. To
set the vector of `covfs`

, use `?setcovfs`

alongside a vector strings of covariance functions in the package
(`?listcov`

).

`setcovfs(om, rep("mat25",3))`

This fixes the dimension of the `outermod`

instance
`om`

to `3`

.

Then we need to give it a set of knot points for each dimension. The
choice of these is still being researched, but choosing points that are
spread out in each dimension that look like our actual data is currently
recommended. You will need to invert a matrix of the size of these knot
points, so it is recommended to keep it small, `<50`

in
general. The function `?setknot`

should be used.

```
= list(seq(0,1,by=0.025),
knotlist seq(0,1,by=0.025),
seq(0,1,by=0.025))
setknot(om, knotlist)
```

The hyperparameters can be set directly through our
`outermod`

object.

```
gethyp(om)
#> inpt1.scale inpt2.scale inpt3.scale
#> 0 0 0
$updatehyp(c(-0.5,-0.5,-0.5))
omgethyp(om)
#> inpt1.scale inpt2.scale inpt3.scale
#> -0.5 -0.5 -0.5
```

An instance of the class `outerbase`

is the equivalent of
a basis matrix with fast computation methods included. It is also
created with `new`

, but it also requires a reference to an
`outermod`

instance and a specific set of prediction points
`x`

.

```
= new(outerbase,
ob # an outermod (reference only)
om, # an input matrix x)
```

This builds a set of basis functions for each dimension, which
sometimes just look like polynomials. They are not quite polynomials,
and some covariance functions give different shapes. The call
`outerbase$getbase`

will allow you to access the basis
functions for each dimension. This is mostly useful for plotting.

```
= ob$getbase(1)
basis_func matplot(x[,1],basis_func[,1:4],
type='l', ylab="func", xlab="first dim")
```

`outermod`

and `outerbase`

are meant to be used
in conjunction with each other. One key ingredient is the
`outermod$selectterms`

function, which allows you to pick
products of basis functions that best represent the current
`outermod`

response.

```
= 60
p = om$selectterms(p) # 60 by 3 matrix
terms head(terms)
#> [,1] [,2] [,3]
#> [1,] 0 0 0
#> [2,] 0 0 1
#> [3,] 0 1 0
#> [4,] 1 0 0
#> [5,] 1 0 1
#> [6,] 0 1 1
```

`outermod$getvar`

returns the vector of variances
associated with these coefficients at these `terms`

.

`= as.vector(om$getvar(terms)) covcoeff `

The specific basis matrix can be formed by getting the basis matrix
at these selected `terms`

. `outerbase$getmat`

will
give a short cut to building this matrix.

```
= ob$getmat(terms)
basismat
= 5
termno = ob$getbase(1)[,terms[termno,1]+1]*
basevec $getbase(2)[,terms[termno,2]+1]*
ob$getbase(3)[,terms[termno,3]+1]
ob
cbind(basevec[1:5],basismat[1:5,5]) # expect equal
#> [,1] [,2]
#> [1,] 1.68450167 1.68450167
#> [2,] -0.30349318 -0.30349318
#> [3,] -0.08489499 -0.08489499
#> [4,] 0.93906322 0.93906322
#> [5,] -0.89735638 -0.89735638
```

This package leverages the insight that Gaussian processes are no
more than linear combinations of basis functions with random
coefficients. This viewpoint is often called the *feature space*
view of Gaussian processes.

To see this, not that if you take the `covcoeff`

and
`basismat`

together, the correlation function is very well
approximated through the following manipulation.

```
= basismat%*%diag(covcoeff)%*%t(basismat)
cormatob
print(round(cormattot[1:5,1:5],3)) # typical gp
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 1.000 0.848 0.515 0.690 0.553
#> [2,] 0.848 1.000 0.661 0.736 0.786
#> [3,] 0.515 0.661 1.000 0.921 0.917
#> [4,] 0.690 0.736 0.921 1.000 0.857
#> [5,] 0.553 0.786 0.917 0.857 1.000
print(round(cormatob[1:5,1:5],3)) # outerbase
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 1.000 0.848 0.515 0.690 0.553
#> [2,] 0.848 1.000 0.661 0.736 0.786
#> [3,] 0.515 0.661 1.000 0.921 0.917
#> [4,] 0.690 0.736 0.921 1.000 0.857
#> [5,] 0.553 0.786 0.917 0.857 1.000
```

They means that we can leverage Bayesian
linear regression to do prediction. This will require assuming that
there is some `noisevar`

, which is also called the nugget in
the Gaussian process literature.

```
= 10^(-4)
noisevar #posterior precision matrix of coefficients
= solve(1/noisevar*t(basismat)%*%basismat+
postcov 1/sigma2hat*diag(1/covcoeff))
#posterior mean of coefficients
= postcov%*%(1/noisevar*t(basismat)%*%y) coeffest
```

Consider predicting at some new `xtest`

to examine if the
inference works the same between the traditional Gaussian process and
the feature space approximation.

```
= new(outerbase,
obtest # same outermod
om, # new input matrix
xtest)
= obtest$getmat(terms) basistest
```

The predictions are nearly equivalent.

```
= basistest%*%coeffest
yhatob plot(yhat, ytest, main="typical gp",
xlab="prediction", ylab="actual")
plot(yhatob, ytest, main = "outerbase equiv.",
xlab="prediction", ylab="actual")
```

The histograms of residuals show similar matching.

```
hist(ytest-yhat, main="typical gp",
xlab="test residuals")
hist(ytest-yhatob, main="outerbase equiv.",
xlab="test residuals")
```

The standardized residuals, which account for the variance, also show similar matching.

```
= basistest%*%postcov%*%t(basistest)
varpredob hist((ytest-yhat)/sqrt(diag(varpred)), main="typical gp",
xlab="standarized test residuals")
hist((ytest-yhatob)/sqrt(diag(varpredob)), main="outerbase equiv.",
xlab="standarized test residuals")
```