Robust Function-on-Function Regression

Functional linear regression is a widely used approach to model functional responses with respect to functional inputs. However, classical functional linear regression models can be severely affected by outliers. We therefore introduce a Fisher-consistent robust functional linear regression model that is able to effectively fit data in the presence of outliers. The model is built using robust functional principal component and least squares regression estimators. The performance of the functional linear regression model depends on the number of principal components used. We therefore introduce a consistent robust model selection procedure to choose the number of principal components. Our robust functional linear regression model can be used alongside an outlier detection procedure to effectively identify abnormal functional responses. A simulation study shows our method is able to effectively capture the regression behaviour in the presence of outliers, and is able to find the outliers with high accuracy. We demonstrate the usefulness of our method on jet engine sensor data. We identify outliers that would not be found if the functional responses were modelled independently of the functional input, or using non-robust methods.


Introduction
Functional Linear Regression (FLR) in the function-on-function case (Ramsay and Dalzell, 1991) is a widely used technique for modelling functional responses with respect to functional inputs.The FLR model is able to capture complex dependency structures as it uses information across time (Morris, 2015).However classical FLR models can be severely affected by outliers as we will demonstrate via a simulation study in Section 6.We therefore develop a robust FLR (RFLR) model, which is able to effectively fit the data in the presence of outliers.The model is built using the robust Functional Principal Component model by Bali et al. (2011) and the multivariate Least Trimmed Squares (MLTS) estimator by Agulló et al. (2008).The RFLR model can be used to identify abnormal functional responses, i.e. samples in which the functional behaviour between the predictor and response curves deviate from normal.
Our study of robust FLR is motivated by the need to identify normal relationships in jet engine sensor data when we expect outliers to be present.The data is collected during Pass-Off tests, which are performed on an engine before deployment.In a Pass-Off test a human controller performs manoeuvres, which can be defined as various engine accelerations and decelerations starting and ending at a set idle speed.During the test, data is captured by sensors measuring engine speed, pressure, temperature and vibration in different parts of the engine.One of the key manoeuvres in a Pass-Off test is the Vibration Survey (VS).In this manoeuvre the engine is accelerated slowly to a certain speed then slowly decelerated.
We have 199 VS datasets, which include the turbine pressure ratio (TPR) that measures the engine speed, and various temperature features including the turbine gas temperature (TGT).In Figure 1 we have plots of the TPR and TGT for the first 30 VS manoeuvres.
To anonymise the data we have transformed the time index onto the interval [0, 1] and the sensor measurements to the range [0, 100].
The VS manoeuvres are performed by a human controller, which causes variability in the TPR curves as can be seen in Figure 1.This variability will naturally affect the TGT curves and may mask the unusual behaviour produced by the engine.Our results in Section 7 support this claim, and show that using a direct outlier detection on the engine temperature curves fails to identify meaningful outliers.Instead this approach picks up curves produced by unusual TPR speed profiles.We therefore require a method of detecting outliers in the presence of the controller induced variability.We expect that the relationship between the engine speed and engine temperature for different VS manoeuvres should be the same irrespective of the way the manoeuvre is performed.For example given a certain engine acceleration we would expect a certain temperature response.If however the response differs from expectation this could be indicative of an engine issue.In Section 5 we will show how RFLR can be used for outlier detection, which is later applied to the jet engine data in Section 7, to identify abnormal behaviour.
The paper is organised as follows.In Section 2 we outline the classical FLR model.In Section 3, we will outline robust Functional Data Analysis (FDA) techniques to obtain a robust FLR model.We also introduce a robust model selection procedure.In Section 4 we prove consistency results for the robust FLR model and the robust model selection procedure.In Section 5, we describe an outlier detection method, which acts on the residuals of the robust FLR model.In Section 6 we perform a simulation study to illustrate the model fitting and outlier detection capabilities of the robust model.In Section 7 we apply the robust model on the engine data and highlight unusual observations that can not be detected by using outlier detection directly on the temperature curves.Finally in Section 8 we provide a conclusion.

