# Damped-sin example: toy function in 1-dimensional x

dampsin <- function(x, a = 10, b = 0.9, c = 1.5) {
  y <- a * sin(4 * pi * x^b) * exp(-c * x)
  return(y)
}

# Could generate the training data using the above function,
# but we will usually read the data like this.
x <- read.table("dampsin_x.txt", header = TRUE)
dim(x)
head(x, n = 3)
y <- read.table("dampsin_y.txt", header = TRUE)
dim(y)
head(y, n = 3)

# Need to install GaSP first
library(GaSP)

# Constant mean, power-exponential correlation function
const_powexp <- Fit(x = x, y = y, reg_model = ~ 1, 
                    cor_family = "PowerExponential", random_error = FALSE)
# Inspect parameter estimates
const_powexp$beta
const_powexp$cor_par
const_powexp$sp_var

# Plot true function (don't know it in practice!) and training data
x_pred <- as.matrix(seq(0, 1, by = 0.01), ncol = 1)
head(x_pred, n = 3)
y_true <- dampsin(x_pred)
plot(x_pred[, 1], y_true, xlab = "x", ylab = "y", type = "l", lty = 1)
points(x$x, y$y, pch = 20, cex = 2.5)

# Predict at the fine grid and add predictions to plot
dampsin_pred <- Predict(const_powexp, x_pred, generate_coefficients = FALSE)
# Doesn't run because GaSP does a lot of error checking.  Fix the error!
dampsin_pred <- Predict(const_powexp, x_pred, generate_coefficients = FALSE)

# Add predictions to plot
lines(x_pred, dampsin_pred$y_pred$Pred, lty = 2, col = "red")
# Add pointwise approximate 95% confidence intervals
lines(x_pred, dampsin_pred$y_pred$Pred - 1.96 * dampsin_pred$y_pred$SE, lty = 3,
      col = "red")
lines(x_pred, dampsin_pred$y_pred$Pred + 1.96 * dampsin_pred$y_pred$SE, lty = 3,
      col = "red")
legend(0.65, 8, legend = c("True function", "Prediction", "Approx. 95% CI"), 
       lty = c(1, 2, 3), col = c("Black", "red", "red"))