{'fraction_non_alphanumeric': 0.04856543841006347, 'fraction_numerical': 0.017777650133547776, 'mean_word_length': 3.727766463000679, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 28, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Switching nonparametric regression models for multi-curve data

We develop and apply an approach for analyzing multi-curve data where each curve is driven by a latent state process. The state at any particular point determines a smooth function, forcing the individual curve to switch from one function to another. Thus each curve follows what we call a switching nonparametric regression model. We develop an EM algorithm to estimate the model parameters. We also obtain standard errors for the parameter estimates of the state process. We consider several types of state processes: independent and identically distributed, independent but depending on a covariate and Markov. Simulation studies show the frequentist properties of our estimates. We apply our methods to a data set of a building's power usage.


Introduction
We develop and apply a method for analyzing multi-curve data where each curve follows a switching nonparametric regression model (De Souza and Heckman, 2014). That is, each curve, over its domain, switches among J unobserved states with each state determining a function. The main goal is to estimate the function corresponding to each state and the parameters of the latent process, along with some measure of accuracy.
We are motivated by the problem of calculating a building's "typical curve" of energy consumption, that is, its expected energy consumption as a function of time and other variables (e.g., weather conditions). Such knowledge allows building managers to compare the building's real-time performance to its "typical" performance which is useful, for instance, for assessing the impact of improvements on a building's energy efficiency. The data set we analyze was provided by PulseEnergy, now part of EnerNOC (www.enernoc.com).
To understand our methodological approach, compare the plots in Figure 1. Figure 1a shows hourly power usage during the months of June and July 2009 in an office building. On some days (holidays and weekends) energy usage is close to zero. We observe that on some business days the energy usage is very high, approximately twice as much as on the other days. This high power consumption occurs on warm days, when the cooling system (also called the chiller) of the building was probably on. Figure 1b presents the building daytime power usage from 9am to 4pm for 44 business days in June and July 2009. Several types of curves can be observed: one type corresponds to days when the cooling system was probably on and another type when the cooling system was off. We also observe that on some days the chiller turned on in the middle of the day. On one day the chiller went on, off and then on again.
Brown, Barrington-Leigh and Brown (2012) consider the data in Figure 1a using a very computer intensive method. They find the "typical curve" by applying a local constant kernel smoother over an extremely large number of data points, and thus, their contribution to the analysis is mainly on improving computational efficiency. They do not consider the special structure we see in Figure 1b. One shortfall of their smoothing method is that they do not model the abrupt changes in level of energy consumption, and thus their approach may oversmooth these changes. Since these changes are real features of the data, they should be modelled explicitly to better understand power usage. Our method exploits the structure of Figure 1b and differs from the approach proposed by Brown, Barrington-Leigh and Brown (2012) in two important ways: by treating each business day as a replicate; and by modelling abrupt changes in the building's power usage as arising from two functions, one function giving power usage when the chiller is off, the other function giving power usage when the chiller is on. The condition "chiller on"/"off" at any particular time is not recorded by the automatic monitoring system. Thus, it can only be inferred from the data, and so the state of the chiller forms a latent process.
De Souza and Heckman (2014) present the case where there is a single realization, a single curve switching among J functions. In that paper, we consider two models for the latent process: one where the states are independent and identically distributed, the other where the sequence of states forms a Markov chain. In addition to estimating all parameters and functions, we derive standard errors for the parameters of the latent process. In the present paper, we extend our 2014 approach into the realm of functional data analysis (Ramsay and Silverman, 2005): we consider the case when there are N curves, called replicates, with each replicate switching among J functions. This is the first work to consider the mixture of multiple functions in functional data analysis. We also consider a third type of latent state process, where the state depends on a time-varying covariate. In our application, the covariate is temperature recorded at a weather station several kilometres from the building. Preliminary data analysis indicates this dependence can be modelled via logistic regression.
Several authors have considered the single realization case from a Bayesian perspective with the smooth functions modeled as realizations of Gaussian processes. See, for instance, Tresp (2001), Rasmussen and Ghahramani (2002) and Ou and Martin (2008). The paper of Ou and Martin (2008) also contains a Bayesian analysis of the replicate case. These papers are discussed in more detail in De Souza and Heckman (2014) and contain methodology that can, in principle, lead to estimation of all J functions and the latent variable process parameters. However, unlike our work, the focus is on the estimation of just one function -the mixture, that is, a weighted average of the J functions.
In a more recent related work, Langrock et al. (2017) consider generalized additive models with a time component, where the predictor is subject to regime changes controlled by an underlying Markov process. The parameter estimates are obtained by a numerical maximum penalized likelihood approach. The authors focus on a single realization case and do not consider the replicate case.
This paper is organized as follows. In Section 2 we provide an overview of the proposed methodology. The solution to the estimation problem is described in Section 3. Some of the calculations are similar to those that appear in De Souza and Heckman (2014); these calculations are given in the Supplementary Material. In Section 4 we present the results of a simulation study. An application of the proposed methodology to a building's power usage data is presented in Section 5. Some discussion is provided in Section 6.
The computing code and the data are available as supplementary material for possible use by interested readers.

Overview of the proposed methodology
We consider a data set with N replicates where replicate k contains n observations y 1k , . . . , y nk and evaluation points x 1 , . . . , x n , which for simplicity are the same across replicates. Observation y ik depends on x i according to a hidden (unobserved) state z ik with possible state values in {1, . . . , J}. If z ik = j the expected response of y ik is f j (x i ). In this work, we assume the replicates are all generated from just one set of functions f 1 , . . . , f J , a reasonable assumption for the power usage data presented in Figure 1b and described in Section 1. We consider three types of hidden states, those that are independent and identically distributed, those that follow a Markov structure and those that are independent but with distribution depending on some covariate(s).
In principal, the xs can differ in value and number across replicates. To proceed, we need only to modify our notation and calculations, since we will model each f j as a linear combination of B-spline basis functions. However, in our Markov state process model, a conceptual challenge arises in interpreting transition probabilities when the xs vary from replicate to replicate.
Our notation is as follows.
• Observed data: x = (x 1 , . . . , x n ) T , fixed across replicates; covariate vectors v 1k , . . . , v nk ; responses y R = (y T 1 , . . . , y T N ) T , where y k = (y 1k , . . . , y nk ) T . • Hidden states: We assume that z 1 , . . . , z N , are independent. Given the hidden states z k , y k = f z k (x) + k , where 1 , . . . , N , are independent and k has a multivariate normal distribution with mean equal to the 0-vector and covariance matrix V, possibly depending on z k . That is, k ∼ M V N (0, V). Therefore, y 1 , . . . , y N are independent and, given the hidden states z k , y k ∼ M V N (f z k (x), V). Our model can be considered a functional data model. In usual functional data modeling, when there is no switching regression, the observations from the kth replicate, y 1k , . . . , y nk , are generated from a single realization of a stochastic process (see, for instance, James, Hastie andSugar, 2000, andYao, Müller andWang, 2005). In our case, for the kth replicate, the observations arise from J stochastic process realizations, f 1k , . . . , f Jk , one for each possible state. The distribution of the kth replicate of the jth stochastic process satisfies E(f jk (x)) = f j (x) with the covariance between f jk (x) and f jk (x * ) generating the covariance matrix V. Thus, V induces a dependence among the observations of the kth realization.
We let γ be the set containing f 1 (x), . . . , f J (x) and the parameters in V. We assume that the distribution of each z k is governed by a parameter vector α. Section 2.1 presents our different choices of V and α.
Our goal is to estimate θ ≡ {α, γ}, along with standard errors or some measure of accuracy for the parameters in α. Similar to De Souza and Heckman (2014) we obtain the parameter estimates by maximizing where p(y k |θ) is the likelihood function based on the observed data from the kth replicate and P is a roughness penalty on the f j s. The exact form of P (f 1 , . . . , f J , λ 1 , . . . , λ J ) is chosen by the user. For our work, we set since the integrated squared second derivative of a function is a common form of roughness penalty (Wahba, 1990). The λ j s are the smoothing parameters, governing the weight of the penalty term. As in De Souza and Heckman (2014) one could also take a Bayesian approach by maximizing (1) with P arising from placing a Gaussian process prior on the f j s. The form of log p(y k |θ) is very complicated, since it involves the distribution of the latent states z k . Therefore, we apply an Expectation-Maximization (EM) algorithm (Dempster, Laird and Rubin, 1977) to maximize (1). We can show (see, for instance, Cappé, Moulines andRydén, 2005 andMcLachlan andKrishnan, 2008) that our EM algorithm generates a sequence of estimates, θ (c) , c ≥ 1, satisfying l(θ (c+1) ) ≥ l(θ (c) ). One could also perform a numerical likelihood maximization as described in MacDonald (2014) and Zucchini, MacDonald and Langrock (2016).
As in the single realization case presented in De Souza and Heckman (2014) we use again the results of Louis (1982) to obtain standard errors for the estimates of the parameters of the latent state process. When the hidden states, z 1k , . . . , z nk , are independent and identically distributed (iid ) we consider J ≥ 2 possible state values. For z 1k , . . . , z nk following a Markov structure we restrict the possible number of states to J = 2. We also obtain standard errors for the intercept and slope parameters for the case where J = 2 and z 1k , . . . , z nk are independent with the distribution of z ik depending on only one covariate. See Section 2 of the Supplementary Material for more details.

Choices of V and α
We consider five models for the covariance of the residual error, V: unrestricted, diagonal with either V = σ 2 I or with entry (i, i) depending on the latent state, and two generated from a "random intercept" covariance structure: a homogeneous random intercept model and a non-homogeneous random intercept model with variability of the intercept depending on the value of the latent state. We usually use V z to denote models where the variability depends on the latent state. However, sometimes we omit the subscript z when referring to a general V. The unknown parameters in V are clear for our first two models. For the third model, the parameters in V z are σ 2 1 , . . . , σ 2 J . For V to follow a homogeneous random intercept model, let where δ k and e ik are independent for all i = 1, . . . , n and k = 1, . . . , N , and the δ k s are iid N (0, τ 2 ) and the e ik s are iid N (0, σ 2 ). Then V will depend on only two parameters and can be written as where I is an n × n identity matrix, 1 is an n-vector of ones and d = τ 2 /σ 2 . Our data analysis (Section 5) requires the more complex covariance structure of a non-homogenous random intercept model, where the variance of the random intercept depends on the state. We define this model for the simple case, where there are J = 2 states. We assume that y ik = f z ik (x i ) + z ik , ik , where 1,ik = δ k + e ik when z ik = 1 and 2,ik = δ k + ϑ k + e ik when z ik = 2. In addition, δ k , ϑ k , and e ik are independent for i = 1, . . . , n, and k = 1, . . . , N , with δ k s iid N (0, τ 2 1 ), ϑ k s iid N (0, τ 2 2 ) and e ik s iid N (0, σ 2 ). Therefore, the covariance matrix for the non-homogeneous random intercept model is given by where d j = τ 2 j /σ 2 and 1 z k is an n-vector with ith entry I(z ik = 2). In our model α is the vector containing the parameters governing the distribution of the hidden states. If z 1k , . . . , z nk are iid, then α is of length J with jth component equal to p(z ik = j|α) ≡ p j . If z 1k , . . . , z nk follow a Markov structure, that is, if p(z ik |z (i−1)k , . . . , z 1k , α) = p(z ik |z (i−1)k , α), i = 2, . . . n, then the parameter vector α consists of the initial probabilities, π j = p(z ik = 1|α), and the transition probabilities, a lj = p(z ik = j|z (i−1)k = l, α), j, l = 1, . . . , J. Note that the transition probabilities do not depend on i or k.
In the case where z 1k , . . . , z nk are independent, with the distribution of z ik depending on a vector of covariates In this case α contains all the regression coefficient vectors β 2 , . . . , β J .

Parameter estimation
Here we present the proposed EM algorithm to obtain the estimates of the parameters in θ. In the M-step, we take the same approach as De Souza and Heckman (2014) and model each f j as a linear combination of K known cubic B-spline basis functions, so that The smoothing parameters, λ 1 , . . . , λ J , can be chosen by a data driven method or subjectively by visual inspection. In Section 3.3, we propose and justify a leave-one-curve-out cross-validation criterion to find the optimal λ j s for the case when V is diagonal and use this method in our application. In our application, when V is based on the nonhomogeneous random intercept model, we choose the smoothing parameters via a "brute force" leave-one-curve-out method, assuming that λ 1 = λ 2 = λ. We use a weighted cross-validation criterion where the weights reflect the uncertainty of the hidden states (see Section 3 of the Supplementary Material). In all of our simulation studies, to reduce computation time, we pre-choose the λ j s by examining a few data sets and visually ensuring that the estimated functions have the same smoothness and shape as the true curves.
Let p(y R , z R |θ) be the joint distribution of the observed and latent data given θ, also called the complete data distribution. The application of the EM algorithm to the replicate case is similar to that of the one realization case considered in De Souza and Heckman (2014), which is based on writing In what follows we present a summary of the E and M steps. See Section 1 of the Supplementary Material for details.
In the E-step we calculate In the M-step, we want to find θ (c+1) that maximizes S(θ, θ (c) ) ≡ Q(θ, θ (c) ) + P (f 1 , . . . , f J , λ 1 , . . . , λ J ) with respect to θ, or at least satisfies S(θ (c+1) , θ (c) ) ≥ S(θ (c) , θ (c) ). Let s be an n-vector of possible hidden states, i.e., each entry of s is in {1, 2, . . . , J}, and let whose value depends on the model assumed for the hidden states. Therefore, disregarding the constant terms, we maximize , and the parameters in V}. Note that θ (c) is fixed and thus so are the p k (s) (c) s. We also consider the smoothing parameters, λ 1 , . . . , λ J , to be fixed. We apply a natural extension of the EM approach, the Expectation-Conditional Maximization (ECM) algorithm (Meng and Rubin, 1993), to obtain the parameter updates θ (c+1) .

