## 1 Introduction

#### Case study and motivation

Prehistoric rock art paintings (see, for example, Figure 1) are exposed to environmental elements, which can accelerate their degradation, increasing the risk of losing such a valuable cultural piece from past societies. A part from and in addition to many other factors, exposure to sunlight can have adverse effects on these systems due to thermal and photochemical degradation of the historic materials [Diez-Herrero_2009].

This is a surface effect and the degree of color change depends on the chemical composition of the pigment. In this study we focused on the study and documentation of the degree of color changing/fading of paintings, patinas and host rock which is crucial for the conservation of these systems [Cassar_2001, del_hoyo_2015].

Materials with higher light sensitivity usually experience a rapid color change during the early stages of exposure, followed by a slower rate after maximum fading has occurred, assuming total disappearance of the substance that produces the color (chromophore) at this second stage of the fading [feller1986determination, giles_1965, giles_1968, johnston1984kinetics]. Thus, color fading is expected to be non-decreasing as a function of time and saturate in the long term.

Microfading Spectrometry (MFS) is a method for assessing light sensitivity color (spectral) variations of cultural heritage objects. This method allows for real time monitoring of spectral changes of a cultural heritage object by undertaking in situ lightfastness measurements on the surface under study [Whitmore_1999, Whitmore_2000, del_hoyo_2010]. MFS technique provides for each measured point on the surface of rock art paintings a time-series of observations that represent potential color fading along time to exposition to light [Whitmore_1999, del_hoyo_2010]. Thus, MFS measurements can be seen as observations of an underlying spatio-temporal stochastic process.

The MFS instrument is very sensitive to movement and glossy surface effects, causing extremely large fluctuations in the measurements. Furthermore, they can be easily contaminated by changes in the lighting conditions when the measurements are performed in outside environments. These large fluctuations and possible systematic noise effects in the observations can cause that models do not satisfy those properties of monotonicity and long-term stabilization of color fading over time.

Existing lightfastness studies on these systems have been limited to analyze a few measured points on the surface of the rock art paintings because these measurements are largely time-consuming and difficult to materialize under harsh conditions. So, in the paper we propose an interpolation procedure in order to extend the analysis to others unobserved locations on the surface of rock art paintings. A complete map of estimated color fading data for the entire surface under study is an important and useful information in order to achieve further successful conservation actions.

#### Background methods

Functional data usually refers to independent realizations of a functional random variable that takes values in a continuous space. Time-series of observations (e.g. color-fading time-series functions) might be the most common case of functional data in 1D,

, but spatially distributed observations can also be seen as functional data in a 2D space, , or spatio-temporal observations can also be seen as functional data in a 3D space, .In order to construct a model useful for making predictions of new functional data as function of new values of the variables in the input space, the process must be considered as an structured process with dependence among observations.

Correlated functional models consider the observed functional data as non-independent functions. A popular approach for correlated space-temporal functional data consists in three-way (spatial (2D) and temporal (1D)) penalized splines models [wood2003thin] with different basis constructions based on Kronecker products [currie2006generalized, lee2011p] or additive basis components [kneib2006structured]. aguilera2017prediction propose a mixture of functional regression model for functional response and penalized spline spatial regression. However, in general, these models become highly parameterized, hardly interpretable, and difficult to fit specially in a Bayesian framework and using sampling methods.

Another powerful approach consists in considering the space-time structured observations as stochastic realizations of a Gaussian process prior with a spatio-temporal covariance function. Gaussian processes (GPs) [Neal_1999, Rasmussen_2006] are a natural and flexible non-parametric prior models for -dimensional (e.g. space and time) functions with multivariate predictors (input variables) in each dimension. In a separable form, the space-time covariance function is a result of the interaction of the two independent processes, space and time [banerjee_2014]. In a non-separable form, the covariance function models space-time interaction [cressie1999classes, de2002nonseparable]. The drawback of GP models is their expensive computational demand in covariance matrix inversion that is in general a operation, with the number of spatial locations and the number of time points.

