### 7.3. Cautionary note on convergence for covariance parameter estimation

The first fitted GP model of Section 7.2 uses the default number 100 of iteration limit. Details on its optimization process can be obtained by removing both occurrences of `trace = FALSE`.

```
R> m2 <- fgpm(sIn = Xs train, fIn = Xf train, sOut = Y train)
```

** Presampling...

** Optimising hyperparameters...

```
iter   10 value 1178.686098
iter   20 value 1169.785515
iter   30 value 1169.358788
iter   40 value 1168.584676
iter   50 value 1168.117087
iter   60 value 1167.926415
iter   70 value 1167.854377
iter   80 value 1167.782569
iter   90 value 1167.236280
iter  100 value 1167.049624
final  value 1167.048204
stopped after 101 iterations
```

The function value is the negated log-likelihood

** Hyperparameters done!

We can also display a convergence code, which is directly provided by the `optim()` function.

```
R> m2@convergence
```

[1] 1

Here the code 1 means that the iteration limit has been reached, which was also pointed out when `m2` was created with the sentence `stopped after 101 iterations`.

In this case, even though the prediction performance of `m1` (created with the same iteration limit as `m2`) was satisfying in Section 7.2, the optimization process did not completely converge. To see this, we can refit the GP model with a larger iteration limit.

```
R> m3 <- fgpm(sIn = Xs train, fIn = Xf train, sOut = Y train,
+ trace = TRUE, control.optim = list(trace = TRUE, maxit = 500))
```

** Presampling...

** Optimising hyperparameters...
### 7.3. Cautionary note on convergence for covariance parameter estimation

The first fitted GP model of Section 7.2 uses the default number 100 of iteration limit. Details on its optimization process can be obtained by removing both occurrences of `trace = FALSE`.

```
R> m2 <- fgpm(sIn = Xs train, fIn = Xf train, sOut = Y train)
```

** Presampling...

** Optimising hyperparameters...

```
iter   10 value 1178.686098
iter   20 value 1169.785515
iter   30 value 1169.358788
iter   40 value 1168.584676
iter   50 value 1168.117087
iter   60 value 1167.926415
iter   70 value 1167.854377
iter   80 value 1167.782569
iter   90 value 1167.236280
iter  100 value 1167.049624
final  value 1167.048204
stopped after 101 iterations
```

The function value is the negated log-likelihood

** Hyperparameters done!

We can also display a convergence code, which is directly provided by the `optim()` function.

```
R> m2@convergence
```

[1] 1

Here the code 1 means that the iteration limit has been reached, which was also pointed out when `m2` was created with the sentence `stopped after 101 iterations`.

In this case, even though the prediction performance of `m1` (created with the same iteration limit as `m2`) was satisfying in Section 7.2, the optimization process did not completely converge. To see this, we can refit the GP model with a larger iteration limit.

```
R> m3 <- fgpm(sIn = Xs train, fIn = Xf train, sOut = Y train,
+ trace = TRUE, control.optim = list(trace = TRUE, maxit = 500))
```

** Presampling...

** Optimising hyperparameters...
iter 10 value 1174.852660
iter 20 value 1173.392364
iter 30 value 1172.160197
iter 40 value 1171.341091
iter 50 value 1170.522756
iter 60 value 1169.983111
iter 70 value 1169.712522
iter 80 value 1167.982626
iter 90 value 1167.150706
iter 100 value 1167.089213
iter 110 value 1167.085926
iter 120 value 1167.083327
iter 130 value 1167.077044
iter 140 value 1166.886057
iter 150 value 1166.884338
iter 160 value 1166.883324
iter 170 value 1166.883100
iter 180 value 1166.882840
iter 190 value 1166.862427
iter 200 value 1166.855279
iter 210 value 1166.849313
iter 220 value 1166.846073
final value 1166.845922
converged

The function value is the negated log-likelihood

** Hyperparameters done!

Now it is displayed that the optimization process did converge. Furthermore, the final negated log-likelihood is smaller (better) than before. The convergence can also be confirmed by displaying the convergence code.

R> m3@convergence

[1] 0

Here the value 0 is obtained, indicating successful completion of the call to the `optim()` function. Hence, the default iteration limit may be too small to reach convergence, in complex settings. This can be monitored simply from the display given when calling `fgpm`, or from the `convergence` slot of a `fgpm` object. The user has the possibility to increase this iteration limit, by using the `control.optim` argument.

We conclude by illustrating the benefit of the `extend` argument of the `update` function, discussed in Section 3.4. First, let us re-estimate the hyperparameters of the `fgpm` model `m1` of Section 7.2, without benefiting from `extend = TRUE`, that is, by starting the likelihood optimization from scratch.
iter 10 value 1174.852660
iter 20 value 1173.392364
iter 30 value 1172.160197
iter 40 value 1171.341091
iter 50 value 1170.522756
iter 60 value 1169.983111
iter 70 value 1169.712522
iter 80 value 1167.982626
iter 90 value 1167.150706
iter 100 value 1167.089213
iter 110 value 1167.085926
iter 120 value 1167.083327
iter 130 value 1167.077044
iter 140 value 1166.886057
iter 150 value 1166.884338
iter 160 value 1166.883324
iter 170 value 1166.883100
iter 180 value 1166.882840
iter 190 value 1166.862427
iter 200 value 1166.855279
iter 210 value 1166.849313
iter 220 value 1166.846073
final value 1166.845922
converged

The function value is the negated log-likelihood

** Hyperparameters done!

Now it is displayed that the optimization process did converge. Furthermore, the final negated log-likelihood is smaller (better) than before. The convergence can also be confirmed by displaying the convergence code.

R> m3@convergence

[1] 0

Here the value 0 is obtained, indicating successful completion of the call to the `optim()` function. Hence, the default iteration limit may be too small to reach convergence, in complex settings. This can be monitored simply from the display given when calling `fgpm`, or from the `convergence` slot of a `fgpm` object. The user has the possibility to increase this iteration limit, by using the `control.optim` argument.

We conclude by illustrating the benefit of the `extend` argument of the `update` function, discussed in Section 3.4. First, let us re-estimate the hyperparameters of the `fgpm` model `m1` of Section 7.2, without benefiting from `extend = TRUE`, that is, by starting the likelihood optimization from scratch.
R> m4 <- update(m1, var.re = TRUE, ls_s.re = TRUE, ls_f.re = TRUE,
+ extend = FALSE, trace = FALSE)

iter 10 value 1192.597779
iter 20 value 1184.610629
iter 30 value 1176.694606
iter 40 value 1176.082474
iter 50 value 1174.242786
iter 60 value 1169.484368
iter 70 value 1168.365643
iter 80 value 1167.828297
iter 90 value 1167.824714
iter 100 value 1167.623696
final value 1167.621387
stopped after 101 iterations

The function value is the negated log-likelihood

As for the creation of m2 above, since we use the default number 100 of iteration limit, the
iteration limit is reached before convergence. Next, using extend = TRUE, we can restart the
likelihood optimization from the non-converged set of hyperparameter estimates of m1.

R> m5 <- update(m1, var.re = TRUE, ls_s.re = TRUE, ls_f.re = TRUE,
+ extend = TRUE, trace = FALSE)

iter 10 value 1166.974505
iter 20 value 1166.964709
iter 30 value 1166.961533
iter 40 value 1166.913197
iter 50 value 1166.906878
iter 60 value 1166.883106
iter 70 value 1166.845932
final value 1166.845921
converged

The function value is the negated log-likelihood

We see that about 70 additional iterations are sufficient for convergence (as opposed to when
starting the likelihood optimization from scratch, in which case 100 iterations are not sufficient
for convergence).

8. Closing discussion