Classical Functional Data Analysis
In this section we give a brief summary of the FDA tools that we will later apply in our model.In the following sections we will use the vector space L 2 (I) which is the Hilbert space of square integrable functions on the compact interval I with the inner product f, g = I f (t)g(t)dt for functions f, g ∈ L 2 (I).
We will define X(t), Y (t) to be univariate stochastic processes defined on I, with mean functions µ X (t) and µ Y (t), and covariance functions C X (s, t) = cov{X(s), X(t)} and C Y (s, t) = cov{Y (s), Y (t)} for all s, t ∈ I.We shall define x(t) = [x 1 (t), ..., x n (t)] and y(t) = [y 1 (t), ..., y n (t)] to be n independent and identically distributed realisations of X(t) and Y (t) respectively.In practice we observe x i (t) and y i (t) at discrete time points.We shall assume for simplicity of exposition that observations are made at equally spaced time points t 1 , ..., t T .
We will outline Functional Linear Regression and Functional Principal Component Analysis with respect to the underlying functions x(t), y(t).In Section 2.3 we need to use the discretely observed data to define a suitable model selection criterion.

Functional Linear Regression
In this section we will introduce the classical FLR model (Ramsay and Dalzell, 1991).In FLR we model the relationship between predictor x i (t) and response y i (t) as: where α(t) is the intercept function, β(s, t) is the regression function and i (t) is the error process.For a fixed t, we can think of β(s, t) as the relative weight placed on x i (s) to predict y i (t).As in Chiou et al. (2016) we will assume the mean functions µ X (t) = 0 and µ Y (t) = 0 which thereby means α(t) = 0.This is a reasonable assumption as in practice we can calculate the mean functions µ X (t) and µ Y (t) efficiently for dense data and then pre-process the data by subtracting µ X (t) and µ Y (t) from the observed curves.
FLR in the function-on-function case can be modelled parametrically (Yao et al., 2005;Chiou et al., 2016) or nonparametrically (Ferraty et al., 2012;Ivanescu et al., 2015;Scheipl et al., 2015).We use a parametric approach which models the regression matrix in terms of pre-defined basis functions.
We will represent x i (t) and y i (t) in terms of (M, K) pre-chosen basis functions φ X j (t), φ Y j (t) respectively: where z im , w ik ∈ R.
We define φ .., z iM ] and . We will then model the regression surface using a double basis expansion (Ramsay and Silverman, 2005): for an M × K regression matrix B. We can then write: (2.3) Letting i (t) = q i φ Y (t) we can reduce Equation (2.3) to: This parametrisation of the residual function is also used by Chiou et al. (2016).We can then estimate B using standard multivariate regression methods typically assuming Gaussian q i .

Functional Principal Component Analysis
In this section we describe Functional Principal Component Analysis (FPCA), which we will use to build data-driven basis functions φ X (t) and φ Y (t) for x i (t) and y i (t), respectively.
These basis functions give effective, low-dimensional representations and will be used in the Functional Linear Regression model described in Section 2.1.
Functional Principal Component Analysis (FPCA) is a method of finding dominant modes of variance for functional data.These dominant modes of variance are called the Functional Principal Components (FPCs).FPCA is also used as a dimensionality reduction tool, as a set of observed curves can be effectively approximated by a linear combination of a small set of FPCs.
The FPCs, φ X m (t) for m = 1, 2, ..., are the eigenfunctions of the covariance function C X (s, t) with eigenvalues λ X m .Note that the eigenfunctions are ordered by the respective eigenvalues.The Karhunen-Loéve theorem (Shang, 2014) shows that x i (t) can be decomposed as where the principal component score The scores z im are realisations from a random variable ξ X m .
We can define the M -truncation as which gives the minimal residual error: To choose M we will use an information criterion outlined in Section 2.3.An analogous procedure is used to find K eigenfunctions φ Y k (t) for y(t).