M-step via an ECM algorithm
The steps of the ECM algorithm are summarized as follows.
1. Hold V and the parameters in α fixed and maximize S * with respect to . That is, maximize the sum of (4) and (5). 2. Hold f 1 (x), . . . , f J (x) and the parameters in α fixed and maximize (4) with respect to the parameters in V, obtaining V (c+1) . 3. Hold f 1 (x), . . . , f J (x) and V fixed and maximize (6) with respect to the parameters in α, obtaining α (c+1) .
The results for Steps 1, 2 and 3 are given below. Details can be found in Section 1.2 of the Supplementary Material. Step We propose a method to update the f j (x)s that is straightforward and yields To do this, let 1 j,s be the n-vector with ith element equal to 1 if s i = j, 0 else. Let I s be the n by nJ matrix, I s = [ diag(1 1,s ) | · · · | diag(1 J,s )]. Then we easily see that f s (x) = I s f . Recall that f j (x) = Bφ j . Let B * be the nJ × KJ block diagonal matrix with each block equal to B and let φ be the JK-vector φ Combining these calculations we see that, to find the f j s that maximize the sum of (4) and (5), we must maximize, as a function of φ, This expression is quadratic in φ and is easily maximized in closed form. Let φ (c+1) be this maximizing φ when we set V = V (c) . So we let f (c+1) = B * φ (c+1) .
For a model with V s ≡ V, with no dependence on the state vector s and no restrictions on the form of V, we show that V (c+1) is . Note that if the values of z k were non-random and known, then p k (s) (c) is a delta function and so V is similar to the sample covariance matrix of the y k s. When V s ≡ V follows a homogeneous random intercept model we update the parameter estimates of the restricted V in (2) as follows. Let σ 2 (c+1) bê with σ 2 replaced by σ 2 (c+1) . Therefore τ 2 (c+1) = d (c+1) × σ 2 (c+1) .

