BFDA: A Matlab Toolbox for Bayesian Functional Data Analysis
JJSS
Journal of Statistical Software
MMMMMM YYYY, Volume VV, Issue II. doi: 10.18637/jss.v000.i00
BFDA : A
MATLAB
Toolbox for Bayesian FunctionalData Analysis
Jingjing Yang
University of Michigan
Peng Ren
Suntrust Banks, Inc.
Abstract
We provide a
MATLAB toolbox,
BFDA , that implements a Bayesian hierarchicalmodel to smooth multiple functional data with the assumptions of the same underlyingGaussian process distribution, a Gaussian process prior for the mean function, and anInverse-Wishart process prior for the covariance function. This model-based approachcan borrow strength from all functional data to increase the smoothing accuracy, as wellas estimate the mean-covariance functions simultaneously. An option of approximatingthe Bayesian inference process using cubic B-spline basis functions is integrated in
BFDA , which allows for efficiently dealing with high-dimensional functional data.Examples of using
BFDA in various scenarios and conducting follow-up functionalregression are provided. The advantages of
BFDA include: (1) Simultaneously smoothsmultiple functional data and estimates the mean-covariance functions in anonparametric way; (2) flexibly deals with sparse and high-dimensional functional datawith stationary and nonstationary covariance functions, and without the requirement ofcommon observation grids; (3) provides accurately smoothed functional data forfollow-up analysis.
Keywords : functional data analysis, Bayesian hierarchical model, Gaussian process,Inverse-Wishart process,
MATLAB .
1. Introduction
Since Ramsay and Dalzell (1991) first coined the term “functional data analysis” (FDA) foranalyzing data that are realizations of a continuous function, many statistical methods andtools have been proposed for FDA. For examples, Graves et al. (2010) provided both R package fda (Ramsay et al. MATLAB package fdaM (Ramsay 2014) for standardfunctional data analysis (Ramsay and Silverman 2002, 2005); Febrero-Bande and de laFuente (2012) provided R package fda.usc for implementing nonparametric functional data a r X i v : . [ s t a t . O T ] F e b BFDA : Bayesian Functional Data Analysis analysis methods (Vieu and Ferraty 2006) with fda (Ramsay et al. et al. (2005a,b) developed the key technique of Functional Principal Component Analysis (FPCA)for analyzing sparse functional data, accompanied by the
MATLAB package
PACE (Yao et al.
WinBUGS (Sturtz et al. R package GPFDA for applying the Bayesian nonparametric Gaussian process (GP)regression models (Shi and Choi 2011). However, the smoothing step that constructsfunctions from noisy discrete data has been neglected by most of the existing FDA methodsand tools, where functional representations are integrated in the analyzing models. On theother hand, most of the existing smoothing methods process functional samples individually(e.g., cubic smoothing spline (CSS) and kernel smoothing (Green and Silverman 1993;Ramsay and Silverman 2005)), which is likely to induce bias when functional samples are ofthe same distribution.Here, we provide a
MATLAB toolbox
BFDA for simultaneously smoothing multiplefunctional observations from the same distribution and estimating the underlyingmean-covariance functions, using a nonparametric Bayesian Hierarchical Model (BHM) withGaussian-Wishart processes (Yang et al. et al.
BFDA is flexible for analyzing sparse and dense functional data without therequirement of common observation grids, suitable for analyzing functional data with bothstationary and nonstationary covariance functions, and efficient for high-dimensionalfunctional data using a Bayesian framework with Approximations by Basis Functions(BABF) (Yang et al.
BFDA provides options for implementing thestandard Bayesian GP regression method, conducting Bayesian principal componentanalysis, and using the fdaM package for follow-up FDA.In the following context, we first review the BHM and BABF methods in Section 2, and thenprovide examples using
BFDA with simulated data in Section 3. In Section 4, we comparethe functional linear regression results by fdaM using the smoothed data by
BFDA and CSS.Last, we conclude with a discussion in Section 5. Details of input options and outputs, as wellas example
MATLAB scripts of generating the simulation results in this paper, are providedin the Appendix.
2. Methods overview
Consider functional data { X i ( t ); t ∈ T , i = 1 , , · · · , n } that are generated from the samestochastic process with independent measurement errors. In order to simultaneously smoothall noisy observations and estimate mean-covariance functions, Yang et al. (2016) proposedthe following Bayesian Hierarchical Model (BHM) with Gaussian-Wishart processes: X i ( t ) = Z i ( t ) + (cid:15) i ( t ); Z i ( · ) ∼ GP ( µ Z ( · ) , Σ Z ( · , · )) , (cid:15) i ( · ) ∼ N (0 , σ (cid:15) ); (1) µ Z ( · ) | Σ Z ( · , · ) ∼ GP (cid:18) µ ( · ) , c Σ Z ( · , · ) (cid:19) , Σ Z ( · , · ) ∼ IW P ( δ, σ s A ( · , · )) , σ (cid:15) ∼ IG ( a (cid:15) , b (cid:15) ); ournal of Statistical Software σ s ∼ IG ( a s , b s );where { Z i ( t ); i = 1 , · · · , n } denotes the underlying true functional data following the sameGP distribution with mean function µ Z ( · ) and covariance function Σ Z ( · , · ), IW P denotesthe Inverse-Wishart process (IWP) prior (Dawid 1981) for the covariance function, IG denotes the Inverse-Gamma prior, and ( µ ( · ) , c, δ, A ( · , · ) , a (cid:15) , b (cid:15) , a s , b s ) are hyper-parametersto be determined. The IWP prior on Σ Z ( · , · ) models the covariance functionnonparametrically and therefore makes the BHM method suitable for analyzing functionaldata with unknown stationary and nonstationary covariance structures. The hyperparameter σ s provides the flexibility of estimating the scale of the covariance structure inthe IWP prior from the data.For the hyper-parameter setup, we take µ ( · ) as the smoothed empirical mean estimate, c as1 for the same data variation on the mean function, δ as 5 for noninformative prior on thevariance function, and determine ( a (cid:15) , b (cid:15) , a s , b s ) by matching the hyper-prior moments withthe empirical estimates. In addition, A ( · , · ) can be taken as the Mat´ern correlation kernel foranalyzing functional data with stationary covariance (default in BFDA ),Matern cor ( d ; ρ, ν ) = 1Γ( ν )2 ν − (cid:18) √ ν dρ (cid:19) ν K ν (cid:18) √ ν dρ (cid:19) , d ≥ , ρ > , ν > , where d denotes the distance between two grid points, ρ is the scale parameter, ν is theorder of smoothness, and K ν ( · ) is the modified Bessel function of the second kind; while asan smoothed empirical covariance estimate for analyzing functional data with nonstationarycovariance.Although the BHM is constructed with infinite-dimensional Gaussian-Wishart processes,practical posterior inference will be conducted in a finite manner, e.g., on the observationgrids { t i } , pooled grid t = ∪ ni =1 t i , or other evaluation grids. For accommodating uncommonobservation grids, especially sparsely observed data, BHM evaluates data functions andmean-covariance functions on the pooled grid, while estimating the unobserved functionaldata by conditioning on the observations (similarly as PACE ).Denoting X i ( t i ) by X t i , Z i ( t i ) by Z t i , µ ( t ) by µ , µ Z ( t ) by µ Z , Σ Z ( t , t ) by Σ Z , and A ( t , t )by A , BHM conducts Bayesian inference for ( { Z i ( t ) } , µ Z , Σ Z , σ (cid:15) , σ s ) by the Monte CarloMarkov Chain (MCMC) algorithm (essentially a Gibbs sampler (Geman and Geman 1984))as follows (refer to Yang et al. (2016) for details):Step 0: Set initial values. Set hyper-parameters ( c, µ , ν, ρ, a (cid:15) , b (cid:15) , a s , b s ). Take ( µ , σ (cid:15) )as the empirical estimates, { Z t i } as the raw data { X t i } , and Σ Z as an identity matrix.Step 1: Conditioning on { X t i } and ( µ Z , Σ Z ), update { Z i ( t ) } from the correspondingconditional multivariate normal (MN) distributions;Step 2: Conditioning on { X t i } and { Z t i } , update the noise variance σ (cid:15) from theconditional Inverse-Gamma (IG) distribution;Step 3: Conditional on { Z i ( t ) } and Σ Z , update µ Z from the conditional MNdistribution;Step 4: Conditioning on { Z i ( t ) } and µ Z , update Σ Z from the conditionalInverse-Wishart (IW) distribution; BFDA : Bayesian Functional Data Analysis
Step 5: Conditioning on Σ Z , Update σ s from the conditional Gamma distribution.Specifically, the averages of posterior samples of { Z i ( t ) , µ Z , Σ Z } are taken as estimates forfunctional signals and mean-covariance functions.In addition, BFDA uses existing
MATLAB package mcmcdiag (S¨arkk¨a and Aki 2014) todiagnosis the MCMC convergence by potential scale reduction factor (PSRF) (Gelman andRubin 1992), and implements the method proposed by Yuan and Johnson (2012) with thepivotal discrepancy measures (PDM) of standardized residuals for the goodness-of-fit diagnosisof the assumed model.
Because BHM (Yang et al. O ( np m ) with n samples, p pooled-grid points, and m MCMC iterations, the method encounters computationalbottleneck for analyzing functional data with large pooled-grid dimension p . To resolve thiscomputational bottleneck issue with high-dimensional data, BFDA enables the option ofusing the BABF method proposed by Yang et al. (2015), in which the posterior inference inBHM is conducted with approximations using basis functions. Here, we briefly outline theBABF method.First select a working grid based on data density, τ = ( τ , τ , · · · , τ L ) (cid:62) ⊂ T , L << p , thenapproximate { Z i ( τ ) } by a system of basis functions (e.g., cubic B-splines). Let B ( · ) = [ b ( · ) , b ( · ) , · · · , b K ( · )] denote K selected basis functions with coefficients ζ i = ( ζ i , ζ i , · · · , ζ iK ) (cid:62) , then Z i ( τ ) = (cid:80) Kk =1 ζ ik b k ( τ ) = B ( τ ) ζ i . Generally, K can beselected as L , then ζ i = B ( τ ) − Z i ( τ ), a linear transformation of Z i ( τ ). Note that even if B ( τ ) is singular or non-square, ζ i can still be written as a linear transformation of Z i ( τ ),e.g., using the generalized inverse (James 1978) of B ( τ ) .Because ζ i is a linear transformation of Z i ( τ ) that follows a MN distribution under theassumptions in Equation 1, the induced Bayesian hierarchical model for { ζ i } is derived as ζ i ∼ M N ( µ ζ , Σ ζ ); µ ζ = B ( τ ) − µ Z ( τ ); Σ ζ = B ( τ ) − Σ Z ( τ , τ ) B ( τ ) −(cid:62) . (2)Further, from the assumed priors of ( µ Z ( · ) , Σ Z ( · , · )) in Equation 1, with Ψ( τ , τ ) = σ s A ( τ , τ ),the following priors of ( µ ζ , Σ ζ ) are also induced: µ ζ | Σ ζ ∼ M N (cid:0) B ( τ ) − µ ( τ ) , c Σ ζ (cid:1) ; (3) Σ ζ ∼ IW ( δ, B ( τ ) − Ψ( τ , τ ) B ( τ ) −(cid:62) ) . Then, the BABF approach has computation complexity O ( nK m ) with the following MCMCalgorithm (refer to Yang et al. (2015) for details):Step 0: Set initial values similarly as in BHM. Set hyper-parameters( c, µ , ν, ρ, a (cid:15) , b (cid:15) , a s , b s ). Take ( µ Z ( τ ) , Σ Z ( τ , τ ) , σ (cid:15) ) as empirical estimates, inducingthe initial values for ( µ ζ , Σ ζ ) by Equation 2.Step 1: Conditioning on { X t i } and ( µ ζ , Σ ζ , σ (cid:15) ), update { ζ i } from the conditional MNdistribution.Step 2: Conditioning on { ζ i } , update µ ζ and Σ ζ respectively from the conditional MNand IW distributions. ournal of Statistical Software { ζ i } , µ ζ , Σ ζ ), approximate { Z t i , µ Z ( t i ), Σ Z ( t i , t i ),Σ Z ( τ , t i ) , Σ Z ( t i , τ ) , Σ Z ( τ , τ ) } by Z t i = B ( t i ) ζ i , µ Z ( t i ) = B ( t i ) µ ζ , Σ Z ( t i , t i ) = B ( t i )Σ ζ B ( t i ) (cid:62) , Σ Z ( τ , t i ) (cid:62) = Σ Z ( t i , τ ) = B ( t i )Σ ζ B ( τ ) (cid:62) , Σ Z ( τ , τ ) = B ( τ )Σ ζ B ( τ ) (cid:62) . Step 4: Conditioning on { Z t i } and { X t i } , update σ (cid:15) by the conditional IG distribution;Step 5: Conditioning on Σ Z ( τ , τ ), update σ s by the conditional Gamma distribution.As a result, the posterior estimates of ( { ζ i } , µ ζ , Σ ζ ) are given by the averages of the MCMCsamples, which are then used to estimate { Z t i , µ Z ( t i ), Σ Z ( t i , t i ) } by the approximationformulas in Step 3. BFDA uses the existing
MATLAB package bspline (Hunyadi 2010) to construct B-spline basisfunctions, using the optimal knot sequence for interpolation at the working grid τ . Theoptimal knot sequence is generated by the MATLAB function optknt (Gaffney and Powell1976; Micchelli et al. et al. (2015) instructed selecting τ torepresent data densities ( L maybe selected by grid search with test data), such as takingthe (cid:16) L +1 , · · · , L × L +1 (cid:17) percentiles of the pooled grid, or the equally-spaced grid for evenlydistributed data.
3. Examples with simulated data
In this Section, we provide examples of using
BFDA to analyze simulated functional datawith stationary and nonstationary covariance functions, common and uncommon (sparse)observation grids, as well as random observation grids. The simulation data used for theexample results were generated with n =30, p =40, au =0, bu = π/ s = √ r =2, nu =3.5, rho =0.5, dense =0.6, and pgrid as the equally spaced grid over (0 , π/
2) with length 40.
BFDA provides the convenience of generating simulated functional data from the same GPwith mean function µ ( t ) = 3 sin (4 t ), stationary covariance function s Matern cor ( d ; ρ, ν ), andnoises ∼ N (0 , ( s/r ) ) by > GausFD_cgrid = sim_gfd(pgrid, n, s, r, nu, rho, dense, cgrid, stat); where pgrid denotes the pooled grid, n denotes the number of functional samples, r denotesthe signal to noise ratio (i.e., the ratio between the signal and noise standard deviations), rho denotes the Mat´ern scale parameter, nu denotes the Mat´ern order of smoothness. Here, cgrid is a boolean indicator that controls the output as either common-grid data on pgrid ( cgrid =1) or uncommon-grid data with a randomly selected proportion (given by dense )of the full data on pgrid ( cgrid =0). In addition, stat =1 specifies simulating stationarydata from GP (3 sin (4 t ) , s Matern cor ( d ; ρ, ν )), while stat =0, specifies simulating data from anonlinearly transformed GP with mean function µ ( t ) = 3( t +0 . sin (4 t / ) and nonstationarycovariance function Σ( t, t (cid:48) ) = s ( t + 0 . t (cid:48) + 0 . cor ( | t / − t (cid:48) / | ; ρ, ν ) . BFDA : Bayesian Functional Data Analysis
Let p denote the length of pgrid . The output GausFD_cgrid is a data structure consisted witha cell of true data (Xtrue cell × n ), a cell of noisy data (Xraw cell × n ), a cell of observationgrids (Tcell × n ), a true covariance matrix on pgrid (Cov true p × p ), and a true mean vectoron pgrid (Mean true × p ). Common grids
First, we need to setup the required parameter structure by function setOptions_bfda . Forexample, to analyze functional data with common observation grids and stationary covariancefunction by BHM, the structure param can be set as > param = setOptions_bfda(’smethod’, ’bhm’, ’cgrid’, 1, ’mat’, 1, ’M’, 10000,’Burnin’, 2000, ’w’, 1, ’ws’, 1); where each parameter is followed by its value, and unspecified parameters are taken asdefault values (Appendix A.1.). Specifically, smethod=’bhm’ denotes using the BHMmethod; cgrid=1 denotes the analyzed data are of common-grid; mat=1 denotes taking A ( · , · ) in Equation 1 as the Mat´ern correlation function; M=10000 denotes the number ofMCMC iterations;
Burnin=2000 denotes the number of MCMC burn-ins; w=1 and ws=1 areused to tune the Gamma priors for σ (cid:15) and σ s .With both param and GausFD_cgrid , we can then call the main function
BFDA() by > [out_cgrid, param] = BFDA(GausFD_cgrid.Xraw_cell, GausFD_cgrid.Tcell,param); for smoothing and estimating the common-grid functional data by BHM. The outputstructure out_cgrid contains smoothed estimates for the signals ( out_cgrid.Z ), meanfunction ( out_cgrid.mu ), covariance function ( out_cgrid.Sigma ), and other parameters inEquation 1, along with the corresponding 95% point-wise credible intervals (Appendix A.1.).The output argument param is the updated parameter structure. Uncommon grids
To apply BHM on stationary functional data of uncommon-grid, e.g.,
GausFD_ucgrid generated by > GausFD_ucgrid = sim_gfd(pgrid, n, s, r, nu, rho, dense, 0, stat); the main function
BFDA can be called by > param_uc = setOptions_bfda(’smethod’, ’bhm’, ’cgrid’, 0, ’mat’, 1, ’M’,10000, ’Burnin’, 2000, ’pace’, 1, ’ws’, 0.1);> [out_ucgrid, param_uc] = BFDA(GausFD_ucgrid.Xraw_cell, GausFD_ucgrid.Tcell,param_uc); where cgrid is set as 0 in param_uc . Example results
In Figure 1(a, b), we show that the smoothed signals by BHM (blue solid) are close to thetruth (red dashed), and the coverage probabilities of the 95% point-wise credible intervals ournal of Statistical Software > .
95, for both scenarios with common and uncommon grids. In addition,the nonparametric mean estimates by BHM (blue solid curves in Figure 1(c, d)) are alsosmooth and close to the truth (red dashed), while the corresponding 95% point-wise credibleintervals (blue dotted) have coverage probabilities > .
9. In addition, we show that theBayesian nonparametric covariance estimates in Figure 2(a, b) are clearly smoother than thesample covariance estimate by using the raw common-grid data in Figure 2(c), and close tothe true stationary covariance in Figure 2(d). Importantly, although 40% information is lostfor the uncommon-grid scenario, BHM still produces good smoothing and estimation results. (a)
SampleTruthBHMBHM 95% CI (b) (c) (d)
Figure 1: Results of analyzing
Stationary functional data by BHM: (a) two sample signalestimates with common grids; (b) two sample signal estimates with uncommon grids; (c)mean estimate with common grids; (d) mean estimate with uncommon grids; along with 95%pointwise credible intervals (blue dots).
BFDA : Bayesian Functional Data Analysis
022 24 (a) (b) (c) (d)
Figure 2: Covariance estimates for
Stationary functional data: (a) BHM estimate withcommon grids; (b) BHM estimate with uncommon grids; (c) sample estimate with rawcommon-grid data; (d) true covariance surface. ournal of Statistical Software Common grids
To apply BHM on functional data with nonstationary covariance function and common grids,e.g.,
GausFD_cgrid_ns generated by > GausFD_cgrid_ns = sim_gfd(pgrid, n, s, r, nu, rho, dense, cgrid, 0); the main function
BFDA() can be called by > param_ns = setOptions_bfda(’smethod’, ’bhm’, ’cgrid’, 1, ’mat’, 0, ’M’,10000, ’Burnin’, 2000, ’pace’, 1, ’ws’, 0.01);> [out_cgrid_ns, param_ns] = BFDA(GausFD_cgrid_ns.Xraw_cell,GausFD_cgrid_ns.Tcell, param_ns);
Here, A ( · , · ) in (2.1) is set as the empirical smooth covariance estimate ( mat =0) that is givenby PACE (Yao et al. pace =1 (default), or by the sample covariance estimateusing CSS smoothed data with pace =0.
Uncommon grids
If the nonstationary functional data are collected on uncommon (sparse) grids, e.g.,
GausFD_ucgrid_ns generated by > GausFD_ucgrid_ns = sim_gfd(pgrid, n, s, r, nu, rho, dense, 0, 0); we only need to set cgrid =0, mat =0 in the parameter structure for the common-grid scenarioand then call the main function
BFDA() by > param_uc_ns = setOptions_bfda(’smethod’, ’bhm’, ’cgrid’, 0, ’mat’, 0, ’M’,10000, ’Burnin’, 2000, ’pace’, 1, ’ws’, 0.01);> [out_ucgrid_ns, param_uc_ns ] = BFDA(GausFD_ucgrid_ns.Xraw_cell,GausFD_ucgrid_ns.Tcell, param_uc_ns); where cgrid is set as 0 in param_uc_ns . Example results
Similarly, as shown in Figures 3 and 4, the BHM estimates of signals and mean-covariancefunctions are close to the truth. Specifically, the 95% pointwise credible intervals of theBHM signal estimates have coverage probabilities > .
95. Although BHM overestimated thecovariance, BHM captured the major covariance structure and produced a smoothed estimate.The magnitude of the BHM estimate can be tuned by ws , where a smaller ws will relativelyshrink the magnitude of BHM covariance estimate. We suggest users to tune this parameteraccording to the magnitude of sample covariance estimate.0 BFDA : Bayesian Functional Data Analysis (a)
SampleTruthBHMBHM 95% CI (b) (c) (d)
Figure 3: Results of analyzing
Nonstationary functional data: (a) two sample signalestimates with common grids; (b) two sample signal estimates with uncommon grids; (c)mean estimate with common grids; (d) mean estimate with uncommon grids; along with 95%pointwise credible intervals (blue dots). (a)
130 10 0 0210 220 (b)
130 10 00210 220 (c)
130 10 0 0210 220 (d)
130 10 0
Figure 4: Covariance estimates for
Nonstationary functional data: (a) BHM estimatewith common grids; (b) BHM estimate with uncommon grids; (c) sample estimate with rawcommon-grid data; (d) true covariance surface. ournal of Statistical Software BFDA also provides the convenience to simulate stationary and nonstationary functional datawith random observation grids from the same GPs as in Section 3.1. For example, a structureof functional data
GausFD_rgrid , with n independent observations and p random grids perobservation (uniformly generated from [ au , bu ], can be generated by > GausFD_rgrid = sim_gfd_rgrid(n, p, au, bu, s, r, nu, rho, stat); where stat=1 specifies simulating from the stationary GP, while stat=0 specifies simulatingfrom the nonstationary GP. Stationary data
To analyze stationary functional data by BABF, simply call the main function
BFDA() by: > param_rgrid = setOptions_bfda(’smethod’, ’babf’, ’cgrid’, 0, ’mat’, 1, ’M’,10000, ’Burnin’, 2000, ’m’, m, ’eval_grid’, pgrid, ’ws’, 1);> [out_rgrid, param_rgrid]= BFDA(GausFD_rgrid.Xraw_cell, GausFD_rgrid.Tcell,param_rgrid); where the working grid τ will be set as the equally spaced quantiles of the pooled grid bydefault, with length m (when tau is not initialized in param_rgrid ). Nonstationary data
For nonstationary functional data, e.g.,
GausFD_rgrid_ns generated by > GausFD_rgrid_ns = sim_gfd_rgrid(n, p, au, bu, s, r, nu, rho, 0); we can call the main function
BFDA() by > param_rgrid_ns = setOptions_bfda(’smethod’, ’babf’, ’cgrid’, 0, ’mat’, 0,’M’, 10000, ’Burnin’, 2000, ’m’, m, ’eval_grid’, pgrid, ’ws’, 0.05);> [out_rgrid_ns, param_rgrid_ns] = BFDA(GausFD_rgrid_ns.Xraw_cell,GausFD_rgrid_ns.Tcell, param_rgrid_ns); where mat is set as 0 in param_rgrid_ns . Example results
With random observation grids, the BABF method can efficiently analyze both stationary andnonstationary functional data, producing smooth estimates for signals and mean-covariancefunctions that are close to the truth (Figures 5, 6). Specifically, the 95% pointwise credibleintervals of signal estimates have coverage probabilities > . BFDA : Bayesian Functional Data Analysis (a)
SampleTruthBABFBABF 95% CI (b) (c) (d)
Figure 5: Results of analyzing functional data with
Random grids by BABF: (a) two samplesignal estimates with stationary data; (b) two sample signal estimates with nonstationarydata; (c) mean estimate with stationary data; (d) mean estimate with nonstationary data;along with 95% pointwise credible intervals (blue dots).
022 24 (a) (b)
20 1 10 0022 24 (c) (d)
20 1 10 0
Figure 6: Covariance estimates for functional data with
Random grids : (a) BABF estimatewith stationary data; (b) BABF estimate with nonstationary data; (c) true covariance surfacefor stationary data; (d) true covariance surface for nonstationary data. ournal of Statistical Software
4. Functional linear regression
We expect follow-up FDA results will be improved by using the accurately smoothedfunctional data produced by
BFDA . Specifically, we show examples of functional linearregression in this Section, considering the following two models, Y = β + (cid:90) X ( t ) (cid:62) β ( t ) dt + (cid:15) , (4) Y ( t ) = β ( t ) + X ( t ) (cid:62) β ( t ) + (cid:15) ( t ) ; (5)where • Y in Equation 4 denotes a n × Y ( t ) = ( y ( t ) , · · · , y n ( t )) (cid:62) in Equation 5 denotes a vector of functional responses; • X ( t ) denotes a n × q design matrix of q functional independent variables; • β ( t ) denotes a q × • β and β ( t ) denote the intercept terms; • (cid:15) and (cid:15) ( t ) denote the error terms.Note that X ( t ) and β ( t ) can also denote nonfunctional covariates and coefficients, becausenonfunctional variables can be thought as constant functions of t . To evaluate the improvement of regression results using the smoothed data by
BFDA , we firstsimulated 30 raw stationary GP trajectories { X i ( t i ) } with random grids uniformly generatedfrom [0 , π/ sim_gfd_rgrid(30, 40, 0, 1.5708, 2.2361, 2, 3.5, 0.5, 1) . Thenwe simulated the scalar responses by Y i = (cid:82) . X i ( t ) t dt + (cid:15) , and the functional responsesby Y i ( t ) = X i ( t ) t + (cid:15) , with (cid:15) ∼ N (0 , fRegress() in fdaM requires inputs of functionaldata with common grids, we interpolated the simulated true data, smoothed data by BABFwith BFDA , and the raw data with noises on the equally spaced common grid (with length 40)over [0 , π/ csaps() with the suggestedoptimal smoothing parameter 1). As a result, the interpolated signals from the raw data areequivalent to the individually smoothed ones by CSS (one example curve is shown in Figure7(a)).With the smoothed data by BABF and CSS, we respectively fitted the functional linearmodels (Equations 4 and 5) using 20 randomly chosen signals, and then tested the predictionresults using the remains. We replicated this fitting process for 100 times, and evaluated theperformance by the average mean square errors (MSEs) of the fitted and predicted responses.
For the case with scalar responses, although the fitted coefficient functions using bothsmoothed data by BABF and CSS are close to the truth (Figure 7(b, c)), with coverage4
BFDA : Bayesian Functional Data Analysis probabilities > .
95 for the 95% confidence intervals, the average MSEs of the fitted andpredicted responses from 100 replications are smaller for using the BABF smoothed datathan the ones using the CSS smoothed data (0.311 vs. 0.388 for fitted responses, 0.497vs. 0.677 for predicted responses, as shown in Table 1). Figure 8 shows the results of anexample replication of this fitting and predicting process with scalar responses. t -4-3-2-10123456 x ( t ) (a) Raw dataBABFCSSTruth (b) (c)
Figure 7: (a) Example estimates of X i ( t ); (b) the estimate of β ( t ) using the smoothed databy BABF; (c) the estimate of β ( t ) using the smoothed data by CSS; along with the truth inthe cyan dotted lines. -5 -4 -3 -2 -1 0 1 2 Truth -5-4-3-2-10123 F i tt ed (a) BABFCSS -5 -4 -3 -2 -1 0 1 2 3
Truth -4-3-2-101234 P r ed i c t ed (b) BABFCSS
Figure 8: (a) Fitted vs. true scalar responses; (b) predicted vs. true scalar responses.
For the case with functional responses, we can see that the fitted intercept term using theBABF smoothed data is closer to the truth (constant 0) with narrower 95% confidence intervalthan the one using the CSS smoothed data (Figure 9(a, c)). In addition, the coefficientfunction using BABF smoothed data has narrower 95% confidence interval but higher coverageprobability (Figure 9(b, d)). Consequently, both fitted and predicted functional responsesusing the BABF smoothed data have smaller point-wise MSEs out of 100 replications, 0.417vs. 1.190 for fitted functional responses, 0.464 vs. 1.354 for predicted functional responses(Table 1). Figure 10 shows the results of an example replication of this fitting and predictingprocess with functional responses. ournal of Statistical Software
Y Y ( t ) Y Y ( t ) Fitted 0.311 (0.061) 0.417 (0.049) 0.388 (0.074) 1.190 (0.186)Predicted 0.497 (0.289) 0.464 (0.112) 0.677 (0.435) 1.354 (0.419)Table 1: Average MSEs of the fitted and predicted responses for 100 replicates, along withthe standard deviations of these MSEs among 100 replicates in the parentheses, for scalarresponses Y and functional responses Y ( t ) . (a) (c) (b) (d) Figure 9: (a) Intercept function β ( t ) using the BABF smoothed data; (b) Intercept function β ( t ) using the CSS smoothed data; (c) coefficient function β ( t ) using the BABF smootheddata; (d) coefficient function β ( t ) using the CSS smoothed data; along with 95% confidenceintervals and true coefficient functions in the cyan dotted lines. t -3-2-10123456 (a) BABFCSSTruth t -5-4-3-2-101 (b) BABFCSSTruth
Figure 10: (a) Example fitted functional responses; (b) example predicted functionalresponses; with true signals in the cyan dotted lines.6
BFDA : Bayesian Functional Data Analysis
5. Discussion
The
MATLAB tool
BFDA presented in this paper can simultaneously smooth multiplefunctional observations and estimate the mean-covariance functions, assuming the functionaldata are from the same GP. The smoothed data by
BFDA are shown to be more accuratethan the conventional individual smoothing methods, thus improving follow-up analysisresults. The advantages of
BFDA include: • Simultaneously smoothing multiple functional samples and estimating mean-covariancefunctions in a nonparametric way; • Flexibly handling functional data with stationary and nonstationary covariancefunctions, common or uncommon (sparse) observation grids; • Efficiently dealing with high-dimensional functional data by the BABF method.
BFDA is suitable for analyzing data that can be roughly assumed as from the same GPdistribution. We recommend using the BHM method for low-dimensional functional datawith common grids or sparse functional data, and using the BABF method forhigh-dimensional data with dense grids (including both common and random grids). Inaddition, we recommend using the Mat´ern function as the prior covariance structure foranalyzing functional data with stationary covariance functions, while using the empiricalcovariance estimate (e.g., the estimate by
PACE is recommended) for analyzing functionaldata with nonstationary covariance functions.The follow-up functional data analysis can be conducted using the existing softwares (e.g., fdaM in MATLAB , fda in R ). Examples are provided in BFDA about using fdaM with thesmoothed data by
BFDA , showing improved regression results than using the individuallysmoothed data. Details about the inputs and outputs of
BFDA , and part of the example
MATLAB scripts are provided in the Appendices. The
BFDA tool and example scripts arefreely available at https://github.com/yjingj/BFDA . We will continue integrating moreoptions of basis functions, Bayesian functional data regression using GPs, and functionalclassification into
BFDA . References
Crainiceanu CM, Goldsmith AJ (2010). “Bayesian Functional Data Analysis Using
WinBUGS .” Journal of Statistical Software , (11). doi:10.18637/jss.v032.i11 . URL .Dawid AP (1981). “Some Matrix-variate Distribution Theory: Notational Considerations anda Bayesian Application.” Biometrika , (1), 265–274.de Boor C (1977). “Computational Aspects of Optimal Recovery.” In CA Micchelli, TJ Rivlin(eds.), Optimal Estimation in Approximation Theory , pp. 69–91. Springer US, Boston,MA. ISBN 978-1-4684-2388-4. doi:10.1007/978-1-4684-2388-4_3 . URL http://dx.doi.org/10.1007/978-1-4684-2388-4_3 . ournal of Statistical Software R Package fda.usc .” Journal of Statistical Software , (1), 1–28. ISSN1548-7660. doi:10.18637/jss.v051.i04 . URL .Gaffney PW, Powell MJD (1976). Optimal interpolation . Springer Berlin Heidelberg.Gelman A, Rubin DB (1992). “Inference from Iterative Simulation Using Multiple Sequences.”
Statistical Science , (4), 457–472. doi:10.1214/ss/1177011136 . URL http://dx.doi.org/10.1214/ss/1177011136 .Geman S, Geman D (1984). “Stochastic relaxation, Gibbs distributions, and the Bayesianrestoration of images.” Pattern Analysis and Machine Intelligence, IEEE Transactions on , , 721–741.Graves S, Hooker G, Ramsay J (2010). Functional Data Analysis with R and MATLAB .Springer-Verlag, New York.Green PJ, Silverman BW (1993).
Nonparametric regression and generalized linear models: aroughness penalty approach . CRC Press.Hunyadi L (2010). bspline : Draw, manipulate and reconstruct B-splines.
MATLAB package,URL .James M (1978). “The Generalised Inverse.”
The Mathematical Gazette , pp. 109–114.Micchelli CA, Rivlin TJ, Winograd S (1976). “The optimal recovery of smooth functions.”
Numerische Mathematik , (2), 191–200.Ramsay JO (2014). fdaM : Functional Data Analysis . MATLAB package, URL .Ramsay JO, Dalzell C (1991). “Some tools for functional data analysis.”
Journal of theRoyal Statistical Society B (Methodological) , (3), 539–572. ISSN 00359246. URL .Ramsay JO, Silverman BW (2002). Applied Functional Data Analysis: Methods and CaseStudies , volume 77. Springer-Verlag, New York.Ramsay JO, Silverman BW (2005).
Functional data analysis . Springer Series in Statistics,2nd edition. Springer-Verlag, New York.Ramsay JO, Wickham H, Graves S, Hooker G (2014). fda : Functional Data Analysis . R package version 2.4.4, URL https://CRAN.R-project.org/package=fda .S¨arkk¨a S, Aki V (2014). “MCMC Diagnostics for MATLAB .” http: // becs. aalto. fi/ en/research/ bayes/ mcmcdiag/ .Shi JQ, Cheng Y (2014). GPFDA : Apply Gaussian Process in Functional data analysis . R package version 2.2, URL https://CRAN.R-project.org/package=GPFDA .Shi JQ, Choi T (2011). Gaussian Process Regression Analysis for Functional Data . CRCPress, Boca Raton, FL.8Sturtz S, Ligges U, Gelman A (2005). “R2WinBUGS: A Package for Running
WinBUGS from R .” Journal of Statistical Software , (1), 1–16. ISSN 1548-7660. doi:10.18637/jss.v012.i03 . URL .Vieu P, Ferraty F (2006). Nonparametric Functional Data Analysis: Theory and Practice .Springer Science & Business Media.Yang J, Cox DD, Lee JS, Ren P, Choi T (2015). “Efficient Bayesian hierarchical functionaldata analysis with basis function approximations using Gaussian-Wishart processes.” eprintarXiv:1512.07568 .Yang J, Zhu H, Choi T, Cox DD (2016). “Smoothing and Meanˆa ˘A¸SCovariance Estimation ofFunctional Data with a Bayesian Hierarchical Model.”
Bayesian Analysis , (3), 649–670. doi:10.1214/15-BA967 . URL http://dx.doi.org/10.1214/15-BA967 .Yao F, M¨uller HG, Wang JL (2005a). “Functional Data Analysis for Sparse LongitudinalData.” Journal of the American Statistical Association , (470), 577–590.Yao F, M¨uller HG, Wang JL (2005b). “Functional Linear Regression Analysis for LongitudinalData.” The Annals of Statistics , (6), 2873–2903.Yao F, M¨uller HG, Wang JL (2015). PACE
Package for Functional Data Analysis andEmpirical Dynamics (
MATLAB ) . Version 2.17, URL .Yuan Y, Johnson VE (2012). “Goodness-of-Fit Diagnostics for Bayesian Hierarchical Models.” Biometrics , (1), 156–164. ISSN 1541-0420. doi:10.1111/j.1541-0420.2011.01668.x .URL http://dx.doi.org/10.1111/j.1541-0420.2011.01668.x . Appendices
A. Inputs and outputs
A.1. Input variables
The main function
BFDA() has three input arguments: • A cell contains all functional data; • A cell contains all grids on which the functional data are collected; • A parameter structure outputted by function setOptions_bfda() , containing allrequired parameters:9 – smethod , specifying the method used for analyzing the functional data. Defaultvalue is ’babf’ for BABF method with basis function approximation; otherchoices are ’bhm’ for BHM method without basis function approximation, ’bgp’ for standard Bayesian GP regression, ’bfpca’ for Bayesian principal componentsanalysis; – Burnin , the number of burn-ins for the MCMC algorithm. Default value is ; – M , the number of iterations for the MCMC algorithm. Default value is ; – cgrid , set as if the functional data are observed on a common-grid, otherwiseset as for uncommon or random grids. Default value is ; – Sigma_est , estimated smooth covariance matrix from previous analysis. Defaultis empty and will be estimated by
PACE or sample estimate from individuallysmoothed data; – mu_est , estimated smooth functional mean from previous analysis. Default isempty and will be set as the smoothed sample mean; – mat , set as to use the Mate´rn covariance function as prior structure for stationaryfunctional data; set as to use the empirical covariance estimate Sigma_est as theprior structure for nonstationary functional data. Default value is ; – nu , order of smoothness for the Mate´rn covariance function. Default is empty andwill be estimated based on Sigma_est ; – delta , shape parameter δ of the IWP. Default is for a non-informative prior; – c , determining the prior covariance for functional mean. Default is ; – w, ws , determining the prior gamma distributions for σ (cid:15) and σ s . Defaults are w =1, ws =0.1. The parameter ws should be tuned for a proper magnitude of theposterior covariance estimate; – pace , if Sigma_est and mu_est are empty, set pace =1 to obtain
Sigma_est and mu_est by PACE , and set pace =0 to use the empirical estimates from theindividually smoothed data by CSS. Default is ; – m, tau , working grid tau is only required for ’babf’ method. Default is emptyand will be set up as the (0 : m − : 100) percentiles of the pooled observation gridwith length m ; – eval_grid , evaluation grid for all functional estimates, only required for ’babf’ methods; – lamb_min, lamb_max, lamb_step , determining the smoothing parametercandidates for general cross validation of the CSS method. Defaults are lamb_min =0.9, lamb_max =0.99, lamb_step =0.01; – a, b , hyper parameters for the gamma distributions in ’bgp’ , and ’bfpca’ . – resid_thin , determine the MCMC thinning steps of the residuals that are usedto test the goodness-of-fit of the model. Default is . A.2. Output variables
The main function
BFDA() has two output arguments, one structure outputted by the specifiedmethod, and the other parameter structure as specified by setOptions_bfda() containingupdated parameter values.0Output structure with smethod = ’bhm’ : • Z, Z_CL, Z_UL , smoothed functional data, lower and upper 95% credible intervals; • Sigma, Sigma_CL, Sigma_UL , functional covariance estimate, lower and upper 95%credible intervals; • Sigma_SE , the empirical covariance estimate by using the smoothed data Z ; • mu, mu_CI , functional mean estimate, 95% credible intervals; • rn, rn_CI , estimate and 95% credible interval for the noise precision; • rs, rs_CI , estimate and 95% credible interval for σ s ; • rho, nu , estimated parameter values for the Mat´ern function; • residuals , MCMC samples of the residuals that are used to test the goodness-of-fit; • pmin_vec , p-values for testing the goodness-of-fit for all functional samples. P-value > < p-value < < smethod = ’babf’ has the following variables that are differentfrom the ones with smethod = ’bhm’ : • Zt , smoothed functional data on the observation grids; • Z_cgrid, Z_cgrid_CL, Z_cgrid_UL , smoothed functional data on the evaluation grid eval_grid , along with lower and upper 95% credible intervals; • Sigma_cgrid, Sigma_cgrid_CL, Sigma_cgrid_UL , functional covariance estimate onthe evaluation grid eval_grid , along with lower and upper 95% credible intervals; • mu_cgrid, mu_cgrid_CI , functional mean estimate on the evaluation grid eval_grid ,along with 95% credible intervals; • Zeta, Zeta_CL, Zeta_UL , estimates for the coefficients of basis functions, along withlower and upper 95% credible intervals; • Sigma_zeta_SE , the empirical covariance estimate with the estimated
Zeta ; • Sigma_zeta, Sigma_zeta_CL, Sigma_zeta_UL , covariance estimate for the coefficientsof basis functions, along with lower and upper 95% credible intervals; • mu_zeta, mu_zeta_CI , mean estimate for the coefficients of basis functions, along with95% credible intervals; • Btau , the basis function evaluations on the working grid tau ; • BT , the basis function evaluations on the observation grids; • Sigma_tau , functional covariance estimate on the working grid tau ;1 • mu_tau , functional mean estimate on the working grid tau ; • optknots , the optimal knots selected by optknt() for evaluations on the working grid tau . B. Example
MATLAB scripts for using
BFDA %% -------- Add pathes of the required MATLAB packages --------% BFDA, bspline, fdaM, mcmcdiag, PACE% replace pwd by the directory of your MATLAB packagesaddpath(genpath(cat(2, pwd, ' /BFDA ' )))addpath(genpath(cat(2, pwd, ' /bspline ' )))addpath(genpath(cat(2, pwd, ' /fdaM ' )))addpath(genpath(cat(2, pwd, ' /mcmcdiag ' )))addpath(genpath(cat(2, pwd, ' /PACErelease2.11 ' )))%% -------- Set up parameters for simulation --------n = 30; % Number of functional samplesp = 40; % Number of pooled grid points, or evaluated grid pointss = sqrt(5); % Standard deviation of functional observationsr = 2; % Signal to noise ratiorho = 1/2; % Scale parameter in the Matern functionnu = 3.5; % Order parameter in the Matern functionpgrid = (0 : (pi/2)/(p-1) : (pi/2)); % Pooled griddense = 0.6; % Proportion of observations on the pooled gridau = 0; bu = pi/2; % Function domainm = 20; % Number of working grid pointsstat = 1; % Specify stationary datacgrid = 1; % Specify common observation grid%% -------- Analyzing stationary functional data with common grid --------% Generate simulated data from GP(3sin(4t), s^2 Matern_cor(d; rho, nu))% with noises from N(0, (s/r)^2)GausFD_cgrid = sim_gfd(pgrid, n, s, r, nu, rho, dense, cgrid, stat);% setup parameters for BFDA% run with BHMparam = setOptions_bfda( ' smethod ' , ' bhm ' , ' cgrid ' , 1, ' mat ' , 1, ... ' M ' , 10000, ' Burnin ' , 2000, ' w ' , 1, ' ws ' , 1);% run with Bayesian Functional PCA% param = setOptions_bfda( ' smethod ' , ' bfpca ' , ' M ' , 50, ' Burnin ' , 20); % run with standard Bayesian Gaussian Process model% param = setOptions_bfda( ' smethod ' , ' bgp ' , ' mat ' , 1, ...% ' M ' , 50, ' Burnin ' , 20);% run with Cubic Smoothing Splines% param = setOptions_bfda( ' smethod ' , ' css ' , ' mat ' , 1, ' M ' , ...% 50, ' Burnin ' , 20, ' pace ' , 0);% call BFDA[out_cgrid, param ] = ...BFDA(GausFD_cgrid.Xraw_cell, GausFD_cgrid.Tcell, param);%% -------- Analyzing stationary functional data with uncommon grid --------GausFD_ucgrid = sim_gfd(pgrid, n, s, r, nu, rho, dense, 0, stat);param_uc = setOptions_bfda( ' smethod ' , ' bhm ' , ' cgrid ' , 0, ' mat ' , 1, ' M ' ,...10000, ' Burnin ' , 2000, ' pace ' , 1, ' ws ' , 0.1);[out_ucgrid, param_uc] = ...BFDA(GausFD_ucgrid.Xraw_cell, GausFD_ucgrid.Tcell, param_uc);%% -------- Analyzing non-stationary functional data with common grid --------GausFD_cgrid_ns = sim_gfd(pgrid, n, s, r, nu, rho, dense, cgrid, 0);param_ns = setOptions_bfda( ' smethod ' , ' bhm ' , ' cgrid ' , 1, ' mat ' , 0, ' M ' ,...10000, ' Burnin ' , 2000, ' pace ' , 1, ' ws ' , 0.01);[out_cgrid_ns, param_ns] = ...BFDA(GausFD_cgrid_ns.Xraw_cell, GausFD_cgrid_ns.Tcell, param_ns);%% -------- Analyzing non-stationary functional data with uncommon grid --------GausFD_ucgrid_ns = sim_gfd(pgrid, n, s, r, nu, rho, dense, 0, 0);param_uc_ns = setOptions_bfda( ' smethod ' , ' bhm ' , ' cgrid ' , 0, ' mat ' , 0, ... ' M ' , 10000, ' Burnin ' , 2000, ' pace ' , 1, ' ws ' , 0.01);[out_ucgrid_ns, param_uc_ns ] = ...BFDA(GausFD_ucgrid_ns.Xraw_cell, GausFD_ucgrid_ns.Tcell, param_uc_ns);%% -------- Analyzing stationary functional data with random grids --------GausFD_rgrid = sim_gfd_rgrid(n, p, au, bu, s, r, nu, rho, stat);param_rgrid = setOptions_bfda( ' smethod ' , ' babf ' , ' cgrid ' , 0, ' mat ' , 1, ... ' M ' , 10000, ' Burnin ' , 2000, ' m ' , m, ' eval_grid ' , pgrid, ' ws ' , 1, ... ' trange ' , [au, bu]);% call BFDA[out_rgrid, param_rgrid]= ...BFDA(GausFD_rgrid.Xraw_cell, GausFD_rgrid.Tcell, param_rgrid);%% -------- Analyzing nonstationary functional data with random grids--------GausFD_rgrid_ns = sim_gfd_rgrid(n, p, au, bu, s, r, nu, rho, 0);param_rgrid_ns = setOptions_bfda( ' smethod ' , ' babf ' , ' cgrid ' , 0, ' mat ' , ...0, ' M ' , 10000, ' Burnin ' , 2000, ' m ' , m, ' eval_grid ' , pgrid, ' ws ' , 0.05, ... ' trange ' , [au, bu]);% call BFDA[out_rgrid_ns, param_rgrid_ns] = ...BFDA(GausFD_rgrid_ns.Xraw_cell, GausFD_rgrid_ns.Tcell, param_rgrid_ns);%% -------- Calculate RMSE (root mean square error) --------display( ' RMSE of the estimated stationary covariance ' )rmse(out_cgrid.Sigma_SE, GausFD_cgrid.Cov_true)display( ' RMSE of the estimated functional data ' )Xtrue_mat = reshape(cell2mat(GausFD_cgrid.Xtrue_cell), p, n);rmse(out_cgrid.Z, Xtrue_mat)display( ' RMSE of the estimated non-stationary covariance ' )Ctrue_ns = cov_ns(pgrid, sf, nu, rho);% calculate the true non-stationary covariance matrixrmse(out_cgrid_ns.Sigma_SE, Ctrue_ns)%% -------- Save simulated data sets and BFDA results --------save( ' ./Examples/Data/Simu_Data.mat ' , ' GausFD_cgrid ' , ' GausFD_ucgrid ' , ... ' GausFD_cgrid_ns ' , ' GausFD_ucgrid_ns ' , ... ' GausFD_rgrid ' , ' GausFD_rgrid_ns ' )save( ' ./Examples/Data/Simu_Output.mat ' , ' out_cgrid ' , ' out_ucgrid ' , ... ' out_cgrid_ns ' , ' out_ucgrid_ns ' , ... ' out_rgrid ' , ' out_rgrid_ns ' ) C. Example
MATLAB scripts for functional regression using fdaM %% -------- Add fdaM path and load the functional data -------- % Replace pwd by the directory of your MATLAB packagesaddpath(genpath(cat(2, pwd, ' /fdaM ' )));load( ' ./Examples/Data/Simu_Data.mat ' );load( ' ./Examples/Data/Simu_Output.mat ' );%% -------- Set sample sizes, training data set, and test data set --------n = 30; % Number of functional curvesp = 40; % Number of pooled grid points, or evaluated grid pointsau = 0; bu = pi/2; % Domain of tpgrid = (au : (bu)/(p-1) : bu); % Pooled gridtrange = [au, bu];sampind = sort(randsample(1:n,20,false)) ;samptest = find(~ismember(1:n, sampind));n_train = length(sampind); n_test = length(samptest);cgrid = 0;Xtrue = zeros(p, n);Xraw = zeros(p, n);Xsmooth = zeros(p, n);if cgrid% Functional observations with common gridsXtrue = reshape(cell2mat(GausFD_cgrid.Xtrue_cell), p, n);Xsmooth = out_cgrid.Z(:, sampind);Xraw = reshape(cell2mat(GausFD_cgrid.Xraw_cell), p, n);else% Functional observations with random gridsfor i = 1:nxi = GausFD_rgrid.Xtrue_cell{i};xrawi = GausFD_rgrid.Xraw_cell{i};ti = GausFD_rgrid.Tcell{i};%interpolate by cubic smoothing splineh = mean(diff(ti));Xtrue(:, i) = csaps(ti, xi, 1/(1 + h^3/6), pgrid);Xraw(:, i) = csaps(ti, xrawi, 1/(1 + h^3/6), pgrid);zi = out_rgrid.Zt{i};Xsmooth(:, i) = csaps(ti, zi, 1/(1 + h^3/6), pgrid);end% Xsmooth = out_rgrid.Z_cgrid;% Xsmooth(1, :) = Xsmooth(2, :) * 0.9;endXtrain = Xsmooth(:, sampind);Xtest = Xsmooth(:, samptest); Xraw_train = Xraw(:, sampind);Xraw_test = Xraw(:, samptest);rmse(Xtrue, Xsmooth)rmse(Xtrue, Xraw)%% -------- Generate response variables --------betamat = (pgrid ' ) .^ 2 ;% -------- Scalar responesdeltat = pgrid(2)-pgrid(1);Avec_true = deltat.*(Xtrue ' *betamat - ...0.5.*(Xtrue(1, :) ' .*betamat(1) + Xtrue(p,:) ' .*betamat(p)) );Avec = Avec_true + normrnd(0, 1, n, 1);Avec_train = Avec(sampind);Avec_test = Avec(samptest);Avec_train_true = Avec_true(sampind);Avec_test_true = Avec_true(samptest);% -------- Functional responsesymat_true = Xtrue .* repmat(betamat, 1, n) ;ymat = ymat_true + normrnd(0, 1, p, n);ymat_train = ymat(:, sampind);ymat_test = ymat(:, samptest);ymat_train_true = ymat_true(:, sampind);ymat_test_true = ymat_true(:, samptest);%% -------- Set up functional data structure xfd, yfd for fdaM ----xnbasis = 20;xbasis = create_bspline_basis(trange, xnbasis, 4);xfd_true = smooth_basis(pgrid, Xtrue, xbasis);xfd = smooth_basis(pgrid, Xtrain, xbasis);xfd_raw = smooth_basis(pgrid, Xraw_train, xbasis);[yfd_samp, df, gcv, beta, SSE, penmat, y2cMap, argvals, y] = ...smooth_basis(pgrid, ymat_train, xbasis);%% -------- Set up the curvature penalty operator -------conbasis = create_constant_basis(trange); % create a constant basiswfd = fd([0, 1], conbasis);wfdcell = fd2cell(wfd);curvLfd = Lfd(2, wfdcell); % Set up xfdcellxfdcell = cell(1, 2);xfdcell{1} = fd(ones(1, n_train), conbasis);xfdcell{2} = xfd;xfd_raw_cell = cell(1, 2);xfd_raw_cell{1} = fd(ones(1, n_train), conbasis);xfd_raw_cell{2} = xfd_raw;% Set up betacell for scalar responsesbetafd0 = fd(0, conbasis);bnbasis = 10;betabasis = create_bspline_basis(trange, bnbasis, 4);betafd1 = fd(zeros(bnbasis, 1), betabasis);betacell_vecy = cell(1, 2);betacell_vecy{1} = fdPar(betafd0);betacell_vecy{2} = fdPar(betafd1, curvLfd, 0);% Set up betacell, yfd_par for functional responsesbetacell_fdy = cell(1, 2);betacell_fdy{1} = fdPar(betafd1, curvLfd, 0);betacell_fdy{2} = fdPar(betafd1, curvLfd, 0);yfd_par = fdPar(yfd_samp, curvLfd, 0);%% ---- Compute cross-validated SSE ' s for a range of smoothing parameters ----%{wt = ones(1, length(sampind));lam = (0:0.1:1);nlam = length(lam);SSE_CV_vecy = zeros(nlam,1);SSE_CV_raw_vecy = zeros(nlam, 1);SSE_CV_fdy = zeros(nlam,1);SSE_CV_raw_fdy = zeros(nlam, 1);for ilam = 1:nlam;lambda_vecy = lam(ilam);betacelli_vecy = betacell_vecy;betacelli_vecy{2} = putlambda(betacell_vecy{2}, lambda_vecy);SSE_CV_vecy(ilam) = fRegress_CV(Avec_train, xfdcell, betacelli_vecy, wt);fprintf( ' Scalar responses, lambda %6.2f: SSE = %10.4f\n ' , ...lam(ilam), SSE_CV_vecy(ilam)); SSE_CV_raw_vecy(ilam) = fRegress_CV(Avec_train, xfd_raw_cell, betacelli_vecy, wt);fprintf( ' Scalar responses, lambda %6.2f: SSE = %10.4f\n ' , ...lam(ilam), SSE_CV_raw_vecy(ilam));betacelli_fdy = betacell_fdy;betacelli_fdy{1} = putlambda(betacell_fdy{1}, lambda_vecy);betacelli_fdy{2} = putlambda(betacell_fdy{2}, lambda_vecy);yfd_par_i = putlambda(yfd_par, lambda_vecy);SSE_CV_fdy(ilam) = fRegress_CV(yfd_par_i, xfdcell, betacelli_fdy, wt);fprintf( ' Functional respones, lambda %6.2f: SSE = %10.4f\n ' , ...lam(ilam), SSE_CV_fdy(ilam));SSE_CV_raw_fdy(ilam) = fRegress_CV(yfd_par_i, xfd_raw_cell, betacelli_fdy, wt);fprintf( ' Functional respones, lambda %6.2f: SSE = %10.4f\n ' , ...lam(ilam), SSE_CV_raw_fdy(ilam));end%}%% -------- Fit the linear model --------lambda = 0.1;wt = ones(1, length(sampind));% --------- Scalar responsesbetacell_vecy{2} = fdPar(betafd1, curvLfd, lambda);fRegressStruct_vecy = fRegress(Avec_train, xfdcell, betacell_vecy, wt);fRegressStruct_raw_vecy = ...fRegress(Avec_train, xfd_raw_cell, betacell_vecy, wt);% Get coefficientsbetaestcell_vecy = fRegressStruct_vecy.betahat;Avec_hat = fRegressStruct_vecy.yhat;intercept_vecy = getcoef(getfd(betaestcell_vecy{1}));disp([ ' Constant term = ' ,num2str(intercept_vecy)])betaestcell_raw_vecy = fRegressStruct_raw_vecy.betahat;Avec_hat_raw = fRegressStruct_raw_vecy.yhat;intercept_raw_vecy = getcoef(getfd(betaestcell_raw_vecy{1}));disp([ ' Constant term = ' ,num2str(intercept_raw_vecy)])display([ ' Scalar reponses: ' , ' fitted mse = ' , ...num2str(mse(Avec_train_true, Avec_hat)), ... ' ; fitted mse_raw = ' ,num2str(mse(Avec_train_true, Avec_hat_raw))]) % Compute Rsquarecovmat = cov([Avec_train, Avec_hat]);Rsqrd = covmat(1,2)^2/(covmat(1,1)*covmat(2,2));disp([ ' R-squared = ' ,num2str(Rsqrd)])covmat_raw = cov([Avec_train, Avec_hat_raw]);Rsqrd_raw = covmat_raw(1,2)^2/(covmat_raw(1,1)*covmat_raw(2,2));disp([ ' raw R-squared = ' ,num2str(Rsqrd_raw)])% Compute sigmaresid_vecy = Avec_train - Avec_hat;SigmaE_vecy = mean(resid_vecy.^2);disp([ ' Scalar responses: SigmaE = ' ,num2str(SigmaE_vecy)])resid_raw_vecy = Avec_train - Avec_hat_raw;SigmaE_raw_vecy = mean(resid_raw_vecy.^2);disp([ ' Scalar responses: Raw SigmaE = ' ,num2str(SigmaE_raw_vecy)])% ---------- Functional responsesbetacell_fdy{1} = fdPar(betafd1, curvLfd, lambda);betacell_fdy{2} = fdPar(betafd1, curvLfd, lambda);yfd_par = fdPar(yfd_samp, curvLfd, lambda);fRegressStruct_fdy = fRegress(yfd_par, xfdcell, betacell_fdy, wt, y2cMap);fRegressStruct_raw_fdy = ...fRegress(yfd_par, xfd_raw_cell, betacell_fdy, wt, y2cMap);betaestcell_fdy = fRegressStruct_fdy.betahat;yfd_hat = fRegressStruct_fdy.yhat;intercept_fdy = eval_fd(pgrid, getfd(betaestcell_fdy{1}));betaestcell_raw_fdy = fRegressStruct_raw_fdy.betahat;yfd_hat_raw = fRegressStruct_raw_fdy.yhat;intercept_raw_fdy = eval_fd(pgrid, getfd(betaestcell_raw_fdy{1}));% MSE of fitted responsesymat_fitted = eval_fd(pgrid, yfd_hat);ymat_fitted_raw = eval_fd(pgrid, yfd_hat_raw);display([ ' mse = ' , num2str(mse(ymat_train_true, ymat_fitted)), ... ' ; mse_raw = ' ,num2str(mse(ymat_train_true, ymat_fitted_raw))])% Compute squared residual correlationresid_fdy = ymat_train_true - ymat_fitted;SigmaE_fdy = cov(resid_fdy ' );resid_raw_fdy = ymat_train_true - ymat_fitted_raw; SigmaE_raw_fdy = cov(resid_raw_fdy ' );%% -------- Recompute the analysis to get confidence limits --------% ------- Scalar responsesstderrStruct_vecy = fRegress_stderr(fRegressStruct_vecy, eye(n_train), SigmaE_vecy);betastderrcell_vecy = stderrStruct_vecy.betastderr;stderrStruct_raw_vecy = ...fRegress_stderr(fRegressStruct_raw_vecy, eye(n_train), SigmaE_raw_vecy);betastderrcell_raw_vecy = stderrStruct_raw_vecy.betastderr;% Constant coefficient standard error:intercept_std_vecy = getcoef(betastderrcell_vecy{1});intercept_ste_raw_vecy = getcoef(betastderrcell_raw_vecy{1});% -------- Functional responsesstderrStruct_fdy = fRegress_stderr(fRegressStruct_fdy, y2cMap, SigmaE_fdy);% fixed a bug in fRegress_stderr.m at line 124:% bstderrfdj = data2fd(bstderrj, tfine, betabasisj); should be% bstderrfdj = data2fd(tfine, bstderrj, betabasisj);betastderrcell_fdy = stderrStruct_fdy.betastderr;stderrStruct_raw_fdy = ...fRegress_stderr(fRegressStruct_raw_fdy, y2cMap, SigmaE_raw_fdy);betastderrcell_raw_fdy = stderrStruct_raw_fdy.betastderr;% Coefficient standard error:intercept_std_fdy = eval_fd(pgrid, betastderrcell_fdy{1});intercept_std_raw_fdy = eval_fd(pgrid, betastderrcell_raw_fdy{1});%% -------- Predict on test data --------%Set up xfd for test dataxfd_test = smooth_basis(pgrid, Xtest, xbasis);xfd_raw_test = smooth_basis(pgrid, Xraw_test, xbasis);% -------- Scalar responsesxfdcell_test = cell(1, 2);xfdcell_test{1} = fd(ones(1, n_test), conbasis);xfdcell_test{2} = xfd_test;xfd_raw_test_cell = cell(1, 2);xfd_raw_test_cell{1} = fd(ones(1, n_test), conbasis);xfd_raw_test_cell{2} = xfd_raw_test;Avec_pred = fRegressPred(xfdcell_test, betaestcell_vecy); Avec_pred_raw = fRegressPred(xfd_raw_test_cell, betaestcell_raw_vecy);display([ ' Scalar responses predict mse = ' , ...num2str(mse(Avec_test_true, Avec_pred)), ... ' ; Scalar responses with raw data predict mse_raw = ' ,...num2str(mse(Avec_test_true, Avec_pred_raw))])% -------- Functional responsesymat_test_pred = ...eval_fd(pgrid, fRegressPred(xfdcell_test, betaestcell_fdy, xbasis));ymat_test_pred_raw = ...eval_fd(pgrid, fRegressPred(xfd_raw_test_cell, betaestcell_raw_fdy, xbasis));display([ ' Functional response prediction mse = ' , ...num2str(mse(ymat_test_true, ymat_test_pred)), ... ' ; Functional responses prediction with Raw data mse_raw = ' ,...num2str(mse(ymat_test_true, ymat_test_pred_raw))]) Affiliation:
Jingjing YangDepartment of BiostatisticsSchool of Public HealthUniversity of Michigan1415 Washington HeightsAnn Arbor, MI, 48109, U.S.A.E-mail: [email protected]; [email protected]
Journal of Statistical Software published by the Foundation for Open Access Statistics
MMMMMM YYYY, Volume VV, Issue II
Submitted: yyyy-mm-dd doi:10.18637/jss.v000.i00