A geostatistical approach for spatio-temporal data is the kriging approach for 1D-functional data [Giraldo_2010], which the 1D-functional data consists in the time-series of the data. The spatial correlation of the time-series is modeled in the covariance function. The dimension of the covariance matrix is and the matrix inversion is a operation, with the number of spatial locations. This approach although requires less computation that the spatio-temporal GPs, has the drawback of being a quite less flexible model in the spatio-temporal structure since the same spatial structure is defined for the whole time-series. A related approach can be found in Baladandayuthapani_2008 where the spatial correlation between the time-series is modeled by defining GPs with a spatial covariance function across the time-series functional coefficients. This construction allows for modeling different spatial structure for the different orders of the coefficients.

Furthermore, inducing monotonicity or gradient to functions can be expressed in terms of the derivative of the function, as either the sign of a derivative for monotonicity or the value of a derivative for gradients. As diferentation is a linear operator, the derivative of both semiparametric models and Gaussian process models can be used as additional observations jointly with the regular observations in order to force the function to fit these properties [rasmussen_2003, shively_2009]. Riihimaki_2011 uses a Gaussian process model and the information of its derivative process to induce monotonicity on the functions using virtual observations of the sign of the derivatives.

On the other hand, monotonicity in semiparametric models can be forced by construction. shively_2009 proposes two approches to obtain monotonic functions, the first using a modified characterization of the smooth monotone functions proposed in ramsay_1998 that allows for unconstrained estimation, and the second using constrained prior distributions for the regression coefficients to ensure the monotonicity. Brezger_2008

induces monotonicity on semiparametric and non-parametric models imposing the restriction that the regression coefficients are ordered, also sometimes known as isotonic regression. These constraints are imposed by specifying truncated prior distributions in order to reject the undesired draws for the parameters in the MCMC sampling.

Reich_2011makes a similar approach imposing order to the regression parameters by means of reparameterizing and constraining these parameters with application to a quantile regression model.

#### Proposed methodology

In this paper, a specific application with the aim of modeling and predicting color fading time-series for new unobserved spatial locations on the surface of rock art paintings is carried out.

We consider spatially correlated functional models, where the functional models correspond to the color fading time-series. The derivative process of the functional models to ensure monotonicity (non-decreasing) and long term stabilization of color fading (long term color fading [del_hoyo_2015]) as a function of time, is taken into account.

The color-fading time-series are modeled as penalized splines-based models. The derivative of a semiparametric model such as a splines model is still a linear model and can be used as a additional observation jointly with the regular observations. Virtual derivative observations of both the values of the partial derivatives and the sign of the partial derivatives are used to induce monotonicity and long-term saturation, respectively.

The spatial correlation of the time-series is modeled by defining multivariate Gaussian process priors distributions over the splines coefficients across time-series. Color fading is mainly related to physicochemical properties of the measured surface. Physicochemical data for all the points on the surface to predict new time-series are hard to obtain. Instead, image color values can be used as input variables to construct and evaluate the correlation among time-series, since these image color variables are related in some way to physicochemical properties of the imaged surface [Malacara_2011]. A multivariate covariance function in a Gaussian process prior allows for using trichromatic image color variables jointly with spatial locations as inputs to evaluate the covariance structure among time-series.

Finally, in order to do model evaluation and comparison, the same model but without derivative information is also fitted. Cross-validation procedures are conducted to compute the posterior predictive checks, the expected log predictive density and the mean square predictive errors in order to do model checking and selection and evaluate the predictive performance.

The rest of the paper is structured as follows. Section 2 describes the available data in detail. Section 3 focuses on the modeling and inference formulation, as well as model checking and model selection procedures. Section 4 describes the results of applying the proposed model to the data set. Section 5 discusses about the results. Finally, Section 6 presents a brief conclusion of the work.

## 2 Data description

The MFS technique allows for real time monitoring of spectral changes of a small area (diameter of 0.5 mm) of a cultural heritage object by exposing it to light. This method was developed by Whitmore_1999, Whitmore_2000 and has been extensively used to study the lightfastness of cultural objects [Ford_2011, Ford_2013, Columbia_2013].