Bayesian Information Criterion for FLR
In the FLR model described in Section 2.1 we need to choose terms M and K. Typically M and K are chosen independently (Yao et al., 2005), however the estimation of β(s, t) also depends on M and K and this should be incorporated into the estimation of these terms.In this section we formulate a Bayesian Information Criterion (BIC) to determine the basis size M and K, similarly to Matsui (2017).
A component of the BIC is the log likelihood, often expressed as a squared error term.
It is tempting to use the squared error resulting from Equation (2.4).However the objective is to fit the data y i so we should use a likelihood of this data instead of a squared error term of basis coefficients.
We have a set of models J = {(M, K)|M = 1, ..., M max , K = 1, ..., K max }, where M max and K max are pre-set maximum number of FPCs that will be considered in the model.Let vector y i be the values of y i (t) evaluated at discrete time points: be the first M principal scores of x i (t) with respect to the FPCs φ X (t) and let φ (K)   be the matrix with (k, i) entry φ Y k (t i ).We assume there exists a true model (M 0 , K 0 ) with associated M 0 × K 0 matrix B M 0 ,K 0 such that where the error i = [ i (t 1 ), ..., i (t T )] is assumed for simplicity to be sampled from N (0, v 2 I T ), where I T is the identity matrix of size T .
For Model (M, K) we define θ M,K = (B M,K , v M,K ) and the prediction ŷM,K We want to identify this true model (M 0 , K 0 ), which we can use to obtain consistent estimates of θ M 0 ,K 0 .
For Model (M, K) we can define the likelihood for sample i as and the log-likelihood l(θ As in Eilers and Marx (1996) BIC where θM,K is the maximum likelihood estimator and the penalty ω(M, K) = M K + 1, in which M K is the number of free parameters in the model and the 1 comes from v. We will To summarise, we estimate the FPCs for X and Y and solve the FLR model for different models (M, K).We then choose model (M * , K * ) n that minimises the BIC criterion.The robust equivalent of this procedure is given in Algorithm 1.

Robust Functional Linear Regression
In Section 2 we have defined the FLR model and have outlined the use of FPCA bases to estimate parameters of the model.In this section we will introduce robust versions of the FDA techniques outlined in Section 2. This will allow us to fit a normality model even in the presence of outliers.We shall also propose a robust BIC procedure for model selection.
We will replace classical FPCA with robust FPCA estimates by Bali et al. (2011) which ensure that outliers do not unduly affect the FPCA estimates.Note that FPCA minimises the residual error given in (2.6).To obtain robust FPCA estimates Bali et al. (2011) minimise a robust scale estimator, using a projection pursuit approach, which iteratively performs a weighted least squares till the estimators stabilise.
We define ỹi (t) = wi φY (t) and assume as in (2.4) that i = qi φY (s).We can now write the robust counterpart of (2.4) as wi = zi B + qi .This is robust as outliers will not be in the subset by definition so shall not affect the model estimation.We will choose a subset of size r = [0.8n].

Robust Bayesian Information Criterion for FLR
The BIC model selection method is known to be non-robust (Machado, 1993).In particular outliers can significantly affect the loglikelihood estimation.We therefore outline a robust BIC (RBIC) model, which, similar to MLTS, maximises over a subset of samples S. RBIC can therefore give good model selection performance in the presence of outliers.
We will define θM,K = ( BM,K , ṽM,K ) as robust estimated parameters for model (M, K) and the robust prediction ỹM, ) T BM,K φ(K) .We define the trimmed likelihood for model (M, K) and set S as We will define S M,K = arg min S∈S l( θM,K , S), where We will denote ( M , K) n = arg min (M,K)∈J RBIC n (M, K), and we will assume that this minimum is unique.
In Algorithm 1 we outline the calculation of the robust FLR model, which incorporates the RBIC procedure.In the algorithm we estimate the model for different values of (M, K) and choose the model with the minimum RBIC value.We consider M = 1, . . ., M max and l = 1, ..., K max where M max , K max are chosen to ensure that 99.99% of the variance in the raw data is captured.

Asymptotic Results
In Section 3 we proposed a Robust FLR model for the function-on-function problem.A minimum criteria for a good model is consistency, i.e. that given an ideal scenario of unlimited data that the estimator will be equal or arbitrarily close to the truth.In this Data: Let (x i , y i ) be mean-corrected time series of length T for i = 1, ..., n.
Algorithm 1: Robust FLR procedure section we shall prove consistency and Fisher-consistency for the robust FLR model.We shall follow a similar approach to Kalogridis and Aelst (2019) who developed a robust FLR model for the scalar-on-function problem.We shall also prove the consistency of the RBIC model selection method outlined in Section 2.3.
Definition 1.Let X 1 , X 2 , ..., X n be a sequence of real-valued random variables.An estimator T n := T (X 1 , X 2 , ..., X n ) of a parameter θ is said to be (asymptotically) consistent if for all > 0 Definition 2. Let X 1 , X 2 , ..., X n be a sequence of real-valued random variables with an associated cumulative distribution function F θ , which depends on an unknown parameter θ.
Let the estimator T n := T (F n ) of a parameter θ, be a function of the empirical distribution function F n .We say this estimator is Fisher-consistent for the parameter θ if Remark 1. Fisher consistency is equivalent to (asymptotic) consistency if the empirical distribution function F n converges pointwise to the true distribution function F θ .This can be shown to be the case for iid real multivariate random variables using the Glivenko-Cantelli theorem (Pollard, 2012).