The maximization in
Step 2 when V s follows the non-homogeneous random intercept model is given in Section 1.2 of the Supplementary Material, for the case of J = 2 states. The ECM algorithm for diagonal V s is given in Section 3.2.
We maximize (6) with respect to the parameters in α, with the calculations depending on the proposed model for the hidden states.
When z 1k , . . . , z nk are iid with p j = p(z ik = j|α) we obtain For Markov z ik s, where the vector α is composed of transition probabilities a lj and initial probabilities π j , we obtain where n s,lj is the number of transitions in s from state l to state j, that is, n s,lj = n i=2 I{s i−1 = l, s i = j}. When z 1k , . . . , z nk are independent with the distribution of z ik depending on some covariate(s), α contains the regression coefficients from our logistic regression model for p(z ik = j|v ik , α) ≡ p j (v ik , α). In this case, (6) which must be maximized numerically, for instance, via a Newton-Raphson method.

ECM algorithm when V is diagonal
Recall that we consider two cases of V diagonal, one with V = σ 2 I and one with V = V z k = diag(σ 2 z 1k , . . . , σ 2 z nk ). We could use the notation and steps of Section 3.1, modifying Step 2 for these types of V. However, it is much easier to re-derive all three steps using the independence of the components of y k in order to rewrite L 1 (γ), and thus S(θ, θ (c) ), in simpler form. We will see below that, instead of the p k (s) (c) s in (4) and (6), we require the simpler The forms of p ik (j) (c) are given in Section 1.3 of the Supplementary Material.
Here, we carry out the calculations of the ECM algorithm for the case that V z k = diag(σ 2 z 1k , . . . , σ 2 z nk ), as they can be easily modified for the case that V = σ 2 I: simply replace σ 2 j by σ 2 .
1. Updating the f j (x)s. Holding the σ 2 j s and the parameters in α fixed and maximizing the sum of (8) and (9) with respect to f j (x) we obtain We let f j (x) (c+1) bef j (x) with σ 2 j in W kj replaced by σ 2(c) j . 2. Updating the σ 2 j s. Holding the f j (x)s and α fixed and maximizing the sum of (7) and (8) with respect to σ 2 j we get . 3. Updating α. Hold the f j (x)s and the σ 2 j s fixed and maximize (10) with respect to the parameters in α. For iid z ik s we obtain For Markov z ik s, we have For z ik s independent with distribution of z ik depending on some covariates we need numerical optimization methods, such as Newton-Raphson, to obtain the coefficient estimates from our logistic regression model for p(z ik = j|v ik , α). So, for example, if there are J = 2 states and the covariate vector is v ik = (1, v ik ) T , we apply a numerical method to obtain β 20 and β 21 that maximize