This article introduced the R package **funGp**, which allows the treatment of regression problems involving multiple scalar and/or functional inputs, and a scalar output, through the fairly general GP model. In addition to being one of the few packages providing functional-input regression under this type of model, **funGp** features the ability to automatically set
R> m4 <- update(m1, var.re = TRUE, ls_s.re = TRUE, ls_f.re = TRUE,
+ extend = FALSE, trace = FALSE)

iter 10 value 1192.597779
iter 20 value 1184.610629
iter 30 value 1176.694606
iter 40 value 1176.082474
iter 50 value 1174.242786
iter 60 value 1169.484368
iter 70 value 1168.365643
iter 80 value 1167.828297
iter 90 value 1167.824714
iter 100 value 1167.623696
final value 1167.621387
stopped after 101 iterations

The function value is the negated log-likelihood

As for the creation of m2 above, since we use the default number 100 of iteration limit, the
iteration limit is reached before convergence. Next, using extend = TRUE, we can restart the
likelihood optimization from the non-converged set of hyperparameter estimates of m1.

R> m5 <- update(m1, var.re = TRUE, ls_s.re = TRUE, ls_f.re = TRUE,
+ extend = TRUE, trace = FALSE)

iter 10 value 1166.974505
iter 20 value 1166.964709
iter 30 value 1166.961533
iter 40 value 1166.913197
iter 50 value 1166.906878
iter 60 value 1166.883106
iter 70 value 1166.845932
final value 1166.845921
converged

The function value is the negated log-likelihood

We see that about 70 additional iterations are sufficient for convergence (as opposed to when
starting the likelihood optimization from scratch, in which case 100 iterations are not sufficient
for convergence).

8. Closing discussion

This article introduced the R package **funGp**, which allows the treatment of regression problems involving multiple scalar and/or functional inputs, and a scalar output, through the fairly general GP model. In addition to being one of the few packages providing functional-input regression under this type of model, **funGp** features the ability to automatically set
the structural parameters of the model for prediction performance. So far, this is a distinctive characteristic rarely found in current open-source regression software. The package is also equipped with multiple plotting functions that enable the quick inspection of the outputs delivered by its main functions: `fgpm()`, `fgpm_factory()`, `predict()` and `simulate()`. Moreover, parallelization has been incorporated as an efficiency booster for the procedures of model construction and structural optimization. The robustness and stability of **funGp** has been validated on numerous analytic examples and on the real life early warning application.