Consistency of the Robust FLR
To prove Fisher-consistency we need to define appropriate probability measures on the predictor X(t), response Y (t) and the residual (t).We will then define conditions by which the robust FPCA and MLTS regression are Fisher-consistent, which will then ensure the Fisher-consistency of β(s, t).We shall also prove consistency of β(s, t) using Remark 1.Following the ideas set by Kalogridis and Aelst (2019), we make 6 assumptions: (C1) X has a finite-dimensional Karhunen-Loéve decomposition, i.e λ X m = 0 for m > M 0 .
(C3) The residual (t) = q φY (t) where q is a Gaussian random variable with mean 0 and covariance matrix Σ.
Let P X be the image measure of X i.e.P X (U ) = P (X ∈ U ) for a Borel set U , and likewise for P Y .We can define the cumulative distribution functions Let F denote the distribution function of (t), which can be defined in the same way as F X and F Y .We can write the functional of the robust estimator β(s, t) as: The functional is Fisher- Conditions C1-C4 are to ensure the FLR problem can be defined by a finite number of terms.Kalogridis and Aelst (2019) show that Conditions C5 and C6 are sufficient for the Fisher-consistency of the robust FPCA estimators by Bali et al. (2011).
Note that x i (t) and y i (t) are defined on a finite number of eigenfunctions, so are defined by finite score vectors.Therefore Corollary 4.1 follows from Lemma 4.1 and Remark 1, which states almost sure convergence of the empirical distribution for iid multivariate random variables.In this case Fisher-consistency is equivalent to consistency.

Consistency of RBIC
We defined RBIC for the FLR problem in Section 3.1.In this section we will prove consistency of RBIC for the FLR problem.We will assume there is a true model, which we previously defined as (M 0 , K 0 ).We can then define overspecified and underspecified models in reference to this true model.We make some assumptions on the behaviour of the likelihood for these two model classes to prove consistency.We also denoted ( M , K) n = min (M,K)∈J RBIC n (M, K), which we will assume is unique.
We will split the candidate models in J into two sets, one is the overspecified models that include the true model J + = {(M, K) ∈ J|M ≥ M 0 and K ≥ K 0 } and underspecified models J − = J c + ∩ J. Recall that r = [αn] for some α ∈ (0, 1), and the likelihood l in (3.2) depends on r terms.
Assumption 1 For (M, K) ∈ J − , there exists some ε M,K > 0 such that This is a reasonable assumption as the underspecified models should give a poorer fit to y i than the true model.
Assumption 2 For (M, K) ∈ J + , there exists some γ M,K > 0 such that This assumption states that the difference in the trimmed loglikelihood is less than a finite γ.The likelihood for the overspecified models and the true model should be close, given the true model is contained within the overspecified models, so the difference in the penalty terms will outweigh the difference in the likelihoods for large enough n.
Note that in Assumption 1 we consider the average difference between the log-likelihoods, whereas in Assumption 2 we look at the total difference.

r
. Using ε M,K from Assumption 1, we can see that −G r < 2ε M,K for sufficiently large r.Using this and Assumption 1 we can show For (M, K) ∈ J + \{(M 0 , K 0 )}, we know that 1 2 (ω(M, K) − ω(M 0 , K 0 )) log(r) > 0 and is monotonically increasing.Therefore there exists N such that for r ≥ N We can show that Note that BIC is a special case of RBIC where r = n, so is also consistent by Theorem 4.1.