Choice of the smoothing parameters when V is diagonal
In principal, we can always compute the smoothing parameters by "leave-onecurve-out" cross-validation. However, for many models, this can be computationally intensive. Fortunately, in the models with V = σ 2 I or V z k = diag(σ 2 z 1k , . . . , σ 2 z nk ), we can shorten calculations by using Theorem 1 below. In this section, we describe our iterative cross-validation procedure, implemented for our data analysis in Section 5.1 for V z k = diag(σ 2 z 1k , . . . , σ 2 z nk ). The steps for V = σ 2 I are the same except withσ 2 replacing theσ 2 j s . In our data analysis we set the initial values, the λ (0) j s, to those that worked well when tested on the data set. We update the λ j s as follows.
1. At iteration i, with λ j = λ (i) j , j = 1, . . . , J, use the ECM algorithm of Section 3.2 to find thep ik (j)s, theσ 2 j s and thef j s. 2. Discard thef j s from Step 1. 3. Let W kj be W kj as defined in (11) but with theσ 2 j s andp ik (j)s replacing the σ 2 j s and p ik (j)s. Treat theσ 2 j s and thep ik (j)s and thus the W kj s as fixed. 4. For j = 1, . . . , J, over a grid of possible λ values, set λ (i+1) j as the value of λ that minimizes the following leave-one-replicate-out cross-validation criterion: wheref (−k) j λ is the function that maximizes 5. Repeat steps 1-4 with λ j = λ (i+1) j , j = 1, . . . , J, until convergence.
We use the final values of the λ j s to obtain all of the parameter estimates from the ECM algorithm as in Section 3.2. Finding λ that minimizes (13) is computationally intensive. Fortunately, we have the following theorem.
Theorem 1 Letf 1λ , . . . ,f Jλ be the maximizers of the sum of (8) and (9), with W kj replaced by W kj . Let H kj be as in (12), but with W kj replaced by W kj . Suppose that I − H kj is invertible and W kj is positive definite, j = 1, . . . , J. Then The proof follows directly from Lemma 2 in the Appendix, which holds in a slightly more general setting.

