Science, Faculty of
Statistics, Department of
DSpace
UBCV
Pedroso Estevam de Souza, Camila
Doctor of Philosophy - PhD
University of British Columbia
In this thesis, we propose a methodology to analyze data arising from a curve that, over its domain, switches among J states. We consider a sequence of response variables, where each response y depends on a covariate x according to an unobserved state z, also called a hidden or latent state. The states form a stochastic process and their possible values are j=1,...,J. If z equals j the expected response of y is one of J unknown smooth functions evaluated at x. We call this model a switching nonparametric regression model. In a Bayesian switching nonparametric regression model the uncertainty about the functions is formulated by modeling the functions as realizations of stochastic processes. In a frequentist switching nonparametric regression model the functions are merely assumed to be smooth. We consider two different data structures: one with N replicates and the other with one single realization. For the hidden states, we consider those that are independent and identically distributed and those that follow a Markov structure. We develop an EM algorithm to estimate the parameters of the latent state process and the functions corresponding to the J states. Standard errors for the parameter estimates of the state process are also obtained. We investigate the frequentist properties of the proposed estimates via simulation studies. Two different applications of the proposed methodology are presented. In the first application we analyze the well-known motorcycle data in an innovative way: treating the data as coming from J>1 simulated accident runs with unobserved run labels. In the second application we analyze daytime power usage on business days in a building treating each day as a replicate and modeling power usage as arising from two functions, one function giving power usage when the cooling system of the building is off, the other function giving power usage when the cooling system is on.
Switching Nonparametric Regression ModelsbyCamila Pedroso Estevam de SouzaBachelor?s Degree in Statistics, University of Campinas, Brazil, 2005Master?s Degree in Statistics, University of Campinas, Brazil, 2008A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFDoctor of PhilosophyinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Statistics)The University Of British Columbia(Vancouver)September 2013c? Camila Pedroso Estevam de Souza, 2013AbstractIn this thesis, we propose a methodology to analyze data arising from a curve that,over its domain, switches among J states. We consider a sequence of responsevariables, where each response y depends on a covariate x according to an unob-served state z, also called a hidden or latent state. The states form a stochasticprocess and their possible values are j = 1, . . . ,J. If z equals j the expected re-sponse of y is one of J unknown smooth functions evaluated at x. We call thismodel a switching nonparametric regression model. In a Bayesian switching non-parametric regression model the uncertainty about the functions is formulated bymodeling the functions as realizations of stochastic processes. In a frequentistswitching nonparametric regression model the functions are merely assumed to besmooth. We consider two different data structures: one with N replicates and theother with one single realization. For the hidden states, we consider those that areindependent and identically distributed and those that follow a Markov structure.We develop an EM algorithm to estimate the parameters of the latent state processand the functions corresponding to the J states. Standard errors for the parameterestimates of the state process are also obtained. We investigate the frequentist prop-erties of the proposed estimates via simulation studies. Two different applicationsof the proposed methodology are presented. In the first application we analyzethe well-known motorcycle data in an innovative way: treating the data as comingfrom J > 1 simulated accident runs with unobserved run labels. In the second ap-plication we analyze daytime power usage on business days in a building treatingeach day as a replicate and modeling power usage as arising from two functions,one function giving power usage when the cooling system of the building is off,the other function giving power usage when the cooling system is on.iiPrefaceThis thesis is an original work of the author, Camila Pedroso Estevam de Souza,under the supervision of Dr. Nancy Heckman.A version of Chapters 2 and 4 together with Section 6.1 was submitted to peerreview in May 2013 and also posted online at arxiv.org. This manuscript is entitled?Switching nonparametric regression models and the motorcycle data revisited? byDe Souza, C. P. E. and Heckman, E. N.. The idea and theory were jointly developedby myself and Dr. Heckman. I conducted all derivations and all computationalwork and the majority of the writing.Davor Cubranic and I developed the R package switchnpreg that is avail-able on CRAN (cran.r-project.org). The package can be used to obtain the param-eter estimates for the one realization case described in Chapter 2. I did the originalcoding and Davor re-structured my code in a modular way and put it in a formappropriate for an R package.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation and data set discussion . . . . . . . . . . . . . . . . . 11.2 Overview of the proposed methodology . . . . . . . . . . . . . . 71.3 One realization case . . . . . . . . . . . . . . . . . . . . . . . . . 91.4 Replicate case . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101.5 Literature review . . . . . . . . . . . . . . . . . . . . . . . . . . 102 Solution for the one realization case . . . . . . . . . . . . . . . . . . 142.1 E-step: general zi?s . . . . . . . . . . . . . . . . . . . . . . . . . 162.1.1 E-step: independent and identically distributed (iid) zi?s . 172.1.2 E-step: Markov zi?s . . . . . . . . . . . . . . . . . . . . . 182.2 M-step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19iv2.3 Standard errors for the parameter estimators of the state process . 232.3.1 The general case . . . . . . . . . . . . . . . . . . . . . . 232.3.2 Standard errors: iid zi?s . . . . . . . . . . . . . . . . . . . 252.3.3 Standard errors: Markov zi?s . . . . . . . . . . . . . . . . 263 Solution for the replicate case . . . . . . . . . . . . . . . . . . . . . . 283.1 E-step details . . . . . . . . . . . . . . . . . . . . . . . . . . . . 303.1.1 E-step: iid zik?s . . . . . . . . . . . . . . . . . . . . . . . 313.1.2 E-step: Markov zik?s . . . . . . . . . . . . . . . . . . . . 323.2 M-step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 333.2.1 ECM algorithm when Vzk = diag(?2z1k , . . . ,?2znk) . . . . . . 413.3 Choice of the smoothing parameters when Vzk = diag(?2z1k , . . . ,?2znk) 443.4 Standard errors for the parameter estimators of the state process . 473.4.1 Standard errors: iid zik?s . . . . . . . . . . . . . . . . . . 483.4.2 Standard errors: Markov zik?s . . . . . . . . . . . . . . . 494 Simulation studies: the one realization case . . . . . . . . . . . . . . 524.1 Simulated data . . . . . . . . . . . . . . . . . . . . . . . . . . . 524.2 Initial values . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 534.3 Choice of the smoothing parameters, the ? j?s . . . . . . . . . . . 554.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 565 Simulation studies: the replicate case . . . . . . . . . . . . . . . . . 695.1 Simulated data . . . . . . . . . . . . . . . . . . . . . . . . . . . 695.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 716 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 806.1 Motorcycle data . . . . . . . . . . . . . . . . . . . . . . . . . . . 806.1.1 Data set background . . . . . . . . . . . . . . . . . . . . 806.1.2 Data analysis . . . . . . . . . . . . . . . . . . . . . . . . 826.2 Building Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . 876.2.1 Results: uncorrelated observations . . . . . . . . . . . . . 876.2.2 Results: correlated observations . . . . . . . . . . . . . . 88v7 Discussion and future work . . . . . . . . . . . . . . . . . . . . . . . 98Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102A The proposed EM algorithm: why it works . . . . . . . . . . . . . . 105B The ECM algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . 108B.1 Formal definition . . . . . . . . . . . . . . . . . . . . . . . . . . 108B.2 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109C Calculation details: standard errors . . . . . . . . . . . . . . . . . . 111D Calculation details: ECM algorithm . . . . . . . . . . . . . . . . . . 114D.1 Step 1 ECM algorithm . . . . . . . . . . . . . . . . . . . . . . . 114D.2 Step 2 ECM algorithm . . . . . . . . . . . . . . . . . . . . . . . 120viList of TablesTable 4.1 Summary of the simulation studies. . . . . . . . . . . . . . . . 59Table 4.2 Number of simulated data sets with values of p?(zi = 1|yi)> 0.2when zi = 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . 59Table 4.3 Number of simulated data sets with values of p?(zi = 2|yi)> 0.2when zi = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . 59Table 4.4 Estimates of ?2 (true value = 5?10?5). . . . . . . . . . . . . 60Table 4.5 Estimates of the parameters of the z process. . . . . . . . . . . 60Table 5.1 Estimates of ?2 (true value = 1?10?5). . . . . . . . . . . . . 72Table 5.2 Estimates of ?2 (true value = 1?10?4). . . . . . . . . . . . . 72Table 5.3 Estimates of the parameters of the z process. . . . . . . . . . . 72Table 6.1 Results for J = 3 with corresponding fitted curves in Figures6.1a and 6.2a. . . . . . . . . . . . . . . . . . . . . . . . . . . 84Table 6.2 Results for J = 4 with corresponding fitted curves in Figures6.1b and 6.2b. . . . . . . . . . . . . . . . . . . . . . . . . . . 84Table 6.3 Results for iid zik?s and V=?2I with corresponding fitted curvesin Figures 6.3. . . . . . . . . . . . . . . . . . . . . . . . . . . 90Table 6.4 Results for iid zik?s and Vzk = diag(?2z1k , . . . ,?2znk) with corre-sponding fitted curves in Figures 6.4. . . . . . . . . . . . . . . 90Table 6.5 Results for Markov zik?s and V = ?2I with corresponding fittedcurves in Figures 6.5. . . . . . . . . . . . . . . . . . . . . . . 90Table 6.6 Results for Markov zik?s and Vzk = diag(?2z1k , . . . ,?2znk) with cor-responding fitted curves in Figures 6.6. . . . . . . . . . . . . . 90viiTable 6.7 Results for iid zik?s and V generated by a random interceptmodel with ?1 = ?2 = 0.1 and corresponding curves in Figures6.7. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91Table 6.8 Results for Markov zik?s and V generated by a random interceptmodel with ?1 = ?2 = 0.1 and corresponding curves in Figures6.8. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91viiiList of FiguresFigure 1.1 Power usage in June and July 2009 in a building monitored byPulse Energy. . . . . . . . . . . . . . . . . . . . . . . . . . . 4Figure 1.2 Daytime power usage from 9am to 4pm on business days inJune and July 2009 in the same building as Figure 1.1. Eachcurve corresponds to a different day. . . . . . . . . . . . . . . 5Figure 1.3 Motorcycle data. Head acceleration in g (one g ? 9.8m/s2)versus the time in milliseconds after impact. . . . . . . . . . . 6Figure 4.1 Simulation 1 (iid zi?s). Example of simulated data along withthe initial and final estimates of f1 and f2 obtained using in (a)the Bayesian approach and in (b) the penalized log-likelihoodapproach. The red dots correspond to z = 2 and the black dotsto z = 1. The solid curves correspond to the true functions f1and f2. The green curves are the initial estimates of f1 andf2 obtained using the function estimate method. The dashedcurves are the final estimates of f1 and f2. . . . . . . . . . . . 61Figure 4.2 Simulation 2 (Markov zi?s with a12 = 0.3 and a21 = 0.4). Ex-ample of simulated data along with the initial and final esti-mates of f1 and f2 as in Figure 4.1. . . . . . . . . . . . . . . 62Figure 4.3 Simulation 3 (Markov zi?s with a12 = 0.1 and a21 = 0.2). Ex-ample of simulated data along with the initial and final esti-mates of f1 and f2 as in Figure 4.1. . . . . . . . . . . . . . . 62ixFigure 4.4 Simulation 1 (iid zi?s). EMSE of both initial (dashed curves)and final (solid curves) estimates of f1 and f2 using the Bayesianand the penalized log-likelihood approaches. Plots (a) and (b)show the results for f1 and (c) and (d) the results for f2. . . . . 63Figure 4.5 Simulation 3 (Markov zi?s with a12 = 0.1 and a21 = 0.2). EMSEof both initial (dashed curves) and final (solid curves) estimatesof f1 and f2 as in Figure 4.4. . . . . . . . . . . . . . . . . . . 64Figure 4.6 Simulation 1 (iid zi?s). EMSE of the final estimates of (a) f1and (b) f2 using the Bayesian (dashed curves) and the penal-ized log-likelihood (solid curves) approaches. . . . . . . . . . 65Figure 4.7 Simulation 3 (Markov zi?s with a12 = 0.1 and a21 = 0.2). EMSEof the final estimates of (a) f1 and (b) f2 as in Figure 4.6. . . . 65Figure 4.8 Box plot of the estimates of ?2 for each simulation study usingthe Bayesian and penalized log-likelihood approaches. Thedashed lines correspond to the true value 5?10?5. . . . . . . 66Figure 4.9 Simulation 1 (iid zi?s). Box plot of the estimates of p1 obtainedusing the Bayesian and penalized log-likelihood approaches.The dashed line corresponds to the true value 0.7. . . . . . . . 67Figure 4.10 Simulation 2 (Markov zi?s with a12 = 0.3 and a21 = 0.4). Boxplots of the estimates of (a) a12 and (b) a21 obtained using theBayesian and the penalized log-likelihood approaches. Thedashed lines correspond to the true values 0.3 and 0.4, respec-tively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68Figure 4.11 Simulation 3 (Markov zi?s with a12 = 0.1 and a21 = 0.2). Boxplots of the estimates of (a) a12 and (b) a21 obtained using theBayesian and the penalized log-likelihood approaches. Thedashed lines correspond to the true values 0.1 and 0.2, respec-tively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68xFigure 5.1 Simulation 1 (iid zik?s with p1 = 0.5). Example of simulateddata along with f?1(x) and f?2(x). The gray dashed curves cor-respond to 20 out of the 100 generated replicates. The blackand red solid curves correspond to the true functions f1 andf2, respectively, evaluated only at x. The black and red dashedcurves correspond to f?1(x) and f?2(x), respectively. . . . . . . 73Figure 5.2 Simulation 1 (iid zik?s with p1 = 0.5). Example of four sim-ulated replicates from data set shown in Figure 5.1. The reddots correspond to z = 2 and the black dots to z = 1. The blackand red curves correspond to the true functions f1 and f2, re-spectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74Figure 5.3 Simulation 2 (Markov zik?s with a12 = 0.3 and a21 = 0.4). Ex-ample of simulated data along with f?1(x) and f?2(x) as in Figure5.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75Figure 5.4 Simulation 2 (Markov zik?s with a12 = 0.3 and a21 = 0.4). Ex-ample of four simulated replicates from data set shown in Fig-ure 5.3 as in Figure 5.2. . . . . . . . . . . . . . . . . . . . . . 76Figure 5.5 EMSE of f?1(x) (black curves) and f?2(x) (red curves) for (a)Simulation 1 and (b) Simulation 2. . . . . . . . . . . . . . . . 77Figure 5.6 Simulation 1 (iid zik?s). Box plots of the estimates of (a) ?2 and(b) ?2. The dashed lines correspond to the true values 10?5 and10?4, respectively. . . . . . . . . . . . . . . . . . . . . . . . 77Figure 5.7 Simulation 2 (Markov zik?s with a12 = 0.3 and a21 = 0.4). Boxplots of the estimates of (a) ?2 and (b) ?2 as in Figure 5.6. . . 78Figure 5.8 Simulation 1 (iid zik?s with p1 = 0.5). Box plot of the estimatesof p1. The dashed line corresponds to the true value 0.5. . . . 78Figure 5.9 Simulation 2 (Markov zik?s with a12 = 0.3 and a21 = 0.4). Boxplot of the estimates of pi1. The dashed line corresponds to thetrue value 0.5. . . . . . . . . . . . . . . . . . . . . . . . . . . 79Figure 5.10 Simulation 2 (Markov zik?s with a12 = 0.3 and a21 = 0.4). Boxplots of the estimates of (a) a12 and (b) a21. The dashed linescorrespond to the true values 0.3 and 0.4, respectively. . . . . 79xiFigure 6.1 Motorcycle data. Final function estimates (solid curves) ob-tained using the penalized log-likelihood approach for (a) J =3 and (b) J = 4. The gray dashed curves correspond to the ini-tial function estimates, which are constant shifts of a smooth-ing spline fit to all the data. . . . . . . . . . . . . . . . . . . . 85Figure 6.2 Motorcycle data. Final function estimates (solid curves) ob-tained using the Bayesian approach for (a) J = 3 and (b) J = 4.The gray dashed curves correspond to the initial function esti-mates, which are constant shifts of a smoothing spline fit to allthe data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86Figure 6.3 Building daytime power usage. Fitted function estimates (solidcurves) assuming iid zik?s and V=?2I. The gray dashed curvescorrespond to the replicates. The red and black dashed curvesare the initial function estimates. The colors red and black cor-respond to the condition chiller on and off, respectively. . . . 92Figure 6.4 Building daytime power usage. Fitted function estimates as-suming iid zik?s and Vzk = diag(?2z1k , . . . ,?2znk). Components ofthe plot are as in Figure 6.3. . . . . . . . . . . . . . . . . . . 93Figure 6.5 Building daytime power usage. Fitted function estimates as-suming Markov zik?s and V = ?2I. Components of the plot areas in Figure 6.3. . . . . . . . . . . . . . . . . . . . . . . . . . 94Figure 6.6 Building daytime power usage. Fitted function estimates as-suming Markov zik?s and Vzk = diag(?2z1k , . . . ,?2znk). Compo-nents of the plot are as in Figure 6.3. The green curve corre-sponds to the replicate where there is a transition from chilleron to off. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95Figure 6.7 Building daytime power usage. Fitted function estimates as-suming iid zik?s and V generated by a random intercept model.Components of the plot are as in Figure 6.3. . . . . . . . . . . 96Figure 6.8 Building daytime power usage. Fitted function estimates as-suming Markov zik?s and V generated by a random interceptmodel. Components of the plot are as in Figure 6.3. . . . . . . 97xiiAcknowledgmentsI would like to express my sincere gratitude to my supervisor Professor NancyHeckman for all her support, guidance and constant encouragement throughoutthe course of my PhD program. Professor Heckman is a truly great mentor andhas always provided me with constructive advice and valuable suggestions. Withher enthusiasm Professor Heckman has also encouraged me to be engaged in theStatistical community by attending and presenting my research work in workshopsand conferences. I am so proud and honored to be one of her students.I am thankful to Professors Paul Gustafson and Alexandre Bouchard-Co?te?,members of my supervisory committee, for pointing out several interesting aspectsof my research and also providing valuable advice. I also want to thank ProfessorJean Opsomer for acting as the external examiner, and Professors Ruben H. Zamarand Jane Z. Wang for serving as the university examiners.I would like to acknowledge Professors Ronaldo Dias and Nancy Garcia of theUniversity of Campinas, Brazil, for having encouraged me to advance my studiestowards my PhD at the University of British Columbia.I want to extend my gratitude to the faculty, staff and fellow graduate studentsof the department of Statistics for the support and great environment.I wish to thank Pulse Energy personnel for the time I spent there during myMitacs-Accelerate research internship and for providing the building data that mo-tivated my thesis work.I would like to express my heartfelt thanks to my beloved parents, Nivaldo andDeliria, and brother Clayton. Lastly and most importantly, I would like to thankmy husband Marcos for his love and unconditional support and our son Nicolas,who was born during my PhD and brought so much joy to our lives.xiiiDedicationTo my parents Nivaldo and DeliriaxivChapter 1Introduction1.1 Motivation and data set discussionThe motivation for my thesis work stems from my four-month Mitacs-Accelerateresearch internship with Pulse Energy in the year of 2009. Pulse Energy is a BritishColumbia based company founded in 2006 with the goal of using modern technol-ogy to reduce energy waste and improve operational efficiency in commercial andinstitutional buildings. The company developed an energy management softwarecalled PulseTM, that allows building managers and users to monitor energy con-sumption online. As a result, managers and occupants can take action in orderto improve the energy efficiency of their building. More information about PulseEnergy can be found on their website: www.pulseenergy.com.One of the company?s main goals is to calculate a building?s ?typical curve?of energy consumption, that is, its expected energy consumption as a function oftime and weather conditions. During my internship I worked with Pulse Energypersonnel providing advice on their existing methodology, which was very com-puter intensive and basically consisted of smoothing over an enormous number ofpoints. After the internship, I started exploiting simplifying structure in the data.I want to use a structure that is flexible, that does not not limit the shapes of theenergy consumption data.To more fully understand my approach to the problem, compare the plots inFigures 1.1 and 1.2. Figure 1.1 shows hourly power usage during the months of1June and July 2009 in a building monitored by Pulse Energy. On some days (hol-idays and weekends) energy usage is very close to zero. We can observe that onsome business days the energy usage is very high, approximately twice as much ason the other days. After some investigation we found out that those business dayswith high power consumption correspond to warm days, when the cooling system(also called the chiller) of the building was probably on. Figure 1.2 presents thebuilding daytime power usage from 9am to 4pm on business days in June and July2009. Several types of curves can be observed: one type corresponds to days whenthe cooling system was probably on and another type when the cooling system wasoff. We can also observe that on some days the chiller turned on in the middle ofthe day. On one day the chiller went on, off and then on again.Pulse Energy?s method considered the data in Figure 1.1, finding the ?typicalcurve? by averaging over a large number of data points, without considering thespecial structure we see in Figure 1.2. My method exploits the structure of Figure1.2, in two important ways: one, by treating each business day as a replicate andthe other, by modeling the building power usage as arising from two functions, onefunction giving power usage when the chiller is off, the other function giving powerusage when the chiller is on. The condition ?chiller on?/?off? at any particular timecannot be observed directly. It can only be inferred from the data, thus forming alatent process.In general, my method can be applied to data arising from a curve that, over itsdomain, switches among J states. The state at any particular point is determinedby a latent process. The state also determines a function. The main goal is toestimate the function corresponding to the jth state and the parameters of the latentprocess, along with some measure of accuracy. We consider two types of analyses:a single realization case and a replicate case. In the single realization case, weconsider one curve switching among J functions (e.g., Figure 1.1 for J = 3). In thereplicate case, we have N curves, called replicates, switching among J functions(e.g., Figure 1.2 with days as replicates and J = 2).I also apply the proposed methodology to another data set, the well-knownmotorcycle data set. The data set consists of 133 measurements of head accelera-tion taken through time in simulated motorcycle accidents. See Figure 1.3. Anal-yses appearing in the literature (e.g., Silverman, 1985, Ha?rdle and Marron, 1995,2Rasmussen and Ghahramani, 2002, Gijbels and Goderniaux, 2004 and Wood, 2011)treat the data as coming from one simulated accident. However, close examinationof the data suggests the measurements are from J > 1 accidents. In the discussionof Silverman (1985) Prof. A. C. Atkinson wrote ?inspection of [the Figure], partic-ularly in the region 30-40 ms, suggests that the data may not be a single time seriesbut are rather the superposition of, perhaps, three series?. Professor Silverman hadno specific information on this point but replied ?Professor Atkinson is right aboutthe motorcycle data in that they are the superposition of measurements made byseveral instruments but the data I have presented are precisely in the form theywere made available to me?. The data structure will most likely remain unclear asthe original report (Schmidt et al., 1981) seems to be no longer available.We analyze the motorcycle data as coming from J > 1 simulated accident runswith hidden (lost) run labels. In order to recover and estimate the functions corre-sponding to each of the accident runs we approach the motorcycle data in a differ-ent way and apply our proposed methodology to the data treating them as a curvethat switches among J > 1 functions with hidden function labels. We choose J us-ing an ad hoc Akaike information criterion (AIC). To our knowledge this is the firsttime that the motorcycle data are analyzed taking into account that they describemore than one simulated accident.30 10 20 30 40 50 6050100150200Time in days (June and July 2009)Power KW/hFigure 1.1: Power usage in June and July 2009 in a building monitored byPulse Energy.49 10 11 12 13 14 15 1680100120140160180200220Hour of dayPower KW/hFigure 1.2: Daytime power usage from 9am to 4pm on business days in Juneand July 2009 in the same building as Figure 1.1. Each curve corre-sponds to a different day.5lllll lllllllllllll lllllllllllllllllllllllllllllllllllllllllllll lllllllllllllllllllllllllllllllllllllllllllllllllllllllll lll llllllll10 20 30 40 50?100?50050TimeHead accelerationFigure 1.3: Motorcycle data. Head acceleration in g (one g? 9.8m/s2) versusthe time in milliseconds after impact.61.2 Overview of the proposed methodologyIn this section I describe a general version of the proposed model. Then in Sections1.3 and 1.4 I describe the one realization case and the replicate case, respectively.Suppose we have a sequence of response variables, y1, . . . ,yn, where yi dependson a covariate xi according to an unobserved state zi, also called a hidden or latentstate. The possible values of the states are j = 1, . . . ,J. If zi = j the expectedresponse of yi is f j(xi). We call this model a switching nonparametric regressionmodel. The objective of the proposed methodology is to estimate the f j?s as well asto estimate the regression error variance and the parameters governing the distribu-tion of the state process. We also obtain standard errors for the proposed parameterestimates of the state process. We consider two types of hidden states, those that areindependent and identically distributed and those that follow a Markov structure.The Markov structure would usually require x to be time.To obtain the parameter estimates we maximizeobserved data log-likelihood+P( f1, . . . , fJ,?1, . . . ,?J), (1.1)where the exact form of P( f1, . . . , fJ,?1, . . . ,?J) is chosen by the user. We considertwo forms of P( f1, . . . , fJ,?1, . . . ,?J) arising from two approaches:? a frequentist approach called penalized log-likelihood estimation;? a Bayesian approach where the posterior density is maximized.These two approaches are computationally similar and consist of using the Expectation-Maximization (EM) algorithm (McLachlan and Krishnan, 2008) to maximize (1.1).In the frequentist approach we assume that f ??j exists almost everywhere and?[ f ??j (x)]2dx < ? for j = 1, . . . ,J and maximize (1.1) withP( f1, . . . , fJ,?1, . . . ,?J) =?J?j=1? j?[ f ??j (x)]2dx. (1.2)The ? j?s are the so called smoothing parameters as they measure the relative impor-tance of fit to the data, as measured by the log-likelihood of the data, and smooth-ness of the f j?s, as quantified by the penalties?[ f ??j (x)]2dx.7The Bayesian approach requires a prior distribution for the parameters. Forthe f j?s, we consider a Gaussian process regression approach (Rasmussen andWilliams, 2006) with f1, . . . , fJ independent and f j a Gaussian process with meanfunction ? j(x)? 0 and covariance function, K, that depends on a vector of param-eters ? j. We place a non-informative prior on the other parameters. Therefore, inthe Bayesian approach we maximize (1.1) withP( f1, . . . , fJ,?1, . . . ,?J) = ?J2log(2pi)? 12J?j=1log |A(? j)|?12J?j=1f j(x)T A(? j)?1 f j(x), (1.3)where A(? j) is an n?n matrix with entries given by[A(? j)]lm = K(xl,xm;? j),with K known. In our applications and simulations we assume the following co-variance functionK(x, t;U j,s2j) = cov( f j(x), f j(t)) =U j exp[?(x? t)22s2j], (1.4)with parameter vector ? j = (U j,s j)?. These parameters control the amount ofsmoothness of each curve f j and, as in the penalized log-likelihood approach, arecalled smoothing parameters.We can extend our approach to consider ? j not equal to the zero function. Forinstance, we might take ? j(x) as a linear combination of known basis functions. Inthis case, in addition to ? j, we would have another parameter vector for f j, namelythe vector containing the basis functions? coefficients.In both approaches we can choose the values of the ? j?s automatically by cross-validation.We take into account two different data structures: one with N replicates andthe other with one realization. In the first case, we suppose we have N data sets,with the kth data set being yik, xik = xi, i = 1, . . . ,nk = n, where the yik?s are gen-8erated from a sequence of hidden states, zik, i = 1, . . . ,n. One could consider xi?sdepending on k, but for ease we fix them across replicates. In the case with one re-alization, we simply have N = 1. Although the one realization case is a particularcase of the replicates case we decided to treat them separately, as their applica-tions are conceptually different. Also, treating the case N = 1 separately simplifiesnotation and exposition.In Sections 1.3 and 1.4 we explain these data structure scenarios as well asintroduce some notation.1.3 One realization caseIn this first scenario we consider a data set with no replicates in the sense that thereis only one set of n pairs xi and yi.Our notation for this case is? observed data: y = (y1, . . . ,yn)T , x = (x1, . . . ,xn)T ;? hidden states: z = (z1, . . . ,zn)T with zi taking a value in {1, . . . ,J};? n-vector of function evaluations: f j(x)= ( f j(x1), . . . , f j(xn))T for j = 1, . . . ,J.We assume that, given the zi, yi = fzi(xi)+?zi?i, for ?1, . . . ,?n are independentwith mean of zero and variance of one. Therefore, given the f j?s and the zi?s, theyi?s are independent with the mean of yi equal to fzi(xi) and the variance equal to?2zi .Define ? to be the set of parameters defining the distribution of y given z and thef j(x)?s. So, for instance, if the ?i?s are normally distributed, then ? = { f j(x) and ?2jfor j = 1, . . . ,J}. Let ? be vector of parameters that determine the joint distribu-tion of z. If the zi?s are independent and identically distributed, then the parametervector ? is of length J with jth component equal to p(zi = j). If the zi?s follow aMarkov structure, the parameter vector ? consists of initial and transition proba-bilities.We are interested in estimating ? ? {?,?}, along with standard errors or somemeasure of accuracy for the parameters in ? .91.4 Replicate caseIn this second scenario we consider a data set with N replicates. Replicate k con-tains n observations y1k, . . . ,ynk and evaluation points x1, . . . ,xn that are the sameacross replicates. There are two possibilities of how the replicates are generated:1. they are all generated from just one set of functions f1, . . . , fJ or2. each replicate is generated from a different set of functions, with replicate kgenerated from f1k, . . . , fJk.In this thesis, we assume replicates are generated according to the first possibilityas it matches the building data presented in Figure 1.2 and described in Section 1.1.Our notation is? observed data: x=(x1, . . . ,xn)T (fixed across replicates); yR =(yT1 , . . . ,yTN)T ,where yk = (y1k, . . . ,ynk)T for k = 1, . . . ,N;? hidden states: zR = (zT1 , . . . ,zTN)T , where zk = (z1k, . . . ,znk)T , k = 1, . . . ,N,and zik takes a value in {1, . . . ,J};? f j(x) = ( f j(x1), . . . , f j(xn))T for j = 1, . . . ,J.We assume that, given f1, . . . , fJ , the replicates, y1, . . . ,yN , are independent.Let fzk(x) =(fz1k(x1), . . . , fznk(xn))T. Given the hidden states zk, yk = fzk(x) +?k, where ?k ? MV N(0,V), independent of f1, . . . , fJ . The vectors ?1, . . . ,?N areindependent so that Cov(?ik,? jl) = 0 for all k 6= l and for all i and j. Therefore,given the f j(x)?s and the hidden states zk, yk ? MV N( fzk(x);V). Except forone case (Section 3.2.1), we do not allow V to depend on the hidden state vectorzk. This assumption greatly simplifies calculations in the estimation problem. Wealso assume that the latent variable vectors, z1, . . . ,zN , are independent with thedistribution of each zk governed by a parameter vector ? as in the one realizationcase. In this scenario the set of parameters ? contains f1(x), . . . , fJ(x) and theparameters in V, and as before we are interested in estimating ? = {?,?}.1.5 Literature reviewSimilar models for the one realization case have appeared in the machine learn-ing literature, where they are called mixture of Gaussian processes models. See10Tresp (2001), Rasmussen and Ghahramani (2002) and Ou and Martin (2008). Thefocus of these papers is on analyzing data from an on-line process, with the goalbeing prediction of a single process. The process is modeled as a mixture of real-izations of Gaussian processes and, just as in our model, the mixture depends onhidden states. Rasmussen and Ghahramani and Ou and Martin take a hierarchi-cal Bayesian approach, not only placing a Gaussian process prior on the f j?s, butalso placing a prior distribution on all parameters, including those that govern theGaussian process. Tresp does not place a prior distribution on the parameters ofthe Gaussian process, but he does use Gaussian processes to model not only thefunctions themselves but also the latent process and the regression error variance.To analyze data, Tresp uses the Expectation-Maximization (EM) algorithm to max-imize the posterior density of the Gaussian processes given the data. Rasmussenand Ghahramani and Ou and Martin use a Monte Carlo method to estimate theposterior distribution of the unknown functions and parameters.A main difference between our model and the model in these three papers isin the distribution of the zi?s. All three papers begin by assuming that the zi?sare independent, conditional on the parameters governing the latent process. InRasmussen and Ghahramani and in Ou and Martin, the process parameters arep j = P(zi = j), j? 1, which are modeled using a Dirichlet distribution. Rasmussenand Ghahramani use a limit of a Dirichlet distribution in order to model an infinitenumber of possible hidden states, to avoid choosing the number of states. Ouand Martin use a finite number of states to avoid computational complexity. Bothpapers use an ad hoc modification of the distribution of the zi?s to allow zi to dependon xi in a smooth way; see Rasmussen and Ghahramani?s equation (5) and Ou andMartin?s equation (13). However, as remarked by Rasmussen and Ghahramani inthe discussion, the properties of the resulting joint distribution of the zi?s are notcompletely known. Tresp?s distributional assumptions are more straightforward:he assumes that the distribution of zi depends on xi according to a logit modelgoverned by a Gaussian process. None of these papers consider Markov zi?s.While the three papers contain methodology that can, in principle, lead to es-timation of the f j?s and the latent variable process parameters, the papers focuson estimation of just one function - the mixture. Thus the resulting methodologyis a form of variable bandwidth smoothing (see, for instance, Fan and Gijbels,111992). In contrast, our goal is estimation of the individual processes that make upthe mixture and estimation of the parameters governing the hidden state process,along with standard errors. We see the distinction between the goals by consideringthe analysis of the motorcycle data: Rasmussen and Ghahramani present just onefunction to summarize the data. We present J > 1 functions, one for each of the Jsimulated accidents, and we estimate the expected proportion of data points fromeach function, and provide standard errors. We also conduct extensive simulationstudies of the frequentist properties of the estimates.A closely related model is the Gaussian mixture model, used in density esti-mation. In the Gaussian mixture model, we assume that the ith data point comesfrom one of a finite set of Gaussian distributions, determined by the value of alatent variable zi. In his work in this area, Bilmes (1998) considered two modelsfor the latent variables: in the first model the latent variables are independent andidentically distributed and in the second they follow a Markov structure. The latercorresponds to a classic hidden Markov model, which is the same as our approachfor Markov z?s if the f j?s are constant. Bilmes provides a very readable descriptionof how he applies the EM algorithm to estimate the parameters for the Gaussiandistributions as well as for the distribution of the zi?s.These models are not to be confused with the work of Shi et al. (2005) onGaussian process mixtures for regression and Chiou (2012) on functional mixtureprediction. These authors analyze data from M independent curves where the entireith curve is a realization of one of J Gaussian processes, determined by the value ofthe latent variable zi. In contrast, like Tresp, Rasmussen and Ghahramani and Ouand Martin, we consider M = 1 observed curve, which switches among J Gaussianprocesses. In addition, we consider M > 1 curves, where each curve switchesamong J Gaussian processes.We see that the literature contains many similar but distinct models with namescontaining the words Gaussian and mixture. For this reason, we prefer to call ourmodels and those of Tresp (2001), Rasmussen and Ghahramani (2002) and Ou andMartin (2008) switching nonparametric regression models as we feel this is moredescriptive. In a Bayesian switching nonparametric regression model the uncer-tainty about the f j?s is formulated by modeling the f j?s as realizations of stochasticprocesses. In a frequentist switching nonparametric regression model the f j?s are12merely assumed to be smooth. To our knowledge, no one has considered a Markovstructure for the latent variables to estimate multiple functions, nor has anyone usedthe non-Bayesian penalized likelihood approach. Our frequentist approach and ourcalculation of standard errors appear to be new.13Chapter 2Solution for the one realizationcaseIn this section we describe how the EM algorithm can be used to obtain the pa-rameter estimates for the single realization case considering both the Bayesian andthe penalized log-likelihood approaches. The E-step of the algorithm is exactlythe same for both approaches. The M-step differs only in the part involving thecalculation of f? j. In the M-step we restrict our calculations to normally distributederrors. Furthermore, in the M-step, for the penalized log-likelihood case we modeleach f j as a linear combination of K known B-spline basis functions b1, . . . ,bK ,that is,f j(x) =K??=1? j?b?(x), (2.1)with {? j1, . . . ,? jK} the set of unknown parameters determining f j. We justify thisform of f j?s heuristically in Section 2.2 using a result of Silverman (1985) .Recall that for the one realization case ? = {?,?}, where ? is the vector con-taining the parameters of the model assumed for z, and ? = { f j(x) and ?2j for j =1, . . . ,J} are the parameters governing the distribution of yi given the f j(x)?s andthe zi?s. Let log p(y|?) be the log-likelihood based on the observed data. Our goal14is to find ?? that maximizesl(?)? log p(y|?)+P( f1, . . . , fJ,?1, . . . ,?J), (2.2)with P( f1, . . . , fJ,?1, . . . ,?J) as in (1.2) or (1.3). The form of log p(y|?) is verycomplicated, since it involves the distribution of the latent zi?s. Therefore, themaximization of (2.2) with respect to ? is a difficult task. In order to tackle thisproblem it is common to apply numerical methods such as the EM algorithm to ob-tain the parameter estimates (McLachlan and Krishnan, 2008). The EM algorithmis usually used to maximize the likelihood function by generating a sequence of es-timates, ? (c), c ? 1, having the property that log p(y|? (c+1)) ? log p(y|? (c)). TheEM algorithm can also be used to maximize (2.2). We show in Appendix A thatour EM algorithm also generates a sequence of estimates, ? (c), c? 1, satisfyingl(? (c+1))? l(? (c)). (2.3)To define our EM algorithm, let p(y,z|?) be the joint distribution of the ob-served and latent data given ? , also called the complete data distribution. In whatfollows view ? as an argument of a function, not as a random variable. Note that wewrite E? (c)[H(y,z)|y]to denote the conditional expected value of H(y,z) assumingthat the data y and z are generated with parameter vector ? (c).The proposed EM algorithm consists of the following two steps based on writ-ing log p(y,z|?) =L1(?)+L2(?), whereL1(?) = log p(y|z,?) =n?i=1log p(yi|zi, fzi(xi),?2zi)andL2(?) = log p(z|?) = log p(z1, . . . ,zn|?). (2.4)1. Expectation step (E-step): calculateQ(? ,? (c))? E? (c)[log p(y,z|?)|y]= E? (c)(L1(?)|y)+E? (c)(L2(?)|y)15that is, calculate the expected value of the logarithm of p(y,z|?) with respectto the distribution of the latent z given the observed y using ? (c) as the truevalue of ? .2. Maximization step (M-step): letS(? ,? (c)) = Q(? ,? (c))+P( f1, . . . , fJ,?1, . . . ,?J)= E? (c)(L1(?)|y)+P( f1, . . . , fJ,?1, . . . ,?J)+E? (c)(L2(?)|y). (2.5)Use the Expectation-Conditional Maximization (ECM) algorithm to find ? (c+1)that maximizes S(? ,? (c)) with respect to ? or at least does not decreaseS(? ,? (c)) from the current value at ? (c). We show in Appendix A that ifS(? (c+1),? (c)) ? S(? (c),? (c)) then (2.3) holds. The ECM algorithm is de-scribed in details in Appendix B.2.1 E-step: general zi?sSinceL1(?) =n?i=1J?j=1I(zi = j) log p(yi|zi = j, f j(xi),?2j ),E? (c)(L1(?)|y) =n?i=1J?j=1p(c)i j log p(yi|zi = j, f j(xi),?2j )wherep(c)i j = p(zi = j|y,?(c)),whose exact form depends on the model for the zi?s.16In a regression model with normal errors, which we assume in the M-step,E? (c)(L1(?)|y) = ?n2log(2pi)? 12n?i=1J?j=1p(c)i j log?2j (2.6)?12n?i=1J?j=1p(c)i j[yi? f j(xi)]2?2j.In the following sections we calculate E? (c)(L1(?)|y) and E? (c)(L2(?)|y) con-sidering different models for the latent variables. Note that E? (c)(L1(?)|y) dependson the model for the zi?s only through p(c)i j .2.1.1 E-step: independent and identically distributed (iid) zi?sIn this section we calculate p(c)i j and E? (c)(L2(?)|y) assuming that the latent vari-ables z1, . . . ,zn are iid with parameter vector ? = (p1, ..., pJ), where p j = p(zi =j|?) and ?Jj=1 p j = 1.Since the zi?s are iid, we obtainp(c)i j = p(zi = j|yi,?(c))=p(yi|zi = j, f j(xi)(c),?2(c)j )? p(c)j?Jl=1 p(yi|zi = l, fl(xi)(c),?2(c)l )? p(c)l. (2.7)Note that we can easily calculate (2.7) when the regression errors are normallydistributed.We now calculate E? (c)(L2(?)|y). For iid zi?s we can writeL2(?) =n?i=1J?j=1I(zi = j) log p j (2.8)and, therefore,E? (c)(L2(?)|y) =n?i=1J?j=1p(c)i j log p j. (2.9)172.1.2 E-step: Markov zi?sHere we calculate p(c)i j and E? (c)(L2(?)|y) assuming a Markov structure for thelatent variables z1, . . . ,zn. In this case the distribution of the zi?s depends on a vector? composed of transition probabilities and initial probabilities pi j = p(z1 = j|?),j = 1, . . . ,J.Specifically, we assume that1. the ith latent variable depends on past latent variables only via the (i? 1)stlatent variable, i.e., p(zi|zi?1, . . . ,z1,?) = p(zi|zi?1,?);2. the transition probabilities do not depend on i, that is,p(zi = j|zi?1 = l,?) = p(zi+s = j|zi+s?1 = l,?)? al j.To compute p(c)i j when the zi?s are Markov, we use the results of Baum et al.(1970). These authors let? (c)i j = p(y1, . . . ,yi,zi = j|?(c))and?(c)i j = p(yi+1, . . . ,yn|zi = j,?(c)),and show how to calculate these recursively using what they call the forward andbackward procedures, respectively. In the forward procedure ? (c)i j is calculatedrecursively as:1. ? (c)1 j = pi(c)j ? p(y1|z1 = j,? (c));2. ? (c)(i+1) j =[J?l=1? (c)il a(c)l j]? p(yi+1|zi+1 = j,? (c)) i = 1, . . . ,n?1.The backward procedure calculates ?(c)i j recursively as:1. Define ?(c)n j = 1.2. Let ?(c)i j =J?l=1[a(c)jl ? p(yi+1|zi+1 = l,?(c))??(c)(i+1)l], i = n?1, . . . ,1.18Note that because of the Markovian conditional independence? (c)i j ?(c)i j = p(y1, . . . ,yi,zi = j|?(c))? p(yi+1, . . . ,yn|zi = j,? (c))= p(y1, . . . ,yi,zi = j|? (c))? p(yi+1, . . . ,yn|zi = j,y1, . . . ,yi,? (c))= p(y,zi = j|? (c)).Thus, we can calculate p(c)i j viap(c)i j =? (c)i j ?(c)i j?Jl=1 ?(c)il ?(c)il.Now let us consider E? (c)(L2(?)|y). Since for Markov zi?sL2(?) =n?i=2J?l=1J?j=1I(zi?1 = l,zi = j) logal j +J?j=1I(z1 = j) logpi j (2.10)E? (c)(L2(?)|y) =n?i=2J?l=1J?j=1p(c)il j logal j +J?j=1p(c)1 j logpi j, (2.11)where p(c)il j ? p(zi?1 = l,zi = j|y,? (c)). This can be expanded asp(c)il j =p(c)(i?1)l?a(c)l j ? p(yi|zi = j,? (c))??(c)i j?(c)(i?1)l,which we easily calculate using the normality assumption for the regression errors.More details on how to obtain the expressions for p(c)i j and p(c)il j can be found inBilmes (1998), Rabiner (1989) and Cappe? et al. (2005).2.2 M-stepFor the M-step we combine our discussion of iid z?s and Markov z?s. We wantto find ? (c+1) that maximizes S(? ,? (c)) in (2.5) with respect to ? or at least pro-duces a value of S no smaller than S(? (c),? (c)). For normally distributed errors19E? (c)(L1(?)|y) is given by (2.6) and, therefore, we can writeS(? ,? (c)) =?n2log(2pi)? 12n?i=1J?j=1p(c)i j log?2j ?12n?i=1J?j=1p(c)i j[yi? f j(xi)]2?2j(2.12)+ P( f1, . . . , fJ,?1, . . . ,?J) (2.13)+ E? (c)(L2(?)|y). (2.14)We maximize S as a function of ? = { f j(x),?2j , j = 1, . . . ,J and ?}. In themaximization ? (c) is fixed, not depending on ? . This implies that the p(c)i j ?s are alsofixed, since their calculation depends on current parameter estimates in ? (c). Wealso consider the smoothing parameters, ?1, . . . ,?J , to be fixed. This maximizationcannot be done analytically, so we apply a natural extension of the EM approach,the ECM algorithm, to guarantee that S(? (c+1),? (c))? S(? (c),? (c)).Because the expression for E? (c)(L2(?)|y) does not depend on the f j?s or ?2j ?s,the f j?s and ?2j ?s that maximize S are the f j?s and ?2j ?s that maximize (2.12) +(2.13). Therefore, the forms of the maximizing f j?s and ?2j ?s will not depend onthe model for the z?s. Their only dependence on the model for the z?s is via thep(c)i j ?s.Thus, to obtain the vector of parameter estimates, ? (c+1), we apply the ECMalgorithm as follows.1. Hold the ?2j ?s and the parameters in ? fixed and maximize (2.12) plus (2.13)with respect to f j(x). LetW j = diag(p(c)1 j /?2j , . . . , p(c)n j /?2j ). (2.15)For the Bayesian approach this is equivalent to maximizing?12J?j=1(y? f j(x))T W j(y? f j(x))?12J?j=1f j(x)T A(? j)?1 f j(x)20obtainingf? j(x) =(I+A(? j)W j)?1A(? j)W jy.Let f (c+1)j (x) be f? j(x) with ?2j in W j replaced by ?2(c)j .For the penalized log-likelihood approach recall that we assume each f j isa linear combination of K known B-spline basis functions (see (2.1)). Sof j(x) = B? j, where ? j = (? j1, . . . ,? jK)T is the vector of coefficients corre-sponding to f j and B is an n?K matrix with entries Bi? = b?(xi). We justifythis assumption heuristically by noting that the maximization criterion for f jgiven by (2.12) plus (2.13) is similar to (5.1) in Silverman (1985) for whichthe maximizer is shown to be a cubic spline. See also Heckman (2012).Thus, we hold the ?2j ?s and ? fixed and maximize?12J?j=1(y?B? j)T W j(y?B? j)?J?j=1? j?Tj R? j,with respect to ? j yielding?? j = (BT W jB+2? jR)?1BT W jy,where R is a K?K matrix with entries R?? ? =?b???(x)b??? ?(x) dx.Let ? (c+1)j be ?? j with ?2j in W j replaced by ?2(c)j . So we let f(c+1)j (x) =B? (c+1)j .2. Now holding the f j(x)?s and the parameters in ? fixed and maximizing(2.12) with respect to ?2j we get??2j =n?i=1p(c)i j[yi? f j(xi)]2n?i=1p(c)i j. (2.16)21Let ?2(c+1)j be ??2j with f j(xi) replaced by f j(xi)(c+1). If we assume ?2j = ?2for all j, then we find that??2(c+1) = 1nJ?j=1n?i=1p(c)i j[yi? f j(xi)(c+1)]2. (2.17)We can obtain better variance estimates by adjusting the degrees of freedomto account for the estimation of the f j?s. We do that by replacing the denom-inator of (2.16) by ?ni=1 p(c)i j ? trace(D jH j) and the denominator of (2.17) byn??Jj=1 trace(D jH j), where D j = diag(p(c)1 j , . . . , p(c)n j ) and H j is the so calledhat matrix satisfying f? j(x) = H jy. For the Bayesian approach H j = (I+A(? j)W j)?1A(? j)W j and for the penalized approach H j = B(BT W jB +2? jR)?1BT W j. This modification is similar to a weighted version of whatWahba (1983) has proposed for the regular smoothing spline case.3. Now we hold the f j(x)?s and ?2j ?s fixed and maximize S with respect to theparameters in ? . Note that (2.12) and (2.13) do not depend on ? , so to find?(c+1), we maximize E? (c)(L2(?)|y) in line (2.14) as a function of ? . Be-cause the form of E? (c)(L2(?)|y) does depend on the model for the z?s, wemust obtain the estimates of ? for each model separately.For iid zi?s, where p j = p(zi = j|?), using (2.9) and Lagrange multiplierswith the restriction that ?Jj=1 p j = 1, we obtain:p(c+1)j =1nn?i=1p(c)i j . (2.18)For Markov z?s the vector ? is composed of transition probabilities al j andinitial probabilities pi j. We first maximize E? (c)(L2(?)|y) given in (2.11)with respect to al j. Holding pi j fixed and using a Lagrange multiplier with22the constraint ?Jj=1 al j = 1, we get:a(c+1)l j =n?i=2p(c)il jn?i=2p(c)(i?1)l.Now let us maximize (2.11) with respect to pi j. Holding al j fixed and usinga Lagrange multiplier with the restriction that ?Jj=1pi j = 1, we obtain:pi(c+1)j = p(c)1 j .2.3 Standard errors for the parameter estimators of thestate processIn this section we use the results of Louis (1982) to obtain standard errors for theestimates of the parameters of the state process. We consider for iid zi?s J ? 2possible state values. For Markov zi?s we restrict the possible number of states toJ = 2 to reduce calculational complexity.Louis (1982) derived a procedure to obtain the observed information matrixwhen the maximum likelihood estimates are obtained using the EM algorithm. Theprocedure requires the computation of the gradient and of the second derivativematrix of the log-likelihood based on the complete data and can be implementedquite easily within the EM steps.2.3.1 The general caseWe propose estimating the covariance matrix of ?? by using the inverse of an esti-mate of I?(?), the observed information matrix if one knew ? and only estimated? . Specifically, we estimate the covariance of ?? by the plug-in estimate [I??(??)]?1,where ?? and ?? are the parameters estimates obtained from our EM procedure.We can find I? easily in our existing EM algorithm since our complete data log-likelihood splits in two parts: L1(?) that only depends on ? and L2(?) that onlydepends on ? . Note that this method ignores the variability in estimating ? .23We first describe how to find I?(?), the information matrix associated with ???,the maximum likelihood estimate of alpha when ? is known. We then show that ???is equal to ?? , our estimate of ? .Suppose that ? is known and let ??? = ???(?) be the maximum likelihood esti-mator of ? given ? , that is, the maximizer of LI(?) ? log p(y|?), the incompletedata log-likelihood. We can obtain ??? using the EM algorithm and the completedata log-likelihood LC(?)? log p(y,z|?) =L2(?)+C. In this case we can derivethe observed information matrix, I?(?), via a direct application of Louis?s proce-dure. Under some regularity conditions, Louis (1982) shows by straightforwarddifferentiation that L?I(?) = E?(L?C(?)|y), L?I(???) = 0, and that the observed infor-mation matrix is given byI?(?) = ?L??I (?)= E?(?L??2 (?)|y)?E?(L ?2(?)L ?2(?)T |y)+L?I(?)L?I(?)T ,(2.19)whereL2(?) is as in (2.4),L ?2(?) and L?I(?) are the gradient vectors ofL2 and LI ,respectively, and L ??2 (?) and L??I (?) are the associated second derivative matrices.The estimate of I?(?) is I?(???). Note that (2.19) needs to be evaluated only atconvergence of the EM algorithm, where L?I is zero. Then, I?(???) = I?(???(?))contains only the first two terms of (2.19). The inverse of I?(???) is the estimatedvariance-covariance matrix of ??? for the known value of ? .To relate these calculations to those of our EM algorithm, where both ? and? are estimated, first note that our ?? and ?? are the maximizers of the criterionlog p(y|?,?)+P( f1, . . . , fJ,?1, . . . ,?J). If the maximizer is unique then ?? = ???(??).Therefore, it is reasonable to estimate the variance-covariance matrix of ?? by thevariance-covariance matrix of ??? plugging in ? = ?? , that is, by calculating theinverse of I??(??).In the next sections we show how to calculate I??(??) for the different modelsfor the zi?s.242.3.2 Standard errors: iid zi?sTo remove the restriction ?Jj=1 p j = 1, we use the parameters p1, . . . , pJ?1 andrewrite (2.8) asL2(?) =n?i=1{J?1?j=1I(zi = j) log p j + I(zi = J) log(1?J?1?j=1p j)}. (2.20)Let L ?2(?) be the (J? 1)? 1 gradient vector of (2.20) with the jth componentgiven by?L2? p j=n?i=1[I(zi = j)p j?I(zi = J)(1??J?1k=1 pk)]and L ??2 (?) be the (J?1)? (J?1) matrix of second derivatives? 2L2? p2j=?n?i=1[I(zi = j)p2j+I(zi = J)(1??J?1k=1 pk)2],? 2L2? p j pl=?n?i=1[I(zi = J)(1??J?1k=1 pk)2]for j 6= l.Consider the (J? 1)? (J? 1) matrix E?(?L ??2 (?)|y) in (2.19) evaluated at? = ?? with ? = ?? . One can show its jlth entry, for j 6= l, is equal to n/(1??J?1k=1 p?k)and its j jth entry isn?i=1(p?i jp?2j+(1??J?1k=1 p?ik)(1??J?1k=1 p?k)2)= n?(1p? j+11??J?1k=1 p?k),where p?i j is p(c)i j evaluated at convergence of the EM algorithm, where we obtain?? = {??, ??}, that is, p?i j = p(zi = j|yi, ??). Note that the simplification above isobtained using the fact that p? j = ?ni=1 p?i j/n.One can also show that the (J? 1)? (J? 1) matrix E?(L ?2(?)L ?2(?)T |y) in(2.19) evaluated at ?? with ? = ?? has off diagonal elements jl equal ton(1??J?1k=1 p?k)?n?i=1[(p?i jp? j?(1??J?1k=1 p?ik)(1??J?1k=1 p?k))?(p?ilp?l?(1??J?1k=1 p?ik)(1??J?1k=1 p?k))]25and jth diagonal element given byn?(1p? j+11??J?1k=1 p?k)?n?i=1(p?i jp? j?(1??J?1k=1 p?ik)(1??J?1k=1 p?k))2.2.3.3 Standard errors: Markov zi?sFor Markov zi?s we show how to obtain standard errors for the estimates of thetransition probabilities for J = 2 possible state values. We apply Louis?s method tofind standard errors for a?12 and a?21 by first consideringL2 in (2.10) with pi1 and pi2fixed. Abusing notation slightly forL2 by omitting pi1 and pi2, we let a11 = 1?a12and a22 = 1?a21 and writeL2(a12,a21) =n?i=2[I(zi?1 = 1,zi = 2) loga12+ I(zi?1 = 1,zi = 1) log(1?a12)+ I(zi?1 = 2,zi = 1) loga21+ I(zi?1 = 2,zi = 2) log(1?a21)].The required two dimensional gradient vector is given by?????L2?a12?L2?a21????=?????n?i=2[I(zi?1 = 1,zi = 2)a12?(I(zi?1 = 1,zi = 1)(1?a12)]n?i=2[I(zi?1 = 2,zi = 1)a21?I(zi?1 = 2,zi = 2)(1?a21)]?????.The associated 2?2 matrix of second derivatives is diagonal with entries? 2L2?a212=?n?i=2[I(zi?1 = 1,zi = 2)a212+I(zi?1 = 1,zi = 1)(1?a12)2]and? 2L2?a221=?n?i=2[I(zi?1 = 2,zi = 1)a221+I(zi?1 = 2,zi = 2)(1?a21)2].26Thus the 2?2 matrix E?(?L ??2 (a12,a21)|y) with pi1 and pi2 fixed, evaluated ata12 = a?12, a21 = a?21 with ? = ?? is given by?????ni=2 p?(i?1)1a?12(1? a?12)00?ni=2 p?(i?1)2a?21(1? a?21)???? ,where p?(i?1)1 = p(zi?1 = 1|yi?1, ??) and p?(i?1)2 = p(zi?1 = 2|yi?1, ??).Calculating the 2? 2 matrix E?(L ?2(a12,a21)L?2(a12,a21)T |y) is straightfor-ward but tedious, involving sums of expectations of indicator functions. The sum-mands require the calculation of p(zi?1 = r,zi = s|y,?), p(zi?2 = r,zi?1 = s,zi =t|y,?) and p(zi?1 = r,zi = s,zi+? = t,zi+?+1 = u|y,?), with r,s, t and u taking val-ues 1 or 2 and ? a positive integer. These conditional probabilities can be calculatedusing Bayes? Theorem and the Markovian conditional independence of the zi?s. Formore details see Appendix C.27Chapter 3Solution for the replicate caseIn this section we describe how the EM algorithm can be used to obtain the param-eter estimates for the replicate case (see Section 1.4) considering both the Bayesianand the penalized log-likelihood approaches. As in the one realization case the E-step of the algorithm is exactly the same for both approaches. The M-step differsonly in the part involving the calculation of f? j. In the M-step, for the penalizedlog-likelihood we also model each f j as a linear combination of K known B-splinebasis functions as in (2.1).Suppose we have a data set with N replicates, all generated from just one set offunctions f1, . . . fJ . We assume the following, as stated in Section 1.4.? Given the hidden states zk, yk = fzk(x)+ ?k, where ?k ?MV N(0,V), inde-pendent of f1, . . . , fJ . The vectors ?1, . . . ,?N are independent so that Cov(?ik,? jl)=0 for all k 6= l and for all i and j. Therefore, given the f j(x)?s and the hiddenstates zk, yk ? MV N( fzk(x);V), that is,p(yk|zk, f1(x), . . . , fJ(x),V) =(2pi)?n2 |V|?12 exp[?12(yk? fzk(x))T V?1(yk? fzk(x))]. (3.1)? The latent variable vectors, z1, . . . ,zN , are independent with the distributionof each zk governed by a parameter vector ? .28The n? n matrix V represents the covariance structure among observationswithin the same replicate k. We consider two types of V: V unrestricted andV generated by a random intercept model. For the latter we assume that yik =fzik(xi)+ ?ik, where- ?ik = ?k + eik;- {?1, . . . ,?N} and {eik for i = 1, . . . ,n, k = 1, . . . ,N} are independent;- ?k?s are iid N(0,?2);- eik?s are iid N(0,?2).Therefore, the variance-covariance matrix V depends on only two parameters andhas entries given by:Vii? = Cov(?ik,?i?k) ={?2 if i 6= i?Var(?ik) = ?2 +?2 if i = i?.(3.2)In Section 3.2.1 we describe the solution for the case when the covariancematrix is diagonal with entries depending on the hidden states. This covariancematrix is written as Vzk = diag(?2z1k , . . . ,?2znk), so that if zik = j the regression errorvariance is ?2j . Note that this approach ignores the within replicate correlation.We are interested in estimating the set of parameters ? = {?,?}, where ? ={ f1(x), . . . , fJ(x) and the parameters in V} and ? contains the parameters from thez process, that is, p(zk|?) = p(zk|?). In this case we want to find ?? that maximizesl(?)?N?k=1log p(yk|?)+P( f1, . . . , fJ,?1, . . . ,?J), (3.3)where log p(yk|?) is the log-likelihood function based on the observed data fromthe kth replicate.Again we use the EM algorithm to obtain parameter estimates that maximize(3.3). Sections 3.1 and 3.2 present the details of the E and M steps, respectively.Recall that yR = (yT1 , . . . ,yTN)T and zR = (zT1 , . . . ,zTN)T . The application of theEM algorithm to the replicate case is very similar to the one realization case de-scribed in Sections 2.1 and 2.2 and it is again based on writing log p(yR,zR|?) =29L1(?)+L2(?), where now for the replicate caseL1(?) = log p(yR|zR,?) =N?k=1log p(yk|zk, fzk(x),V),andL2(?) = log p(zR|?) =N?k=1log p(zk|?). (3.4)3.1 E-step detailsIn this step we calculate Q(? ,? (c))=E? (c)(log p(yR,zR|?)|yR)=E? (c)(L1(?)|yR)+E? (c)(L2(?)|yR).Let s be an n-vector of possible hidden states, i.e., each entry of s is an integerj = 1,2, . . . , or J. SinceL1(?) =N?k=1?alls?sI(zk = s) log p(yk|zk = s, fs(x),V),E? (c)(L1(?)|yR)) =N?k=1?alls?spk(s)(c) log p(yk|zk = s, fs(x),V)=?Nn2log2pi? N2log |V|?12N?k=1?alls?spk(s)(c)[(yk? fs(x))T V?1(yk? fs(x))],(3.5)wherepk(s)(c) ? p(zk = s|yR,? (c))= p(zk = s|yk,? (c)),whose exact form depends on the model for the hidden states zik?s. It is importantto note that in the building daytime power usage application n is small so that thesum over all possible s?s is not too large.30We calculate E? (c)(L1(?)|yR) and E? (c)(L2(?)|yR) assuming iid zik?s (Section3.1.1) and Markov zik?s (Section 3.1.2). Note that the calculation of E? (c)(L1(?)|yR)depends on the model for the zik?s only via pk(s)(c).3.1.1 E-step: iid zik?sIn this section we calculate pk(s)(c) and E? (c)(L2(?)|yR) assuming that, for k =1, . . . ,N, the latent variables z1k, . . . ,znk are iid with parameter vector ? =(p1, . . . , pJ),where p j = p(zik = j|?) and ?Jj=1 p j = 1.Since the zik?s are iid, we obtainpk(s)(c) =p(yk|zk = s,? (c))? p(zk = s|? (c))p(yk|? (c))=p(yk|zk = s, fs(x)(c),V(c))?n?i=1p(zik = si|?(c))?all u?sp(yk|zk = u, fu(x)(c),V(c))?n?i=1p(zik = ui|?(c))=p(yk|zk = s, fs(x)(c),V(c))?n?i=1p(c)si?all u?sp(yk|zk = u, fu(x)(c),V(c))?n?i=1p(c)ui. (3.6)Note that we can calculate (3.6) considering the normality assumption in (3.1).Now let us calculate E? (c)(L2(?)|y). For iid zik?s we can writeL2(?) =N?k=1?alls?sI(zk = s) log p(zk = s|?)=N?k=1?alls?sI(zk = s)n?i=1log psi=N?k=1?alls?sI(zk = s)n?i=1J?j=1I{si = j} log p j. (3.7)31Therefore,E? (c)(L2(?)|yR) =N?k=1?alls?spk(s)(c)J?j=1ns, j log p j, (3.8)where ns, j = ?ni=1 I{si = j}, i.e., the number of entries in s that are equal to j.Again pk(s)(c) can be calculated using (3.6).3.1.2 E-step: Markov zik?sIn this section we calculate pk(s)(c) and E? (c)(L2(?)|yR) assuming that, for k =1, . . . ,N, the latent variables z1k, . . . ,znk follow a Markov structure with parametervector ? . In this case ? contains the transition probabilities and initial probabilitiespi j = p(zik = 1|?), j = 1, . . . ,J. As in the one realization case we assume that1. the ith latent variable depends on past latent variables only via the (i? 1)stlatent variable, i.e., p(zik|z(i?1)k, . . . ,z1k,?) = p(zik|z(i?1)k,?);2. the transition probabilities do not depend on i, that is,p(zik = j|z(i?1)k = l,?) = p(z(i+s)k = j|z(i+s?1)k = l,?)? al j.For Markov zik?s we obtainpk(s)(c) = p(zk = s|yk,? (c))=p(yk|zk = s,? (c))? p(zk = s|? (c))p(yk|? (c))=p(yk|zk = s, fs(x)(c),V(c))?pi(c)s1n?i=2a(c)si?1si?all u?sp(yk|zk = u, fu(x)(c),V(c))?pi(c)u1n?i=2a(c)ui?1ui. (3.9)32Now let us consider E? (c)(L2(?)|yR). Since for Markov zik?sL2(?) =N?k=1?alls?sI(zk = s) log p(zk = s|?)=N?k=1?alls?sI(zk = s)[logpis1 +n?i=2logasi?1si]=N?k=1?alls?sI(zk = s)[J?j=1I(s1 = j) logpi j +J?j=1J?l=1ns,l j logal j](3.10)E? (c)(L2(?)|yR) =N?k=1?alls?spk(s)(c)J?j=1I(s1 = j) logpi j+N?k=1?alls?spk(s)(c)J?j=1J?l=1ns,l j logal j. (3.11)where ns,l j = ?ni=2 I(si?1 = l,si = j), that is, the number of times we go from statel to state j in the vector s.3.2 M-stepAs in the one realization case we combine our discussion of iid z?s and Markov z?s.In this step we want to find ? (c+1) that maximizesS(? ,? (c)) =C?N2log |V|?12N?k=1?alls?spk(s)(c)[(yk? fs(x))T V?1(yk? fs(x))](3.12)+ P( f1, . . . , fJ,?1, . . . ,?J) (3.13)+ E? (c)(L2(?)|yR) (3.14)with respect to ? or at least produces a value no smaller than S(? (c),? (c)).We maximize S as a function of ? = {?, f1(x), . . . , fJ(x) and the parame-ters in V}. In the maximization ? (c) is fixed, not depending on ? . This im-33plies that the pk(s)(c)?s are also fixed, since their calculation depends on the fixedcurrent parameter estimates in ? (c). We also consider the smoothing parameters,?1, . . . ,?J , to be fixed. As in the one realization case, this maximization cannot bedone analytically, so again we use the ECM algorithm to obtain ? (c+1) satisfyingS(? (c+1),? (c))? S(? (c),? (c)).Because the expression for E? (c)(L2(?)|yR) does not depend on the f j?s or V,the f j?s and V that maximize S are the f j?s and V that maximize (3.12) + (3.13).Therefore, the form of the maximizing f j?s or V will not depend on the model forthe z?s. Their only dependence on the model for the z?s is via the pk(s)(c)?s.The steps of the ECM algorithm can be summarized as follows.1. Hold V and the parameters in ? fixed and maximize (3.12) + (3.13) withrespect to the f j(x)?s, obtaining f1(x)(c+1), . . . , fJ(x)(c+1).2. Hold the f j(x)?s and the parameters in ? fixed and maximize (3.12) withrespect to the parameters in V, obtaining V(c+1).3. Now hold the f j(x)?s and V fixed and maximize (3.14) with respect to theparameters in ? , obtaining ?(c+1).In the following we present the detailed steps of the ECM algorithm.1. Updating the f j(x)?s (any V). Hold V and the parameters in ? fixed andmaximize S(? ,? (c)) with respect to the f j(x)?s. This is equivalent to max-imizing (3.12) + (3.13). So, disregarding the constant terms, we want tomaximizeN?k=1?alls?spk(s)(c)yTk V?1 fs(x)?12N?k=1?alls?spk(s)(c) fs(x)T V?1 fs(x) (3.15)+ P( f1, . . . , fJ,?1, . . . ,?J) (3.16)with respect to f1(x), . . . , fJ(x). We show in Appendix D.1 that (3.15) canbe written in terms of f1(x), . . . , fJ(x) asN?k=1yTk V?1J?j=1D(c)k j f j(x)?12N?k=1J?j=1J?m=1f j(x)T Skj,m fm(x),34where D(c)k j is an n?n diagonal matrix with ith diagonal element given by(D(c)k j )ii = p(zik = j|yk,?(c)) (3.17)and Skj,m is an n?n matrix with entries(Skj,m)ii? = (V?1)ii? p(zik = j,zi?k = m|yk,? (c)). (3.18)Note that when i = i? and j 6= m, (Skj,m)ii = 0. For i = i? and j = m (Skj, j)ii =(V?1)ii p(zik = j|yk,? (c)).We calculate the conditional probabilities required in (3.17) and (3.18) asfollows:p(zik = j|yk,? (c)) = ?s?s:si= jp(c)k (s)andp(zik = j,zi?k = m|yk,? (c)) = ?s?s:si= j,si?=mp(c)k (s).We propose two methods to maximize (3.15) + (3.16). Method 2 is straight-forward and yields an estimate of f = ( f1(x)T , . . . , fJ(x)T )T in a closed form.However, this method does involve big matrices, which can give rise to somecomputational issues. Method 1 involves smaller matrices than Method 2 asf1(x), . . . , fJ(x) are sequentially updated. This method may not necessarilymaximize (3.15) + (3.16), but will not decrease it. The two methods aredescribed in detail as follows.Method 1: sequential updatingFor each j, hold f1(x), . . . , f j?1(x), f j+1(x), . . . , fJ(x) fixed and maximize(3.15) + (3.16) with respect to f j(x).For the Bayesian approachP( f1, . . . , fJ,?1, . . . ,?J) =?12J?j=1f j(x)T A(? j)?1 f j(x)35and we show in Appendix D.1 that f? j(x) is the solution of the followinglinear system:(I+A(? j)N?k=1Skj, j)f? j(x) = (3.19)A(? j)[N?k=1(D(c)k j V?1yk?12 ?m:m6= jSkj,m fm(x))], (3.20)where the matrices D(c)k j , Skj,m are as in (3.17) and (3.18), respectively.We update the f j(x)(c)?s as follows. In (3.19) and (3.20), let V = V(c). Recallthat Skj,m also depends on the entries of V. Now- let f1(x)(c+1) be f?1(x) with f2(x), . . . , fJ(x) = f2(x)(c), . . . , fJ(x)(c);- for j = 2, . . . ,J?1, let f j(x)(c+1) be f? j(x) with fm(x) = fm(x)(c+1) form = 1, . . . , j?1 and fm(x) = fm(x)(c) for m = j+1, . . . ,J;- and let fJ(x)(c+1) be f?J(x) withf1(x), . . . , fJ?1(x) = f1(x)(c+1), . . . , fJ?1(x)(c+1).For the penalized log-likelihood approach recall that we assume f j is a linearcombination of K known basis functions, so f j(x) = B? j, where ? j is theK? 1 vector of coefficients corresponding to f j and B is an n?K matrixwith entries Bi? = b?(xi). Therefore, (3.15) + (3.16) can be written asN?k=1ykV?1J?j=1Dk jB? j?12n?k=1J?j=1J?m=1?Tj BT Skj,mB?m??j=1? j?Tj R? j, (3.21)where R is a K ?K matrix with entries R?? ? =?b???(x)b??? ?(x)dx, and thematrices Dk j, Skj,m are as in (3.17) and (3.18), respectively. So we can nowmaximize (3.21) as a function of ?1, . . . ,?J .For each j, hold ?1, . . . ,? j?1,? j+1, . . . ,?J fixed and maximize (3.21) with re-36spect to ? j. We show in Appendix D.1 that ?? j is the solution of the followinglinear system: [BT( N?k=1Skj, j)B+2? jR]?? j =BT V?1N?k=1Dk jyk?12N?k=1?m:m 6= jBT Skj,mB?m.We update the ? (c)j ?s sequentially, as in the Bayesian approach, and set f j(x)(c+1) =B? (c+1)j .Method 2: simultaneous updatingMaximize (3.15) + (3.16) with respect to the whole nJ-vector f=( f1(x)T , . . . , fJ(x)T )T .For the Bayesian approach we show in Appendix D.1 that f? is the solution ofthe following linear system:[diag(A(?1), . . . ,A(?J))M+ I]f? = diag(A(?1), . . . ,A(?J))v, (3.22)whereM =N?k=1Sk (3.23)is an nJ?nJ matrix withSk =????Sk1,1 Sk1,2 ? ? ? Sk1,J....... . ....SkJ,1 SkJ,2 ? ? ? SkJ,J???? ;andv =N?k=1D(c)Tk V?1yk (3.24)is an nJ-vector with D(c)k =[D(c)k1 ? ? ?D(c)kJ].37Let f(c+1) be f? with V = V(c). Recall again that the matrices Skj,m also dependon the entries of V.For the penalized log-likelihood approach we now maximize (3.21) with re-spect to the whole JK-vector ? = (?T1 , . . . ,?TJ )T . We show in Appendix D.1that ?? is the solution of the following linear system:[B?T MB?+2diag(?1R, . . . ,?JR)]?? = B?T v, (3.25)where B? is the nJ? JK matrix given byB? =??????B 0 ? ? ? 00 B ? ? ? 00 0. . ....0 0 ? ? ? B??????;and M and v are defined as in (3.23) and (3.24), respectively.Let ? (c+1) be ?? with V = V(c). So we let f(c+1) = B?? (c+1).Some advantages of the penalized log-likelihood approach are as follows.? Depending on the number of basis functions, the penalized log-likelihoodapproach may involve smaller matrices than the Bayesian approach.? The matrix B has the nice property of being banded, which may speedcalculations.? We need to estimate only J real-valued smoothing parameters while inthe Bayesian approach each ? j may be a vector of parameters.2. Updating V for unrestricted V. Hold the f j(x)?s and the parameters in ?fixed and maximize (3.12) with respect to V obtainingV? =1NN?k=1?alls?spk(s)(c)(yk? fs(x))(yk? fs(x))T.See Appendix D.2 for the calculation details. Note that if the values of zk areknown then V? is the sample covariance matrix.38Let V(c+1) be V? with fs(x) = fs(x)(c+1).Updating V generated by a random intercept model. Note that we canrewrite V given by (3.2) asV = ?2??????1+ ?2/?2 ?2/?2 ? ? ? ?2/?2?2/?2 1+ ?2/?2 ? ? ? ?2/?2....... . ....?2/?2 ?2/?2 ? ? ? 1+ ?2/?2??????= ?2(I+d11T),where I is an n?n identity matrix, 1 is an n?1 vector of ones and d = ?2/?2.To see how the parameter estimation in this case relates to the simple linearrandom intercept model refer to Section 2.4 of Demidenko (2004).We can calculate the determinant and inverse of V using the so-called dimension-reduction formulas (see Section 2.2.3 of Demidenko, 2004) as follows.|V|= ?2n|I+d11T |= ?2n(1+d1T 1) = ?2n(1+nd),andV?1 =1?2(I?d1+nd11T).Thus, we can rewrite (3.12) asC?N2log(1+nd)+d2?2(1+nd)N?k=1?alls?spk(s)(c)[(yk? fs(x))T 1]2(3.26)?Nn2log?2? 12?2N?k=1?alls?spk(s)(c)(yk? fs(x))T (yk? fs(x)). (3.27)We hold the f j(x)?s, ? and ?2 fixed and maximize (3.26) with respect to d39obtainingd? =1?2Nn2N?k=1?alls?spk(s)(c)[(yk? fs(x))T 1]2?1n. (3.28)Now holding the f j(x)?s and the parameters in ? fixed and maximizing(3.26) + (3.27) with respect to ?2 with d replaced by d? we obtain??2 = 1N(n?1)(N?k=1?alls?spk(s)(c)(yk? fs(x))T (yk? fs(x)) (3.29)?1nN?k=1?alls?spk(s)(c)[(yk? fs(x))T 1]2). (3.30)We update the parameter estimates of the restricted V as follows. Let ?2(c+1) =??2 and d(c+1) be d? with ?2 replaced by ?2(c+1) so that ?2(c+1) = d(c+1)??2(c+1).3. Updating ? (any V). We now hold the f j(x)?s and V fixed and maximizeS(? ,? (c)) with respect to the parameters in ? . Note that (3.12) and (3.13)do not depend on ? , so to find ?(c+1), we maximize E? (c)(L2(?)|yR) in line(3.14) with respect to ? . Because the form of E? (c)(L2(?)|yR) does dependon the model for the z?s, we must obtain the estimates for the parameters in? according to each proposed model.For iid zik?s, where p j = p(zik = j|?), using (3.8) and a Lagrange multiplierwith the restriction that ?Jj=1 p j = 1, we obtain:p(c+1)j =1NnN?k=1?alls?spk(s)(c)ns, j.For Markov zik?s the vector ? is composed of transition probabilities al j andinitial probabilities pi j. We first maximize E? (c)(L2(?)|yR) given in (3.11)with respect to al j. Holding pi j fixed and using a Lagrange multiplier with40the constraint ?Jj=1 al j = 1, we get:a(c+1)l j =N?k=1?alls?spk(s)(c)n?i=2I(si?1 = l,si = j)N?k=1?alls?spk(s)(c)n?i=2I(si?1 = l).Now let us maximize (3.11) with respect to pi j. Holding al j fixed and usinga Lagrange multiplier with the restriction that ?Jj=1pi j = 1, we obtain:pi(c+1)j =1NN?k=1?alls?spk(s)(c)I(s1 = j).3.2.1 ECM algorithm when Vzk = diag(?2z1k , . . . ,?2znk)Suppose that V = Vzk = diag(?2z1k , . . . ,?2znk), that is, given the f j?s and the zik?s,the observations y1k, . . . ,ynk are independent with yik distributed as N( fzik(xi);?2zik).Then most of Steps 1 and 3 of the ECM algorithm in Section 3.2 are still valid butStep 2 must be changed. However, it is easier to re-derive all three steps again asS(? ,? (c)) has a simpler form when the covariance matrix is diagonal than whenobservations are correlated.We can show that for Vzk = diag(?2z1k , . . . ,?2znk), S(? ,?(c)) is given byS(? ,? (c)) =C?12N?k=1n?i=1J?j=1p1k( j)(c) log?2j ?12N?k=1J?j=1(yk? f j(x))T Wk j(yk? f j(x))(3.31)+ P( f1, . . . , fJ,?1, . . . ,?J) (3.32)+ E? (c)(L2(?)|yR), (3.33)where pik( j)(c) = p(zik = j|yk,? (c)) andWk j = diag(p1k( j)(c)/?2j , . . . , pnk( j)(c)/?2j ). (3.34)41Similarly to the one realization case (see Equation 2.7), for the case of iid zik?s weobtainpik( j)(c) = p(zik = j|yk,? (c))= p(zik = j|yik,? (c))=p(yik|zik = j, f j(xi)(c),?2(c)j )? p(c)j?Jl=1 p(yik|zik = l, fl(xi)(c),?2(c)j )? p(c)l. (3.35)For Markov zik?s we calculate pik( j)(c) viapik( j)(c) =?ik( j)(c)?ik( j)(c)?Jl=1 ?ik(l)(c)?ik(l)(c),where ?ik( j)(c) = p(y1k, . . . ,yik,zik = j|? (c)) and ?ik( j)(c) = p(y(i+1)k, . . . ,ynk|zik =j,? (c)) are obtained using the forward and backward procedures described in Sec-tion 2.1.2, respectively.We want to maximize S(? ,? (c)) as a function of ? = { f j(x),?2j , j = 1, . . . ,J,and ?}. Because this maximization cannot be done analytically we apply the ECMalgorithm as follows.1. Updating the f j(x)?s. Hold the ?2j ?s and the parameters in ? fixed and max-imize (3.31) + (3.32) with respect to f j(x).Considering the Bayesian approach we obtainf? j(x) =(I+A(? j)N?k=1Wk j)?1A(? j)N?k=1Wk jyk. (3.36)For the penalized log-likelihood approach we obtainf? j(x) = B(BTN?k=1Wk j B+2? jR)?1BTN?k=1Wk jyk. (3.37)In both approaches we let f j(x)(c+1) be f? j(x) with ?2j in Wk j replaced by?2(c)j .422. Updating the ?2j ?s. Now holding the f j(x)?s and ? fixed and maximizing(3.31) with respect to ?2j we get?? j2 =N?k=1n?i=1pik( j)(c)[yik? f j(xi)]2N?k=1n?i=1pik( j)(c).Let ?2(c+1)j be ??2j with f j(xi) = f j(xi)(c+1). If we assume ?2j = ?2 for j =1, . . . ,J, then Vzk = V = ?2I and we find that?2(c+1) =N?k=1n?i=1J?j=1pik( j)(c)[yik? f j(xi)(c+1)]2Nn.3. Updating ? . Now we hold the f j(x)?s and the ?2j ?s fixed and maximizeS with respect to the parameters in ? . This is equivalent to maximizingE? (c)(L2(?)|yR) in line (3.33) as a function of ? . Because the form ofE? (c)(L2(?)|yR) does depend on the model for the hidden states, we mustobtain the estimates for the parameters in ? according to each proposedmodel. For iid zik?s we obtainp(c+1)j =1NnN?k=1n?i=1pik( j)(c)For Markov zik?s we geta(c+1)l j =N?k=1n?i=2p(z(i?1)k = l,zik = j|yk,? (c))N?k=1n?i=2p(z(i?1)k = l|yk,? (c))andpi(c+1)j =1NN?k=1p1k( j)(c).433.3 Choice of the smoothing parameters whenVzk = diag(?2z1k, . . . ,?2znk)We find the optimal ? j?s iteratively starting with initial values ? (0)j ?s. We set the? (0)j ?s to a value that worked well when tested on the data set.We update the ? j?s as follows.1. With ? j = ? (i)j , j = 1, . . . ,J, use the ECM algorithm of Section 3.2.1 to findthe p?ik( j)?s, the ??2j ?s and the f? j?s.2. Discard the f? j?s from Step 1.3. Treat the ??2j ?s and the p?ik( j)?s as fixed and thus Wk j as in (3.34) as fixed.Over a grid of possible values set ? (i+1)j as the value of ? that minimizes thefollowing ?leave-one-replicate-out? cross-validation criterion:CVj(? ) =N?k=1(yk? f?(?k)j?)T Wk j(yk? f?(?k)j?)(3.38)where f? (?k)j? ? f?(?k)j? (x) is the maximizer ofS(?k)( f j) =?12N?r=1:r 6=k(yr? f j(x))T Wr j(yr? f j(x))+P( f j,? ). (3.39)where for the Bayesian approach P( f j,? ) =?12 f j(x)T A(? ) f j(x) and for thepenalized log-likelihood approach P( f j,? ) =???Tj R? j.4. Repeat 1-3 with ? j = ? (i+1)j , j = 1, . . . ,J, till convergence.We use the final values of the ? j?s to obtain all of the parameter estimates from theECM algorithm as in Step 1.Finding ? that minimizes (3.38) is computationally intensive. Fortunately, wecan prove the following theorem.Theorem 1. The ?leave-one-replicate-out? estimate of f j, f?(?k)j? , satisfiesf? (?k)j? ? yk = (I?Hk j(? ))?1( f? j? ? yk) (3.40)44and, therefore, (3.38) can be expressed asCVj(? ) =N?k=1[(I?Hk j(? ))?1( f? j? ?yk)]TWk j[(I?Hk j(? ))?1( f? j? ?yk)].Theorem 1 is based on the following lemma:Lemma 1. Let f? (?k)j? be the maximizer ofS(?k)( f j) = ?12(f? (?k)j? ? f j(x))T Wk j(f? (?k)j? ? f j(x))?12N?r=1:r 6=k(yr? f j(x))T Wr j(yr? f j(x))+P( f j,? ), (3.41)that is, the maximizer of (3.31) + (3.32) in Section 3.2.1 with yk replaced by f?(?k)j? .Then, f? (?k)j? = f?(?k)j? .Proof. For simplicity let k = 1. We want to show that f? (?1)j? = f?(?1)j? . We knowf? (?1)j? maximizes S(?1)( f j) in (3.39) and, therefore,S(?1)( f? (?1)j? )?S(?1)( f? (?1)j? )? 0.We also know that f? (?1)j? is the maximizer of S(?1)( f j) in (3.41). Thus, S(?1)( f?(?1)j? )?S(?1)( f? (?1)j? )? 0, that is,?12(f? (?1)j? ? f?(?1)j?)T W1 j(f? (?1)j? ? f?(?1)j?)?12N?r=2(yr? f?(?1)j?)T Wr j(yr? f?(?1)j?)+ P( f? (?1)j? ,? )?[?12(f? (?1)j? ? f?(?1)j?)T W1 j(f? (?1)j? ? f?(?1)j?)?12N?r=2(yr? f?(?1)j?)T Wr j(yr? f?(?1)j?)+P( f? (?1)j? ,? )]? 0. (3.42)45We can rewrite (3.42) as?12(f? (?1)j? ? f?(?1)j?)T W1 j(f? (?1)j? ? f?(?1)j?)? S(?1)( f? (?1)j? )?S(?1)( f? (?1)j? )? 0,(3.43)which implies that(f? (?1)j? ? f?(?1)j?)T W1 j(f? (?1)j? ? f?(?1)j?)? 0,and because W1 j is positive definite it also implies f?(?1)j? = f?(?1)j? .Proof of the Theorem 1. The maximizer of (3.31) + (3.32) , f? j? , given by (3.36) or(3.37) can be written asf? j? =N?k=1Hk j(? )yk, (3.44)where for the Bayesian approachHk j(? ) =(I+A(? )N?r=1Wr j)?1A(? )Wk j,and for the penalized log-likelihood approachHk j(? ) = B(BTN?r=1Wr j B+2?R)?1BT Wk j.So, f? (?k)j? , the maximizer of (3.31) + (3.32) with yk replaced by f?(?k)j? , can bewritten asf? (?k)j? =N?r=1:r 6=kHr j(? )yr +Hk j(? ) f? (?k)j?= f? j? ?Hk j(? )yk +Hk j(? ) f?(?k)j? . (3.45)46From Lemma 1 we know f? (?k)j? = f?(?k)j? . Thus,f? (?k)j? = f? j? ?Hk j(? )yk +Hk j(? ) f?(?k)j? .Now subtracting yk on both sides of (3.3) we obtainf? (?k)j? ?yk = (I?Hk j(? ))?1( f? j? ?yk).3.4 Standard errors for the parameter estimators of thestate processAs in the one realization case we use the results of Louis (1982) to obtain standarderrors for the estimates of the parameters of the state process. We consider for iidzik?s J ? 2 possible state values. For Markov zik?s we restrict the possible numberof states to J = 2.As described in details in Section 2.3.1 we estimate the covariance of ?? by theplug-in estimate [I??(??)]?1, where ?? and ?? are the parameters estimates obtainedfrom our EM procedure and I?(?) is the observed information matrix if one knew? and only estimated ? . Using the results from Louis (1982) we can writeI?(?) = E?(?L ??2 (?)|yR)?E?(L ?2(?)L ?2(?)T |yR)+L?I(?)L?I(?)T , (3.46)where L2(?) is as in (3.4), L ?2(?) and L ??2 (?) are the gradient vector and theassociated second derivative matrix of L2, respectively, and L?I(?) is the gradientvector of the incomplete data log-likelihood. Recall from Section 2.3.1 that (3.46)needs to be evaluated only at convergence of the EM algorithm, where L?I is zero.In the next sections we show how to calculate I??(??) for the different modelsfor the zik?s.473.4.1 Standard errors: iid zik?sTo remove the restriction ?Jj=1 p j = 1, we use the parameters p1, . . . , pJ?1 andrewrite (3.7) asL2(?) =N?k=1?alls?sI(zk = s)[J?1?j=1ns, j log p j +ns,J log(1?J?1?j=1p j)], (3.47)where ns, j =?ni=1 I{si = j}. LetL?2(?) be the (J?1)?1 gradient vector of (3.47)with the jth component given by?L2? p j=N?k=1?alls?sI(zk = s)[ns, jp j?ns,J(1??J?1m=1 pm)]and L ??2 (?) be the (J?1)? (J?1) matrix of second derivatives? 2L2? p2j=?N?k=1?alls?sI(zk = s)[ns, jp2j+ns,J(1??J?1m=1 pm)2],? 2L2? p j pl=?N?k=1?alls?sI(zk = s)[ns,J(1??J?1m=1 pm)2]for j 6= l.Consider the (J? 1)? (J? 1) matrix E?(?L ??2 (?)|yR) in (3.46) evaluated at? = ?? with ? = ?? . One can show its jlth entry, for j 6= l, is equal to Nn/(1??J?1m=1 p?m) and its j jth entry isNn?(1p? j+11??J?1m=1 p?m),One can also show that the (J?1)? (J?1) matrix E?(L ?2(?)L ?2(?)T |yR) in(3.46) evaluated at ?? with ? = ?? has off diagonal elements jl equal toN?k=1?alls?sp?k(s)(ns, jp? j?ns,J(1??J?1m=1 p?m))(ns,lp?l?ns,J(1??J?1m=1 p?m))?N?k=1[?alls?sp?k(s)(ns, jp? j?ns,J(1??J?1m=1 p?m))? ?alls?sp?k(s)(ns,lp?l?ns,J(1??J?1m=1 p?m))]48and jth diagonal element given byN?k=1?alls?sp?k(s)(ns, jp? j?ns,J(1??J?1m=1 p?m))2?N?k=1[?alls?sp?k(s)(ns, jp? j?ns,J(1??J?1m=1 p?m))]2,where p?k(s) is pk(s)(c) evaluated at convergence, where we obtain ?? = {??, ??}, thatis, p?k(s) = p(zk = s|yk, ??).3.4.2 Standard errors: Markov zik?sFor Markov zik?s we show how to obtain standard errors for the estimates of theparameters of the z process for J = 2 possible state values. Recall that for the onerealization case we fix pi1 and find standard errors for a?12 and a?21. Because now forthe replicate case we can get much better estimates of pi1 we obtain standard errorsfor p?i1, a?12 and a?21 simultaneously. We apply Louis?s method considering L2 in3.10. We let pi2 = 1?pi1, a11 = 1?a12 and a22 = 1?a21 and writeL2(?) =N?k=1?alls?sI(zk = s)[I(s1 = 1) logpi1 + I(s1 = 2) log(1?pi1)+ ns,11 log(1?a12)+ns,12 loga12+ ns,21 loga21 +ns,22 log(1?a21),]where ns,l j = ?ni=2 I(si?1 = l,si = j).The required three dimensional gradient vector L ?2(?) is given by??????????L2?pi1?L2?a12?L2?a21?????????=???????????N?k=1?alls?sI(zk = s)[I(s1 = 1)pi1?I(s1 = 2)(1?pi1)]N?k=1?alls?sI(zk = s)[ns,12a12?ns,11(1?a12)]N?k=1?alls?sI(zk = s)[ns,21a21?ns,22(1?a21)]???????????.49The associated 3?3 matrix of second derivatives L ??2 (?) is diagonal with en-tries? 2L2?pi21=?N?k=1?alls?sI(zk = s)[I(s1 = 1)pi21+I(s1 = 2)(1?pi1)2],? 2L2?a212=?N?k=1?alls?sI(zk = s)[ns,12a212+ns,11(1?a12)2]and? 2L2?a221=?N?k=1?alls?sI(zk = s)[ns,21a221+ns,22(1?a21)2].Thus the 3?3 matrix E?(?L ??(?)|yR) in (3.46) evaluated at ? = ?? with ? = ??is diagonal with entriesNp?i1(1? p?i1),1a?12(1? a?12)N?k=1?alls?sp?k(s)n?i=2I(si?1 = 1) and1a?21(1? a?21)N?k=1?alls?sp?k(s)n?i=2I(si?1 = 2).The 3?3 matrix E?(L ?2(?)L ?2(?)T |yR) in (3.46) evaluated at ?? with ? = ?? issymmetric with diagonal elementsN?k=1?alls?sp?k(s)(I(s1 = 1)p?i1?I(s1 = 2)(1? p?i1))2?N?k=1[?alls?sp?k(s)(I(s1 = 1)p?i1?I(s1 = 2)(1? p?i1))]2,N?k=1?alls?sp?k(s)(ns,12a?12?ns,11(1? a?12))2?N?k=1[?alls?sp?k(s)(ns,12a?12?ns,11(1? a?12))]2andN?k=1?alls?sp?k(s)(ns,21a?21?ns,22(1? a?21))2?N?k=1[?alls?sp?k(s)(ns,21a?21?ns,22(1? a?21))]2.The [1,2], [1,3] and [2,3] off diagonal elements are50N?k=1?alls?sp?k(s)(I(s1 = 1)p?i1?I(s1 = 2)(1? p?i1))(ns,12a?12?ns,11(1? a?12))?N?k=1[?alls?sp?k(s)(I(s1 = 1)p?i1?I(s1 = 2)(1? p?i1))? ?alls?sp?k(s)(ns,12a?12?ns,11(1? a?12))],N?k=1?alls?sp?k(s)(I(s1 = 1)p?i1?I(s1 = 2)(1? p?i1))(ns,21a?21?ns,22(1? a?21))?N?k=1[?alls?sp?k(s)(I(s1 = 1)p?i1?I(s1 = 2)(1? p?i1))? ?alls?sp?k(s)(ns,21a?21?ns,22(1? a?21))]andN?k=1?alls?sp?k(s)(ns,12a?12?ns,11(1? a?12))(ns,21a?21?ns,22(1? a?21))?N?k=1[?alls?sp?k(s)(ns,12a?12?ns,11(1? a?12))? ?alls?sp?k(s)(ns,21a?21?ns,22(1? a?21))],respectively.51Chapter 4Simulation studies: the onerealization caseWe carry out simulation studies for the one realization case considering that thezi?s can take values 1 or 2 and they can be either iid or follow a Markov structure.The parameters of interest are estimated using both the Bayesian and the penalizedlog-likelihood approaches presented in Chapter 2. For each simulation study 300independent data sets are generated. Note that all tables and figures for the onerealization case are located at the end of this chapter.4.1 Simulated dataWe consider three types of simulation studies according to three different types ofsimulated data. Table 4.1 presents a summary of these simulation studies. In allstudies we use the same vector of evaluation points x and the same true functionsf1 and f2. The vector x = (x1, . . . ,xn)T consists of n = 199 equally spaced points,1,1.5, . . . ,99.5,100. The true functions evaluated at x, that is, the vectors f1(x) andf2(x), are obtained by sampling from a multivariate normal distribution with meanof zero and covariance matrix determined by (1.4). We let the parameter U j in (1.4)be equal to 1/(s j?2pi)so that ? j = s j. We consider ?1 = s1 = 28 and ?2 = s2 = 38for f1(x) and f2(x), respectively.After we obtain the commonly used f1(x) and f2(x) we generate a set of sim-52ulated data as follows.1. Generate the zi?s according to the specified model (iid or Markov).2. Generate the yi?s with common regression error variance ?2 as follows:? if zi = 1, yi = f1(xi)+??i;? if zi = 2, yi = f2(xi)+??i,where ?2 = 5?10?5 and ?i has a N(0,1) distribution.3. Repeat steps 1 and 2 S times obtaining S different data sets.Due to the computational complexity of our approach we consider S = 300. OurEM algorithm is quite fast, however the whole procedure also involves a couple ofcross-validation (CV) steps which makes it more complex and slower.Note again that f1(x) and f2(x) are fixed across all simulated data sets.In our first simulation study we generate data assuming that the zi?s are inde-pendent with p1 = p(z1 = 1) = 0.7. Figure 4.1 shows an example of a data set ofthis type.For the second and third simulation studies we consider Markov zi?s. In thesecond study we use transition probabilities a12 = p(zi = 2|zi?1 = 1) = 0.3 anda21 = p(zi = 1|zi?1 = 2) = 0.4 and in the third study a12 = 0.1 and a21 = 0.2. Forboth studies we use initial probabilities pi1 = pi2 = 0.5. Figures 4.2 and 4.3 showexamples of data sets from the second and third studies, respectively. Compared toFigure 4.1, in Figures 4.2 and 4.3 we observe that the system can stay in one statefor a long time, giving information about just one of the functions during that rangeof x values. This is more pronounced in Figure 4.3 when a12 and a21 are small.4.2 Initial valuesTo analyze the data via our EM algorithm, we need to provide initial values of allof the parameter estimates.We set the initial estimates of all of the parameters governing the distributionof the zi?s to 0.5.53In order to automatically obtain initial values for f1(x) and f2(x) we first assigntemporary values to the latent variables, creating two groups of observations (oneconsisting of all (xi,yi)?s with temporary zi value equal to 1, the other consistingof the remaining (xi,yi)?s). For each group we calculate an initial estimate of thecorresponding f j(x). For the Bayesian approach our initial estimate of f j(x) is itsposterior mean, commonly used in Gaussian regression (Rasmussen and Williams,2006), with initial ?2 = 0.005. For the penalized log-likelihood approach we usesmoothing splines to obtain the initial estimates of the f j?s. The smoothing param-eters for both approaches are chosen by cross-validation (CV). For the Bayesianapproach we also add the restriction that if the ? j obtained by CV is smaller than10 we use a bigger value, in this case 15. Two methods are used to assign thetemporary values to the zi?s.? Function estimate method: we fit one curve, m?(?), to the whole data set usinga standard cubic smoothing spline. If yi ? m?(xi), we set zi = 1, otherwise weset zi = 2.? Residual-based method: we obtain m?(?) as in the function estimate methodand calculate the residuals yi? m?(xi). We then divide the evaluation intervalinto small subintervals. Within each subinterval, we consider all residualscorresponding to xi?s within that subinterval. We use the k-means algorithm(k = 2) to partition these residuals into two groups. We label the group withthe smaller mean as group 1 (zi?s = 1). The k-means algorithm is a clusteringmethod that aims to partition observations into k groups such that the sumof the squared differences between each observation and its assigned groupmean is minimized (Johnson and Wichern, 2008).The function estimate method is used in Simulations 1 and 2, and the residual-based method in Simulation 3. The green curves in Figures 4.1, 4.2 and 4.3 areexamples of initial functions.In Simulation 3, we use the residual-based method because the probabilities ofchanging from one state to another are small, that is, the z process tends to remainin one state for an extended period, as shown in Figure 4.3. During this period, weonly obtain information from one of the f j?s and in this case the function estimate54method fails to produce good initial values. The residual-based method requiresa choice of number of sub-intervals and sub-interval lengths. Our choice of sub-intervals is based on examination of results from a few data sets. The chosensub-intervals are as follows: [1, 34], [34.5, 67.5] and [68, 100].A wide range of reasonable initial estimates of f1 and f2 yields good final esti-mates. For all data sets considered we can always find reasonable initial estimatesby eye. For most data sets the proposed automatic methods work. However, for afew data sets the automatic procedures produce obviously bad initial estimates thatdo not allow the method to recover.To obtain an initial estimate of ?2 we first use the initial function estimatesand the two temporary groups of observations to obtain the sample variance of theresiduals of each group separately adjusting for the correct degrees of freedom. Wethen set the initial estimate of ?2 equal to the pooled variance.4.3 Choice of the smoothing parameters, the ? j?sWe find the optimal ? j?s iteratively starting with initial values ? (0)j , j = 1, . . . ,J.The choice of the ? (0)j ?s is important to obtain good final estimates of the f j?s. Wefound out that very small ? (0)j ?s do not lead to good final estimates as more pointstend to be initially misclassified. Recall that for the Bayesian approach we restrictthe smoothing parameters of the initial function estimates to be greater than 10.So, for the Bayesian approach, in each data set we set the ? (0)j ?s to the values usedto obtain the initial function estimates. For the penalized log-likelihood approach,we use one value of ? (0)j for all data sets: we set ?(0)j to a value that worked wellwhen tested on a couple of simulated data sets.We update the ? j?s as follows.1. With ? j = ? (i)j , j = 1, . . . ,J, use the EM algorithm of Chapter 2 to find thep?i j?s, ??2 and the f? j?s.2. Discard the f? j?s from Step 1.3. Treat ??2 and the p?i j?s as fixed and thus W j as in (2.15) as fixed. For each ?55on a grid G and each j = 1, . . . ,J, calculatef? j? (x) = H j(? )y,where for the Bayesian approach H j(? ) = (I+A(? )W j)?1A(? )W j, and forthe penalized log-likelihood approach H j(? ) = B(BT W jB+2?R)?1BT W j.4. For each j = 1, . . . ,J, set ? (i+1)j as the value in the grid G that minimizes thefollowing cross-validation (CV) criterion:CVj(? ) =1nn?i=1p?i j(yi? f? j? (xi)1? (H j? )ii)2,where (H j? )ii is the ith entry of the diagonal of H j(? ). This CV criterion issimilar to the so-called weighted CV score presented by Silverman (1985).5. Repeat 1-4 with ? j = ? (i+1)j , j = 1, . . . ,J, till convergence.We use the final values of the ? j?s to obtain all of the parameter estimates from theEM algorithm as in Step 1.4.4 ResultsFigures 4.1, 4.2 and 4.3 show the fitted values f?1(x) and f?2(x) (dashed curves) fora simulated data set from each of simulation studies 1, 2 and 3, respectively.We assess the quality of the estimated functions via the pointwise empiricalmean squared error (EMSE). The empirical mean squared error of f? j at a givenpoint xi is calculated as follows:EMSE j(xi) =1SS?s=1[f? sj (xi)? f j(xi)]2, (4.1)where f? sj is the estimate of f j in the sth simulated data set and S is the total numberof simulated data sets, in this case S = 300. We calculate the EMSE for both initialand final estimates of the f j?s. Figures 4.4 and 4.5 present the pointwise EMSEof both initial and final estimates of f1 and f2 using the Bayesian and penalized56log-likelihood approaches for simulation studies 1 and 3, respectively. The plotsfor Simulation 2 are omitted because they are very similar to the ones obtainedfor Simulation 1. In all three simulation studies we observe the presence of edgeeffects, that is, the EMSE j values are higher at the edges than at the middle of theevaluation interval. We also observe that in all simulation studies the final estimatesproduce smaller values of EMSE j than the initial estimates, which indicates thatthe proposed methodology improves the initial naive estimates.Figures 4.6 and 4.7 present the EMSE of just the final estimates of f1 andf2 using the Bayesian and penalized log-likelihood approaches for Simulations 1and 3, respectively. The plots from Simulation 2 are again omitted as they arevery similar to the ones from Simulation 1. When we compare the Bayesian andthe penalized log-likelihood approaches we observe that for Simulations 1 and 2the Bayesian approach produces slightly smaller values of EMSE j for the finalfunction estimates than the penalized log-likelihood approach. In Simulation 3the Bayesian approach produces slightly smaller values of EMSE for f?1 than thepenalized log-likelihood approach except for the right edge. In Simulation 3, thereis no clear winner for estimating f2. In addition, both approaches produce valuesof EMSE for f?2 that are larger than the values obtained in Simulations 1 and 2.To further study the quality of our method, we look at possible cases of misclas-sification. We consider values of p?(zi = z|yi) > 0.2 when the true zi 6= z. Table 4.2presents the number of simulated data sets that have values of p?(zi = 1|yi) > 0.2when zi = 2. Table 4.3 shows the number of simulated data sets with values ofp?(zi = 2|yi) > 0.2 when zi = 1. For Simulations 1 and 2 all misclassifications oc-cur at the edges of the evaluation interval, indicating that the proposed methodsometimes has problems at the edges. The same is true in Simulation 3 for the pe-nalized log-likelihood approach. In Simulation 3, the Bayesian approach leads toproblems not just at the edges: two data sets have large values of p?(zi = 2|yi)> 0.2when zi = 1 at the middle of the evaluation interval when f1 and f2 are closer.Table 4.4 presents the mean and standard deviation of all 300 estimates of ?2for each simulation study considering both the Bayesian and the penalized log-likelihood approaches. All estimates are obtained adjusting the degrees of freedomto account for the estimation of the f j?s. We can observe that the means for theBayesian approach are closer to the true value of ?2 (5? 10?5) than the means57obtained using the penalized log-likelihood approach. The medians (not includedin the table) for the Bayesian approach are also closer to 5?10?5 than the mediansfor the penalized log-likelihood approach. See also Figure 4.8.Table 4.5 contains the mean and the standard deviation of the estimates ofthe parameters of the zi?s for each simulation study considering both the Bayesianand the penalized log-likelihood approaches. Note that the standard deviations ofthe estimates are close to the values of the means of the proposed standard errors(s.e.?s), as desired. Table 4.5 also shows the empirical coverage percentages of botha 90% and a 95% confidence interval. We consider confidence intervals of the formmean of the parameter estimates ? z?/2?proposed s.e.,where z?/2 is the ?/2 quantile of a standard normal distribution with ? = 0.1 and0.05. The empirical coverage percentages for Simulations 1 and 2 are very closeto the true level of the corresponding confidence interval. In Simulation 3 this isnot case. In particular, some of the 90% confidence intervals for a21 are basedon estimates of a21 that are so poor that even the 95% confidence intervals do notcontain the true value of a21. In this case the 95% confidence intervals have roughlythe same empirical coverage as the 90% confidence intervals.Figures 4.9, 4.10 and 4.11 present the box plots of the estimates of the param-eters of the z process for Simulations 1, 2 and 3, respectively.58Table 4.1: Summary of the simulation studies.Sim?ntypeof zi?szi?sparameters?2j ?sf j(x)?MV N(0,A(? j))1 iidp1 = 0.7?21 = ?22 =5?10?5?1 = 28?2 = 38p2 = 0.32Markovpi1 = pi2 = 0.5a12 = 0.3a21 = 0.43pi1 = pi2 = 0.5a12 = 0.1a21 = 0.2Table 4.2: Number of simulated data sets with values of p?(zi = 1|yi) > 0.2when zi = 2.Approach Simulation 1 Simulation 2 Simulation 3Bayes 0? 0? 4Penalized log-likelihood 0? 0? 1? p?(zi = 1|yi) was never larger than 0.04 when the true zi = 2.Table 4.3: Number of simulated data sets with values of p?(zi = 2|yi) > 0.2when zi = 1.Approach Simulation 1 Simulation 2 Simulation 3Bayes 3 3 14Penalized log-likelihood 2 1 1359Table 4.4: Estimates of ?2 (true value = 5?10?5).Simulation approach mean?105 (SD??105)1 Bayesian 4.984 (0.491)PL?? 4.912 (0.498)2 Bayesian 4.982 (0.556)PL 4.916 (0.579)3 Bayesian 4.935 (0.517)PL 4.880 (0.519)?SD = standard deviation, ??PL = penalized log-likelihood.Table 4.5: Estimates of the parameters of the z process.empirical coverageSim?n z parameters approach mean (SD?) mean of s.e.?s 90% 95%1 p1 = 0.7 Bayesian 0.699 (0.032) 0.032 90.7% 95.7%PL?? 0.699 (0.032) 0.032 90.7% 95.7%2 a12 = 0.3 Bayesian 0.300 (0.043) 0.043 90.0% 94.3%PL 0.300 (0.043) 0.043 90.0% 94.3%a21 = 0.4 Bayesian 0.399 (0.053) 0.053 90.3% 95.7%PL 0.399 (0.053) 0.053 90.3% 95.7%3 a12 = 0.1 Bayesian 0.105 (0.025) 0.027 93.0% 97.3%PL 0.105 (0.024) 0.027 93.3% 97.3%a21 = 0.2 Bayesian 0.208 (0.052) 0.050 91.7% 92.7%PL 0.208 (0.052) 0.050 91.7% 92.7%?SD = standard deviation, ??PL = penalized log-likelihood.60lllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll0 20 40 60 80 100?0.20?0.15?0.10?0.050.00xf j(x)(a) [binary] 40 60 80 100?0.20?0.15?0.10?0.050.00xf j(x)(b) Penalized log-likelihoodFigure 4.1: Simulation 1 (iid zi?s). Example of simulated data along with theinitial and final estimates of f1 and f2 obtained using in (a) the Bayesianapproach and in (b) the penalized log-likelihood approach. The red dotscorrespond to z = 2 and the black dots to z = 1. The solid curves cor-respond to the true functions f1 and f2. The green curves are the initialestimates of f1 and f2 obtained using the function estimate method. Thedashed curves are the final estimates of f1 and f2.61lllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll0 20 40 60 80 100?0.20?0.15?0.10?0.050.00xf j(x)(a) [binary] 40 60 80 100?0.20?0.15?0.10?0.050.00xf j(x)(b) Penalized log-likelihoodFigure 4.2: Simulation 2 (Markov zi?s with a12 = 0.3 and a21 = 0.4). Exam-ple of simulated data along with the initial and final estimates of f1 andf2 as in Figure 4.1.lllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll0 20 40 60 80 100?0.20?0.15?0.10?0.050.00xf j(x)(a) [binary] 40 60 80 100?0.20?0.15?0.10?0.050.00xf j(x)(b) Penalized log-likelihoodFigure 4.3: Simulation 3 (Markov zi?s with a12 = 0.1 and a21 = 0.2). Exam-ple of simulated data along with the initial and final estimates of f1 andf2 as in Figure 4.1.620 20 40 60 80 1000.0e+005.0e?061.0e?051.5e?052.0e?05xEMSE(x)(a) f1: Bayesian0 20 40 60 80 100xEMSE(x)1e?052e?053e?05(b) f1: penalized log-likelihood0 20 40 60 80 100xEMSE(x)1e?042e?043e?04(c) f2: Bayesian0 20 40 60 80 1000e+002e?044e?046e?04xEMSE(x)(d) f2: penalized log-likelihoodFigure 4.4: Simulation 1 (iid zi?s). EMSE of both initial (dashed curves) andfinal (solid curves) estimates of f1 and f2 using the Bayesian and thepenalized log-likelihood approaches. Plots (a) and (b) show the resultsfor f1 and (c) and (d) the results for f2.630 20 40 60 80 1000e+001e?042e?043e?044e?04xEMSE(x)(a) f1: Bayesian0 20 40 60 80 1000e+002e?044e?046e?04xEMSE(x)(b) f1: penalized log-likelihood0 20 40 60 80 1000e+001e?042e?043e?04xEMSE(x)(c) f2: Bayesian0 20 40 60 80 1000e+002e?044e?046e?048e?04xEMSE(x)(d) f2: penalized log-likelihoodFigure 4.5: Simulation 3 (Markov zi?s with a12 = 0.1 and a21 = 0.2). EMSEof both initial (dashed curves) and final (solid curves) estimates of f1and f2 as in Figure 4.4.640 20 40 60 80 1000.0e+005.0e?061.0e?051.5e?052.0e?052.5e?053.0e?05xEMSE(x)(a) EMSE f?10 20 40 60 80 1000e+001e?052e?053e?054e?055e?056e?05xEMSE(x)(b) EMSE f?2Figure 4.6: Simulation 1 (iid zi?s). EMSE of the final estimates of (a) f1and (b) f2 using the Bayesian (dashed curves) and the penalized log-likelihood (solid curves) approaches.0 20 40 60 80 1000e+002e?054e?056e?058e?051e?04xEMSE(x)(a) EMSE f?10 20 40 60 80 100xEMSE(x)01e?042e?043e?04(b) EMSE f?2Figure 4.7: Simulation 3 (Markov zi?s with a12 = 0.1 and a21 = 0.2). EMSEof the final estimates of (a) f1 and (b) f2 as in Figure 4.6.65l3.5e?054.0e?054.5e?055.0e?055.5e?056.0e?05Bayesian Penalized log?likelihood(a) Simulation 1llll3.5e?054.5e?055.5e?056.5e?05Bayesian Penalized log?likelihood(b) Simulation 2ll3.5e?054.0e?054.5e?055.0e?055.5e?056.0e?056.5e?05Bayesian Penalized log?likelihood(c) Simulation 3Figure 4.8: Box plot of the estimates of ?2 for each simulation study usingthe Bayesian and penalized log-likelihood approaches. The dashed linescorrespond to the true value 5?10?5.660.650.700.75Bayesian Penalized log?likelihoodFigure 4.9: Simulation 1 (iid zi?s). Box plot of the estimates of p1 ob-tained using the Bayesian and penalized log-likelihood approaches. Thedashed line corresponds to the true value 0.7.67llllll lllll0.200.250.300.350.40Bayesian Penalized log?likelihood(a) Box plot of a?12llll0.300.350.400.450.500.550.60Bayesian Penalized log?likelihood(b) Box plot of a?21Figure 4.10: Simulation 2 (Markov zi?s with a12 = 0.3 and a21 = 0.4). Boxplots of the estimates of (a) a12 and (b) a21 obtained using the Bayesianand the penalized log-likelihood approaches. The dashed lines corre-spond to the true values 0.3 and 0.4, respectively.llll lllll0.040.060.080.100.120.140.160.18Bayesian Penalized log?likelihood(a) Box plot of a?12llllllllll0.100.150.200.250.300.350.40Bayesian Penalized log?likelihood(b) Box plot of a?21Figure 4.11: Simulation 3 (Markov zi?s with a12 = 0.1 and a21 = 0.2). Boxplots of the estimates of (a) a12 and (b) a21 obtained using the Bayesianand the penalized log-likelihood approaches. The dashed lines corre-spond to the true values 0.1 and 0.2, respectively.68Chapter 5Simulation studies: the replicatecaseIn this chapter we present the simulation studies for the replicate case. We considerthat the zik?s can take values 1 or 2 and they can be either iid or follow a Markovstructure. We consider a covariance matrix V generated by a random interceptmodel as in (3.2). The parameters of interest are estimated using the penalizedlog-likelihood approach described in Chapter 3. For each simulation study 300independent data sets are generated, each containing 100 independent replicates.Note that all tables and figures for the replicate case are located at the end of thischapter.5.1 Simulated dataWe consider two types of simulation studies according to two different types ofsimulated data. In the first study we generate data considering iid zik?s and inthe second study we consider Markov zik?s. In all studies we use the same vectorof evaluation points x and the same true functions f1 and f2. The vector x =(x1, . . . ,xn)T consists of n = 10 equally spaced points, 1,12, . . . ,89,100. The truefunction f2 is the same f2 function used in the studies for the one realization caseand f1 is now simply f2?0.1.We generate a set of simulated data as follows.691. Generate the zik?s according to the specified model (iid or Markov).2. Generate the yik?s as follows:? if zik = 1, yik = f1(xi)+?k + eik;? if zik = 2, yik = f2(xi)+?k + eik,where ?k and eik are independent, ?k has a normal N(0,?2) distribution and?ik has a N(0,?2) distribution. We set ?2 = 10?4 and ?2 = 10?5.3. Repeat steps 1 and 2 N times obtaining N replicates.4. Repeat steps 1, 2 and 3 S times obtaining S simulated data sets.We consider S = 300 and N = 100 due to the computational complexity of ourapproach, which in this case involves calculating p(zk = s|yk,? (c)) for all possibles?s. Note again that f1(x) and f2(x) are fixed across all replicates and data sets.In our first simulation study we generate data assuming that the zik?s are inde-pendent with p1 = p(zik = 1) = 0.5. Figure 5.1 shows an example of a data setof this type displaying 20 of its 100 replicates. Figure 5.2 presents four replicatesfrom Figure 5.1 plotted separately.In the second simulation study we consider Markov zik?s with transition prob-abilities a12 = p(zi = 2|zi?1 = 1) = 0.3 and a21 = p(zi = 1|zi?1 = 2) = 0.4 andinitial probabilities pi1 = pi2 = 0.5. Figure 5.3 presents an example of a data setfrom the second simulation study with 20 of its 100 replicates. Figure 5.4 showsfour replicates from Figure 5.3 plotted separately.We analyze the data using the penalized log-likelihood approach and the pro-posed EM algorithm. To apply our EM algorithm, we need to provide initial valuesof all of the parameter estimates. In this case we set the initial values to be the trueparameter values to speed up computation. Note that we did try initial values thatwere different than the true parameter values and the EM algorithm also converged,but it took longer than starting from the truth as expected.We did not choose the values of ?1 and ?2 by a data based method. Rather weconsider the smoothing parameters ?1 and ?2 to be fixed and equal to 10?4 in bothsimulation studies. The value 10?4 was chosen by eye when tested on a couple ofsimulated data sets.705.2 ResultsFigures 5.1 and 5.3 show the fitted values f?1(x) and f?2(x) (dashed curves) for asimulated data set from each of simulation studies 1 and 2, respectively.We assess the quality of the estimated functions via the pointwise empiricalmean squared error (EMSE) as in (4.1). We observe that in both simulation studiesf?1(x) and f?2(x) produce very small values of EMSE. See Figure 5.5.Tables 5.1 and 5.2 present the mean and standard deviation of all 300 estimatesof ?2 and ?2, respectively, for each simulation study. In addition, Figures 5.6and 5.7 show the box plots of the estimates of ?2 and ?2 for Simulations 1 and 2,respectively. We can observe that in both simulation studies we are underestimatingthe values of ?2 and ?2 as in this case we do not adjust the degrees of freedom toaccount for the estimation of the f j?s.Figure 5.8 shows the box plot of all 300 estimates of p1. Figures 5.9 and 5.10present the box plots of the estimates of the initial probability pi1 and transitionprobabilities a12 and a21, respectively. We observe in all cases that the median ofthe estimates is very close to the true parameter value.Table 5.3 contains the mean and the standard deviation of the estimates of theparameters of the latent process for each simulation study. Note that the standarddeviations of the estimates are close to the values of the means of the proposedstandard errors (s.e.?s), as desired. Table 5.3 also shows the empirical coveragepercentages of both a 90% and a 95% confidence interval. As in the one realizationcase, we consider confidence intervals of the formmean of the parameter estimates ? z?/2?proposed s.e.,where z?/2 is the ?/2 quantile of a standard normal distribution with ? = 0.1 and0.05. The empirical coverage percentages for both simulation studies are very closeto the true level of the corresponding confidence interval.71Table 5.1: Estimates of ?2 (true value = 1?10?5).Simulation mean?105 (SD??105)1 0.978 (0.046)2 0.978 (0.045)?SD = standard deviation.Table 5.2: Estimates of ?2 (true value = 1?10?4).Simulation mean?104 (SD??104)1 0.977 (0.152)2 0.977 (0.152)?SD = standard deviation.Table 5.3: Estimates of the parameters of the z process.empirical coverageSim?n z parameters mean (SD?) mean of s.e.?s 90% 95%1 p1 = 0.5 0.499 (0.016) 0.016 90.3% 95.7%2 pi1 = 0.5 0.502 (0.050) 0.050 89.7% 95.7%a12 = 0.3 0.300 (0.021) 0.020 90.0% 94.3%a21 = 0.4 0.401 (0.024) 0.025 89.7% 95.3%?SD = standard deviation.720 20 40 60 80 100?0.15?0.10?0.050.00xf j(x)Figure 5.1: Simulation 1 (iid zik?s with p1 = 0.5). Example of simulated dataalong with f?1(x) and f?2(x). The gray dashed curves correspond to 20out of the 100 generated replicates. The black and red solid curvescorrespond to the true functions f1 and f2, respectively, evaluated onlyat x. The black and red dashed curves correspond to f?1(x) and f?2(x),respectively.73llllllllll0 20 40 60 80 100?0.15?0.10?0.050.00xf j(x)llllllllll0 20 40 60 80 100?0.15?0.10?0.050.00xf j(x)ll llllllll0 20 40 60 80 100?0.15?0.10?0.050.00xf j(x)llllll l lll0 20 40 60 80 100?0.15?0.10?0.050.00xf j(x)Figure 5.2: Simulation 1 (iid zik?s with p1 = 0.5). Example of four simulatedreplicates from data set shown in Figure 5.1. The red dots correspond toz = 2 and the black dots to z = 1. The black and red curves correspondto the true functions f1 and f2, respectively.740 20 40 60 80 100?0.15?0.10?0.050.00xf j(x)Figure 5.3: Simulation 2 (Markov zik?s with a12 = 0.3 and a21 = 0.4). Exam-ple of simulated data along with f?1(x) and f?2(x) as in Figure 5.1.75ll llll llll0 20 40 60 80 100?0.15?0.10?0.050.000.05xf j(x)lllllll lll0 20 40 60 80 100?0.15?0.10?0.050.000.05xf j(x)llllll l lll0 20 40 60 80 100?0.15?0.10?0.050.000.05xf j(x)llllll l lll0 20 40 60 80 100?0.15?0.10?0.050.000.05xf j(x)Figure 5.4: Simulation 2 (Markov zik?s with a12 = 0.3 and a21 = 0.4). Exam-ple of four simulated replicates from data set shown in Figure 5.3 as inFigure 5.2.760 20 40 60 80 1009.0e?071.0e?061.1e?061.2e?061.3e?061.4e?061.5e?06xEMSE(x)(a) Simulation 10 20 40 60 80 1009.0e?071.0e?061.1e?061.2e?061.3e?061.4e?061.5e?06xEMSE(x)(b) Simulation 2Figure 5.5: EMSE of f?1(x) (black curves) and f?2(x) (red curves) for (a) Sim-ulation 1 and (b) Simulation 2.llll8.50e?069.00e?069.50e?061.00e?051.05e?051.10e?05(a) Box plot of ??2l0.000060.000080.000100.000120.00014(b) Box plot of ??2Figure 5.6: Simulation 1 (iid zik?s). Box plots of the estimates of (a) ?2 and(b) ?2. The dashed lines correspond to the true values 10?5 and 10?4,respectively.77lllll8.50e?069.00e?069.50e?061.00e?051.05e?051.10e?05(a) Box plot of ??2l0.000060.000080.000100.000120.00014(b) Box plot of ??2Figure 5.7: Simulation 2 (Markov zik?s with a12 = 0.3 and a21 = 0.4). Boxplots of the estimates of (a) ?2 and (b) ?2 as in Figure 5.6.0.460.480.500.520.54Figure 5.8: Simulation 1 (iid zik?s with p1 = 0.5). Box plot of the estimatesof p1. The dashed line corresponds to the true value 0.5.78ll0.400.450.500.550.600.65Figure 5.9: Simulation 2 (Markov zik?s with a12 = 0.3 and a21 = 0.4). Boxplot of the estimates of pi1. The dashed line corresponds to the true value0.5.l0.260.280.300.320.340.36(a) Box plot of a?12ll0.350.400.45(b) Box plot of a?21Figure 5.10: Simulation 2 (Markov zik?s with a12 = 0.3 and a21 = 0.4). Boxplots of the estimates of (a) a12 and (b) a21. The dashed lines corre-spond to the true values 0.3 and 0.4, respectively.79Chapter 6ApplicationsIn this chapter we apply the proposed methodology to two different data sets: themotorcycle data set (Figure 1.3) and the building data set (Figure 1.2). The ap-plication to the motorcycle data is an example of the one realization case wherewe treat the data as coming from one curve that switches among J > 1 functions,one for each simulated accident. The fit to the building data is an example of thereplicate case, where we have N curves, the replicates, switching between J = 2functions.6.1 Motorcycle data6.1.1 Data set backgroundThe so-called motorcycle data set (Figure 1.3) is a well-known and widely useddata set, especially in the fields of nonparametric regression and machine learning.It consists of n = 133 measurements of head acceleration (in g) taken through time(in milliseconds) after impact in simulated motorcycle accidents. A table contain-ing the raw data can be found in Ha?rdle (1990) and the data are also available inthe software R.The data were collected by Schmidt et al. (1981) and became very popularafter appearing in Silverman (1985). Since then many different methodologieshave been applied to the motorcycle data. A Google search shows that this data set80appears in more than one hundred articles and book chapters.The motorcycle data set seems to be quite popular as an example among re-searchers in areas involving choice of smoothing parameter (e.g., Ha?rdle and Mar-ron, 1995), heterogeneity of the variance (e.g., Silverman, 1985) and estimation ofa regression function (or its first derivative) with jump discontinuities (e.g., Gijbelsand Goderniaux, 2004). It is important to say that all these analyses of the mo-torcycle data have something in common: they all treat the 133 measurements ascoming from one simulated accident. However, as we discussed in the introduc-tion, a close examination of the data suggests the measurements are from J > 1accidents.More recently the motorcycle data set has become a benchmark data set formachine learning techniques involving mixtures of Gaussian processes (e.g., Yangand Ma, 2011, Schiegg et al., 2012 and La?zaro-Gredilla et al., 2012). These tech-niques still treat the data as coming from one simulated accident and use the dif-ferent Gaussian processes only as a mechanism to account for the heterogeneity ofthe variance due to different phases in the data. Therefore, they do not considerthe possibility that the different processes in the mixture can actually correspondto multiple runs of accidents. As an example consider the work done by Schiegget al. (2012). These authors present a machine learning technique called a Markovlogic mixture of Gaussian processes. They apply their proposed methodology tothe motorcycle data in order to fit a single function that takes into account threephases in the data, which they call riding, impact and hitting the ground. Indeed,when, for comparison with their method, Schiegg et al. fit a mixture of three Gaus-sian processes (or experts) as in Tresp (2001) they write ?the experts of the mixtureof Gaussian processes have no specific meaning.?In this section we fit our proposed methodology to the motorcycle data settreating the measurements as coming from J > 1 functions (one for each simulatedaccident) with hidden (unknown) function labels. We choose J using an ad hocAkaike?s information criterion (AIC). Unlike the machine learning literature, weprovide standard errors for the parameters governing the latent switching process.Note that all tables and figures are located in the end of Section 6.1.2.816.1.2 Data analysisWe fit the proposed switching nonparametric regression model to the motorcycledata considering the one realization case described in Chapter 2. We assume thatthe hidden states, which correspond to the unknown accident run labels, are iid.We estimate the model parameters using both the Bayesian and the penalized log-likelihood approaches. The smoothing parameters, the ? j?s, are selected by cross-validation as in Section 4.3. For the Bayesian approach we set the covarianceparameter U j in (1.4) to be fixed atU j =U =?ni=1(yi? y?)2n?1??ni=1(yi? f? (xi))2n? trace(H),where f? is a regular smoothing spline fit to the data and H is the corresponding hat(or smoothing) matrix. So, U j is ?known? and ? j = s j.The data set contains 39 time points with multiple acceleration measurements,which our method does not currently allow. So, for a direct application of ourmethodology to the data, we jitter each of those time points by adding a very smallrandom noise.We fit the model considering J = 2, . . . ,6 functions. We choose J to minimizethe ad hoc AIC?2log p(y|??)+2( J?j=1trace(H j)+number of estimated variances+ J?1),where H j is the hat matrix corresponding to f? j(x).We apply the proposed model considering both equal and different regressionerror variances. However, the choice of J is not so obvious when we use a com-mon regression error variance. Therefore, the results presented here are obtainedconsidering different error variances.For the Bayesian approach the AIC is minimum for a model with J = 4. How-ever, for the penalized log-likelihood it is minimum for a model with J = 3. As bothresults appear sensible, we present the results for both J = 3 and J = 4. Figures 6.1aand 6.1b show the estimated functions obtained using the penalized log-likelihoodapproach for J = 3 and J = 4, respectively, and Figures 6.2a and 6.2b show the82same information using the Bayesian approach.Table 6.1 presents the estimated model parameters when J = 3 for both theBayesian and the penalized log-likelihood approaches. The table also shows thechosen smoothing parameter used in the estimation of each function. We can ob-serve that there is some qualitative agreement between the results from the twoapproaches. For instance, the green curve has the largest variance. The mixingproportion estimates agree, well within the reported standard errors. Although thevalues of the ? j?s are not comparable between the two approaches, we see that forboth methods, the green curve has the smallest value of ? j, indicating that the greencurve is the least smooth curve.Table 6.2 is similar to Table 6.1 and presents the results for J = 4.Based on the obtained results we believe the motorcycle data come from J = 3or 4 different simulated accident runs with different regression error variances. IfJ = 3 then roughly 27% ? 5%, 40% ? 5% and 34% ? 5% of the data come fromsimulated accidents corresponding to the black, red and green curves in Figure6.1a, respectively. If J = 4 then roughly 22% ? 4%, 22% ? 5%, 20% ? 4% and36% ? 5% of the data come from simulated accidents corresponding to the black,red, green and blue curves in Figure 6.2b, respectively.83Table 6.1: Results for J = 3 with corresponding fitted curves in Figures 6.1aand 6.2a.Approach curve ??2j p? j (s.e.) ? jPenalized black 43.054 0.269 (0.047) 0.893log-likelihood red 14.227 0.395 (0.053) 0.890(Fig. 6.1a) green 171.048 0.337 (0.053) 0.167Bayesian black 50.134 0.272 (0.047) 3.901(Fig. 6.2a) red 8.593 0.361 (0.050) 5.005green 184.478 0.367 (0.052) 2.512Table 6.2: Results for J = 4 with corresponding fitted curves in Figures 6.1band 6.2b.Approach curve ??2j p? j (s.e.) ? jPenalized black 50.348 0.211 (0.042) 1.591log-likelihood red 7.623 0.232 (0.050) 1.213(Fig. 6.1b) green 7.245 0.244 (0.051) 0.825blue 135.968 0.313 (0.049) 0.545Bayesian black 48.376 0.223 (0.043) 3.782(Fig. 6.2b) red 5.473 0.221 (0.045) 4.805green 1.106 0.198 (0.042) 5.377blue 145.440 0.358 (0.048) 3.35584llll lllllllllllll lllllllllllllllllllllllllllllllllllllllllllll llllllllllllllllllllllllllllllllllllllllllll lllllllllllll lll llllllll10 20 30 40 50?150?100?50050TimeHead acceleration(a) J = 3llll lllllllllllll lllllllllllllllllllllllllllllllllllllllllllll llllllllllllllllllllllllllllllllllllllllllll lllllllllllll lll llllllll10 20 30 40 50?150?100?50050TimeHead acceleration(b) J = 4Figure 6.1: Motorcycle data. Final function estimates (solid curves) obtainedusing the penalized log-likelihood approach for (a) J = 3 and (b) J =4. The gray dashed curves correspond to the initial function estimates,which are constant shifts of a smoothing spline fit to all the data.85llll lllllllllllll lllllllllllllllllllllllllllllllllllllllllllll llllllllllllllllllllllllllllllllllllllllllll lllllllllllll lll llllllll10 20 30 40 50?150?100?50050TimeHead acceleration(a) J = 3llll lllllllllllll lllllllllllllllllllllllllllllllllllllllllllll llllllllllllllllllllllllllllllllllllllllllll lllllllllllll lll llllllll10 20 30 40 50?150?100?50050TimeHead acceleration(b) J = 4Figure 6.2: Motorcycle data. Final function estimates (solid curves) obtainedusing the Bayesian approach for (a) J = 3 and (b) J = 4. The graydashed curves correspond to the initial function estimates, which areconstant shifts of a smoothing spline fit to all the data.866.2 Building DataThe building data shown in Figure 1.2 of Chapter 1 consist of daytime hourlypower usage of a building from 9am to 4pm (n = 8) on N = 44 business days inJune and July 2009. As already described in Section 1.1 the data present differenttypes of curves: one type corresponds to days with high power consumption whenthe chiller of the building was probably on and another type when the chiller wasprobably off. On some days the chiller turned on in the middle of the day. On oneday the chiller went on, off and then on again.We apply our proposed methodology to the building data treating each day as areplicate and modeling power usage as arising from J = 2 functions, one functiongiving power usage when the chiller is off ( j = 1), and the other function givingpower when the chiller is on ( j = 2). The condition ?chiller on?/?off? cannot beobserved directly and thus forms a latent process.In Section 6.2.1 we present the results assuming the observations within a repli-cate are uncorrelated and in Section 6.2.2 we show the results when we assume theobservations are correlated with covariance matrix V generated by a random in-tercept model. Note that all tables and figures are located in the end of Section6.2.2.6.2.1 Results: uncorrelated observationsIn this section we fit the proposed model to the building data assuming the ob-servations within a replicate are uncorrelated, that is, with diagonal covariancematrix. We consider two models for V: in one, V depends on the hidden states(Vzk = diag(?2z1k , . . . ,?2znk)) and in the other V = ?2I. We consider the penalizedlog-likelihood approach and use the ECM algorithm described in Section 3.2.1 toestimate the model parameters for both iid and Markov zik?s. The smoothing pa-rameters, the ? j?s, are chosen by cross-validation as described in detail in Section3.3.Figures 6.3 and 6.4 present the fitted functions for iid hidden states zik?s whenwe assume V = ?2I and Vzk = diag(?2z1k , . . . ,?2znk), respectively. We can observethat the fitted curves giving power usage when the chiller is off are very similarin the two figures. The estimated curve giving power usage when the chiller is on87obtained assuming V = ?2I is slightly smoother than the one obtained assumingVzk = diag(?2z1k , . . . ,?2znk).Tables 6.3 and 6.4 present the parameter estimates and chosen ? j?s for iid zik?swhen we assume V = ?2I and Vzk = diag(?2z1k , . . . ,?2znk), respectively. We can seethat the estimates of p j = p(zik = j) agree within the reported standard errors. Wealso observe in Table 6.4 that the estimated variance when the chiller is on is muchhigher than when the chiller is off.Figures 6.5 and 6.6 present the fitted curves for Markov zik?s when we assumeV = ?2I and Vzk = diag(?2z1k , . . . ,?2znk), respectively. As in the iid case, the fittedcurve giving power usage when the chiller is on obtained assuming V = ?2I isslightly smoother than the one obtained assuming Vzk = diag(?2z1k , . . . ,?2znk).Tables 6.5 and 6.6 show the estimated model parameters when the zik?s followa Markov structure with V = ?2I and Vzk = diag(?2z1k , . . . ,?2znk), respectively. Thetables also show the chosen smoothing parameters, the ? j?s. We observe in Table6.5 where we assume V =?2I that the estimate of a12 is very small and in Table 6.6that the estimate of a12 is zero. These estimates can be explained by consideringthe green curve in Figure 6.6. This green curve corresponds to the only replicate inthe data set where there is a transition from ?chiller on? to ?off?, and that transitionactually happens gradually, which our model does not allow. When we assumedifferent variances the method incorrectly classifies all observations from the greenreplicate in Figure 6.6 as coming from the condition ?chiller on?, failing to detectthe only transition from ?chiller on? to ?off?, and, therefore, getting zero as anestimate of a12.6.2.2 Results: correlated observationsIn this section we fit the proposed switching nonparametric regression model tothe building data assuming the observations within a replicate are correlated withcovariance matrix V generated by a random intercept model as in (3.2). We con-sider the penalized log-likelihood approach and use the ECM algorithm describedin Section 3.2 to obtain the parameter estimates. We conduct the analysis assum-ing the hidden states zik?s are iid or follow a Markov structure. The smoothingparameters, the ? j?s, are fixed and chosen by eye.88Tables 6.7 and 6.8 present the parameter estimates obtained assuming V isgenerated by a random intercept model for iid and Markov zik?s, respectively. Weobserve that the estimates of the covariance parameters ?2 and ?2 are very similarin the two tables. We also see that the estimates of p1 and p2 in Table 6.7 agreewithin the reported standard errors with the estimates obtained in Tables 6.3 and6.4 where we assume the observations are uncorrelated. Similarly, the estimates ofpi1, pi2, a12 and a21 in Table 6.8 agree with estimates obtained in Table 6.5.Figures 6.7 and 6.8 show the fitted functions when we assume V is generatedby a random intercept model for iid and Markov zik?s, respectively. We can observethat the fitted functions are very similar in the two figures. We can also see prob-lems with the obtained function estimates. The fitted functions corresponding tothe condition ?chiller on? are underestimating power usage and the fitted functionscorresponding to ?chiller off? are overestimating power usage.To obtain better function estimates for the building data one can consider acovariance structure generated by a random intercept with variances depending onthe hidden states. In this case, one would assume that yik = fzik(xi)+ ?zik, ik, where- if zik = 1 then ?1,ik = ?k + eik;- if zik = 2 then ?2,ik = ?k +?k + eik;- ?k, ?k, and eik are independent for i = 1, . . . ,n and k = 1, . . . ,N;- ?k?s are iid N(0,?21 );- ?k?s are iid N(0,?22 );- eik?s are iid N(0,?2).Therefore, the covariance matrix would be given byVzk = ?2I+ ?21 11T + ?21 1zk 1Tzk , (6.1)where 1 is an n-vector of ones and 1Tzk is an n-vector with entries I(zik = 2).Writing the steps of the ECM algorithm and fitting the building data consider-ing the covariance structure above is part of the proposed future work.89Table 6.3: Results for iid zik?s and V = ?2I with corresponding fitted curvesin Figures 6.3.curve (chiller condition, j) ??2 p? j (s.e.) ? jblack (off, j = 1)103.50.665 (0.025) 0.020red (on, j = 2) 0.335 (0.025) 0.078Table 6.4: Results for iid zik?s and Vzk = diag(?2z1k , . . . ,?2znk) with correspond-ing fitted curves in Figures 6.4.curve (chiller condition, j) ??2j p? j (s.e.) ? jblack (off, j = 1) 12.7 0.658 (0.025) 0.073red (on, j = 2) 355.4 0.342 (0.025) 0.006Table 6.5: Results for Markov zik?s and V = ?2I with corresponding fittedcurves in Figures 6.5.a?12(s.e.)a?21(s.e.)curve (chiller condition, j) ??2 p?i j (s.e.) ? jblack (off, j = 1)103.10.705 (0.069) 0.024(0.011)0.00991(0.00986)0.019red (on, j = 2) 0.295 (0.069) 0.084Table 6.6: Results for Markov zik?s and Vzk = diag(?2z1k , . . . ,?2znk) with corre-sponding fitted curves in Figures 6.6.curve (chiller condition, j) ??2j p?i j a?12 a?21 ? jblack (off, j = 1) 12.2 0.6820.015 < 10?160.050red (on, j = 2) 400.2 0.318 0.00690Table 6.7: Results for iid zik?s and V generated by a random intercept modelwith ?1 = ?2 = 0.1 and corresponding curves in Figures 6.7.curve (chiller condition, j) ??2 ??2 p? j (s.e.)black (off, j = 1)14.3 162.40.665 (0.025)red (on, j = 2) 0.335 (0.025)Table 6.8: Results for Markov zik?s and V generated by a random interceptmodel with ?1 = ?2 = 0.1 and corresponding curves in Figures 6.8.a?12(s.e.)a?21(s.e.)curve (chiller condition, j) ??2 ??2 p?i j (s.e.)black (off, j = 1)14.3 161.90.705 (0.069) 0.024(0.011)0.00990(0.00985)red (on, j = 2) 0.295 (0.069)919 10 11 12 13 14 15 1680100120140160180200220Hour of dayKW/hFigure 6.3: Building daytime power usage. Fitted function estimates (solidcurves) assuming iid zik?s and V = ?2I. The gray dashed curves corre-spond to the replicates. The red and black dashed curves are the initialfunction estimates. The colors red and black correspond to the conditionchiller on and off, respectively.929 10 11 12 13 14 15 1680100120140160180200220Hour of dayKW/hFigure 6.4: Building daytime power usage. Fitted function estimates assum-ing iid zik?s and Vzk = diag(?2z1k , . . . ,?2znk). Components of the plot areas in Figure 6.3.939 10 11 12 13 14 15 1680100120140160180200220Hour of dayKW/hFigure 6.5: Building daytime power usage. Fitted function estimates assum-ing Markov zik?s and V = ?2I. Components of the plot are as in Figure6.3.949 10 11 12 13 14 15 1680100120140160180200220Hour of dayKW/hFigure 6.6: Building daytime power usage. Fitted function estimates assum-ing Markov zik?s and Vzk = diag(?2z1k , . . . ,?2znk). Components of the plotare as in Figure 6.3. The green curve corresponds to the replicate wherethere is a transition from chiller on to off.959 10 11 12 13 14 15 1680100120140160180200220Hour of dayKW/hFigure 6.7: Building daytime power usage. Fitted function estimates assum-ing iid zik?s and V generated by a random intercept model. Componentsof the plot are as in Figure 6.3.969 10 11 12 13 14 15 1680100120140160180200220Hour of dayKW/hFigure 6.8: Building daytime power usage. Fitted function estimates assum-ing Markov zik?s and V generated by a random intercept model. Com-ponents of the plot are as in Figure 6.3.97Chapter 7Discussion and future workIn this thesis we proposed a model to analyze data arising from a curve that, overits domain, switches among J states. We called this model a switching nonpara-metric regression model. We considered two different data structures: one with Nreplicates and the other with one single realization. Two different approaches forthe estimation of the parameters were proposed: a Bayesian approach and a fre-quentist approach called penalized log-likelihood approach. For the hidden states,we considered those that are iid and those that follow a Markov structure.Overall our main contributions include the introduction and development of thefrequentist approach to the problem, the parameter estimation for the replicate case,the calculation of standard errors for the parameter estimates of the latent processand the study of the frequentist properties of the proposed estimates via simulation.We presented two different applications of the proposed methodology. In the firstapplication we analyzed the well-known motorcycle data in an innovative way:treating the data as coming from J > 1 simulated accident runs with unobservedrun labels. In the second application we analyzed the building data (Figure 1.2)treating each business day as a replicate and modeling power usage as arising fromtwo functions, one function giving power usage when the cooling system of thebuilding is off, the other function giving power usage when the cooling system ison.The R package switchnpreg was developed and it is available on CRAN(cran.r-project.org). The package can be used to obtain the parameter estimates for98the one realization case described in Chapter 2.Future work on different topics is presented as follows.1. Inclusion of the estimation of the parameters for the replicate case in theR package switchnpreg. The package is coded in a very efficient andmodular way, so that it is easy to make changes or add new things to it.2. Study of a cross-validation (CV) criterion when V is not diagonal. In Sec-tion 3.3 we derived a CV criterion to obtain the optimal smoothing param-eters when the observations within a replicate are uncorrelated, that is, withdiagonal covariance matrix. As future work we would like to study throughsimulations a CV criterion when the observations within a replicate are cor-related, that is, when the covariance matrix V is no longer diagonal. We havealready written down a weighted ?leave-one-replicate-out? CV criterion andfound a short-cut formula for it. We would like to investigate if the criterionworks well by testing it on a number of simulated data sets and on the build-ing data. The proposed CV criterion when V is not diagonal is presented asfollows.We propose to find the optimal ? j?s as the vector ? = (?1, . . . ,?J)T that min-imizes the following ?leave-one-replicate-out? CV criterion:CV (? ) =N?k=1(wk? f?(?k)?)TSk(wk? f?(?k)?), (7.1)where wk = S?1k DkV?1yk and f?(?k)? is the maximizer of (3.15) plus (3.16)leaving the kth replicate out. Similarly to the case when Vzk = diag(?2z1k , . . . ,?2znk) in Section 3.3, we can prove that (7.1) can be expressed by the follow-ing short-cut formula:CV (? ) =N?k=1[(I?Hk(? ))?1(wk? f?? )]TSk[(I?Hk(? ))?1(wk? f?? )],(7.2)where f?? = ?NK=1 Hk(? )wk is the maximizer of (3.15) plus (3.16) given by99(3.22) or (3.25) and Hk(? ) depends on the approach considered. For theBayesian approachHk(? ) =[diag(A(?1), . . . ,A(?J))M+ I]?1diag(A(?1), . . . ,A(?J))Sk,and for the penalized log-likelihood approachHk(? ) = B?[B?T MB?+2diag(?1R, . . . ,?JR)]?1B?T Sk.3. Extension to a covariance structure depending on the hidden states. Thederivation of Steps 1 and 2 of the ECM algorithm must be modified. Step3 remains the same as it depends on the covariance matrix only through thepk(s)(c)?s. In this case, for Steps 1 and 2 we would maximize S(? ,? (c)) onpage 33 with (3.12) replaced byC?12N?k=1?alls?spk(s)(c)[log |Vs|+(yk? fs(x))T V?1s (yk? fs(x))].First we would consider the case when the covariance matrix is generated bya random intercept model with variances depending on the hidden states asin (6.1). Then we would extend the approach to covariance matrices arisingfrom more general random effect models such as a random intercept andslope model.4. Extension of the current work to a latent process z depending on some co-variate(s). So, for example, if there are J = 2 states, we might suppose thatp(zi = 1|vi)? p(vi) follows a logit model with two parameters, that is,logp(vi)1? p(vi)= ?0 +?1vi.So ? = (?0,?1). In Step 3 of the ECM algorithm we must now maximizeE? (c)(L2(?)|y) =n?i=1[? log(1+ e(?0+?1vi))+ p(c)i1 (?0 +?1vi)](7.3)100with respect to ?0 and ?1, keeping the p(c)i1 ?s fixed. We can maximize (7.3)using a numerical optimization method, such as Newton-Raphson. Note thatSteps 1 and 2 of the ECM algorithm remain the same.5. Study and development of other criteria to select the number of states J. Forinstance, we would like to exploit a criterion within the Bayesian frameworksuch as the Bayes factor discussed by Kass and Raftery (1995). We havetried an ad-hoc Bayesian information criterion (BIC) or (Schwarz criterion)to choose J in the motorcycle data application, however, the BIC did notwork well in that case. The BIC is an approximation to the log Bayes factorand, therefore, it might be better in our case to directly use Bayes factor.As future work we would like to determine what sort of Bayes factor wouldwork well in our model.101BibliographyBaum, L., Petrie, T., Soules, G., and Weiss, N. (1970), ?A maximization techniqueoccurring in the statistical analysis of probabilistic functions of Markovchains,? The Annals of Mathematical Statistics, 41, 164?171. ? pages 18Bilmes, J. (1998), ?A gentle tutorial of the EM algorithm and its application toparameter estimation for Gaussian mixture and hidden Markov models,?International Computer Science Institute, 4, 126. ? pages 12, 19Cappe?, O., Moulines, E., and Ryde?n, T. (2005), Inference in Hidden MarkovModels, Springer Verlag. ? pages 19Chiou, J. (2012), ?Dynamical functional prediction and classification, withapplication to traffic flow prediction,? The Annals of Applied Statistics, 6,1588?1614. ? pages 12Demidenko, E. (2004), Mixed Models: Theory and Applications, John Wiley &Sons. ? pages 39Fan, J. and Gijbels, I. (1992), ?Variable bandwidth and local linear regressionsmoothers,? Annals of Statistic, 20, 2008?2036. ? pages 11Gijbels, I. and Goderniaux, A. (2004), ?Bootstrap test for change-points innonparametric regression,? Journal of Nonparametric Statistics, 16, 591?611.? pages 3, 81Ha?rdle, W. (1990), Applied Nonparametric Regression, Cambridge UniversityPress. ? pages 80Ha?rdle, W. and Marron, J. (1995), ?Fast and simple scatterplot smoothing,?Computational Statistics & Data Analysis, 20, 1?17. ? pages 2, 81Heckman, N. (2012), ?Reproducing Kernel Hilbert Spaces made easy,? StatisticsSurveys, 6, 113?141. ? pages 21102Johnson, R. A. and Wichern, D. W. (2008), Applied Multivariate StatisticalAnalysis, 6th Ed., Pearson. ? pages 54, 120Kass, R. E. and Raftery, A. E. (1995), ?Bayes factors,? Journal of the AmericanStatistical Association, 90, 773?795. ? pages 101La?zaro-Gredilla, M., Van Vaerenbergh, S., and Lawrence, N. D. (2012),?Overlapping mixtures of Gaussian processes for the data association problem,?Pattern Recognition, 45, 1386?1395. ? pages 81Louis, T. (1982), ?Finding the observed information matrix when using the EMalgorithm,? Journal of the Royal Statistical Society Series B, 44, 226?233. ?pages 23, 24, 47McLachlan, G. and Krishnan, T. (2008), The EM Algorithm and Extensions, 2ndEd., Wiley New York. ? pages 7, 15, 107, 110Meng, X. and Rubin, D. (1993), ?Maximum likelihood estimation via the ECMalgorithm: a general framework,? Biometrika, 80, 267?278. ? pages 108Ou, X. and Martin, E. (2008), ?Batch process modelling with mixtures ofGaussian processes,? Neural Computing & Applications, 17, 471?479. ?pages 11, 12Rabiner, L. (1989), ?A tutorial on hidden Markov models and selectedapplications in speech recognition,? Proceedings of the IEEE, 77, 257?286. ?pages 19Rasmussen, C. and Ghahramani, Z. (2002), ?Infinite mixtures of Gaussian processexperts,? in Advances in Neural Information Processing Systems 14:Proceedings of the 2001 Conference, The MIT Press, vol. 2, pp. 881?888. ?pages 3, 11, 12Rasmussen, C. E. and Williams, C. K. I. (2006), Gaussian Processes for MachineLearning, The MIT Press. ? pages 8, 54Schiegg, M., Neumann, M., and Kersting, K. (2012), ?Markov logic mixtures ofGaussian processes: towards machines reading regression data,? in Proceedingsof the 15th International Conference on Artificial Intelligence and Statistics,vol. 22. ? pages 81Schmidt, G., Mattern, R., and Schu?ler, F. (1981), ?Biomechanical investigation todetermine physical and traumatological differentiation criteria for themaximum load capacity of head and vertebral column with and without103protective helmet under the effects of impact,? EEC Research Program onBiomechanics of Impacts. Final report Phase III, Project G5, Institut fu?rRechtsmedizin, Universita?t Heidelberg. ? pages 3, 80Shi, J., Murray-Smith, R., and Titterington, D. (2005), ?Hierarchical Gaussianprocess mixtures for regression,? Statistics and Computing, 15, 31?41. ?pages 12Silverman, B. (1985), ?Some aspects of the spline smoothing approach tonon-parametric regression curve fitting,? Journal of the Royal StatisticalSociety Series B, 47, 1?52. ? pages 2, 14, 21, 56, 80, 81Tresp, V. (2001), ?Mixtures of Gaussian processes,? in Advances in NeuralInformation Processing Systems 13: Proceedings of the 2000 Conference, TheMIT Press, pp. 654?660. ? pages 11, 12, 81Wahba, G. (1983), ?Bayesian ?confidence intervals? for the cross-validatedsmoothing spline,? Journal of the Royal Statistical Society Series B, 45,133?150. ? pages 22Wood, S. N. (2011), ?Fast stable restricted maximum likelihood and marginallikelihood estimation of semiparametric generalized linear models,? Journal ofthe Royal Statistical Society Series B, 73, 3?36. ? pages 3Yang, Y. and Ma, J. (2011), ?An efficient EM approach to parameter learning ofthe mixture of gaussian processes,? in Proceedings of the 8th InternationalConference on Advances in Neural Networks, Springer-Verlag, pp. 165?174. ?pages 81104Appendix AThe proposed EM algorithm:why it worksWe use an EM algorithm to find ?? that maximizesl(?) = log p(y|?)+P( f1, . . . , fJ,?1, . . . ,?J), (A.1)where log p(y|?) is the observed data log-likelihood and P( f1, . . . , fJ,?1, . . . ,?J) isas in Section 1.2.The derivation of the EM algorithm can be summarized in two steps:1. l(?) is expressed in terms of the latent data;2. We use this expression to come up with a sequence of ? (c)?s which will notdecrease l(?), that is, with l(? (c+1))? l(? (c))? 0.Let S(? ,? (c)) be as on page 16:S(? ,? (c)) = Q(? ,? (c))+P( f1, . . . , fJ,?1, . . . ,?J),where Q(? ,? (c)) = E? (c)[log p(y,z|?)|y]=?zlog p(y,z|?)p(z|y,? (c)). We nowshow that any sequence {? (c), c ? 1}, that satisfies S(? (c+1),? (c)) ? S(? (c),? (c))also satisfies l(? (c+1))? l(? (c)).105The logarithm of the joint distribution of y and z (also called the complete datadistribution) can be written as p(y,z|?) = p(z|y,?)p(y|?) and using this in (A.1),we getl(?) = log p(y,z|?)? log p(z|y,?)+P( f1, . . . , fJ,?1, . . . ,?J). (A.2)Multiplying both sides of (A.2) by p(z|y,? (c)) and summing over all possiblevalues of z preserves the equality. The left side is unchanged, since l(?) does notdepend on z. Since the penalty term P( f1, . . . , fJ,?1, . . . ,?J) does not depend on zthe right hand side of (A.2) becomes?zlog p(y,z|?)p(z|y,? (c))??zlog p(z|y,?)p(z|y,? (c))+P( f1, . . . , fJ,?1, . . . ,?J),Thus, we can writel(?) = Q(? ,? (c))?H(? ,? (c))+P( f1, . . . , fJ,?1, . . . ,?J), (A.3)whereH(? ,? (c)) =?zlog p(z|y,?)p(z|y,? (c)).Therefore, we see that finding a value ? (c+1) so that l(? (c+1)) ? l(? (c)) is thesame as finding ? (c+1) so thatS(? (c+1),? (c))?S(? (c),? (c))? H(? (c+1),? (c))?H(? (c),? (c)). (A.4)Note that for any ?H(? ,? (c))?H(? (c),? (c)) = E? (c)(log{p(z|y,?)/p(z|y,?(c))}|y)? log[E? (c){p(z|y,?)/p(z|y,?(c))}|y](A.5)= log?zp(z|y,?)p(z|y,? (c))p(z|y,? (c))= log?zp(z|y,?) = log1= 0,106where the inequality in (A.5) is a consequence of Jensen?s inequality and the con-cavity of the logarithmic function.Therefore, we can guarantee that (A.4) holds if we guarantee that S(? (c+1),? (c))?S(? (c),? (c)) ? 0. To do that choose ? (c+1) as any value of ? which does not de-crease S(? ,? (c)) from the current value at ? (c). If ? (c+1) is the maximizing value,then S(? (c+1),? (c))? S(? (c),? (c)).To increase the rate of convergence of the algorithm we recommend taking? (c+1) = argmax?S(? ,? (c)) if possible.More details on the EM algorithm can be found in McLachlan and Krishnan(2008).107Appendix BThe ECM algorithmRecall that for the success of the EM algorithm, we require that S(? (c+1),? (c)) ?S(? (c),? (c)), that is, we do not need to maximize S(? ,? (c)), but merely increaseit. Sometimes the M-step of the EM algorithm is complicated and, therefore, thisalgorithm becomes less attractive. In many cases, however, the M-step becomesrelatively simple if the maximization is undertaken conditional on some of theparameters (or some functions of the parameters). Meng and Rubin (1993) usethis idea in an extension of the EM algorithm, which they call the ECM algorithmfor expectation-conditional maximization algorithm. The ECM algorithm replacesa complicated M-step of the EM algorithm with several computationally simplerconditional maximization steps (CM-steps). In our case, each CM-step maximizesS(? ,? (c)) = Q(? ,? (c)) + P( f1, . . . , fJ,?1, . . . ,?J) or at least does not decrease itfrom the current value at ? (c) subject to constraints on the set of parameters ? ,where the collection of all constraints is such that the maximization is over the fullparameter space of ? .B.1 Formal definitionInstead of updating ? (c) to ? (c+1) directly via one M-step, we replace the M-stepby R > 1 CM-steps. Let ? (c+r/R) denote the value of ? on the rth CM-step ofthe update, where ? (c+r/R) is the value that maximizes S(? ,? (c)) subject to the108constraintgr(?) = gr(? (c+(r?1)/R)).Here, {gr(?),r = 1, . . . ,R} is a set of pre-selected vector-valued functions typ-ically ?choosing? components of ? , for example, g1(?) = ?1. Thus, ? (c+r/R) satis-fiesS(? (c+r/R),? (c))? S(? ,? (c)) for all ? such that gr(?) = gr(? (c+(r?1)/R)). (B.1)In particular S(? (c+r/R),? (c))? S(? (c+(r?1)/R),? (c)).The value of ? on the final CM-step, ? (c+R/R) = ? (c+1), is taken to be the inputfor the next iteration. From (B.1) we haveS(? (c+1),? (c))? S(? (c+(R?1)/R),? (c))? S(? (c+(R?2)/R),? (c))? ?? ? ? S(? (c),? (c)).As noted before in Appendix A, this inequality is a sufficient condition for l(? (c+1))?l(? (c)) to hold. This shows that the ECM algorithm has the same desirable conver-gence properties as the EM algorithm.B.2 ExampleConsider the M-step presented in Section 2.2 for iid zi?s. We want to increase thevalue of S(? ,? (c)) with respect to ? = {p j,?2j , f j(x)} from the value at ? = ? (c).In this case we have a total of R = 3 CM-steps described as follows.? Step r = 1:Let g1(p j,?2j , f j(x)) be a partition of ? , that is, g1(p j,?2j , f j(x))= (p j,?2j )T .Thus, we want to maximize S(? ,? (c)) over ? subject tog1(p j,?2j , f j(x)) = g1(p(c)j ,?2(c)j , f j(x)(c)),i.e., maximize S(? ,? (c)) over f j(x) holding p j = p(c)j and ?2j = ?2(c)j . As aresult we obtain? (c+1/R) = {p(c+1/3)j ,?2(c+1/3)j , f j(x)(c+1/3)}= {p(c)j ,?2(c)j , f? j(x)}, (B.2)109where f? j(x) = H jy. The form of the hat matrix H j depends on the approachconsidered (Bayesian or penalized log-likelihood). See page 21.? Step r = 2:Let g2(p j,?2j , f j) = {p j, f j(x)}. We want to maximize S(? ,? (c)) over ?subject tog2(p j,?2j , f j) = g2(p(c+1/3)j ,?2(c+1/3)j , f(c+1/3)j ),i.e., maximize S(? ,? (c)) over ?2j holding p j = p(c+1/3)j = p(c)j and f j(x) =f j(x)(c+1/3) = f? j(x). As a result we obtain? (c+2/3) = {p(c+2/3)j ,?2(c+2/3)j , f j(x)(c+2/3)}= {p(c)j , ??2j , f? j(x)}, (B.3)where ??2j =n?i=1[p(c)i j (yi? f? j(xi))2]/n?i=1p(c)i j as in (2.16).? Step r = 3:Let g3(p j,?2j , f j) = {?2j , f j(x)}. We want to maximize S(? ,? (c)) over ?subject tog3(p j,?2j , f j) = g3(p(c+2/3)j ,?2(c+2/3)j , f(c+2/3)j ),i.e., maximize S(? ,? (c)) over p j holding ?2j = ?2(c+2/3)j = ??2j and f j(x) =f j(x)(c+2/3) = f? j(x).As a result we obtain? (c+3/3) = ? (c+1) = {p?(c+1)j ,?2(c+1)j , f j(x)(c+1)}= { p? j, ??2j , f? j(x)}, (B.4)where p? j = ?ni=1 p(c)i j /n as in (2.18).More examples and details about the ECM algorithm can be found in McLach-lan and Krishnan (2008).110Appendix CCalculation details: standarderrorsIn this appendix section we provide more details of the calculation of the matrixE?(L ?2(a12,a21)L?2(a12,a21)T |y), which is required to obtain the standard errorsfor a?12 and a?21 in Section 2.3.3.Let p?i,i?(r,s, t,u) = p(zi?1 = r,zi = s,zi??1 = t,zi? = u|y, ??). Straightforwardcalculations show that the 2?2 matrix E?(L ?2(a12,a21)L?2(a12,a21)T |y) evaluatedat convergence of the EM algorithm has off diagonal elements equal to?i6=i?[p?i,i?(1,2,2,1)a?12a?21?p?i,i?(1,2,2,2)a?12(1? a?21)?p?i,i?(1,1,2,1)(1? a?12)a?21+p?i,i?(1,1,2,2)(1? a?12)(1? a?21)],(C.1)first diagonal element equal to?ni=2 p?(i?1)1a?12(1? a?12)+?i6=i?[p?i,i?(1,2,1,2)a?212?2?p?i,i?(1,2,1,1)a?12(1? a?12)+p?i,i?(1,1,1,1)(1? a?12)2](C.2)and second diagonal element equal to?ni=2 p?(i?1)2a?21(1? a?21)+?i6=i?[p?i,i?(2,1,2,1)a?221?2?p?i,i?(2,1,2,2)a?21(1? a?21)+p?i,i?(2,2,2,2)(1? a?21)2]. (C.3)111Note that each sum for i 6= i? in (C.1), (C.2) and (C.3) can be split into smallersums according to the difference between i and i?. For instance, we can split the firstsum in (C.1) into four sums according to the following cases: i = i??1, i = i?+1,i = i??1?? with ?? 1 and i = i?+1+? with ?? 1. Therefore, the first term in(C.1) can be written as1a?12a?21[n?i?=3p(zi??2 = 1,zi??1 = 2,zi? = 1|y, ??)+n?i=3p(zi?2 = 2,zi?1 = 1,zi = 2|y, ??)+n?3??=1n?(?+1)?i=2p(zi?1 = 1,zi = 2,zi+? = 2,zi+?+1 = 1|y, ??)+n?3??=1n?(?+1)?i?=2p(zi??1 = 2,zi? = 1,zi?+? = 1,zi?+?+1 = 2|y, ??)].Let ??i j and ??i j be the forward and backward quantities, respectively, calculatedat convergence and (A?)s,t be the st-th element of the matrix corresponding tothe ?-th power of the 2? 2 matrix A of transition probabilities. The conditionalprobabilities above can be calculated as follows.p(zi?2 = r,zi?1 = s,zi = t|y, ??) =p(zi?2 = r,zi?1 = s,zi = t,y|??)?2r=1?2s=1?2t=1 p(zi?2 = r,zi?1 = s,zi = t,y|??)withp(zi?2 = r,zi?1 = s,zi = t,y|??) = ??it ? p(yi|zi = t, f?t(xi), ??2t )? a?st? p(yi?1|zi?1 = s, f?s(xi), ??2s )? a?rs? ??(i?2)randp(zi?1 = r,zi = s,zi+? = t,zi+?+1 = u|y, ??) =p(zi?1 = r,zi = s,zi+? = t,zi+?+1 = u,y|??)?2r=1?2s=1?2t=1 p(zi?1 = r,zi = s,zi+? = t,zi+?+1 = u,y|??))112withp(zi?1 = r,zi = s,zi+? = t,zi+?+1 = u,y|??) =??(i+?+1)u? p(yi+?+1|zi+?+1 = u, f?u(xi+?+1), ??2u )? a?tu ? p(yi+?|zi+? = t, f?t(xi+?), ??2t )?(A?)s,t ? p(yi|zi = s, f?s(xi), ??2s )? a?rs? ??(i?1)r.113Appendix DCalculation details: ECMalgorithmD.1 Step 1 ECM algorithmIn this section we describe in more detail Step 1 of the ECM algorithm for thereplicate case presented in Section 3.2.Step 1: Hold V and the parameters in ? fixed and maximize (3.12) + (3.13) withrespect to the f j(x)?s. So, disregarding the constant terms, we want to maximizeS? =N?k=1?alls?spk(s)(c)yTk V?1 fs(x)?12N?k=1?alls?spk(s)(c) fs(x)T V?1 fs(x)(D.1)+ P( f1, . . . , fJ,?1, . . . ,?J) (D.2)with respect to f1(x), . . . , fJ(x).We can rewrite (D.1) in terms of f j(x) as follows. The first term in (D.1) can114be written asN?k=1?alls?spk(s)(c)yTk V?1 fs(x) =N?k=1yTk V?1 ?alls?spk(s)(c) fs(x)=N?k=1yTk V?1J?j=1??????????s?s:s1= jpk(s)(c) f j(x1)...?s?s:sn= jpk(s)(c) f j(xn)?????????=N?k=1yTk V?1J?j=1D(c)k j f j(x), (D.3)where D(c)k j is an n?n diagonal matrix with ith diagonal element given by(D(c)k j )ii = p(zik = j|yk,?(c)) = ?s?s:si= jpk(s)(c). (D.4)Now we show that the second term in (D.1) can be written as?12N?k=1?alls?spk(s)(c) fs(x)T V?1 fs(x) == ?12N?k=1?alls?spk(s)(c)n?i=1n?i?=1fsi(xi)(V?1)ii? fsi? (xi?)= ?12N?k=1J?j=1J?m=1????i,i?f j(xi)(V?1)ii? ?s?s:si= j,si?=mp(c)k (s) fm(xi?)???= ?12N?k=1J?j=1J?m=1[?i,i?f j(xi) (V?1)ii? p(zik = j,zi?k = m|yk,? (c)) fm(xi?)]= ?12N?k=1J?j=1J?m=1f j(x)T Skj,m fm(x), (D.5)115where Skj,m is an n?n matrix with entries(Skj,m)ii? = (V?1)ii? p(zik = j,zi?k = m|yk,? (c)). (D.6)Note that when i = i? and j 6= m, (Skj,m)ii = 0. For i = i? and j = m (Skj, j)ii =(V?1)ii p(zik = j|yk,? (c)).Therefore, substituting the first and second terms of (D.1) with (D.3) and (D.5),respectively, we obtainS? =N?k=1yTk V?1J?j=1D(c)k j f j(x)?12N?k=1J?j=1J?m=1f j(x)T Skj,m fm(x) (D.7)+ P( f1, . . . , fJ,?1, . . . ,?J). (D.8)Two proposed methods to maximize S? are as follows. We note that both meth-ods involve finding critical points, that is, taking derivates, setting equal to zero andsolving for the critical point. We have not been able to show that the critical pointsare actually maximizers: to prove this requires showing that certain matrices arepositive definite. However, in our experience with simulations and data analysis,all of the matrices in question have been positive definite.Method 1: sequential updatingCycle through f1(x), . . . , fJ(x), that is, hold f1(x), . . . , f j?1(x), f j+1(x), . . . , fJ(x)fixed and maximize S? with respect to f j(x).To find the maximizing f j(x) for the Bayesian approach, recall thatP( f1, . . . , fJ,?1, . . . ,?J) =?12J?j=1f j(x)T A(? j)?1 f j(x). (D.9)Taking the derivative of (D.7) + (D.9) with respect to f j and setting it equal to zerowe obtain[N?k=1Skj, j +A(? j)?1]f? j(x) =N?k=1(D(c)k j V?1yk?12 ?m:m 6= jSkj,m fm(x)),116which can be written as(I+A(? j)N?k=1Skj, j)f? j(x) = (D.10)A(? j)[N?k=1(D(c)k j V?1yk?12 ?m:m6= jSkj,m fm(x))]. (D.11)Therefore, f? j(x) is the solution of the linear system above. The solution is themaximizer of S? provided the matrix[?Nk=1 Skj, j +A(? j)?1] is positive definite.We update the f j(x)(c)?s as follows. In (D.10) and (D.11), let V = V(c). Recallthat Skj,m also depends on the entries of V. Now- let f1(x)(c+1) be f?1(x) with f2(x), . . . , fJ(x) = f2(x)(c), . . . , fJ(x)(c);- for j = 2, . . . ,J?1, let f j(x)(c+1) be f? j(x) withf1(x), . . . , f j?1(x) = f1(x)(c+1), . . . , f j?1(x)(c+1)andf j+1(x), . . . , fJ(x) = f j+1(x)(c), . . . , fJ(x)(c);- let fJ(x)(c+1) be f?J(x) with f1(x), . . . , fJ?1(x) = f1(x)(c+1), . . . , fJ?1(x)(c+1).To maximize S? with respect to f j(x) recall that, for the penalized log-likelihoodapproach, we assume f j(x) = B? j, where ? j is the K?1 vector of coefficients cor-responding to f j and B is an n?K basis function matrix with entries Bi? = b?(xi).Therefore, (D.7) + (D.8) can be written asN?k=1ykV?1J?j=1Dk jB? j?12n?k=1J?j=1J?m=1?Tj BT Skj,mB?m??j=1? j?Tj R? j, (D.12)where R is a K?K matrix with entries R?? ? =?b???(x)b??? ?(x)dx, and the matricesDk j, Skj,m are as in (D.4) and (D.6), respectively. So we can now maximize (D.12)117as a function of the vector of coefficients, ?1, . . . ,?J .For each j, hold ?1, . . . ,? j?1,? j+1, . . . ,?J fixed and maximize (D.12) with re-spect to ? j. Taking the derivative of (D.12) with respect to ? j and setting it equalto zero we obtain [BT( N?k=1Skj, j)B+2? jR]?? j =BT V?1N?k=1Dk jyk?12N?k=1?m:m6= jBT Skj,mB?m.Therefore, ?? j is the solution of the linear system above. The solution is the maxi-mizer of (D.12) provided[BT(?Nk=1 Skj, j)B+2? jR]is positive definite. We updatethe ? (c)j ?s sequentially, as in the Bayesian approach, and set f j(x)(c+1) = B?(c+1)j .Method 2: simultaneous updatingMaximize S? with respect to the whole nJ-vector f = ( f1(x)T , . . . , fJ(x)T )T .For the Bayesian approach we let? A = diag(A(?1), . . . ,A(?J));? D(c)k =[D(c)k1 , ? ? ? ,D(c)kJ]with D(c)k j as in (D.4);? Sk be the nJ?nJ matrix????Sk1,1 Sk1,2 ? ? ? Sk1,J....... . ....SkJ,1 SkJ,2 ? ? ? SkJ,J????with SKj,m as in (D.6).118Thus, we can rewrite S? given by (D.7) + (D.8) in terms of f as follows.S? =N?k=1yTk V?1D(c)k f?12N?k=1fT Sk f?12fT A?1f.=N?k=1yTk V?1D(c)k f?12fT(N?k=1Sk +A?1)f.? vT f?12fT(M+A?1)f, (D.13)wherev =N?k=1D(c)Tk V?1yk (D.14)is a vector of dimension nJ andM =N?k=1Sk (D.15)is a nJ?nJ matrix. Taking the derivative of (D.13) with respect to f and setting itequal to zero we obtain(M+A?1)f = v,which is equivalent to(AM+ I)f = Av. (D.16)Therefore, f? is the solution of the linear system in (D.16). The solution is themaximizer of (D.13) provided[M+A?1]is positive definite. Now let f(c+1) be f?with V = V(c). Recall that the matrices Sk in M also depend on the entries of V.For the penalized log-likelihood approach we maximize (D.12) with respectto the whole JK-vector ? = (?T1 , . . . ,?TJ )T . We can write (D.12) in terms of ? asfollows.N?k=1yTk V?1D(c)k B?? ? 12?T B?T(N?k=1Sk)B?? ??T diag(?1R, . . . ,?JR)? , (D.17)119where B? is the nJ? JK matrix given byB? =??????B 0 ? ? ? 00 B ? ? ? 00 0. . ....0 0 ? ? ? B??????.Taking the derivative of (D.17) with respect to ? we obtain[B?T MB?+2diag(?1R, . . . ,?JR)]?? =B?T v.Therefore, ?? is the solution of the linear system above with v and M defined as in(D.14) and (D.15), respectively. The solution is the maximizer of (D.17) provided[B?T MB?+ 2diag(?1R, . . . ,?JR)]is positive definite. Let ? (c+1) be ?? with V =V(c). So we let f(c+1) = B?? (c+1).D.2 Step 2 ECM algorithmIn this section we provide more details of the calculations of Step 2 of the ECMalgorithm described in Section 3.2 for unrestricted V.Step 2: Hold the f j(x)?s and the parameters in ? fixed and maximize (3.12) inSection 3.2 with respect to V. Using standard results from matrix algebra we canrewrite (3.12) asN2log |V?1|?12N?k=1?alls?spk(s)(c)tr[V?1(yk? fs(x))(yk? fs(x))T]. (D.18)Now by straightforward application of Results 4.10 and 4.11 of Johnson and Wich-ern (2008) with B =N?k=1?alls?spk(s)(c)(yk? fs(x))(yk? fs(x))Twe obtainV? =1NN?k=1?alls?spk(s)(c)(yk? fs(x))(yk? fs(x))T.120Let V(c+1) be V? with fs(x) = fs(x)(c+1).121
Thesis/Dissertation
eng
Statistics
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Switching nonparametric regression models
Text