In this practical case we seek to evaluate and document the degree of color fading over time and space on rock art paintings caused by direct solar irradiation. Each measured point on the surface gives rise to a time-series that represents color fading (Color differences/Delta E* value, Malacara_2011) over time. The MFS measurements have a duration of 10 minutes. The sampling frequency will be once per minute, such that the resulting time-series will consist of 11 () time points . Due to these measurements are largely time-consuming and difficult to materialize, specially in these rock art systems, the MFS dataset available consisted of 13 observed locations on the surface () only. Figure 1-left shows their pixel locations on a color image of the study area; each location incorporates color fading time-series of observations (Figure 1-right). The rock art paintings study area is located in the Cova Remígia rock-shelter, Castellón (Spain).

The number of available input variables to evaluate the covariance between spatial locations is 5 (), which are the three image color variables, Hue (), Saturation () and Intensity (), and the two spatial coordinates and . Table 1

contains summary statistics for these input variables, which have been rescaled by dividing by their standard deviations. The

and spatial coordinates were divided by their common standard deviation.Each time-series is modeled as a spline-based function. The order or number of spline knots of the spline-based functions model is set to 3 () and placed uniformly through the time points variable.

## 3 Correlated time-series with derivative information

#### Penalized splines models

We consider a continuous stochastic process based on time-series of observations observed at spatial locations on a surface. Thus, denotes the th value of the time-series at the th location, with representing the time points and representing the spatial locations. The time-varying dimension of the data is modeled using basis functions models, so the

-vector of values of a time-series

,, is considered to follow a Normal distribution depending on an underlying penalized quadratic-splines-based function

and noise variance

. This can be expressed as:(1) |

In the previous equation, is the linear part of the model for the time-series , where is the matrix of linear functional values with the th row , and are the linear coefficients for the time-series . On the other hand, is the non-linear or splines-based part of the model for time-series , where are the splines coefficients for time-series , and the order of the splines model and number of knots, and is the matrix of penalized quadratic-spline functions values:

(2) |

where is the matrix of the quadratic-spline functions values, and the matrix of penalization of Crainiceanu_2005. An element of the matrix is

(3) |

and of the matrix is , with and are the th and th pre-fixed knots corresponding to the th and th spline function, respectively.

The observational model for the matrix of observations takes the form

with

(4) | |||

(5) |

In the previous equations is the matrix of linear coefficients with the th column , is the matrix of splines coefficients with th column . and are the matrices of linear and splines function values, respectively, and and are the th rows of the matrices and , respectively.

#### Correlating the splines models

In order to establish a correlation structure among time-series , the matrix of spline coefficients is considered as a realization of a continuous stochastic process, modeled as a (zero-mean) Gaussian process prior.

With the aim of simplify the covariance structure, null covariances between splines coefficients belonging to different spline knots can be considered, that is equivalent to define independent Gaussian process priors, one for each th row of matrix , i.e.,

for . is the covariance function for the vector of coefficients

is the standard deviation of the Gaussian process which controls the overall scale or magnitude of the range of values of the vector . Thus, specific covariance structure for each th vector of splines coefficients can be specified.Furthermore, in this work we will simplify even more the covariance structure considering the same spatial structure for every vector of splines coefficients , i.e.,

(6) |

for , and is a common covariance function for every vector of coefficients .

The covariance matrix is computed by a squared exponential covariance function [Rasmussen_2006], dependent on the matrix of input variables, where is the number of inputs variables, and the vector of lengthscale hyperparameters ,

where . The vector contains the lengthscale parameters for the input variables. For the input variables , and correspond the parameters , and , respectively. The spatial coordinates input variables and are sharing the lengthscale parameter (), such that the covariance function depends on the (Euclidean) distance between spatial coordinates. Hyperparameter controls the smoothness of the covariance function or rate of decay of the correlation in the direction of the th predictor, so that, the larger the the smoother the correlation function and the smoother the posterior functions for ; Although the scale of is dependent on the scale of the input data along the dimension . The squared exponential covariance function prior is a stationary function with respect to the input variables and implies independence among the contributing effects of the different input variables.