Simulation study
We carry out a simulation study under three different designs considering that the hidden states, the z ik s, can take values 1 or 2. For each design 300 independent data sets are generated, each with N = 100 replicates. In design 1, z 1k , . . . , z nk are iid and V follows the homogeneous random intercept model as in (2). In design 2, z 1k , . . . , z nk follow a Markov structure and V also follows the homogeneous random intercept model. In design 3, z 1k , . . . , z nk are independent with the distribution of z ik depending on a univariate covariate, v ik . In this third design, we take V = σ 2 I. To study all three designs we use the same vector of evaluation points x and the same true functions f 1 and f 2 . The vector x = (x 1 , . . . , x n ) T consists of n = 10 equally spaced points, 1, 12, 23, . . . , 89, 100. The true function f 2 is the same we used in the simulation study presented in De Souza and Heckman (2014). The true function f 1 is simply f 2 − 0.1. In the third study, for each simulated data set, we generate v ik , k = 1, . . . , n, i = 1, . . . , N . Figure 2 contains example data sets generated from each of the three designs.
For Designs 1 and 2 we generate each simulated data set as follows.
2. Generate the y ik s according to the homogeneous random intercept model of Section 2.1 with τ 2 = 10 −4 and σ 2 = 10 −5 . 3. Repeat steps 1 and 2 N = 100 times to obtain a data set of 100 replicates.
For Design 3 we generate each simulated data set as follows.
We analyze the simulated data under each design using the proposed EM algorithm. We set initial parameter values to the true parameter values to speed up computation. We did try initial values that were different than the true parameter values and the EM algorithm also converged, but it took longer than when starting from the truth, as expected.
The values of λ 1 and λ 2 are fixed and equal to 10 −4 in the study of all designs. We choose this value by examining a few simulated data sets and a range of lambda values. We find that the results of these preliminary analyses are not sensitive to the choice of smoothing parameter over a wide range of lambda values.