Outlier Detection
There is a rich literature of outlier detection methods for functional data (FD).There are functional depth based methods such as the thresholding approach by Febrero-Bande et al. ( 2008) and the functional boxplot by Sun and Genton (2011).Alternatively we can use methods based on outlyingness measures such as Arribas-Gil and Romo (2014), and Dai and Genton (2018).For multivariate FD there exist outlier detection methods such as Rousseeuw et al. (2018) and Hubert et al. (2015).These methods do not model the dependency between the functional response and functional input, and may therefore miss important outliers.This will be shown in the simulation study in Section 6. RFLR can model this dependency structure, which can improve the detection of outliers.We therefore suggest an outlier detection algorithm which uses RFLR to model the dependency structure.
Using residuals from the model we can apply standard outlier detection approaches.The outliers in the residuals will be samples that are not well explained by the model which fits the majority of the curves.
For an outlier we expect the residual curve r i (t) = y i (t) − ỹi (t) to deviate in behaviour from the other residuals.Traditionally, we would use the integrated square error to identify outliers.However using a functional depth approach (Febrero-Bande et al., 2008) is more effective in identifying outliers in functional data, in particular shape outliers that are not unusual if viewed at each time point but are abnormal across the entire trajectory.The approach assigns a depth value to samples r i (t).Samples with small depth values lie far away from the other samples.
We will use the h-modal depth (Cuevas et al., 2007) to rank samples r i .For a given kernel G h (typically Gaussian with bandwidth h), the h-modal depth of r i with respect to r = {r 1 , ..., r n } is given by: (5.1) The h-modal depth has two useful properties.First, it uses a distance metric therefore samples further away from the centre will be given a smaller depth value.Second, in the case of multiple "normal" types behaviour, the h-modal depth works effectively as it doesn't assume there is one centre.
In the algorithm we need to choose the bandwidth h and a threshold C to identify outliers.The bandwidth h is taken to be the 15th percentile of the empirical distribution of {||r i − r j ||, i, j = 1, ..., n} (Febrero-Bande et al., 2008).The threshold C is chosen such that P (D(r i |r, h) ≤ C) = δ, where δ is a pre-chosen percentile.To estimate the threshold C they use a bootstrapping approach, which estimates a value of C for different random sets of samples and then aggregates these estimates.We describe the outlier detection algorithm in Algorithm 2.
1. Use Algorithm 1 to obtain φY k (t), zm and B. 2. Calculate residual curves r i (t). 3. Estimate bandwidth h. 4. For each r i (t) calculate D(r i |r, h). 5. Estimate C for given percentile δ. 6.If D(r i |r, h) < C sample i is an outlier.

Simulation Study
In this section we will provide a simulation study to investigate the finite sample properties of RBIC and robust FLR (RFLR) in comparison to BIC and classical FLR (CFLR).In the simulation study we will generate data using a FLR process and corrupt a certain number of samples, which will be the outliers.The outliers have been designed to be undetectable, if the response curves are considered independently of the predictor curves.Therefore standard functional data outlier detection algorithms such as those discussed in Section 5 will perform poorly.
The main motivation for the RFLR model is to obtain good model fitting in the presence of outliers.In this simulation study we compare the fitting error (FE) given in (6.1), for the non-outlier samples using the robust model, which uses RFLR and RBIC with the classical approach using CFLR and BIC.We define the indicator variable u i = 1 if sample i is an outlier and 0 otherwise.Letting ŷi (t) be the estimation of y i (t) and given that proportion a of the samples have been contaminated then FE is given by: Next we compare the outlier detection capabilities of robust and classical approaches using the receiver operating characteristic (ROC) curve to determine the sensitivity/specificity trade-off for different thresholds.If we have perfect outlier detection for all thresholds then the area under the curve (AUC) of the ROC curve would be 1.We can therefore use the AUC value as a measure of outlier detection accuracy regardless of threshold.
FPCA is performed by taking the principal components of a 200 cubic B-spline representation of each of the predictor and response curves (Ramsay and Silverman, 2005).
The robust FPCA approach outlined in Section 3 is performed using the CR algorithm proposed by Croux and Ruiz-Gazen (1996) on the same B-spline coefficients.The MLTS estimator is calculated using the heuristic given by Agulló et al. (2008) using different trimming proportions (1 − α) for α ∈ [0, 1].