#### Function value constraints (zero-order constraints)

The modeling procedure has to ensure that the function values are zero at the subset of starting time points, that is, the time-series must start in zero. The function values in (5) evaluated at the subset results

where and correspond to matrices and evaluated at the subset of observations, respectively.

This property can be specified by using virtual observations equal to zero at these points, , and the Dirac Delta function as an observational model for these observations,

(7) | |||||

While the rest of observations are considered to be contaminated with Gaussian noise ,

(8) | |||||

where denote the dataset without the subset of observations, . And denotes the latent function values without the subset of observations, .

#### Splines functional models with derivatives

Color degradation is expected to be non-decreasing as a function of time and stabilize in the long term, that is, the time-series must be monotonically non-decreasing and stabilize eventually. These properties can be expressed in terms of the first order partial derivative of the penalized splines model for each time-series. The partial derivative of the penalized splines function of time-series , in (1), with respect to the time input variable takes the form

(9) |

where is the partial derivative with respect to time of matrix of penalized spline function values (2),

In the previous equation, the element of the partial derivative of the matrix is the partial derivative of the element (3),

#### Derivative constraints (first-order constraints)

The derivative of function in (5) is:

(10) |

In order to impose a saturation constraint for long term stabilization, virtual observations of the values of partial derivatives with respect to the time input variable equal to zero are considered. Most of the rock art painting systems analyzed so far show stabilization at or before the 10 minutes of MFS monitoring measurements. Therefore, virtual observations of the partial derivatives equal to zero are considered at the subset of ending points of the set of time-series,

where and denotes the partial derivative and function values, respectively, at the subset of points. A Dirac Delta function is considered as an observational model for these observations,

(11) | |||||

where is the partial derivative of the functional models in (9) at the subset points .

On the other hand, the time-series are guaranteed to be non-decreasing as a function of time when the partial derivative of the spline function is positive. So, virtual observations of the sign of the partial derivatives equal to one at the subset of desired time points where to induce monotonicity are considered,

The probit function can be used as a likelihood for the signs of the partial derivatives [Riihimaki_2011],

(12) |

The function

is the standard Normal cumulative distribution function and

is a parameter controlling the strictness of the constraint. When approaches zero (), function approaches the step function. We set the parameter equal to .#### Likelihood function

The joint likelihood of both regular observations (noise observations (8) and noise-free observations (7)) and derivative observations for monotonicity (12) and for long-term stabilization (11), given the parameters , and , the function values and , and the derivatives function values , results:

(13) | |||||

#### Bayesian inference

Bayesian inference is based on the posterior joint distribution of parameters given the data, which is proportional to the product of the likelihood and priors:

where is the likelihood of the model (13) and is the Gaussian process prior for the parameters (6). The vector contains the hyperparameters and . We set normal prior distributions for the linear parameters , , and positive half-normal prior distributions for the hyperparameters , , and ,

, and gamma distributions for the hyperparameters

, for all . These correspond to weakly informative prior distributions based on prior knowledge about the magnitude of the parameters.Predictive inference of new output values and for a new sequence of input values can be computed by sampling from the posterior joint predictive distribution .

To posterior distribution of interest

is in general intractable. Hence, to estimate both the parameter posterior distribution and the posterior predictive distribution for this model, simulation and/or distributional approximations methods must be used. Simulation methods based on Markov chain Monte Carlo

[brooks_2011] are general sampling methods to obtain samples from the joint posterior distribution.The proposed model considers the same spatial covariance structure between splines coefficients belonging to different order of the splines knots. This requires computation expense in covariance matrix inversion, where denotes the number of spatial locations.

For large data sets where iterative simulation algorithms are too slow, modal and distributional approximation methods can be as efficient and approximate alternatives [Gelman_2013].

### 3.1 Model checking, predictive performance and model selection

Common procedures of checking normality and tendencies on the fitted residuals can be used. However, the posterior predictive checks, which are also known as the

leave-one-out probability integral transformation