Results
The three plots in Figure 2 show the fitted valuesf 1 (x) andf 2 (x) (dashed curves) for a data set generated from each simulation design.
We assess the quality of the estimated functions via the pointwise empirical mean squared error (EMSE) as in De Souza and Heckman (2014). For all designŝ f 1 (x) andf 2 (x) produce very small values of EMSE (< 2×10 −6 ). However, when generating data according to Design 3, the EMSE values forf 1 (x) are larger than for Designs 1 and 2.
We observe that in all cases we are slightly underestimating the values of the variance parameters. This may be due to the challenges of correctly adjusting the degrees of freedom in the estimates, in order to account for the estimation of the f j s. Recall that, in Designs 1 and 2, the error variance satisfies 10 5 σ 2 = 1 and in Design 3, 10 5 σ 2 = 5. The averages of our estimates of 10 5 σ 2 (with standard errors) under Designs 1, 2 and 3 are, respectively, 0.978 (0.046), 0.977 (0.045) and 4.919 (0.238). In Designs 1 and 2, we have an additional variance parameter, namely, the variance of the random effect intercept, with 10 4 τ 2 = 1. In these cases, the averages of 10 4 τ 2 are equal to 0.977 with standard deviations equal to 0.152. Table 1 contains the mean and the standard deviation of the estimates of the parameters of the latent process under each simulation design, along with the averages of our proposed standard errors (SEs). Note that the standard deviations of the estimates are close to the values of the means of the proposed SEs, as desired. Table 1 also shows the empirical coverage percentages of a 90% and a 95% confidence interval. We consider confidence intervals of the form "mean of the parameter estimates ±z α/2 × proposed SE", where z α/2 is the α/2 quantile of a standard normal distribution with α = 0.1 and 0.05. The empirical coverage percentages under all three simulation designs are very close to the true level of the corresponding confidence interval.

Analysis of the power usage data
The data shown in Figure 1b consist of daytime hourly power usage of a building from 9am to 4pm (n = 8 observations in a day) on N = 44 business days in June and July 2009. For the same days and hours we also have available the temperature at a local weather station. We apply our proposed methodology to these data treating each day as a replicate and modelling power usage as arising from J = 2 functions, one function giving power usage when the chiller is off (j = 1), and the other function giving power usage when the chiller is on (j = 2). In Section 5.1 we present the results assuming the covariance matrix V is diagonal and in Section 5.2 we present the results when we assume V is generated by the non-homogeneous random intercept model as in (3).

