{'fraction_non_alphanumeric': 0.04208910352768094, 'fraction_numerical': 0.016630952761697636, 'mean_word_length': 4.565681737273051, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Robust joint modeling of sparsely observed paired functional data

A reduced-rank mixed effects model is developed for robust modeling of sparsely observed paired functional data. In this model, the curves for each functional variable are summarized using a few functional principal components, and the association of the two functional variables is modeled through the association of the principal component scores. Multivariate scale mixture of normal distributions is used to model the principal component scores and the measurement errors in order to handle outlying observations and achieve robust inference. The mean functions and principal component functions are modeled using splines and roughness penalties are applied to avoid overfitting. An EM algorithm is developed for computation of model fitting and prediction. A simulation study shows that the proposed method outperforms an existing method which is not designed for robust estimation. The effectiveness of the proposed method is illustrated in an application of fitting multi-band light curves of Type Ia supernovae.


Introduction
This paper is motivated by a need in astrostatistics to empirically quantify the shapes of Type Ia supernova light curves and to quantitatively relate the shape parameters with the intrinsic properties of Type Ia supernovae (SNeIa). SNeIa are "standardardizable" candles in cosmology, widely used to measure the expansion rate of the Universe. The evolution of SNeIa can be described by light curves which are measurements of brightness as a function of timeby measuring the photon flux through different astronomical filters. One of the remarkable features of SNeIa is their homogeneity nature and thereby, the luminosity evolution dispersion is low. However, SNeIa consist of a considerable amount of peculiarities with differing light curve shapes, adding more heterogeneity and complexity to the light curve modeling. Peculiar SNeIa include SN1991bg-like ones which are known to have faster decline rate after the maximum brightness, to show less or no sign of shoulder in band and to lack the secondary peak in band compared with normal SNeIa (Filippenko, Richmond, Branch, Gaskell, Herbst, Ford, Treffers, Matheson, Ho, Dey et al., 1992). The SN1991bg-like SNeIa are manually removed from the training sample in some existing light curve fitting models, e.g. Guy, Astier, Baumont, Hardin, Pain, Regnault, Basa, Carlberg, Conley, Fabbro et al. (2007), to mitigate the disturbing effect. He, Wang and Huang (2018) pointed out that SNeIa light curve observations can be viewed as observed functional data and developed a functional principal component analysis (FPCA) model to fit a collection of well-observed SNeIa samples accumulated in the past years. The principal component functions are useful for fitting SNeIa light curves in future supernova surveys. The goal of this work is to extend the previous work in two directions: 1. To jointly model observations from different filter bands. 2. To develop a robust method that is resistant to outlying observations. For simplicity of presentation, we focus on jointly modeling light curve data from two filter bands. Extension to multi-band is straightforward with some notational complications. Our solution approach is built upon Zhou, Huang and Carroll (2008), where a reduced-rank mixed effects model is constructed to model paired functional data, splines are used to model the principal component functions, and principal component scores are modeled as normally distributed random effects.
Due to the high dimensionality of functional data, the outliers of functional data have many different facets. They can be displayed as outlying measurements at some points, or as outlying shapes of the entire function. Lange, Little and Taylor (1989) suggested that the t distribution provides a powerful tool for handling outliers in a wide range of settings due to its robust property. Inspired by this work, we propose to use the t distribution to model the principle component scores and the measurement errors to obtain robustness of the methodology. More generally, our framework incorporates a broader class of distributions for modeling the principal component scores, i.e., scale mixture of normal (SMN) distributions (Andrews and Mallows, 1974), for which the t distribution is a special case. By treating the scale parameter in the hierarchical representation of SMN as a latent variable, we develop an EM algorithm to estimate parameters in our robust reduced-rank mixed effects model for paired functional data.
Several methods have been developed for robust modeling of single curve functional data. In a discussion paper, Locantore, Marron, Simpson, Tripoli, Zhang and Cohen (1999) proposed an approach that projects the data onto a sphere or an ellipse around a robust estimate of the center of the data and then performs the usual PCA on the projected data. See also the ten lively discussions of the paper and the rejoinder on improving the methodology and various aspects of robust functional data analysis. Gervini (2008) extended the work of Locantore et al. (1999), introduced the concepts of functional median and functional spherical principal components (PC), and studied the robustness properties of the approach. Hyndman and Ullah (2007), Shang (2009), andBali, Boente, Tyler andWang (2011) proposed a projection pursuit (PP) approach for robust functional PCA. Gervini (2009) used splines to model the principal component weight functions, and modeled the functional PC scores using heavy-tailed distribution to achieve robustness. Sawant, Billor and Shin (2012) used robust functional PCA for functional outlier detection. Boente and Salibian-Barrera (2015) proposed a new class of estimators for principal components based on robust scale estimators. While these works have focused on functional data of single curves, this paper develops a method for robust modeling of paired curves, which is new in functional data analysis. Our joint modeling approach allows strength borrowing through modeling correlation of the paired curves and thereby improves statistical efficiency. Treating the scale parameter of SMN as a latent variable in an EM algorithm for robust estimation is also a novel approach in the functional data literature. Moreover, with the exception of Gervini (2009), all robust methods mentioned above for single curve functional data require that entire functions are observed, therefore they are not directly applicable to sparsely observed functional data.
There is a related literature of functional canonical correlation analysis that can be used to study the correlation of paired functional data. In particular, by extending Leurgans, Moyeed and Silverman (1993), Boente and Kudraszow (2022) developed a robust smoothed canonical correlation analysis method for functional data. But this method assumes that the whole functions are observed and does not address the issue of only having sparse observations of the functions. In contrast, our method has the advantage of being able to estimate individual functions or curves with a small number of observations by borrowing strength across curves.
The rest of the paper is organized as follows. §2 gives specifics of the proposed model. §3 develops the estimation method and computation algorithm, and addresses model selection issues. §4 presents results from a simulation study.
§5 applies the proposed method with the method of Zhou et al. (2008) using a real Type Ia supernova dataset.

Reduced-rank mixed effects models with scale mixture of normals
Let ( ) and ( ) denote the measurements of two separate functional variables at time . Zhou et al. (2008) proposed the following reduced-rank mixed-effects (RRME) model where ( ) and ( ) are the mean curves for the two functional variables, ( ) = ( 1 ( ), 2 ( ), … , ( )) T and = ( 1 ( ), 2 ( ), … , ( )) T are vectors of principal component functions, and are subject specific vectors of principal component scores and are the random effects, ( ) and ( ) are subject specific measurement errors. In Zhou et al. (2008), the principle component scores , and the error terms ( ), ( ), are all assumed to follow normal distributions. The normality assumption makes the model fitting of (1) not resistant to the influence of outlying data points.
To introduce robustness in the RRME model, we propose to replace the normal distribution by a distribution with fat tails. The scale mixture of normal (SMN) distribution introduced by Andrews and Mallows (1974) will serve our purpose. We say that ∼ SMN( , ; ) with position parameter ∈ , scale parameter > 0 and mixing distribution , if it can be written as = + −1∕2 0 , where 0 ∼  (0, ), and is a positive random variable with distribution . Alternatively, has the following hierarchical representation: When  is a point mass distribution concentrated at point 1, then SMN( , ; ) reduces to  ( , ). When  =  = ( ∕2, ∕2), then follows a generalized t distribution with degrees of freedom , mean and scale parameter (Theodossiou, 1998). When  =  = ( , 1), follows the slash distribution with degrees of freedom (Osorio, 2016). Both the generalized t distribution and the slash distribution have fatter tails than the normal distribution, making them suitable to model outlying observations. When → ∞, both ( ∕2, ∕2) and ( , 1) will converge to the point mass at 1, and it follows that both the generalized t and slash distributions converge to a normal distribution.
To estimate unknown functions in the robust reduced-rank mixed effects model, we follow the same idea as Zhou et al. (2008) to represent them as spline functions and reduce the problem to estimation of spline coefficients. Specifically, we represent the mean curves ( ), ( ) and elements of ( ) and ( ), as a member of the same space of spline functions with dimension . The basis of the spline space, denoted by ( ), is chosen to be orthonormal; that is, the elements of ( ) = { 1 ( ), … , ( )} T satisfy ∫ ( ) ( ) = , or collectively, ∫ ( ) ( ) T = . The construction of an orthonormal basis using B-splines follows the procedure described in Appendix 1 of Zhou et al. (2008). Let and be -dimensional vectors of spline coefficients such that Let and be respectively × and × matrices of spline coefficients such that In general, the basis expansions in (5) and (6) are only approximations. Assuming the functions are smooth, the approximations can be very good if a sufficiently large is used (De Boor, 1978).
In practice, the functions are observed at discrete points (referred to as times) and the observation times can be different for different functions. Let 1 , … , and 1 , … , denote the observation times of and respectively. We refer to ( ), = 1, … , and ( ), = 1, … , as observations. Write Plugging (5) and (6) into (1), we see that the robust reduced-rank mixed-effects model implies the following model for observed functional data, It follows from (3) and (4) that To make the model identifiable, we require that T = and T = , then we have

Penalized likelihood estimation
In order to have enough flexibility of modeling arbitrary smooth functions, the number of spline basis functions needs to be sufficiently large. Following a common strategy in the nonparametric function estimation literature, e.g., Eilers and Marx (1996), we employ roughness penalties in a penalized likelihood formulation to prevent overfitting.
Let (Ξ; , ) denote the contribution to the likelihood from subject , where Ξ denote collectively unknown parameters , , , , 2 , 2 , , , , and from the distribution . The likelihood for all observations is . The method of penalized log-likelihood minimizes the criterion We propose that the roughness penalty take the form of integrated squared second derivatives, where , , , and are four penalty parameters. Applying the basis expansions (5) and (6), the penalty term can be written as where = ∫ ′′ ( ) ′′ ( ) T is the penalty matrix, and and are, respectively, the -th columns of and . Our use of four penalty parameters gives the flexibility of allowing different amounts of smoothing for the mean curves and principal component curves.
Direct minimization of the penalized log-likelihood (9) is challenging due to the complication of the likelihood evaluation. We use the latent variable representation of the model given in (7) and (8) and then apply the EM algorithm (Dempster, Laird and Rubin, 1977) to minimize the objective function (9). Treating the principle component scores , and the latent scale variables as missing data, we obtain the following complete data likelihood function where , , are the conditional density functions of , , and ( T , T ), respectively; is the density function of the latent variable . With an irrelevant constant ignored, it follows that Given the current estimate of the parameter set at step , denoted as Ξ ( ) , the EM algorithm maximizes with respect to the parameter set Ξ the following penalized conditional expected log-likelihood to obtain an update Ξ ( +1) . The algorithm iterates until convergence is reached. Under some regularity conditions, the algorithm is guaranteed to converge to a local minimizer of the penalized log-likelihood (9) (Wu, 1983). Details of the EM algorithm is given in the next section.

The EM algorithm 3.2.1. E step
In the E step, we need to obtain the penalized conditional expected log-likelihood (12). Based on (11), we need only to evaluate the conditional expectation of , log , , , given observations and , for = 1, … , . First, we know that the conditional density of given the observation and is where ( ) is the marginal density of the latent variable , and ( , | ) is the conditional density function of the -th subject which follows a multivariate normal distribution In the E step, we need to compute the conditional expectation of and log , denoted aŝ We give below the expressions of conditional expectations of and log in the RRME-t and RRME-slash model.
• RRME-slash model. In this model, the latent variable follows ( , 1). Following equation (13), the conditional density of is Note that the conditional distribution is a truncated Gamma distribution. We obtain the conditional expectation of and log as following: is the integral of a truncated Gamma distribution (Moore, 1982;Osorio, 2016).
Next, we evaluate the conditional expectation of and . Conditional on , the cross-covariance matrix of the vector of random effects and the vector of data is cov , Therefore, we obtain that where Therefore, Note that̄ and̄ are independent with . Using (14) and the law of iterated expectation, we obtain that Therefore,

M step
In the M step, we update the parameters by minimizing the penalized expected log-likelihood obtained in the E step. Since the parameters are well separated in the objective function, we update the estimates of the parameters sequentially in the following order: (1) 2 and 2 , (2) and , (3) , , and , (4) . Details of the updating formula are given below. The parameters appeared on the right hand side of equations are fixed at the values obtained from the previous iteration of the EM algorithm.
1. We update the estimates of 2 and 2 with the estimates of the other parameters fixed, 2. We update the estimates of and with the estimates of other parameters fixed, 3. To estimate , and , we need minimize the penalized expected log-likelihood with the orthonormal constraint of , . Firstly, to obtain the estimate of , we need only to minimize the objective function Denote * and * as the minimizers of the penalized expected log-likelihood, we update the columns of * and * sequentially. Write * = ( * 1 , * 2 , … , * ) and * = ( * 1 , … , * ). The update formulas arê * = ∑ Note that * and * may not be orthonormal and̂ ,̂ may not be diagonal with decreasing values. For identifiablity, we perform the following singular value decompositions to get the estimateŝ ,̂ and̂ such that̂ wherê and̂ are orthogonal matrices,̂ and̂ are diagonal matrices with decreasing diagonal elements. Moreover, [log ( )| , ; Ξ ( ) ] to update the estimate of , where Ξ ( ) is the current estimates of the parameters. With RRME-t model, we obtain that There is no closed-form solution of̂ . We suggest using an iterative optimization method, such as the secant method, to obtain̂ . With RRME-slash model, we have that The update formula iŝ

Model Selection
Specification of B-Splines. The number of knots and the positions of the knots are not crucial in many applications as long as sufficiently many knots are used to cover the data range, since the roughness penalty helps regularize the estimation and prevent overfitting (Eilers and Marx, 1996). We found that using 10-20 knots is often sufficient.
Choice of Penalty Parameters. When the number of principal components is fixed, the -fold ( = 5 or 10) within-subject cross-validation (CV) is used to select the penalty parameters, using the mean squared error as the metric to measure goodness-of-fit. The downhill simplex method of Nelder and Mead (1965) is used to search the optimal penalty parameters.
Number of Principal Components. The within-subject cross-valiation can also be used to select the number of principal components. Specifically, for each pair of ( , ), let̂ ,̂ ,̂ ,̂ be the penalty parameters selected by minimizing the within-subject cross-validation, and denote the corresponding CV value as ( , ;̂ ,̂ ,̂ ,̂ ). We may choose the number of principal component as Our experience suggests that (̂ ,̂ ) may be too large than necessary. Inspired by the idea of Cattell's scree test (Cattell, 1966), we propose to choose the smallest , whose CV value is not greater than the smallest CV value by a factor of . To be specific, denote the set

Simulation
In §4.1, we present the simulation setup with 3 scenarios. Then we applied RREM-t, RRME-slash and RRMEnormal models on each scenario. The fitting performance and the result of selecting the number of principle component will be shown in §4.2.

Simulation Setup
In this simulation, we generate the simulation data based on equations (1), (3) and (4) We consider two levels for the scale parameters: 2 = 2 = 0.04, 2 = 2 = 0.25. We simulate the data in the following four scenario setups. The first three scenarios consider datasets containing outliers, and the last scenario considers datasets without outliers.
• Scenario 3: Let the latent variable = 1. Note that all random variables follow normal distributions. First, we generate data from this model. Next, we sample 1% observations from ( ) and ( ), = 1, … , , = 1, … , , and with probability 1∕2 add or subtract a random number between 8-10 to each of the selected observations. The selected observations now become outliers to the data set. See an example in Figure 1. The data are generated from an RRME-normal with some artificial outliers added. This is a scenario that none of the models is correct.
• Scenario 4: The latent variable ≡ 1. Note that this model is an RRME-normal model.
For each scenario, with each specific degrees of freedom where applicable, combined with each level of scale parameters, we set up the simulation with = 100 subjects; for each subject, the number of observed time points is generated from 1+ (15, 0.9); the first time point is set at 0 and the rest time points are generated independently from the uniform distribution (0, 1). Figure 2 shows examples of the observations from three selected pairs of curves in each of the four scenarios.

Result
To check the performance of our proposed model on fitting the mean curves, principal component curves and individual curves, we ran each simulation setup = 500 times. The simulated datasets were fitted using RRMEt, RRME-slash and RRME-normal model separately. We used cubic B-splines with 10 interior knots for the basis functions. The penalty parameters were selected using the 10-fold within-subject CV and the downhill simplex method. We fixed the number of principal components ( , ) = (2, 2) in the model fitting.
The integrated squared errors (ISE) between a fitted curveĥ and the corresponding true curve ℎ is defined by We use the mean ISE (MISE) to assess the quality of the fitted curves. To be specific, for mean curves and principal component curves, MISE between fitted curves {ĥ } =1 and the corresponding real curves {ℎ } =1 is defined by For individual curves, the MISE between fitted curves {ĥ } and the corresponding real curves {ℎ }, = 1, … , , = 1, … , , is defined as following, Table 1 presents the MISEs for the mean curves, the principal component curves as well as the individual curves, under the four simulation scenarios. In general, when 2 = 2 = 0.04, the MISEs are smaller than that when 2 = 2 = 0.25. The performance of RRME-t model and RRME-slash model are similar. When the data were generated from scenario 1 and 2, compared with the RRME-normal model, the RRME-t and RRME-slash model have similar MISEs on the mean curves, but have smaller MISEs on the principal component curves and individual curves. When the data were generated from the RRME-t with = 2 or RRME-slash model with = 1, the MISEs for using RRME-t or RRME-slash model are smaller than that using RRME-normal model by 86%-91% on principal component curves, and by 23%-41% on individual curves. However, as the degrees of freedom increases, the MISEs using RRME-t or RRME-slash model get closer to that using RRME-normal model. When the data were generated from Scenario 3, the RRME-normal model with outliers, compared with MISEs using RRME-normal model, the MISEs using RRME-t or RRME-slash model are smaller by 11%-21% on mean curves, by 15%-64% on principal component curves, and similar MISEs on individual curves. When the data were generated from Scenario 4, the MISE for RRME-T, RRME-slash and RRME-normal are comparable, which means that the robust models perform well when there is no outlier.
From Table 1, we observe that performance of the proposed robust method is not sensitive to the choice of the t or the slash distribution, i.e., using either of the two distributions gives similar results. This phenomenon is also observed in our extensive simulation studies not reported here. In practice, one can use either the t or the slash distribution, but we would recommend using the t distribution because it is familiar to statisticians and the calculation is simpler.
To evaluate the proposed method on the selection of the number of principal components, we used the datasets simulated from Scenario 1 with degrees of freedom = 2, 5 and Scenario 2 with = 1, 2, with the two levels of scale parameters. Datasets simulated from Scenario 1 were fitted with RRME-t model and datasets simulated from Scenario 2 were fitted with RRME-slash model. For all model fitting, we used cubic b-splines with 10 interior knots for the basis functions. The penalty parameters were selected using 10-fold CV and the simplex method. The numbers of principal components were selected from {( , ) ∶ 1 ≤ ≤ 3, 1 ≤ ≤ 3} with = 0, 0.01, 0.05. From Table 2, we can see that for the two Scenarios, the levels of the scale parameters have little effect on the rate of correct choice of the number of principal components except when = 0 and = 5 or = 2. For the same model with the same scale parameter, the rate of correct choice increases as the degrees of freedom increases. As shown in Table 2, we can see that the rate of correct choice increases dramatically as the CV value relaxes by a factor . We suggest to use = 0.01 or 0.05 in practice.

Type Ia Supernova Light Curve Example
In this section, we applied RRME-t and RRME-normal model to a data set consisting of 102 Type Ia supernova multi-bands light curves studied by He et al. (2018). In this data set, about 9% SNeIa are identified as SN1991bg-like supernovae and about 82% SNeIa are identified as normal SNeIa. The rest of the data consist of one unlabeled Type Ia supernova, three super-Chandrasekhar SNeIa or "SNe Iax" whose photometric characteristics are similar as SN1991bglike SNeIa and five overluminous subtypes with different photometric features from normal or SN1991bg-like SNeIa. Since light curves of normal SNeIa have similar shapes and differ from those of SN1991bg-like supernovae, we can consider SN1991bg-like supernovae are the outliers.
Each Type Ia supernova has several light curves corresponding to different astronomical filters or passbands. We focused on and band in this study. Essential astronomical corrections were performed and all light curves were aligned by setting the peak magnitude to be zero and the Julian date of maxima to be zero. The redshift ( ) effect were removed by dividing Julian dates by a factor of (1 + ), transforming Julian date to phase. We treat each light curve as a functional observation whose explanatory variable is phase and response variable is magnitude. Note that brighter object has smaller magnitude. RRME-t and RRME-normal model were utilized to fit and band light curves. To model the curves on the time interval [−10, 50], we used cubic B-splines with 15 knots equally placed on [−5, 45] and two extra knots on each side of this interval, −10, −7.5 and 47.5, 50, which is consistent with the knots selected in He et al. (2018). The penalty parameters were selected by 10-fold within-subject cross-validation (CV) taking mean square error (MSE) as the metric and the downhill simplex method. Table 3 presents the CV values of different numbers of principal components for the RRME-t model and the RRME-normal model. Based on the selection criterion introduced in §3.3, for both models, we selected 2 principal components for R band and 3 principal components for I band. In addition, we would like to point out that the selection result remains the same for = 0, and 0.01.
Some example fits are shown in Figure 3, where one can compare the quality of the fits by the RRME-t and RRME-normal model for three normal SNeIa: SN2004eo, SN2005iq, SN2006ac and one SN1991bg-like supernova SN2002fb. The observed data of SN2002fb illustrate the characteristics of SN1991bg-like SNeIa, namely the absence of a shoulder in band and the missing of the secondary maxima in band which are supposed to occur about 20 -35 days past maximum light for normal SNeIa. Both the RRME-t and RRME-normal models provide reasonably good fits for these observed data. The RRME-t model provides better fits for three normal SNeIa, especially in band where the secondary peak feature possessed by three normal SNeIa are not fully captured by RRME-normal model. However, RRME-normal model performs pooly on the three normal SNeIa, which reflects the strong leverage of the outlier under the normal assumption. Generally speaking, the within-subject CV values of RRME-t model are smaller than RRME-normal model for both bands. The 10-fold within-subject CV values of RRME-t model of R band and I band are 5.49 and 10.25, respectively. The 10-fold within-subject CV values of RRME-normal model of R band and I band are 5.53 and 14.22, respectively. The RRME-t model shows significant improvement over the RRME-normal model in fitting I band light curve. Figure 4 shows the estimated mean functions and the effects of the principal component curves for and band obtained with the RRME-t model and RRME-normal model. The RRME-normal model provides similar results with the ones derived by He et al. (2018) in terms of the mean curves and the first two principal component curves in and bands. The mean curves estimated by the RRME-t and RRME-normal models are similar in and band, respectively. Furthermore, for R band, RRME-t and RRME-normal model have similar first and second principal components, where the first principal components control the overall stretch of the light curve which is an important feature for SNeIa photometric study, and the second principal components contribute slightly to the magnitude around the second peak. For I band, RRME-t and RRME-normal model produce similar first principal components that adjust the light curve shape after the maximum brightness. The second principal components from both models control the brightening rate before the maximum; however, the second principal component of RRME-t model focuses on the increasing rate and the declining rate before and after the second peak which adapts to the intrinsic variability of the second peak locations for normal SNeIa, whereas the counterpart of RRME-normal model mainly adjusts the overall magnitude around the secondary peak. The third principal component curve of RRME-t model provides further adjustment to the magnitude around the second peak but the contribution from the counterpart of RRME-normal model is negligible. With a fixed relation between the brightening rate before the maximum brighteness and the curve shape around second peak, the principal component curves of RRME-normal model presents less flexibility in precisely modeling the shape of the light curves, especially at around 20 days post the maximum, explaining the inferior fitting results of three normal SNeIa in I bands in Figure 3. Figure 2: Observations from three selected pairs of curves generated for each scenario are presented, RRME-T(2), RRMEslash(1), RRME-normal and contaminated RRME-normal models from top to bottom, with 2 = 2 = 0.04. Different color represents observations from different functions for each scenario. The unobserved full functions are represented by the dashed lines. The two snow points in the third row are contaminated point outliers.
Y Z RRME−t(2) RRME−slash (1) Contaminated RRME−normal   Table 2 Selection of the number of principal components in the simulation study of §4. Data are generated from the RRME-t or RRME-slash models with the true number of components ( , ) = (2, 2), different levels of scale parameters and , and degrees of freedom . Reported are the percentages of correct choice of the number of principal components ( , ) = (2, 2) from different based on 500 simulation runs.  Table 3 Type Ia supernova light curve example. The CV values of the RRME-t and RRME-normal models for different number of principal components. Each reported number equals to the actual number multiplied by 1000.  Figure 3: Four specific supernova stars, SN2004eo, SN2005iq, SN2006ac and SN2002fb. Each column shows the observations and the corresponding fitted curves. The two rows correspond to the R band and the I band, respectively. The black dots are observations; the red solid lines are fitted curves using the RRME-normal model, and the blue dashed lines are fitted curves using the RRME-t model. The horizontal axis is phase and the vertical axis is aligned magnitude. Phase Magnitude RRME−normal RRME−t Figure 4: The estimated mean functions and effects of functional principal components for the type Ia supernova light curve example. The first two columns show the results for R band while the last two columns show the results of I band when applying the RRME-t and RRME-normal models. The -th row of each column represents the effect of -th principal component (PC). The solid lines are estimated mean curves and the "+" and "-" points represent the mean curve ± 2 × standard deviation of principal score × PC function.