(LOO-PIT), can be used as a more rigorous procedure in order to guarantee good model performance and ensure that the model is compatible with the observed data. They are based on computing the probability of new predictions to be lower or equal to their corresponding actual observations following a leave-one-out cross-validation procedure [gelfand_1992, Gelman_2013]:where is the subset of observation indices of new predictions in the cross-validation. are the new observations from the predictive distribution at the subset , and are the actual observations at this subset

. The similarity or provenance of these probabilities from standard uniform distributions endorses these probabilities with the desirable property of having the same interpretation across models, which implies good fit to the data and good prediction

[Bayarri_2000]. Using simulation methods for estimating and predicting a Bayesian model, computing the probability of a predicted value being minor the observed one is straightforward through the collection of simulated values.The expected log predictive density (ELPD) evaluates, by averaging over all the steps in the leave-one-out cross-validation procedure, how far new data is from the model while taking the posterior uncertainties into account. It is based on the log-density of new data given the model [vehtari_2012]:

where denotes the cardinality of the subset of observation indices in the cross-validation. denotes the dataset without the subset of observations, .

The mean square predictive error (MSE) also evaluates, by averaging over all the steps in the leave-one-out cross-validation procedure, how far new data is from the model by using the distance (error) between the actual observation () and the predictive mean ():

Following a leave one observation out cross-validation scheme (CV1), the subset of observation indices will be just a single observation , . The LOO-PIT, ELPD, and MSE will be computed following CV1. The LOO-PIT will essentially be useful for model checking, ensuring that the model is compatible with the data. The ELPD and MSE will evaluate the predictive performance of individual observations .

The end goal of this work is to predict complete color-fading time-series at new unobserved locations. In order to do that, a leave one location out cross-validation scheme (CV2) can be performed, where the subset of observation indices will be a complete time-series of a specific spatial location , . The statistics ELPD and MSE will be computed following CV2. Plots of predicted new time-series superimposed to their corresponding actual observations will be shown in order to visually evaluate the predictive performance. Model selection can be done by comparing the predictive performance between models using the ELPD and MSE statistics. The best model is who maximizes the ELPD and/or minimizes the MSE.

## 4 Experimental results

The posterior distributions and predictive distributions have been estimated by Hamiltonian Monte Carlo sampling methods [neal_2011] using the Stan software [Carpenter_2017]. Three simulation chains with different initial values have been launched. The convergence of the simulation chains was evaluated by the split-Rhat convergence diagnosis [Gelman_2013] and the effective sample size of the chains. A value of 1 in the split-Rhat convergence statistic indicates convergence of the simulation chain, although conventionally accepted values of convergence would be between 1 and 1.1. In our case, a split-Rhat value lower than 1.05 has been obtained for all parameter simulation chains.

The estimated marginal posterior distributions for parameters are shown in Figure 2. is the overall scale of the Gaussian process prior for the vectors of splines coefficients . is the noise of the observations, and are the lengthscale parameters associated with the input variables of the squared exponential covariance function in the Gaussian process priors. As mentioned above, the spatial coordinates input variables and are sharing the lengthscale parameter (), such that the covariance function depends on the (Euclidean) distance between spatial coordinates. The lengthscale parameters , and are associated to the input variables , and , respectively.

Figure 3

shows the means and 95% pointwise credible intervals of predictive distributions,

, evaluated over the observed data and plotted versus the time points and superimposed to the observed MFS data . Additionally, the means of the predictive distribution of the model without the inclusion of the derivatives, , are also plotted for comparison. In Figure 4, the means and 95% pointwise credible intervals of the posterior distributions of the functions and , and , are plotted.Figure 5 shows the frequency histograms of the posterior predictive checks (LOO-PIT) by following the cross-validation scheme CV1, for both models, with and without derivatives.

Figure 6 shows predictive distributions of the proposed model with derivatives, , and of the model without derivatives, , for new time-series by following the cross-validation scheme CV2. In this way, predictive performance of complete new time series for unobserved locations is evaluated. The actual data of the time-series are also plotted to visually evaluate the predictions.