Scenarios
We will generate samples x(t) using a FPCA based model with mean function µ X (t) = −10(t − 0.5) 2 + 2 for t ∈ [0, 1] and eigenfunctions: The principal scores are sampled from Gaussian distributions with mean 0 and variances 40, 10 and 1 for the eigenfunctions respectively.Note that we do not create any outliers in the FPCA decompositions of the predictor curves.We generate 400 predictor curves x 1 (t), ..., x 400 (t), which are observed at T = 500 equidistant points in the interval [0, 1].
The samples y(t) will have eigenfunctions: and mean function µ Y (t) = 60 exp(−(t − 1) 2 ).We will generate where B will have random entries between [−3, 3].We generate non-outlier curves: where the residual function i (t) = q i φ Y (t) + d i where q i and d i are sampled iid from N (0, 0.1).We will consider three cases when the proportion of outliers are a = 0.1, 0.2 and 0.3.
In Scenario 1 outliers will be generated by replacing B with B 1 = B + R where R has random entries sampled from N (0, 0.5) giving β 1 (s, t) = φ X (s) T B 1 φ Y (t).Outliers y i (t) are given by In Scenario 2 we generate outliers by adding a random B-spline function p(t) defined on an interval of length 1/10.Letting , l] for l ∼ N (2, 1), then the outliers y i (t) are given by Note that the outliers in Scenario 1 affect the regression function across the entire interval whereas the outliers in Scenario 2 only affect a small interval of the curves.
In Figure 2 we have a plot of the predictor curves x i (t) and response curves y i (t) with outliers from Scenario 1 and Scenario 2. The figure shows the outliers are masked by the variability in the curves and therefore cannot by identified using standard outlier detection algorithms.To make the outliers clearer we have plotted the residuals of the response curves using the true regression function and mean functions.In the bottom row of Figure 2 we can see that the outliers in Scenario 2 are localised to a fixed interval whereas in Scenario 1 the outliers affect the response curve at all time points.
The RFLR model depends on the proportion of trimming α.To investigate the effect of the trimming we will consider trimming proportions α = 0.1, 0.2 and 0.3.We shall also investigate the performance using BIC and RBIC with fixed trimmed sample size of r = [0.8n].
We sample 400 predictor and response curve datasets and generate classical and robust models to calculate the average FE (6.1).In Tables 1 and 2 we present the results for Scenario 1 and 2 respectively.The CFLR model gives a smaller FE value in the case of no-outliers a = 0, however the robust model still gives good model fits.If we compare the FE using BIC and RBIC, we can see that BIC gives better model choices when a = 0.This is due to BIC using all the data and in particular using samples in the tails of the distribution.In the presence of outliers the robust model outperforms the classical model, and as expected the difference in FE increases as the number of outliers increases.We should also note that RBIC is giving better model choices than BIC when outliers are present.Next, we can see using trimming proportion α = 0.1 we obtain significantly large FE values when a = 0.3.However the FE values for α = 0.2 and 0.3 are very similar in the (1)     i (t) for Scenario 2. The residual curves are generated using the true regression function and mean functions.In each scenario there are 5 outliers each in a distinctive colour.The predictors curves x i (t) are identical for both scenarios, and the response curves look very similar due to mean and functional components masking the outliers.However the residuals are clearly distinctive.Average AUC values over 100 replications for Scenario 2, using proportion of outliers a= 0.1, 0.2 and 0.3.Using Direct compared to classic FPCA with BIC, and using robust FPCA with RBIC and trimming levels α = 0.1, 0.2 and 0.3.points on the B-spline representations to be our inputs x i (t) and y i (t).
We will be applying the outlier detection algorithm described in Algorithm 2, which uses RFLR.We will compare these outliers with those detected on the temperature curves directly and using CFLR and BIC in Algorithm 2. We can look at the residuals curves to determine if the outliers do indeed look abnormal.In particular we want to show that using functional regression we are able to determine outliers that would otherwise be missed by investigating the temperature curves directly.
Using the depth based outlier detection (Direct) (Febrero-Bande et al., 2008) directly on the temperature curves (with a default threshold of δ = 0.01), we obtain the outliers in Table 5.We can see that the outliers in the TPR are the same as the outliers in the temperature features.This suggests the outliers being identified are arising from the controller induced variability.We therefore need to model the dependency between the control feature (TPR) and the temperature features.
We applied the outlier detection algorithm given in Algorithm 2 using CFLR and BIC with threshold δ = 0.01.The outliers identified are given in Table 5.The residuals curves are shown in Figure 5, with the outliers coloured in blue.It is not clear from this plot that the outliers are truly different from the other data.
Lastly we applied Algorithm 2 using RFLR and RBIC with threshold δ = 0.01.The outlier samples are given in Table 5 for each temperature feature.In Figure 6 we have the residual curves using RFLR.We can see that the RFLR model fits the majority of the temperature curves well.The outliers that are picked up clearly look abnormal, with significant deviations from the general behaviour.The RFLR model is therefore able to identify interesting behaviour, which may otherwise have been undetected.Engineers have informed us that Sample 24 comes from an engine in which they detected damaged hardware.All the other outliers in the RFLR column of Table 5 were also noted to come from engines that displayed odd behaviour during the Pass-Off test.This is not the case for the outliers reported in the CFLR column.
In Figure 4 we have a plot of the temperature parameters with the outliers identified using the curves directly in green, those using the RFLR model in red and those detected by both in purple.We can see that the outliers from the RFLR model do not necessarily appear as abnormal if we look at the temperature curves directly.Sample 106 is identified as an outlier by multiple temperature features and also when the depth based outlier detection is used on the temperature curves directly.Comparing the outliers identified using a classical approach, we can see Sample 24 is identified as an outlier multiple times using the classical and robust approaches.However most of the outliers from the classical approaches differ with the outliers identified using the robust approach.We can also see that the outliers using the RFLR are significantly more distinctive than the outliers using CFLR.