Results: diagonal V
In this section we consider two models for V: V = σ 2 I and V = V z k = diag(σ 2 z 1k , . . . , σ 2 z nk ). We use the ECM algorithm described in Section 3.2 to estimate the model parameters considering iid z ik s, Markov z ik s and z ik s that are independent with distribution depending on temperature. The smoothing parameters, the λ j s, are chosen by cross-validation as described in detail in Section 3.3.
Figures 3a and 3b present the fitted functions for iid hidden states z ik s when we assume V = σ 2 I and V z k , respectively. We can observe that the fitted curves are very similar in the two figures. The estimated curve giving power usage when the chiller is on, obtained assuming V = σ 2 I, is slightly smoother than the one obtained assuming V z k . Table 2 presents the parameter estimates and chosen λ j s. We can see that the estimates of p j = p(z ik = j) from the two models for V agree within the reported standard errors. We also observe in the lower half of the Table that the estimated variance when the chiller is on is much higher than when the chiller is off. Figures 3c and 3d present the fitted curves for Markov z ik s when we assume V = σ 2 I and V z k , respectively. As in the iid case, the fitted curve giving power usage when the chiller is on obtained assuming V = σ 2 I is slightly smoother than the one obtained assuming V z k . Table 3 provides information on the estimated model parameters and the chosen smoothing parameters. As in the iid case, the estimated variance when the chiller is on is much higher than when the chiller is off. We observe that the estimates of a 21 , the transition probability from "chiller on" to "chiller off", are very small or equal to zero. Any estimate of a 21 is expected to be small, as there is only one replicate in the data set where we observe this transition. The estimate of zero is reasonable when we assume different variances;â 21 is zero because the transition happens gradually, which our model does not allow, and the method incorrectly classifies all observations as coming from the condition "chiller on", failing to detect the transition. This replicate is the green curve in Figure 3d. Figure 4a presents the fitted curves when we assume the z ik s are independent with the distribution of z ik depending on temperature via the following logistic regression model: log p(chiller on | temperature) p(chiller off | temperature) = β 0 + β 1 temperature. Table 4 shows the corresponding estimated model parameters assuming V z k along with the chosen smoothing parameters, the λ j s. We observe in Table 4 that the standard error forβ 1 is very small and by considering a confidence interval of the formβ 1 ± 1.96 × SE(β 1 ) we conclude that the coefficient β 1 is statistically significant.

Results: correlated observations generated by the non-homogeneous random intercept model
In the analyses of Section 5.1, we see that the variability in energy consumption when the chiller is on is higher than when the chiller is off. Thus, models such as V = σ 2 I or V following the homogeneous random intercept model may not be appropriate. Therefore, to model this heterogeneity in variance and the correlation between observations, we fit the proposed switching nonparametric regression model to the power usage data assuming the covariance matrix V is generated by the non-homogeneous random intercept model as in (3). We use the ECM algorithm described in Section 3 and in Section 1.2 of the Supplementary Material to obtain the parameter estimates. We conduct the analysis assuming the hidden states z ik s are iid. We assume that λ 1 = λ 2 = λ and choose the smoothing parameters via a "brute force" leave-one-curve-out method over a grid of possible values of λ (see Table 1 of Supplementary Material). Table 5 presents the parameter estimates. We observe that the estimates of p 1 and p 2 in Table 5 agree within the reported standard errors with the estimates obtained in Table 2 where we assume the observations are uncorrelated. Figure  4b shows the corresponding fitted curves. We can observe that the fitted function corresponding to the condition "chiller on" is lower than that in Figures 3a to 4a. The non-homogeneous random intercept model appears to "explain" days of high power usage by a larger variability of the "chiller on" random intercept. Thus the replicates with very high power usage have less of an impact on the final fitted "chiller on" curve.

Discussion
We have introduced a method for the analysis of data arising from random samples of a process with a complex structure. The structure depends on a latent state process where each state corresponds to a true smooth regression function. The estimation techniques and standard error calculations were developed for several specific cases of state processes and error covariances. We have considered restrictive covariance structures, save for the case where V is completely unrestricted. While the covariance models we consider may not capture all of the dependencies in a data set, our techniques and ideas should carry over to more complex time series modelling of the error process. For instance, we can model more complicated covariance structures via random regression approaches, such as with B-spline basis functions or with lines that have random slopes in addition to random intercepts. Similarly, we can use our methods to consider more complex models for the latent process, such as a Markov model with covariatedependent transition probabilities. Further useful extensions might incorporate a dependence among replicates; for instance, in studying energy consumption of several buildings, one would want to incorporate a random "building" effect.
Theorem 1 is based on the following lemmas, which frame the problem for fixed j and fixed λ (so these are dropped in notation) and with general matrices W r , r = 1, . . . , N . Lemma 1 holds for general penalties, while Lemma 2 places further restrictions, restrictions that hold in our setting. Throughout, we assume that all maximizers exist.

Proof of Lemma 2
Note thatf ( * k) , as defined in Lemma 1, is the maximizer of S with y k replaced byf (−k) . By the assumption of the form of the maximizer of S,f ( * k) (x) can be written asf From Lemma 1 we knowf (−k) (x) =f ( * k) (x). Thus, Now subtracting y k from both sides of this equation, we obtain  Table 2 Data analysis results for iid z ik s for V = σ 2 I and Vz k = diag(σ 2 z 1k , . . . , σ 2 z nk ), with corresponding fitted curves in Figures 3a and 3b, respectively. Note that for V = σ 2 I the estimateσ 2 does not depend on j and, therefore, its value appears in the middle row.