Table 2 shows the ELPD and MSE statistics computed by following the two different cross-validation schemes, CV1 and CV2, and for both models, i.e. with and without derivatives.

In order to make spatial continuous maps of color fading estimates, predictions of color fading time-series, , have been computed for all the spatial pixel locations of the rock art painting image (Figure 1). In Figure 7, six images representing the spatial distribution and their evolution over time of those color fading estimates are shown. The images correspond to the time points , , , and , respectively. Color changes higher that aproximately 3.5 (deltaE*) are considered to be perceptible changes for the human eye [Malacara_2011].

## 5 Discussion

The input variables , , and the spatial coordinates (, ) have been previously standardized. The fact that the lengthscale parameters and , corresponding to the variables and , respectively, are relatively small (the mode of is 1.1 and the mode of is 0.8) indicates that the function is non-linear with those variables or that the rate of decay of the correlation is high so that variations in the input variables imply a quick decrease in the correlations allowing for the non-linear effects. The variables and spatial coordinates have larger lengthscales (modes of 20.0 and 8.9, respectively) so that the function depends on and spatial distances in a less non-linear way. Especially, the variable with a lengthscale of 20.0 contributes with a constant effect of one to the correlation, and implies a constant function to the Gaussian process, so being an irrelevant variable to the functions.

Modeling and smoothing over the temporal dimension of the data has been achieved using penalized splines of order 3, as can be seen in Figure 3.

The use of a noise-free pseudo observation model at the starting time points, , forces the predictive distributions to be zero at those starting points (Figures 3 and 6).

The inclusion of observations of the sign of the partial derivatives for all the data points, , has induced monotonicity (non-decreasing) through the whole time-series. In the same way, the inclusion of observations of the value of the partial derivative equal to zero at the ending time points, , and a noise-free pseudo-observation model for these observations, has induced to reach a stationary state at the ending of the time-series. These constraints have to be considered either in the observed and predicting data points.

In the Figure 4-right, the derivative function values of the process for the observed data are plotted. Derivatives are always positive, approaching zero at the ending time points, which means the time-series are always non-decreasing and stabilizing in the long term. Derivative functions are linear, as can be seen in those plots, since functions have been modeled as quadratic splines-based functions.

Monotonicity and long term stabilization properties of the time-series were not ensured using the model without derivative information (Figures 4-left and 6), hence, the proposed model with derivative information yields a better fit for the functions dynamics, improving their interpretability. In this sense, the analysis of the color fading time-series using a model without derivative information could not be properly done because the temporal degradation, specially at the last time points, is unrealistic.

The LOO probability integral transformation (LOO-PIT) for both models, with and without derivatives, have been computed following the leave-one-out cross-validation scheme CV1. The frequency histograms of these LOO-PIT shows similarities to uniform distributions for both models, which means good model performances and that the models are compatible with the observed data.

The cross-validation scheme CV2 based on leaving a whole time-series out of the training data, has been carried out in order to evaluate the prediction performance of complete new time-series at new locations. Figure 6 shows the predicted time-series at new locations following the cross-validation scheme CV2, using both models, with and without derivatives. Better predictive performance is appreciated for the model with derivatives. Predictions are closer to the actual data and predictive intervals are narrower, because monotonicity and saturation constraints improves and reduces the credible intervals of the predictions. The model without derivatives shows decreasing patterns on the functions which are not consistent with the prior knowledge. Credible intervals, especially at the last part of the splines-based functions (time-series), are getting much bigger. Here, it can be seen as imposing the monotonicity and saturation constraints on the splines functions improve considerably the credible intervals.

The ELPD and MSE statistics are mostly better for the model with derivatives, maximizing the ELPD and minimizing the MSE (Table 2), in both cross-validation schemes. Therefore, the results of these statistics confirm that the model with derivatives is closer to new data, either in terms of the expected log-density or the mean error.

Prediction sensitivity due to the short set of data available have been found. It can be seen in Figure 5 poor predictive performance when compared to the observed data. This is due to the high sensitivite of the model to leaving some data out.