Conclusion
There exist a number of functional regression models for functional inputs and responses, however these methods are not robust to outliers.We have introduced a robust FLR model that is able to produce good model fits in the presence of outliers.Alongside the robust

Figure 1 :
Figure 1: Plots of the first 30 TPR and TGT time series.
robust estimate of the regression matrix B, we will use the Multivariate Least Trimmed Squares (MLTS) estimator by Agulló et al. (2008), to mitigate the affect of outliers with respect to the regression relationship.Given α ∈ [0, 1] we can define r = [αn] as the α proportion of samples rounded to the nearest integer, and the set S = {S ⊂ {1, ..., n}, |S| = r}.The objective of MLTS is to find a subset S such that S = arg min S∈S i∈S || wi − zi B|| 2 .

Figure 2 :
Figure 2: Left: Plots of the predictor curves x i (t), response curves y for Scenario 1. Right: Plots of the predictor curves x i (t), response curves y

Figure 3 :
Figure 3: ROC curve for one instance of Scenario 1 and 2 with the proportion of outlier a= and proportion trimmed α = 0.2.
. There are a number of temperature features measured within an engine including the TGT, discussed previously.In addition we have four other temperature readings T25, T30, TCAR and TCAF, from sensors measuring temperature in different parts of the engine.All the temperature features are shown in Figure4.The TCAR is particularly interesting as it has two distinct curve behaviours.It is also worth noting that the temperature values are distinctively higher at the end of the manoeuvre than at the beginning even though the engine speeds are the same.This highlights the trajectorydependent behaviour that we seek to model.The VS manoeuvres time series are of similar length.To standardise we have fitted a B-spline basis of 400 basis functions to each to ensure the time series are well approximated.Then we have taken 1000 equally spaced

Figure 4 :
Figure4: Plots of the TPR, T25, T30, TGT, TCAR and TCAF time series with outliers using robust FLR in red; those using the curves directly in green and those for both in purple.

Figure 5 :
Figure 5: Plots of the residuals of the T25, T30, TGT, TCAR and TCAF with outliers using classical FLR in blue.

Figure 6 :
Figure 6: Plots of the residuals of the T25, T30, TGT, TCAR and TCAF with outliers using robust FLR in red.

Table 1 :
Average fitting errors (FE) for 100 replications for Scenario 1, using classic FPCA and robust FPCA with different amount of trimming in the MLTS estimator and using models selected by BIC and RBIC.

Table 2 :
Average fitting errors (FE) for 100 replications for Scenario 2, using classic FPCA and robust FPCA with different amount of trimming in the MLTS estimator and using models selected by BIC and RBIC.