inst/doc/co2.R

In GPFDA: Gaussian Process for Functional Data Analysis

## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
        echo = TRUE, results = 'hold', warning=F, cache=F, 
  #dev = 'pdf', 
  message=F, 
  fig.width=5, fig.height=5,
  tidy.opts=list(width.cutoff=75), tidy=FALSE
)
old <- options(scipen = 1, digits = 4)

## -----------------------------------------------------------------------------
library(GPFDA)

## -----------------------------------------------------------------------------
co2 <- datasets::co2

## -----------------------------------------------------------------------------
y <- co2[1:60]
n <- length(y)
input <- 5*(1:n)/n

## -----------------------------------------------------------------------------
inputNew <- 5*seq(1/n, 1, length.out=1000) # input test for prediction

## ---- warnings=F, results=F, fig.width=8, fig.height=5------------------------
set.seed(789)
fit1 <- gpr(input=input, response=y, Cov='pow.ex', gamma=1, meanModel='t',
                      nInitCandidates=50, trace=2)

## -----------------------------------------------------------------------------
sapply(fit1$hyper, exp)

## ---- warnings=F, results=F, fig.width=8, fig.height=5------------------------
plot(fit1, main="exponential kernel - fitted model")

## ---- warnings=F, results=F, fig.width=8, fig.height=5------------------------
pred1 <- gprPredict(train = fit1, inputNew=inputNew)
plot(pred1, main="exponential kernel - predictions")

## -----------------------------------------------------------------------------
cov.custom <- function(hyper, input, inputNew=NULL){
  hyper <- lapply(hyper, exp)
  if(is.null(inputNew)){
    inputNew <- as.matrix(input)
  }else{
    inputNew <- as.matrix(inputNew)
  }
  input <- as.matrix(input)
  A1 <- distMat(input=input, inputNew=inputNew, A=as.matrix(hyper$custom.w), 
                        power=2)
  sept <- outer(input[,1], inputNew[,1], "-")
  A2 <- hyper$custom.u*(sin(pi*sept))^2
  customCov <- hyper$custom.v*exp(-A1-A2)
  return(customCov)
}

## -----------------------------------------------------------------------------
Dloglik.custom.w <- function(hyper, input, AlphaQ){
  A1 <- distMatSq(input=input, A=as.matrix(exp(hyper$custom.w)), power=2)
  Dcov <- - cov.custom(hyper, input)*A1
  out <- 0.5*sum(diag(AlphaQ%*%Dcov))
  return(out)
}

Dloglik.custom.u <- function(hyper, input, AlphaQ){
  sept <- outer(input[,1], input[,1], "-")
  A2 <- exp(hyper$custom.u)*(sin(pi*sept))^2
  Dcov <- - cov.custom(hyper, input)*A2
  out <- 0.5*sum(diag(AlphaQ%*%Dcov))
  return(out)
}

Dloglik.custom.v <- function(hyper, input, AlphaQ){
  Dcov <- cov.custom(hyper, input)
  out <- 0.5*sum(diag(AlphaQ%*%Dcov))
  return(out)
}

## -----------------------------------------------------------------------------
D2custom.w <- function(hyper, input, inv.Q, Alpha.Q){
  Cov <- cov.custom(hyper, input)
  A1 <- distMatSq(input=input, A=as.matrix(exp(hyper$custom.w)), power=2)
  D1cov <- -Cov*A1
  D2cov <- -(D1cov + Cov)*A1
  D2c.w <- D2(D1cov, D2cov, inv.Q, Alpha.Q)
  return(D2c.w)
}

D2custom.u <- function(hyper, input, inv.Q, Alpha.Q){
  Cov <- cov.custom(hyper, input)
  sept <- outer(input[,1], input[,1], "-")
  A2 <- exp(hyper$custom.u)*(sin(pi*sept))^2
  D1cov <- - Cov*A2
  D2cov <- -(D1cov + Cov)*A2
  D2c.u <- D2(D1cov, D2cov, inv.Q, Alpha.Q)
  return(D2c.u)
}

D2custom.v <- function(hyper, input, inv.Q, Alpha.Q){
  out <- cov.custom(hyper, input)
  return(out)
}

## -----------------------------------------------------------------------------
diag.custom <- function(hyper, input){
  Qstar <- rep(exp(hyper$custom.v), nrow(input))
  return(Qstar)
}

## ---- message=F, results=F, fig.width=8, fig.height=5-------------------------
fit2 <- gpr(input=input, response=y, Cov='custom',
            NewHyper=c('custom.w','custom.u','custom.v'), gamma=2, meanModel='t',
            nInitCandidates=50, trace=2, useGradient = F)

## -----------------------------------------------------------------------------
sapply(fit2$hyper, exp)

## ---- message=F, fig.width=8, fig.height=5------------------------------------
plot(fit2, main="customised kernel - fitted model")

## ---- message=F, fig.width=8, fig.height=5------------------------------------
pred2 <- gprPredict(train = fit2, inputNew=inputNew)
plot(pred2, main="customised kernel - predictions")

## ---- include = FALSE---------------------------------------------------------
options(old)