We envisage the extension of the package in various directions, including the addition of complementary levels of the structural parameters already available (e.g., other projection methods, distance functions and kernel functions), incorporation of additional structural parameters (e.g., type of transformation function for the output) and implementation of alternative structural optimization methods (e.g., brute-force search and Genetic Algorithms). Our code implementations were made with these potential extensions in mind and thus, the package is architECTured under the principle of modularity for easy extension. Further potential extensions, both in terms of methodology and implementation, are possible. It would be valuable to allow for functional outputs and for multi-dimensional functional inputs, for instance two-dimensional maps numerically given by matrices. For multi-dimensional functional inputs, tensor product splines may be considered. As discussed in Section 3, an extension of the **funGp** package from simple to universal Kriging could be investigated, hence allowing to include a trend in the GP model. Note also that the automatic structural configuration procedure of Section 5 currently aims at returning a “best” structural configuration. Since close-to-optimal predictive performances could be achieved by very different structural configurations, an interesting extension of the **funGp** package could be to carry out Bayesian model averaging of multiple structural configurations, in the spirit of Zhang and Taflanidis (2019). Sections 4.5 and 4.6 of Garnett (2022) also discuss Bayesian priors on structural configurations for GPs, together with efficient methods for exploring the posterior. In relation with this, note that the Bayesian framework on the covariance function could enable to have credible intervals that take the uncertainty on the covariance parameters into account. Another extension perspective of the **funGp** package is to allow for covariance computations for functional inputs based on $L^2$ norms applied directly to the functions rather than their projections, and weighted by functional correlation lengths. This would provide new distances, sharing commonalities with (6), but where $\theta_f^{(k)}$ would be a function. The **funGp** package is available from GitHub (https://github.com/djbetancourt-gh/funGp) and CRAN (https://CRAN.R-project.org/package=funGp). We encourage the community to make contributions in the proposed topics or any other found relevant.

## Acknowledgments

The **funGp** package was conceived within the RISCOPE project, funded by the French Agence Nationale de la Recherche (grant ANR-16-CE04-0011). We thank the ANR for its support. We are also grateful to Fabrice Gamboa (IMT, ANITI) for his helpful recommendations on GP modeling and to Juliette Garcia (ENAC) for her assistance on the stabilization of the ant colony algorithm for structural parameter optimization. We are very grateful to two anonymous reviewers and an action editor for constructive comments and suggestions that lead to an improvement of the manuscript.
the structural parameters of the model for prediction performance. So far, this is a distinctive characteristic rarely found in current open-source regression software. The package is also equipped with multiple plotting functions that enable the quick inspection of the outputs delivered by its main functions: `fgpm()`, `fgpm_factory()`, `predict()` and `simulate()`. Moreover, parallelization has been incorporated as an efficiency booster for the procedures of model construction and structural optimization. The robustness and stability of **funGp** has been validated on numerous analytic examples and on the real life early warning application.

We envisage the extension of the package in various directions, including the addition of complementary levels of the structural parameters already available (e.g., other projection methods, distance functions and kernel functions), incorporation of additional structural parameters (e.g., type of transformation function for the output) and implementation of alternative structural optimization methods (e.g., brute-force search and Genetic Algorithms). Our code implementations were made with these potential extensions in mind and thus, the package is architECTured under the principle of modularity for easy extension. Further potential extensions, both in terms of methodology and implementation, are possible. It would be valuable to allow for functional outputs and for multi-dimensional functional inputs, for instance two-dimensional maps numerically given by matrices. For multi-dimensional functional inputs, tensor product splines may be considered. As discussed in Section 3, an extension of the **funGp** package from simple to universal Kriging could be investigated, hence allowing to include a trend in the GP model. Note also that the automatic structural configuration procedure of Section 5 currently aims at returning a “best” structural configuration. Since close-to-optimal predictive performances could be achieved by very different structural configurations, an interesting extension of the **funGp** package could be to carry out Bayesian model averaging of multiple structural configurations, in the spirit of Zhang and Taflanidis (2019). Sections 4.5 and 4.6 of Garnett (2022) also discuss Bayesian priors on structural configurations for GPs, together with efficient methods for exploring the posterior. In relation with this, note that the Bayesian framework on the covariance function could enable to have credible intervals that take the uncertainty on the covariance parameters into account. Another extension perspective of the **funGp** package is to allow for covariance computations for functional inputs based on $L^2$ norms applied directly to the functions rather than their projections, and weighted by functional correlation lengths. This would provide new distances, sharing commonalities with (6), but where $\theta_f^{(k)}$ would be a function. The **funGp** package is available from GitHub (https://github.com/djbetancourt-gh/funGp) and CRAN (https://CRAN.R-project.org/package=funGp). We encourage the community to make contributions in the proposed topics or any other found relevant.

## Acknowledgments

The **funGp** package was conceived within the RISCOPE project, funded by the French Agence Nationale de la Recherche (grant ANR-16-CE04-0011). We thank the ANR for its support. We are also grateful to Fabrice Gamboa (IMT, ANITI) for his helpful recommendations on GP modeling and to Juliette Garcia (ENAC) for her assistance on the stabilization of the ant colony algorithm for structural parameter optimization. We are very grateful to two anonymous reviewers and an action editor for constructive comments and suggestions that lead to an improvement of the manuscript.
References

Allard D, Senoussi R, Porcu E (2016). “Anisotropy Models for Spatial Data.” *Mathematical Geosciences*, **48**(3), 305–328. doi:10.1007/s11004-015-9594-x.

Bachoc F (2013). “Cross Validation and Maximum Likelihood Estimations of Hyper-Parameters of Gaussian Processes with Model Misspecification.” *Computational Statistics & Data Analysis*, **66**, 55–69. doi:10.1016/j.csda.2013.03.016.

Bachoc F, Ammar K, Martinez JM (2016). “Improvement of Code Behavior in a Design of Experiments by Metamodeling.” *Nuclear Science and Engineering*, **183**(3), 387–406. doi:10.13182/nse15-108.

Ben Salem M, Bachoc F, Roustant O, Gamboa F, Tomaso L (2019). “Gaussian Process-Based Dimension Reduction for Goal-Oriented Sequential Design.” *SIAM/ASA Journal on Uncertainty Quantification*, **7**(4), 1369–1397. doi:10.1137/18m1167930.

Betancourt J (2020). *Functional-Input Metamodeling: An Application to Coastal Flood Early Warning*. Ph.D. thesis, Université Toulouse 3-Paul Sabatier. URL https://theses.hal.science/tel-03117879/.

Betancourt J, Bachoc F, Klein T, Gamboa F (2020a). “Ant Colony Based Model Selection for Functional-Input Gaussian Process Regression.” *Research Report hal-02532713*, Institut de Mathématiques de Toulouse & Institut Universitaire de France, UMR 5219, Université de Toulouse, CNRS, UPS IMT. Deliverable: D3.b, Reference: WP3.2, URL https://hal.archives-ouvertes.fr/hal-02532713.

Betancourt J, Bachoc F, Klein T, Idier D, Pedreros R, Rohmer J (2020b). “Gaussian Process Metamodeling of Functional-Input Code for Coastal Flood Hazard Assessment.” *Reliability Engineering & System Safety*, p. 106870. doi:10.1016/j.ress.2020.106870.

Betancourt J, Bachoc F, Klein T, Rohmer J (2024). **funGp**: Gaussian Process Models for Scalar and Functional Inputs. R package version 1.0.0, URL https://CRAN.R-project.org/package=funGp.

Blum C (2005). “Ant Colony Optimization: Introduction and Recent Trends.” *Physics of Life Reviews*, **2**(4), 353–373. doi:10.1016/j.plrev.2005.10.001.

Bonabeau E, Dorigo M, Theraulaz G (1999). *Swarm Intelligence: From Natural to Artificial Systems*. Oxford University Press.

Brockhaus S, Rügamer D, Greven S (2020). “Boosting Functional Regression Models with FDboost.” *Journal of Statistical Software*, **94**(10), 1–50. doi:10.18637/jss.v094.i10.

Byrd RH, Lu P, Nocedal J, Zhu C (1995). “A Limited Memory Algorithm for Bound Constrained Optimization.” *SIAM Journal on Scientific Computing*, **16**(5), 1190–1208. doi:10.1137/0916069.

Chen T, Dai B, Wang R, Liu D (2014). “Gaussian-Process-Based Real-Time Ground Segmentation for Autonomous Land Vehicles.” *Journal of Intelligent & Robotic Systems*, **76**(3-4), 563–582. doi:10.1007/s10846-013-9889-4.
References

Allard D, Senoussi R, Porcu E (2016). “Anisotropy Models for Spatial Data.” *Mathematical Geosciences*, **48**(3), 305–328. doi:10.1007/s11004-015-9594-x.

Bachoc F (2013). “Cross Validation and Maximum Likelihood Estimations of Hyper-Parameters of Gaussian Processes with Model Misspecification.” *Computational Statistics & Data Analysis*, **66**, 55–69. doi:10.1016/j.csda.2013.03.016.

Bachoc F, Ammar K, Martinez JM (2016). “Improvement of Code Behavior in a Design of Experiments by Metamodeling.” *Nuclear Science and Engineering*, **183**(3), 387–406. doi:10.13182/nse15-108.

Ben Salem M, Bachoc F, Roustant O, Gamboa F, Tomaso L (2019). “Gaussian Process-Based Dimension Reduction for Goal-Oriented Sequential Design.” *SIAM/ASA Journal on Uncertainty Quantification*, **7**(4), 1369–1397. doi:10.1137/18m1167930.

Betancourt J (2020). *Functional-Input Metamodeling: An Application to Coastal Flood Early Warning*. Ph.D. thesis, Université Toulouse 3-Paul Sabatier. URL https://theses.hal.science/tel-03117879/.

Betancourt J, Bachoc F, Klein T, Gamboa F (2020a). “Ant Colony Based Model Selection for Functional-Input Gaussian Process Regression.” *Research Report hal-02532713*, Institut de Mathématiques de Toulouse & Institut Universitaire de France, UMR 5219, Université de Toulouse, CNRS, UPS IMT. Deliverable: D3.b, Reference: WP3.2, URL https://hal.archives-ouvertes.fr/hal-02532713.

Betancourt J, Bachoc F, Klein T, Idier D, Pedreros R, Rohmer J (2020b). “Gaussian Process Metamodeling of Functional-Input Code for Coastal Flood Hazard Assessment.” *Reliability Engineering & System Safety*, p. 106870. doi:10.1016/j.ress.2020.106870.

Betancourt J, Bachoc F, Klein T, Rohmer J (2024). **funGp**: Gaussian Process Models for Scalar and Functional Inputs. R package version 1.0.0, URL https://CRAN.R-project.org/package=funGp.

Blum C (2005). “Ant Colony Optimization: Introduction and Recent Trends.” *Physics of Life Reviews*, **2**(4), 353–373. doi:10.1016/j.plrev.2005.10.001.

Bonabeau E, Dorigo M, Theraulaz G (1999). *Swarm Intelligence: From Natural to Artificial Systems*. Oxford University Press.

Brockhaus S, Rügamer D, Greven S (2020). “Boosting Functional Regression Models with FDboost.” *Journal of Statistical Software*, **94**(10), 1–50. doi:10.18637/jss.v094.i10.

Byrd RH, Lu P, Nocedal J, Zhu C (1995). “A Limited Memory Algorithm for Bound Constrained Optimization.” *SIAM Journal on Scientific Computing*, **16**(5), 1190–1208. doi:10.1137/0916069.

Chen T, Dai B, Wang R, Liu D (2014). “Gaussian-Process-Based Real-Time Ground Segmentation for Autonomous Land Vehicles.” *Journal of Intelligent & Robotic Systems*, **76**(3-4), 563–582. doi:10.1007/s10846-013-9889-4.
Cheng Y, Shi JQ (2016). **flars**: Functional LARS. R package version 1.0, URL https://CRAN.R-project.org/package=flars.

Chiles JP, Delfiner P (2009). *Geostatistics: Modeling Spatial Uncertainty*, volume 497. John Wiley & Sons.

Cressie N (1990). "The Origins of Kriging." *Mathematical Geology*, **22**(3), 239-252. doi: **10.1007/bf00889887**.

Dancik GM, Dorman KS (2008). “mlegp: Statistical Analysis for Computer Models of Biological Systems Using R.” *Bioinformatics*, **24**(17), 1966–1967. doi: **10.1093/bioinformatics/btn329**.

De Boor C (1978). *A Practical Guide to Splines*. Springer-Verlag. doi: **10.1007/978-1-4612-6333-3**.

De G Matthews AG, Van der Wilk M, Nickson T, Fujii K, Boukouvalas A, León-Villagrá P, Ghahramani Z, Hensman J (2017). “GPflow: A Gaussian Process Library Using Tensor-Flow.” *The Journal of Machine Learning Research*, **18**(40), 1–6.

Deville Y, Ginsbourger D, Roustant O (2021). **kergp**: Gaussian Process Laboratory. R package version 0.5.5, URL https://CRAN.R-project.org/package=kergp.

Dorigo M, Gambardella LM (1997). “Ant Colony System: A Cooperative Learning Approach to the Traveling Salesman Problem.” *IEEE Transactions on Evolutionary Computation*, **1**(1), 53–66. doi: **10.1109/4235.585892**.

Drucker H, Burges CJC, Kaufman L, Smola AJ, Vapnik V (1997). “Support Vector Regression Machines.” In *Advances in Neural Information Processing Systems*, pp. 155–161.

Dubrule O (1983). “Cross Validation of Kriging in a Unique Neighborhood.” *Journal of the International Association for Mathematical Geology*, **15**(6), 687–699. doi: **10.1007/bf01033232**.

Febrero-Bande M, Oviedo de la Fuente M (2012). “Statistical Computing in Functional Data Analysis: The R Package fda.usc.” *Journal of Statistical Software*, **51**(4), 1–28. doi: **10.18637/jss.v051.i04**.

Ferraty F, Vieu P (2006). “Reference Manual for Implementing NonParametric Functional Data Analysis (NPFDA).” *Companion Manual of the Book: Nonparametric Functional Data Analysis*, Theory and Practice, Springer-Verlag, New York.

Fonseca DJ, Navaresse DO, Moynihan GP (2003). “Simulation Metamodeling through Artificial Neural Networks.” *Engineering Applications of Artificial Intelligence*, **16**(3), 177–183. doi: **10.1016/S0952-1976(03)00043-5**.

Forrester A, Sobester A, Keane A (2008). *Engineering Design via Surrogate Modelling: A Practical Guide*. John Wiley & Sons.

Friedman JH (1991). “Multivariate Adaptive Regression Splines.” *The Annals of Statistics*, **19**(1), 1–67. doi: **10.1214/aos/1176347963**.
Cheng Y, Shi JQ (2016). **flars**: Functional LARS. R package version 1.0, URL https://CRAN.R-project.org/package=flars.

Chiles JP, Delfiner P (2009). *Geostatistics: Modeling Spatial Uncertainty*, volume 497. John Wiley & Sons.

Cressie N (1990). "The Origins of Kriging." *Mathematical Geology*, **22**(3), 239-252. doi: 10.1007/bf00889887.

Dancik GM, Dorman KS (2008). “mlegp: Statistical Analysis for Computer Models of Biological Systems Using R.” *Bioinformatics*, **24**(17), 1966–1967. doi:10.1093/bioinformatics/btn329.

De Boor C (1978). *A Practical Guide to Splines*. Springer-Verlag. doi:10.1007/978-1-4612-6333-3.

De G Matthews AG, Van der Wilk M, Nickson T, Fujii K, Boukouvalas A, León-Villagrá P, Ghahramani Z, Hensman J (2017). “GPflow: A Gaussian Process Library Using Tensor-Flow.” *The Journal of Machine Learning Research*, **18**(40), 1–6.

Deville Y, Ginsbourger D, Roustant O (2021). **kergp**: Gaussian Process Laboratory. R package version 0.5.5, URL https://CRAN.R-project.org/package=kergp.

Dorigo M, Gambardella LM (1997). “Ant Colony System: A Cooperative Learning Approach to the Traveling Salesman Problem.” *IEEE Transactions on Evolutionary Computation*, **1**(1), 53–66. doi:10.1109/4235.585892.

Drucker H, Burges CJC, Kaufman L, Smola AJ, Vapnik V (1997). “Support Vector Regression Machines.” In *Advances in Neural Information Processing Systems*, pp. 155–161.

Dubrule O (1983). “Cross Validation of Kriging in a Unique Neighborhood.” *Journal of the International Association for Mathematical Geology*, **15**(6), 687–699. doi:10.1007/bf01033232.

Febrero-Bande M, Oviedo de la Fuente M (2012). “Statistical Computing in Functional Data Analysis: The R Package **fda.usc**.” *Journal of Statistical Software*, **51**(4), 1–28. doi: 10.18637/jss.v051.i04.

Ferraty F, Vieu P (2006). “Reference Manual for Implementing NonParametric Functional Data Analysis (NPFDA).” Companion Manual of the Book: Nonparametric Functional Data Analysis, Theory and Practice, Springer-Verlag, New York.

Fonseca DJ, Navaresse DO, Moynihan GP (2003). “Simulation Metamodeling through Artificial Neural Networks.” *Engineering Applications of Artificial Intelligence*, **16**(3), 177–183. doi:10.1016/s0952-1976(03)00043-5.

Forrester A, Sobester A, Keane A (2008). *Engineering Design via Surrogate Modelling: A Practical Guide*. John Wiley & Sons.

Friedman JH (1991). “Multivariate Adaptive Regression Splines.” *The Annals of Statistics*, **19**(1), 1–67. doi:10.1214/aos/1176347963.
Gajardo A, Bhattacharjee S, Carroll C, Chen Y, Dai X, Fan J, Hadjipantelis PZ, Han K, Ji H, Zhu C, Müller HG, Wang JL (2022). **fdapace**: Functional Data Analysis and Empirical Dynamics. R package version 0.5.9, URL https://CRAN.R-project.org/package=fdapace.

García-Portugués E, Álvarez-Liébana J (2022). **goffda**: Goodness-of-Fit Tests for Functional Data. R package version 0.1.1, URL https://CRAN.R-project.org/package=goffda.

Gardner J, Pleiss G, Weinberger KQ, Bindel D, Wilson AG (2018). “Gpytorch: Blackbox Matrix-Matrix Gaussian Process Inference with GPU Acceleration.” *Advances in Neural Information Processing Systems*, **31**. doi:10.26226/m.64c26777632e9539aa87d541.

Garnett R (2022). *Bayesian Optimization*. Cambridge University Press.

Goldsmith J, Scheipl F, Huang L, Wrobel J, Di C, Gellar J, Harezlak J, McLean MW, Swihart B, Xiao L, Crainiceanu C, Reiss PT (2023). **refund**: Regression with Functional Data. R package version 0.1-32, URL https://CRAN.R-project.org/package=refund.

GPy Authors (2012). “GPy: A Gaussian Process Framework in Python.” URL https://github.com/SheffieldML/GPy.

Gu M, Palomo J, Berger J (2022). **RobustGaSP**: Robust Gaussian Stochastic Process Emulation. R package version 0.6.5, URL https://CRAN.R-project.org/package=RobustGaSP.

Idier D, Aurouet A, Bachoc F, Baills A, Betancourt J, Gamboa F, Klein T, López-Lopera AF, Pedreros R, Rohmer J, Thibault A (2021). “A User-Oriented Local Coastal Flooding Early Warning System Using Metamodelling Techniques.” *Journal of Marine Science and Engineering*, **9**(11), 1191. doi:10.3390/jmse9111191.

Iooss B, Marrel A (2019). “Advanced Methodology for Uncertainty Propagation in Computer Experiments with Large Number of Inputs.” *Nuclear Technology*, **205**(12), 1588–1606. doi:10.1080/00295450.2019.1573617.

Jolliffe I (2002). *Principal Component Analysis*. 2nd edition. Springer-Verlag. doi:10.1007/978-1-4757-1904-8.

Lee Y, Park JS (2017). “Model Selection Algorithm in Gaussian Process Regression for Computer Experiments.” *Communications for Statistical Applications and Methods*, **24**(4), 383–396. doi:10.5351/csam.2017.24.4.383.

Li Y, Soleimani H, Zohal M (2019). “An Improved Ant Colony Optimization Algorithm for the Multi-Depot Green Vehicle Routing Problem with Multiple Objectives.” *Journal of Cleaner Production*, **227**, 1161–1172. doi:10.1016/j.jclepro.2019.03.185.

Liu X, Guillas S (2017). “Dimension Reduction for Gaussian Process Emulation: An Application to the Influence of Bathymetry on Tsunami Heights.” *SIAM/ASA Journal on Uncertainty Quantification*, **5**(1), 787–812. doi:10.1137/16m1090648.

MacDonald B, Ranjan P, Chipman H (2015). “GPfit: An R Package for Fitting a Gaussian Process Model to Deterministic Simulator Outputs.” *Journal of Statistical Software*, **64**(12), 1–23. doi:10.18637/jss.v064.i12.
Gajardo A, Bhattacharjee S, Carroll C, Chen Y, Dai X, Fan J, Hadjipantelis PZ, Han K, Ji H, Zhu C, Müller HG, Wang JL (2022). **fdapace**: Functional Data Analysis and Empirical Dynamics. R package version 0.5.9, URL https://CRAN.R-project.org/package=fdapace.

García-Portugués E, Álvarez-Liébana J (2022). **goffda**: Goodness-of-Fit Tests for Functional Data. R package version 0.1.1, URL https://CRAN.R-project.org/package=goffda.

Gardner J, Pleiss G, Weinberger KQ, Bindel D, Wilson AG (2018). “Gpytorch: Blackbox Matrix-Matrix Gaussian Process Inference with GPU Acceleration.” *Advances in Neural Information Processing Systems*, **31**. doi:10.26226/m.64c26777632e9539aa87d541.

Garnett R (2022). *Bayesian Optimization*. Cambridge University Press.

Goldsmith J, Scheipl F, Huang L, Wrobel J, Di C, Gellar J, Harezlak J, McLean MW, Swihart B, Xiao L, Crainiceanu C, Reiss PT (2023). **refund**: Regression with Functional Data. R package version 0.1-32, URL https://CRAN.R-project.org/package=refund.

GPy Authors (2012). “GPy: A Gaussian Process Framework in Python.” URL https://github.com/SheffieldML/GPy.

Gu M, Palomo J, Berger J (2022). **RobustGaSP**: Robust Gaussian Stochastic Process Emulation. R package version 0.6.5, URL https://CRAN.R-project.org/package=RobustGaSP.

Idier D, Aurouet A, Bachoc F, Baills A, Betancourt J, Gamboa F, Klein T, López-Lopera AF, Pedreros R, Rohmer J, Thibault A (2021). “A User-Oriented Local Coastal Flooding Early Warning System Using Metamodelling Techniques.” *Journal of Marine Science and Engineering*, **9**(11), 1191. doi:10.3390/jmse9111191.

Iooss B, Marrel A (2019). “Advanced Methodology for Uncertainty Propagation in Computer Experiments with Large Number of Inputs.” *Nuclear Technology*, **205**(12), 1588–1606. doi:10.1080/00295450.2019.1573617.

Jolliffe I (2002). *Principal Component Analysis*. 2nd edition. Springer-Verlag. doi:10.1007/978-1-4757-1904-8.

Lee Y, Park JS (2017). “Model Selection Algorithm in Gaussian Process Regression for Computer Experiments.” *Communications for Statistical Applications and Methods*, **24**(4), 383–396. doi:10.5351/csam.2017.24.4.383.

Li Y, Soleimani H, Zohal M (2019). “An Improved Ant Colony Optimization Algorithm for the Multi-Depot Green Vehicle Routing Problem with Multiple Objectives.” *Journal of Cleaner Production*, **227**, 1161–1172. doi:10.1016/j.jclepro.2019.03.185.

Liu X, Guillas S (2017). “Dimension Reduction for Gaussian Process Emulation: An Application to the Influence of Bathymetry on Tsunami Heights.” *SIAM/ASA Journal on Uncertainty Quantification*, **5**(1), 787–812. doi:10.1137/16m1090648.

MacDonald B, Ranjan P, Chipman H (2015). “GPfit: An R Package for Fitting a Gaussian Process Model to Deterministic Simulator Outputs.” *Journal of Statistical Software*, **64**(12), 1–23. doi:10.18637/jss.v064.i12.
Markussen B (2023). **fdaMixed**: Functional Data Analysis in a Mixed Model Framework. R package version 0.6.1, URL https://CRAN.R-project.org/package=fdaMixed.

Marrel A, Iooss B, Van Dorpe F, Volkova E (2008). "An Efficient Methodology for Modeling Complex Computer Codes with Gaussian Processes." *Computational Statistics & Data Analysis*, **52**(10), 4731-4744. doi:10.1016/j.csda.2008.03.026.

Mori H, Kurata E (2008). "Application of Gaussian Process to Wind Speed Forecasting for Wind Power Generation." In *2008 IEEE International Conference on Sustainable Energy Technologies*, pp. 956-959. IEEE.

Muehlenstaedt T, Fruth J, Roustant O (2017). "Computer Experiments with Functional Inputs and Scalar Outputs by a Norm-Based Approach." *Statistics and Computing*, **27**(4), 1083-1097. doi:10.1007/s11222-016-9672-z.

Nanty S, Helbert C, Marrel A, Pérot N, Prieur C (2016). "Sampling, Metamodeling, and Sensitivity Analysis of Numerical Simulators with Functional Stochastic Inputs." *SIAM/ASA Journal on Uncertainty Quantification*, **4**(1), 636-659. doi:10.1137/15m1033319.

Neufeld A, Heggeseth B (2019). **splinetree**: Longitudinal Regression Trees and Forests. R package version 0.2.0, URL https://CRAN.R-project.org/package=splinetree.

Nilsson J, De Jong S, Smilde AK (1997). "Multiway Calibration in 3D QSAR." *Journal of Chemometrics*, **11**(6), 511-524. doi:10.1002/(sici)1099-128x(199711/12)11:6<511::aid-cem488>3.0.co;2-w.

Oakley J, O'Hagan A (2002). "Bayesian Inference for the Uncertainty Distribution of Computer Model Outputs." *Biometrika*, **89**(4), 769-784. doi:10.1093/biomet/89.4.769.

Østergård T, Jensen RL, Maagaard SE (2018). "A Comparison of Six Metamodeling Techniques Applied to Building Performance Simulations." *Applied Energy*, **211**, 89-103. doi: 10.1016/j.apenergy.2017.10.102.

Papaioannou I, Ehre M, Straub D (2019). "PLS-Based Adaptation for Efficient PCE Representation in High Dimensions." *Journal of Computational Physics*, **387**, 186-204. doi: 10.1016/j.jcp.2019.02.046.

Pedregosa F, Varoquaux G, Gramfort A, Michel V, Thirion B, Grisel O, Blondel M, Prettenhofer P, Weiss R, Dubourg V, Vanderplas J, Passos A, Cournapeau D, Brucher M, Perrot M, Duchesnay E (2011). "scikit-learn: Machine Learning in Python." *Journal of Machine Learning Research*, **12**(85), 2825-2830.

Ramsay JO, Graves S, Hooker G (2023). **fda**: Functional Data Analysis. R package version 6.1.4, URL https://CRAN.R-project.org/package=fda.

Ramsay JO, Silverman BW (2007). *Applied Functional Data Analysis: Methods and Case Studies*. Springer-Verlag. doi:10.1007/b98886.

Rasmussen CE, Nickisch H (2010). "Gaussian Processes for Machine Learning (GPML) Toolbox." *The Journal of Machine Learning Research*, **11**, 3011-3015.
Markussen B (2023). **fdaMixed**: Functional Data Analysis in a Mixed Model Framework. R package version 0.6.1, URL https://CRAN.R-project.org/package=fdaMixed.

Marrel A, Iooss B, Van Dorpe F, Volkova E (2008). "An Efficient Methodology for Modeling Complex Computer Codes with Gaussian Processes." *Computational Statistics & Data Analysis*, **52**(10), 4731-4744. doi:10.1016/j.csda.2008.03.026.

Mori H, Kurata E (2008). "Application of Gaussian Process to Wind Speed Forecasting for Wind Power Generation." In *2008 IEEE International Conference on Sustainable Energy Technologies*, pp. 956-959. IEEE.

Muehlenstaedt T, Fruth J, Roustant O (2017). "Computer Experiments with Functional Inputs and Scalar Outputs by a Norm-Based Approach." *Statistics and Computing*, **27**(4), 1083-1097. doi:10.1007/s11222-016-9672-z.

Nanty S, Helbert C, Marrel A, Pérot N, Prieur C (2016). "Sampling, Metamodeling, and Sensitivity Analysis of Numerical Simulators with Functional Stochastic Inputs." *SIAM/ASA Journal on Uncertainty Quantification*, **4**(1), 636-659. doi:10.1137/15m1033319.

Neufeld A, Heggeseth B (2019). **splinetree**: Longitudinal Regression Trees and Forests. R package version 0.2.0, URL https://CRAN.R-project.org/package=splinetree.

Nilsson J, De Jong S, Smilde AK (1997). "Multiway Calibration in 3D QSAR." *Journal of Chemometrics*, **11**(6), 511-524. doi:10.1002/(sici)1099-128x(199711/12)11:6<511::aid-cem488>3.0.co;2-w.

Oakley J, O'Hagan A (2002). "Bayesian Inference for the Uncertainty Distribution of Computer Model Outputs." *Biometrika*, **89**(4), 769-784. doi:10.1093/biomet/89.4.769.

Østergård T, Jensen RL, Maagaard SE (2018). "A Comparison of Six Metamodeling Techniques Applied to Building Performance Simulations." *Applied Energy*, **211**, 89-103. doi:10.1016/j.apenergy.2017.10.102.

Papaioannou I, Ehre M, Straub D (2019). "PLS-Based Adaptation for Efficient PCE Representation in High Dimensions." *Journal of Computational Physics*, **387**, 186-204. doi:10.1016/j.jcp.2019.02.046.

Pedregosa F, Varoquaux G, Gramfort A, Michel V, Thirion B, Grisel O, Blondel M, Prettenhofer P, Weiss R, Dubourg V, Vanderplas J, Passos A, Cournapeau D, Brucher M, Perrot M, Duchesnay E (2011). "scikit-learn: Machine Learning in Python." *Journal of Machine Learning Research*, **12**(85), 2825-2830.

Ramsay JO, Graves S, Hooker G (2023). **fda**: Functional Data Analysis. R package version 6.1.4, URL https://CRAN.R-project.org/package=fda.

Ramsay JO, Silverman BW (2007). *Applied Functional Data Analysis: Methods and Case Studies*. Springer-Verlag. doi:10.1007/b98886.

Rasmussen CE, Nickisch H (2010). "Gaussian Processes for Machine Learning (GPML) Toolbox." *The Journal of Machine Learning Research*, **11**, 3011-3015.
Rasmussen CE, Williams CKI (2006). *Gaussian Processes for Machine Learning*, volume 2.
MIT Press, Cambridge. doi:10.7551/mitpress/3206.001.0001.

R Core Team (2023). *R: A Language and Environment for Statistical Computing*. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.

Rohmer J, Idier D (2012). “A Meta-Modelling Strategy to Identify the Critical Offshore Conditions for Coastal Flooding.” *Natural Hazards & Earth System Sciences*, **12**(9). doi: 10.5194/nhess-12-2943-2012.

Roustant O, Ginsbourger D, Deville Y (2012). “DiceKriging, DiceOptim: Two R Packages for the Analysis of Computer Experiments by Kriging-Based Metamodeling and Optimization.” *Journal of Statistical Software*, **51**(1), 1–55. doi:10.18637/jss.v051.i01.

Sacks J, Welch WJ, Mitchell TJ, Wynn HP (1989). “Design and Analysis of Computer Experiments.” *Statistical Science*, **4**, 409–423. doi:10.1214/ss/1177012413.

Saltelli A, Tarantola S, Campolongo F, Ratto M (2004). *Sensitivity Analysis in Practice: A Guide to Assessing Scientific Models*, volume 1. John Wiley & Sons.

Santner TJ, Williams BJ, Notz WI, Williams BJ (2003). *The Design and Analysis of Computer Experiments*. Springer-Verlag. doi:10.1007/978-1-4757-3799-8.

Savitsky TD (2016). “Bayesian Nonparametric Mixture Estimation for Time-Indexed Functional Data in R.” *Journal of Statistical Software*, **72**(2), 1–34. doi:10.18637/jss.v072.i02.

Shi JQ, Cheng Y, Konzen E (2023). **GPFDA: Gaussian Process for Functional Data Analysis.**
R package version 3.1.3, URL https://CRAN.R-project.org/package=GPFDA.

Simpson TW, Poplinski JD, Koch PN, Allen JK (2001). “Metamodels for Computer-Based Engineering Design: Survey and Recommendations.” *Engineering with Computers*, **17**(2), 129–150. doi:10.1007/p100007198.

Singh SS, Singh K, Kumar A, Biswas B (2020). “ACO-IM: Maximizing Influence in Social Networks Using Ant Colony Optimization.” *Soft Computing*, **24**(13), 10181–10203. doi: 10.1007/s00500-019-04533-y.

Snelson E, Ghahramani Z (2005). “Sparse Gaussian Processes Using Pseudo-Inputs.” *Advances in Neural Information Processing Systems*, **18**. doi:10.7551/mitpress/7503.003. 0041.

Stein ML (1999). *Interpolation of Spatial Data: Some Theory for Kriging*. Springer-Verlag.
doi:10.1007/978-1-4612-1494-6.

Tsiliki G, Munteanu C, Seoane J, Fernandez-Lozano C, Sarimveis H, Willighagen E (2015). “RRegrs: An R Package for Computer-Aided Model Selection with Multiple Regression Models.” *Journal of Cheminformatics*, **7**(1), 46. doi:10.1186/s13321-015-0094-2.

Vanhatalo J, Riihimäki J, Hartikainen J, Jylänki P, Tolvanen V, Vehtari A (2012). “Bayesian Modeling with Gaussian Processes Using the GPstuff Toolbox.” *arXiv 1206.5754*, arXiv.org E-Print Archive. doi:10.48550/arXiv.1206.5754.
Rasmussen CE, Williams CKI (2006). *Gaussian Processes for Machine Learning*, volume 2.
MIT Press, Cambridge. doi:10.7551/mitpress/3206.001.0001.

R Core Team (2023). *R: A Language and Environment for Statistical Computing*. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.

Rohmer J, Idier D (2012). “A Meta-Modelling Strategy to Identify the Critical Offshore Conditions for Coastal Flooding.” *Natural Hazards & Earth System Sciences*, **12**(9). doi: 10.5194/nhess-12-2943-2012.

Roustant O, Ginsbourger D, Deville Y (2012). “DiceKriging, DiceOptim: Two R Packages for the Analysis of Computer Experiments by Kriging-Based Metamodeling and Optimization.” *Journal of Statistical Software*, **51**(1), 1–55. doi:10.18637/jss.v051.i01.

Sacks J, Welch WJ, Mitchell TJ, Wynn HP (1989). “Design and Analysis of Computer Experiments.” *Statistical Science*, **4**, 409–423. doi:10.1214/ss/1177012413.

Saltelli A, Tarantola S, Campolongo F, Ratto M (2004). *Sensitivity Analysis in Practice: A Guide to Assessing Scientific Models*, volume 1. John Wiley & Sons.

Santner TJ, Williams BJ, Notz WI, Williams BJ (2003). *The Design and Analysis of Computer Experiments*. Springer-Verlag. doi:10.1007/978-1-4757-3799-8.

Savitsky TD (2016). “Bayesian Nonparametric Mixture Estimation for Time-Indexed Functional Data in R.” *Journal of Statistical Software*, **72**(2), 1–34. doi:10.18637/jss.v072.i02.

Shi JQ, Cheng Y, Konzen E (2023). **GPFDA: Gaussian Process for Functional Data Analysis.**
R package version 3.1.3, URL https://CRAN.R-project.org/package=GPFDA.

Simpson TW, Poplinski JD, Koch PN, Allen JK (2001). “Metamodels for Computer-Based Engineering Design: Survey and Recommendations.” *Engineering with Computers*, **17**(2), 129–150. doi:10.1007/p100007198.

Singh SS, Singh K, Kumar A, Biswas B (2020). “ACO-IM: Maximizing Influence in Social Networks Using Ant Colony Optimization.” *Soft Computing*, **24**(13), 10181–10203. doi: 10.1007/s00500-019-04533-y.

Snelson E, Ghahramani Z (2005). “Sparse Gaussian Processes Using Pseudo-Inputs.” *Advances in Neural Information Processing Systems*, **18**. doi:10.7551/mitpress/7503.003. 0041.

Stein ML (1999). *Interpolation of Spatial Data: Some Theory for Kriging*. Springer-Verlag.
doi:10.1007/978-1-4612-1494-6.

Tsiliki G, Munteanu C, Seoane J, Fernandez-Lozano C, Sarimveis H, Willighagen E (2015). “RRegrs: An R Package for Computer-Aided Model Selection with Multiple Regression Models.” *Journal of Cheminformatics*, **7**(1), 46. doi:10.1186/s13321-015-0094-2.

Vanhatalo J, Riihimäki J, Hartikainen J, Jylänki P, Tolvanen V, Vehtari A (2012). “Bayesian Modeling with Gaussian Processes Using the GPstuff Toolbox.” *arXiv 1206.5754*, arXiv.org E-Print Archive. doi:10.48550/arXiv.1206.5754.
Wilson A, Nickisch H (2015). “Kernel Interpolation for Scalable Structured Gaussian Processes (KISS-GP).” In *International Conference on Machine Learning*, pp. 1775–1784. PMLR.

Yuan J, Nian V (2018). “Ship Energy Consumption Prediction with Gaussian Process Meta-model.” *Energy Procedia*, **152**, 655–660. doi:10.1016/j.egypro.2018.09.226.

Zhang H, Wang Y (2010). “Kriging and Cross-Validation for Massive Spatial Data.” *Environmetrics: The Official Journal of the International Environmetrics Society*, **21**(3-4), 290–304. doi:10.1002/env.1023.

Zhang J, Taflanidis AA (2019). “Bayesian Model Averaging for Kriging Regression Structure Selection.” *Probabilistic Engineering Mechanics*, **56**, 58–70. doi:10.1016/j.probengmech.2019.02.002.
Wilson A, Nickisch H (2015). “Kernel Interpolation for Scalable Structured Gaussian Processes (KISS-GP).” In *International Conference on Machine Learning*, pp. 1775–1784. PMLR.

Yuan J, Nian V (2018). “Ship Energy Consumption Prediction with Gaussian Process Meta-model.” *Energy Procedia*, **152**, 655–660. doi:10.1016/j.egypro.2018.09.226.

Zhang H, Wang Y (2010). “Kriging and Cross-Validation for Massive Spatial Data.” *Environmetrics: The Official Journal of the International Environmetrics Society*, **21**(3-4), 290–304. doi:10.1002/env.1023.

Zhang J, Taflanidis AA (2019). “Bayesian Model Averaging for Kriging Regression Structure Selection.” *Probabilistic Engineering Mechanics*, **56**, 58–70. doi:10.1016/j.probengmech.2019.02.002.
## A. funGp at a glance

The **funGp** package allows the easy construction of GP regression models with special treatment of functional inputs, an aspect scarcely treated by other R packages. The package is constituted by two major modules which are described in the paragraphs that follow.

**Regression.** This module is responsible for the core functionalities of the package: (i) Creation and diagnostic of GP regression models, (ii) prediction at unobserved input coordinates, (iii) conditional simulation, and (iv) model updating. These procedures are implemented in the five functions listed in Table 2. All the listed functions are set up with default arguments that allow an effortless first interaction with the package. However, high degree of customization is enabled to let the users get greater control over the outputs and extend the scope of the analysis as they become familiar with the overall functioning of the package. As an example, the `fgpm()` function allows the customization of the so-called structural parameters and the definition of custom hyperparameters, among other features.

<table><caption>Table 2: Base **funGp** functionalities related to the class "fgpm".</caption><thead><tr><td>Function</td><td>Type</td><td>Description</td></tr></thead><tbody><tr><td><code>fgpm()</code></td><td>Function</td><td>Creation of GP regression models. Outputs an object of S4 class <code>"fgpm"</code> representing the fitted model</td></tr><tr><td><code>predict()</code></td><td>Method</td><td>Prediction at new input points. Delivers the estimated output values and corresponding uncertainty at specified points</td></tr><tr><td><code>simulate()</code></td><td>Method</td><td>Conditional simulation. Delivers the simulated output values and corresponding uncertainty at specified points</td></tr><tr><td><code>update()</code></td><td>Method</td><td>Efficient model update. Allows to add, delete or substitute data, and also to replace or re-estimate hyperparameters</td></tr><tr><td><code>plot()</code></td><td>Method</td><td>Diagnostic plot for GP model; receives an <code>fgpm</code> object as main argument</td></tr></tbody></table>

**Structural optimization.** The set of functions listed in Table 2 could constitute by themselves a regression package. However, **funGp** goes further and incorporates the automatic calibration of structural parameters of the model, a functionality widely neglected by other regression tools. Structural optimization enables the efficient exploration of the space of possible model configurations, seeking for optimal predictive performance. Diagnostic plots are provided for a quick evaluation of the selected model and the evolution of the optimization. The main procedures of this module are implemented in the three functions listed in Table 3.

## B. Manual prediction

As explained in Section 3.2, GP predictions are accessible in **funGp** through the `predict()` method. However, as a verification/inspection procedure, it is also possible to recreate the different data structures necessary for making the predictions manually based on the equations (2) and (3) of Gaussian conditioning. The recovery of these data structures is explained below, departing from an hypothetical model `m1` created with `fgpm()`.
## A. funGp at a glance

The **funGp** package allows the easy construction of GP regression models with special treatment of functional inputs, an aspect scarcely treated by other R packages. The package is constituted by two major modules which are described in the paragraphs that follow.

**Regression.** This module is responsible for the core functionalities of the package: (i) Creation and diagnostic of GP regression models, (ii) prediction at unobserved input coordinates, (iii) conditional simulation, and (iv) model updating. These procedures are implemented in the five functions listed in Table 2. All the listed functions are set up with default arguments that allow an effortless first interaction with the package. However, high degree of customization is enabled to let the users get greater control over the outputs and extend the scope of the analysis as they become familiar with the overall functioning of the package. As an example, the `fgpm()` function allows the customization of the so-called structural parameters and the definition of custom hyperparameters, among other features.

<table><caption>Table 2: Base **funGp** functionalities related to the class "fgpm".</caption><thead><tr><td>Function</td><td>Type</td><td>Description</td></tr></thead><tbody><tr><td><code>fgpm()</code></td><td>Function</td><td>Creation of GP regression models. Outputs an object of S4 class <code>"fgpm"</code> representing the fitted model</td></tr><tr><td><code>predict()</code></td><td>Method</td><td>Prediction at new input points. Delivers the estimated output values and corresponding uncertainty at specified points</td></tr><tr><td><code>simulate()</code></td><td>Method</td><td>Conditional simulation. Delivers the simulated output values and corresponding uncertainty at specified points</td></tr><tr><td><code>update()</code></td><td>Method</td><td>Efficient model update. Allows to add, delete or substitute data, and also to replace or re-estimate hyperparameters</td></tr><tr><td><code>plot()</code></td><td>Method</td><td>Diagnostic plot for GP model; receives an <code>fgpm</code> object as main argument</td></tr></tbody></table>

**Structural optimization.** The set of functions listed in Table 2 could constitute by themselves a regression package. However, **funGp** goes further and incorporates the automatic calibration of structural parameters of the model, a functionality widely neglected by other regression tools. Structural optimization enables the efficient exploration of the space of possible model configurations, seeking for optimal predictive performance. Diagnostic plots are provided for a quick evaluation of the selected model and the evolution of the optimization. The main procedures of this module are implemented in the three functions listed in Table 3.

## B. Manual prediction

As explained in Section 3.2, GP predictions are accessible in **funGp** through the `predict()` method. However, as a verification/inspection procedure, it is also possible to recreate the different data structures necessary for making the predictions manually based on the equations (2) and (3) of Gaussian conditioning. The recovery of these data structures is explained below, departing from an hypothetical model `m1` created with `fgpm()`.
<table><caption>Table 3: Complementary **funGp** functionnalities related to the class "Xfgpm" enabling structural optimization.</caption><thead><tr><th>Function</th><th>Type</th><th>Description</th></tr></thead><tbody><tr><td>fgpm_factory()</td><td>Function</td><td>Structural optimization of GP models. Outputs an object of S4 class "Xfgpm" containing the optimized model along with records from the optimization process</td></tr><tr><td>plot()</td><td>Method</td><td>which = "diag": Information on the absolute and relative quality of the selected model</td></tr><tr><td>plot()</td><td>Method</td><td>which = "evol": Verification plot showing the evolution in the quality of the explored structural configurations along the optimization</td></tr><tr><td>summary()</td><td>Method</td><td>Lists the evaluated structural configurations with their prediction performances</td></tr></tbody></table>

* Vector of observed output values $y$: this vector is stored at the @sOut slot of the fgpm object associated with the model.

R> y <- m1@sOut

* Unconditional $n \times n_*$ training-prediction cross-covariance matrix $C((X, F), (X_*, F_*))$ and $n_* \times n_*$ prediction auto-covariance matrix $C((X_*, F_*), (X_*, F_*))$: these matrices can be obtained from a full prediction generated by predict(), as shown below.

R> preds <- predict(m1, sIn.pr = sIn.pr, fIn.pr = fIn.pr, detail = "full")

R> K.tp <- preds$K.tp

R> K.pp <- preds$K.pp

* Unconditional $n \times n$ training auto-covariance matrix $C((X, F), (X, F))$: this matrix can be regenerated from the lower triangular matrix of its LU decomposition, stored at the @preMats$L$ slot of the model.

R> K.tt <- tcrossprod(m1@preMats$L)

Once the y, K.tp, K.pp and K.tt data structures have been recovered, the predicted output values and corresponding confidence intervals can be reproduced based on Gaussian conditioning as follows (with the two first lines providing the predictive means and variances and the next two providing the limits of the 95% confidence intervals).

R> y.pr <- t(K.tp) %*% solve(K.tt) %*% y

R> v.pr <- diag(K.pp - t(K.tp) %*% solve(K.tt) %*% K.tp)

R> ll <- y.pr - 1.96 * sqrt(v.pr)

R> ul <- y.pr + 1.96 * sqrt(v.pr)

We remark that this way of computing GP predictions is computationally inefficient, and is only recommended for purposes of verification and examination of the involved data structures. Predictions and simulations performed through the predict() and simulate() functions in **funGp** are highly optimized for computational efficiency.
<table><caption>Table 3: Complementary **funGp** functionnalities related to the class "Xfgpm" enabling structural optimization.</caption><thead><tr><th>Function</th><th>Type</th><th>Description</th></tr></thead><tbody><tr><td>fgpm_factory()</td><td>Function</td><td>Structural optimization of GP models. Outputs an object of S4 class "Xfgpm" containing the optimized model along with records from the optimization process</td></tr><tr><td>plot()</td><td>Method</td><td>which = "diag": Information on the absolute and relative quality of the selected model</td></tr><tr><td>plot()</td><td>Method</td><td>which = "evol": Verification plot showing the evolution in the quality of the explored structural configurations along the optimization</td></tr><tr><td>summary()</td><td>Method</td><td>Lists the evaluated structural configurations with their prediction performances</td></tr></tbody></table>

* Vector of observed output values $y$: this vector is stored at the @sOut slot of the fgpm object associated with the model.

R> y <- m1@sOut

* Unconditional $n \times n_*$ training-prediction cross-covariance matrix $C((X, F), (X_*, F_*))$ and $n_* \times n_*$ prediction auto-covariance matrix $C((X_*, F_*), (X_*, F_*))$: these matrices can be obtained from a full prediction generated by predict(), as shown below.

R> preds <- predict(m1, sIn.pr = sIn.pr, fIn.pr = fIn.pr, detail = "full")

R> K.tp <- preds$K.tp

R> K.pp <- preds$K.pp

* Unconditional $n \times n$ training auto-covariance matrix $C((X, F), (X, F))$: this matrix can be regenerated from the lower triangular matrix of its LU decomposition, stored at the @preMats$L$ slot of the model.

R> K.tt <- tcrossprod(m1@preMats$L)

Once the y, K.tp, K.pp and K.tt data structures have been recovered, the predicted output values and corresponding confidence intervals can be reproduced based on Gaussian conditioning as follows (with the two first lines providing the predictive means and variances and the next two providing the limits of the 95% confidence intervals).

R> y.pr <- t(K.tp) %*% solve(K.tt) %*% y

R> v.pr <- diag(K.pp - t(K.tp) %*% solve(K.tt) %*% K.tp)

R> ll <- y.pr - 1.96 * sqrt(v.pr)

R> ul <- y.pr + 1.96 * sqrt(v.pr)

We remark that this way of computing GP predictions is computationally inefficient, and is only recommended for purposes of verification and examination of the involved data structures. Predictions and simulations performed through the predict() and simulate() functions in **funGp** are highly optimized for computational efficiency.
Affiliation:

José Betancourt, François Bachoc, Thierry Klein

Institut de Mathématiques de Toulouse (IMT), Université de Toulouse

École Nationale de l'Aviation Civile (ENAC)

118 Route de Narbonne

31062 Toulouse Cedex 9, France

E-mail: fungp.rpack@gmail.com, francois.bachoc@math.univ-toulouse.fr,
thierry.klein@math.univ-toulouse.fr

URL: https://www.linkedin.com/in/djbetancourt,
https://www.math.univ-toulouse.fr/~fbachoc/
https://perso.math.univ-toulouse.fr/klein

Déborah Idier, Jérémy Rohmer

Bureau de Recherches Géologiques et Minières (BRGM)

3, av. Claude Guillemin

45060 Orléans Cedex 2, France

E-mail: D.Idier@brgm.fr, J.Rohmer@brgm.fr

URL: https://www.researchgate.net/profile/Deborah_Idier2
https://www.researchgate.net/profile/Jeremy_Rohmer

Yves Deville

Alpestat

E-mail: deville.yves@alpestat.com
Affiliation:

José Betancourt, François Bachoc, Thierry Klein

Institut de Mathématiques de Toulouse (IMT), Université de Toulouse

École Nationale de l'Aviation Civile (ENAC)

118 Route de Narbonne

31062 Toulouse Cedex 9, France

E-mail: fungp.rpack@gmail.com, francois.bachoc@math.univ-toulouse.fr,
thierry.klein@math.univ-toulouse.fr

URL: https://www.linkedin.com/in/djbetancourt,
https://www.math.univ-toulouse.fr/~fbachoc/
https://perso.math.univ-toulouse.fr/klein

Déborah Idier, Jérémy Rohmer

Bureau de Recherches Géologiques et Minières (BRGM)

3, av. Claude Guillemin

45060 Orléans Cedex 2, France

E-mail: D.Idier@brgm.fr, J.Rohmer@brgm.fr

URL: https://www.researchgate.net/profile/Deborah_Idier2
https://www.researchgate.net/profile/Jeremy_Rohmer

Yves Deville

Alpestat

E-mail: deville.yves@alpestat.com