The spatial distribution and evolution over time, of color fading estimates for an entire area is shown in Figure 7. The spatial distribution over time of color fading values seems to be quite unvarying. This was expected since the time-series adjusted in the different locations have similar patterns, smooth, monotonic increasing and tending to stabilize in the long term (Figure 4). Color fading values, specially when they are low values, like in this case (the maximun is of 7.64), are not worth being converted to the RGB color space and plotted as an image, because the color changes will be not visible in a RGB image. The best way is plotting color deltaE* values (Figure 7). The science of colorimetry argues that deltaE* values higher that aproximately 3.5 would be perceptible for the human eye looking at the real object.

The actual equivalency of the time points () used in MFS measurements in years depends on the hours and intensity of sunlight that affects the paintings on a changing everyday basis. Without proper monitoring of light, this equivalency is difficult to obtain. Although this aspect of the research was not considered in the current study, future work will include an evaluation of the location and geographical orientation of the paintings together with long-term monitoring of light and UV radiation with the aim of estimating the dose acting upon the paintings in years.

## 6 Conclusion

Color is an important aspect of documentation and conservation of historic and cultural heritage, such as rock art paintings. The MFS data are largely time-consuming and difficult to materialize, specially in these rock art systems, so an interpolation procedure is needed in order to make predictions of MFS data on unobserved locations. Furthermore, these measurements in these systems are contaminated with large fluctuations, so the consideration of constraints in the modeling in order to overcome possible modeling issues that may arise due to these large fluctuations is highly encouraged.

We have formulated a model for correlating the MFS time-series of observations, with the addition of constraints related to their derivatives in order to fit some desired properties and minimize the effects of largely fluctuations in the original observations. Penalized splines of order 3 has been used for modeling the time-series. The time-series have been correlated through the spline coefficients using Gaussian processes models with multivariate predictors. In the practical case, we have been able to include colorimetric variables and the spatial coordinates variables in the covariance function, and demonstrated the colorimetric values being useful to correlate the MFS data. The contribution of the Euclidean distance (through the spatial coordinates variables) on this covariance structure is found to be quite weak.

Considering noise-free pseudo observations at the starting time points has allowed the functions to start at zero. The use of first order constraints related to the derivative of the functions have allowed the functions to fulfill the properties of being non-decreasing and stabilizing in the long term, and demonstrated better either in terms of predictive performance or application-specific interpretability.

The computation requirements to invert the covariance matrix in the proposed correlated functional models have been reduced substantially compared to using a Gaussian process models with a spatio-temporal covariance function. The proposed general correlated functional models framework, where a complete covariance structure among splines coefficients is considered, requires computational demand in covariance matrix inversion. Whereas spatio-temporal Gaussian processes require . Notice that the number of knots in spline-based models is usually much lower than the number of time points (). In case a null covariance is considered between splines coefficients belonging to different spline knots, the computational expenses of the proposed model becomes . And furthermore, if the same spatial structure is considered for the splines coefficients belonging to different spline knots, as we implemented in the present paper, the computational expenses of the proposed model becomes .

Reliable color fading estimates evolution maps can be elaborated by means of using the proposed model with derivative information in comparison to the model without derivative information.

Finally, multivariate covariance metrics and zero and first order constraints might be very hard to implement outside of a Bayesian framework and Gaussian process models. Gaussian process is flexible enough which allowed us to properly model this complex covariance structure of the time-series dependent on different input covariates. The Bayesian framework has allowed us to jointly use normal distributed observations with probit distributed observations of the sign of the partial derivatives, allowing to fulfill the determinants on the behavior of the functions.

## Acknowledgments

The authors gratefully acknowledge the support from the Spanish Ministerio de Economía y Competitividad to the projects HAR2014-59873-R and MTM2016-77501-P, as well as the grant PPIC-2014-001 funded by Consejería de Educación, Cultura y Deportes (JCCM) and FEDER. The authors want to express their gratitude to the General Directorate of Culture and Heritage, Conselleria d’Educació, Investigació, Cultura i Esport, Generalitat Valenciana for letting us access and carry out research at the archaeological site.

Comments

There are no comments yet.