Review on Energy Resilience

A classification-based approach to monitoring the safety of dynamic systems

Shengtong Zhong a, Helge Langseth a,n, Thomas Dyhre Nielsen b

a Department of Computer and Information Science, The Norwegian University of Science and Technology, Trondheim, Norway b Department of Computer Science, Aalborg University, Aalborg, Denmark

a r t i c l e i n f o

Article history: Received 5 December 2012 Received in revised form 16 May 2013 Accepted 31 July 2013 Available online 8 August 2013

Keywords: Monitoring Anomaly detection Dynamic classification Hidden Markov models Approximate inference

a b s t r a c t

Monitoring a complex process often involves keeping an eye on hundreds or thousands of sensors to determine whether or not the process is stable. We have been working with dynamic data from an oil production facility in the North sea, where unstable situations should be identified as soon as possible. Motivated by this problem setting, we propose a general model for classification in dynamic domains, and exemplify its use by showing how it can be employed for activity detection. We construct our model by using well known statistical techniques as building-blocks, and evaluate each step in the model- building process empirically. Exact inference in the proposed model is intractable, so in this paper we experiment with an approximate inference scheme.

& 2013 Elsevier Ltd. All rights reserved.

1. Introduction

A typical task for the risk and reliability engineer is to monitor the status of a dynamic system, like, e.g., a chemical process. Doing so will often mean tending to a large number of sensors, each of them updating their readings on a regular basis. Real-life processes have their own natural dynamics when everything is running according to plan; “outliers” may on the other hand be seen as indications that the process is leaving its stable state, and thereby becoming more dangerous. Thus, the engineer would like to know if the system is unstable in order to ensure that the proper corrective actions are implemented as soon as the system becomes unsafe. Unfortunately, it may be difficult to measure the status of the system directly, and one will typically only have access to indirect status indicators, which need to be analyzed and com- bined in a statistical model. Formally, detecting the instantaneous status of a system described by a collection Y ¼ fY1; Y2; …; Yng of random variables is identical to classification, where an object described by a value assignment y ¼ fy1; y2; …; yng is mapped to one of a set of possible labels (or classes). The labels for an object are represented by a class variable C, and are denoted spðCÞ. We will focus on real-valued attributes in this paper, meaning that yARn. In a probabilistic framework, it is well-known that the

optimal classifier will label an object y by the class label ĉ, where

ĉ ¼ arg min c AspðCÞ

∑ c′AspðCÞ

Lðc; c′ÞPðc′jyÞ ð1Þ

and Lðc; c′Þ is the loss-function encoding the cost of mis- classification. Learning a classifier therefore amounts to estimating the probability distribution PðC ¼ cjyÞ.

The engineer may not only want to assess the instantaneous status of a system, but rather to detect if the system is about to become unstable (that is, to predict future problems). This would give a system operator the chance to implement countermeasures before anyone is exposed to an increased level of risk. Classifiers that fail to take the dynamic aspect of a process into account will not be able to make accurate predictions, and will therefore not be able to recognize a problem under development. In dynamic classification, the task is to assign a class label to an object at each time step. To support the classification, objects are characterized by a new observation at each time step as well. We use Yt ¼ fYt1; Yt2; …; Ytng to denote the random variables describing the object at time t, where yt ¼ fyt1; yt2; …; ytng is a specific value assignment to these variables. The collected observations from time t¼1 and up to time t is denoted as y1:t. The set of possible labels (or classes) for the time series at time t is represented by a class variable Ct and denoted spðCtÞ. With the observations y1:t from time step 1 to t, the optimal classifier will label y1:t by the class label ĉt at time t, where

ĉ t ¼ arg min

ct AspðCtÞ ∑

c′AspðCtÞ Lðct; c′ÞPðc′jy1:tÞ;

confer also Eq. (1).

Contents lists available at ScienceDirect

journal homepage: www.elsevier.com/locate/ress

Reliability Engineering and System Safety

0951-8320/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ress.2013.07.016

n Corresponding author. Tel.: +47 73596488; fax: +47 73594466. E-mail addresses: [email protected] (S. Zhong),

[email protected] (H. Langseth), [email protected] (T.D. Nielsen).

Reliability Engineering and System Safety 121 (2014) 61–71

In a risk and reliability setting, the desire to build efficient statistical models that are flexible yet easy to understand for domain experts has led to reduced focus on traditional frame- works like fault trees. On the other hand, the Bayesian network (BN) framework [28,17] has received increased attention from the community over the last decade [22], partly because BNs have proven to be an attractive alternative to classical reliability formalisms, see e.g., [33,18]. BNs have also been used extensively for classification [8,21,35].

The dynamic Bayesian network framework [12] supports the specification of dynamic processes, and has already found numer- ous applications in reliability engineering, see, e.g., [20,27]. A simple instance of this framework is the hidden Markov model (HMM), which has also been considered for classification purposes [15,31,6]; to this end the “hidden” node in the HMM is used as the classification node, and the attributes at time t are assumed to be independent of those at time t+1 given the class label at either of the two points in time. Further simplification can be obtained by assuming that all attributes at one time step are conditionally independent given the class label at that time step; the resulting model by [26] is known as a dynamic naïve Bayes (dNB) classifier. The dNB models can be efficiently estimated from data due to the relatively small number of parameters required to specify them.

To the best of our knowledge, there has been no systematic investigation into the properties of probabilistic classifiers and their applicability to real-life dynamic data. In this paper we will take a step in that direction by examining the underlying assump- tions of some well-known probabilistic classifiers and their natural extensions to dynamic domains. We do so by carefully linking our analysis back to a real-life dataset, and the result of this analysis is a classification model, which can be used to, e.g., help prevent unwanted events by automatically analyzing a data stream and raise an alert if the process is entering an unstable state. For the discussions to be concrete, we will tie the model development to the task of activity recognition in offshore oil drilling; this is further described in Section 2. In Section 3 we give a general overview of the dynamic classification scheme, and we also propose a specific classification model called a dynamic latent classification model (abbreviated to dLCM). Next, we look at inference and learning in dLCMs (Section 4), before reporting on their classification accuracy in Section 5. Finally, in Section 6 we conclude and give directions for future research.

2. The domain and the dataset

Offshore oil drilling is a complex process, potentially with major risks to the safety of the operators involved (see, e.g., [34]). Further, the drilling process in itself is extremely expensive, leading to a focus on cost efficient operation, including high demands wrt. the reliability of the equipment employed. This has again resulted in a plethora of data being collected – either for real-time analysis of the state of the ongoing operations or to enable investigations after an event has occurred. We will consider one such dataset from an oil production installation in the North Sea. Data, consisting of 62 variables, is captured every 5 s. The data is monitored in real time by experienced engineers, who have a number of tasks to perform ranging from understanding the situation on the platform, in order to avoid a number of either dangerous or costly situations, to optimization of the drilling operation. The variables that are collected in this dataset cover measurements taken both topside (like flow rates) and down-hole (like, for instance, the gamma rate).

The overall drilling process can be broken down into a series of activities that are performed iteratively as the depth of the well increases. Recognizing which activity is performed at a given point

in time is called activity recognition, and is the focus of the present paper. Out of the 62 attributes that are collected, domain experts have selected the following 9 attributes as the most important for activity detection: Depth Bit Measured, Depth Hole Measured, Block Position, Hookload, Weight On Bit, Revolutions Per Minute, Torque, Mud Flow In, and Standpipe Pressure.

In the Wellsite Information Transfer Specification (WITS), a total of 34 different activities with associated activity codes are defined. Each activity has its separate purpose and consists of a set of actions. Out of the 34 different drilling activities in total, only a handful are really important to recognize. The important activities in our analysis, which roughly correspond to those that constitute most of the total well drilling time, are described next

WITS2 – Drilling: The activity occurs when the well is gaining depth by crushing rock at the bottom of the hole and removing the crushed pieces (cuttings) out of the well- bore. Thus, the drill string is rotating during this activity, and mud is circulated at low speed to transport out the cuttings. The activity is interrupted by other activities, but continues until the well reaches the reservoir and production of oil may commence.

WITS3 – Connection: This activity involves changing the length of the drill-string, by either adding or removing pieces of drill-pipe.

WITS8 – Tripping in: This is the act of running the drill string into the well hole.

WITS9 – Tripping out: Tripping out means pulling the drill string out of the well bore.

It what follows, the remaining activities will collectively be grouped under the label Others.

Knowing which activity is performed at any point in time is important in several contexts: firstly, the operation of an offshore installation can be monitored by groups of experts located else- where (typically in on-shore control-rooms). These experts are shielded from the offshore-operation in that they only observe visualizations of streams of data. Important aggregations, like which activity is performed, helps them better understand the situation on-site.

Secondly, operators are consistently looking for more cost- efficient ways of drilling, and the sequencing of activities during an operation is important for hunting down potential time-sinks.

Thirdly, some undesired events can only happen during specific activities, and knowing the current activity is therefore of high importance. For instance, apparent early warnings of undesired events can be given more credence if that event can actually occur during the current activity and no weight if the event is impos- sible. From a safety perspective, this allows for a better early warning system with a lower rate of false alarms.

Finally, it is worth mentioning that activity recognition is a task that also finds applications in areas as diverse as health care [32] and video analysis [23]. In this paper we develop a model for dynamic classification and exemplify the process in the oil drilling domain, but other safety and reliability applications of the devel- oped model are readily available.

3. From static to dynamic Bayesian classifiers

In this section we develop a general framework for performing dynamic classification. The framework will be specified incremen- tally by examining its expressivity relative to the oil production data. In Section 5 we further justify the framework by setting up an empirical study using the oil production data. In the study we analyze the accuracy results for the sequence of models that are

S. Zhong et al. / Reliability Engineering and System Safety 121 (2014) 61–7162

being considered in this section and which lead to the proposed modeling framework.

3.1. Static classifiers

Standard (static) classifiers like NB [5] or TAN [8] assume that the class variable and attributes at different time points are independent given the model. This independence assumption is clearly violated in many domains and, in particular, in domains that specify a process evolving over time. To validate the indepen- dence assumptions in practice, we can for instance, compare the marginal distribution of the class variable with the conditional distribution of the class variable given its value at the previous time step. From the results using the oil production data, we see a considerable correlation between the class variable of consecutive time slices. In particular, if the system was in the drilling activity at time t�1, the probability of being in the drilling activity also at time t changes from 0.325 (static classifier) to 0.997 (dynamic classifier). The reason for this dramatic difference is that the system tends to remain in the drilling activity as soon as drilling has commenced, an effect that the static classifier is unable to represent. One way to capture this dependence is to explicitly take the dynamics of the process into account, i.e., to look at dynamic classifiers.

3.2. A simple dynamic classifier

The temporal dynamics of the class variable can be described using, e.g., a first order Markov model, where PðCtjC1:t�1Þ ¼ PðCtjCt�1Þ for all t. By combining this temporal model with the class-conditional observation model for the attributes we have the well-known hidden Markov model (HMM) [29]. The HMM model is described by a prior distribution over the class variable PðC0Þ, a conditional observation distribution PðYtjCtÞ, and transition prob- abilities for the class variable PðCtjCt�1Þ; we assume that the model is stationary, i.e., PðYtjCtÞ ¼ PðYsjCsÞ and PðCtjCt�1Þ ¼ PðCsjCs�1Þ, for all s; tZ1. HMMs have previously been used in reliability contexts. For example, Smyth [31] considers fault detection in dynamical systems, Durand and Gaudoin [6] applies HMMs for modeling the failure and debugging process of software, and Zamalieva et al. [36] use HMMs for online labeling of event sequences wrt. failure and non-failure scenarios.

With a continuous observation vector, the typical way of modeling the conditional distribution is to use a class- conditional multivariate Gaussian distribution with mean μc and covariance matrix Σc, i.e., YjfC ¼ cg � Nðμc; ΣcÞ [10]. Unfortunately, learning a full covariance matrix involves estimating a number of parameters that is quadratic in the number of attributes, which may result in over-fitting when data is scarce compared to the number of free parameters. One approach to alleviate this problem is to introduce additional independence assumptions about the domain being modeled. Specifically, by assuming that all variables are independent given the class variable, we will at each time step have a NB model defined by a diagonal covariance matrix, thus requiring only O(n) parameters to be learned, where n ¼ jYj is the number of attributes in the model. This structure corresponds to

the dNB model for dynamic domains. A graphical representation of the resulting independence assumptions can be seen in Fig. 1 in the form of a 2TBN [25].

As for the (static) NB model, the independence assumptions encoded in the dNB model are often violated in real-world settings. For example, if we consider the measured flow of drilling fluid going into the well (Mud Flow In) and the observed pressure in the well (Stand Pipe Pressure), and plot their values conditioned on the class variable (activities tripping in and tripping out), it is evident that there is a conditional correlation between the two attributes given the class, see Fig. 2; similar results are also obtained when considering other pairs of attributes.

3.3. Modeling dependence between attributes

There are several approaches to model attribute dependence. For example, Friedman et al. [9] propose an extension of the TAN model [8] to facilitate continuous domains. In the TAN framework, each attribute is allowed to have at most one parent besides the class variable. As an alternative, [21] present the latent classifica- tion model (LCM), which can be seen as combining the NB model with a factor analysis model [7].

An LCM offers a natural extension of the NB model by introdu- cing continuous latent variables Z ¼ ðZ1; …; ZkÞ as children of the class variable C and parents of all the attributes Y ¼ ðY1; …; YnÞ. The latent variables and the attributes work as a factor analyzer focusing on modeling the correlation structure among the attributes.

Following the approach by [21], we introduce latent variables to encode conditional dependencies among the attributes. Speci- fically, for each time step t we have the vector Zt ¼ ðZt1; …; ZtkÞ of

Fig. 1. Attributes are assumed to be conditionally independent given the class variable (equivalent structure for HMM) with n¼4.

Fig. 3. In each time step, the conditional dependencies between the attributes are encoded by the latent variables ðZt1; Zt2Þ.

0 0.5 1 1.5 2 2.5 3 x 107

0

0.005

0.01

0.015

0.02

0.025

0.03

0.035

Correlation between pair of attributes

Mud Flow In

S ta

nd P

ip e

P re

ss ur

e

Fig. 2. Scatter plot of the Mud Flow Int (x-axis) and the Stand Pipe Pressuret (y-axis) for the two classes in the oil production data (black “+” is the tripping in class and grey “o” is the tripping out class). The conditional correlation between the two attributes is evident.

S. Zhong et al. / Reliability Engineering and System Safety 121 (2014) 61–71 63

latent variables that appear as children of the class variable and parents of all the attributes (see Fig. 3). The latent variable Zt is assigned a multivariate Gaussian distribution conditional on the class variable and the attribute vector Y is also assumed to be a multivariate Gaussian distribution conditional on the latent vari- ables

ZtjfCt ¼ ctg � Nðμct ; Σct Þ;

YtjfZt ¼ ztg � NðLzt; ΘÞ;

where Σct and Θ are diagonal matrixes and L is the transition matrix; note that the stationarity assumption is encoded in the model.

In this model, the latent variables capture the dependencies between the attributes. They are conditionally independent given the class but marginally dependent. Furthermore, the same map- ping, L, from the latent space to the attribute space is used for all classes, and hence, the relation between the class and the attributes is conveyed by the latent variables only.

The model in Fig. 3 assumes that the attributes in different time slices are independent given the class variable. This assumption implies that the temporal dynamics is captured at the class level only. When the state specification of the class variable is coarse, then this assumption will rarely hold (obviously, the finer the granularity of the state specification of the class variable, the more appropriate this assumption will be). For the oil production data, this assumption does not hold as we can see in Fig. 4, which shows the conditional correlation of the Stand Pipe Pressure attribute in successive time slices in both tripping in and tripping out activities.

We propose to address this apparent short-coming by model- ing the dynamics of the system at the level of the latent variables, which semantically can be seen as a compact representation of the “true” state of the system. At this abstraction level, the modeling approach is related to that of Kohda and Cui [20] who propose a factorial HMM, where latent/unobserved variables are used to model temporal dynamics in safety monitoring systems. To be more specific, we encode the state specific dynamics by assuming that the latent variable vector Zt follows a linear multivariate Gaussian distribution conditioned on Zt�1:

ZtjfZt�1 ¼ zt�1; Ct ¼ ctg � NðAct zt�1; Σct Þ

where Act encodes the class conditional transition dynamics for the latent variables. A graphical representation of the model is given in Fig. 5, and will be referred to as a dynamic latent classification model (dLCM). Observe that conditional on the class

variables, the state specific model dynamics is related to a factorized Kalman filter model.

3.4. Modeling non-linear systems

One of the main assumptions in the model above is that there is a linear mapping from the latent variables to the attribute space as well as a linear mapping from the latent variables in a given time slice to the latent variables in the succeeding time slice (i.e., that the state specific dynamics are linear). When the class variable takes a fixed value, the model is equivalent to a linear dynamical system (LDS) [1], also known as a linear state-space model.

Given sufficient dimension of the latent space, an LDS can represent any complex real-world process, although the computa- tional cost can make the LDS model infeasible in practice if the target process exhibits a complicated behavior (see Appendix C). In order to reduce the computational cost while maintaining the representational power, we introduce a discrete mixture variable M for each time slice as done by [21] for static domains (see Fig. 6). A related dynamic model is the switching state-space model (SSSM), which was proposed to combine discrete and continuous dynamics [13]. The representational power and computational efficiency of the SSSM have been well demonstrated [13,2]. The model we propose differs from the SSSM not only by the introduction of a class variable, but also by our model using the discrete class variables to carry the dynamics over time (whereas this is achieved by the latent mixture node for the SSSM). SSSMs focus on modeling non-linear real-world processes with one single system state (class), whereas our model is intended to capture the non-linearity of multiple system states (classes) at the same time.

We call our model a dynamic latent classification model (dLCM), and note that we for each time slice can regard the model as a combination of a naïve Bayes model and a mixture of factor analyzers. In this case, the mixture variable follows a multinomial distribution conditioned on the class variable, and the attributes Yt follow a multivariate Gaussian distribution conditioned on the latent variables and the discrete mixture variable, i.e.

MtjfCt ¼ ctg � PðMtjCt ¼ ctÞ; YtjfZt ¼ zt; Mt ¼ mtg � NðLmt zt; Θmt Þ;

where 1rmt rjspðMÞj.

4. Learning and inference

In what follows we discuss algorithms for performing inference and learning in the proposed models.

4.1. Inference

By inference we mean to calculate the conditional distribution over some variables of interest given observations of others, for instance calculating Pðctjy1:tÞ. As the number of time-slices for which data has been collected may be quite large, we are looking for an efficient way of calculating these probability distributions.

As our model is structurally similar to the linear dynamical systems model [1], we will find inspiration for our inference scheme from the inference algorithms associated with those models [30,1]. Therefore, we will consider two phases of inference, the forward (or filtering) phase and the backward (smoothing) phase. The results of the forward and backward calculations are given in the next subsections; further details can be found in Appendix A.

0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0

0.005

0.01

0.015

0.02

0.025

0.03

0.035

Fig. 4. Scatter plot of the Stand Pipe Pressure (x-axis is value of Stand Pipe Pressure at time t and y-axis is value of Stand Pipe Pressure at time t+1) for the two classes in the oil production data (black “+” is tripping in class and gray “o” is tripping out class). The conditional correlation of this attribute over time is evident.

S. Zhong et al. / Reliability Engineering and System Safety 121 (2014) 61–7164

4.1.1. Filtering using forward recursion The goal of the filtering phase is to quantify the uncertainty

over the state of the system at time t given the observations we have up to and including time t. Primarily, the variable of interest is the class variable Ct, but simultaneously, the latent variables ðZt; MtÞ also convey information, and are therefore also of rele- vance. We thus calculate the probability distribution pðct; zt; mtjy1:tÞ during filtering, and this is done recursively in t. This means that pðct; zt; mtjy1:tÞ is found using the related results from the previous time-step, pðct�1; zt�1; mt�1jy1:t�1Þ, combined with information about how likely the new observation yt is and how the system evolves over time (see also Appendix A)

pðzt; mt; ct; jy1:tÞppðytjzt; mtÞpðmtjctÞ

� ∑ ct�1

pðctjct�1Þ Z zt�1

pðztjzt�1; ctÞ ∑ mt�1

pðzt�1; mt�1; ct�1jy1:t�1Þ dzt�1:

ð2Þ However, by examining the inference rule above we see that

exact filtering is intractable (scaling exponentially with t, see also [24,2]) because neither the class variables nor the mixture variables are observed: at time t¼1, pðz1; m1; c1jy1Þ is built up by a single Gaussian. However, at time-step t¼2, due to the summation over the class variable C1 and mixture variable M1 in Eq. (2), pðz2; m2; c2jy1:2Þ will contain a mixture of jspðCÞj � jspðMÞj Gaussians; the model contains a mixture of jspðCÞj2 � jspðMÞj2 Gaussians at t¼3, and so on. To control this explosion in computational complexity, we will resort to a Gaussian collapse method [3,2]. The Gaussian collapse guarantees that the distribution pðzt; mt; ct; jy1:tÞ is repre- sented by a single Gaussian at any time-step t; details are given in Appendix A.

4.1.2. Smoothing using the backward recursion Similar to the forward pass, the backward pass also relies on a

recursive computation. Note that where the traditional forward– backward algorithm would make use of the backward-phase to calculate pðytþ1:Tjzt; mt; ctÞ [1], we rather follow [2] and calculate pðzt; mt; ctjy1:TÞ instead. The idea is to compute pðzt; mt; ctjy1:TÞ from the corresponding result of the previous recursive step,

pðztþ1; mtþ1; ctþ1jy1:TÞ. As the calculations are a bit involved they are not presented here, but can be found in A.2.

There are two approximations involved when completing the backward recursion: firstly, similar to [2], we approximate pðztþ1jmt:tþ1; ct:tþ1; y1:TÞ by pðztþ1jmtþ1; ctþ1; y1:TÞ. Secondly, we observe that pðzt; mt; ctjy1:TÞ will become a mixture of Gaussians, and that the number of components increases exponentially in T�t (see Eq. (A.3)). We therefore resort to the same solution strategy as for the forward phase, and approximate pðzt; mt; ctjy1:TÞ by a single Gaussian at each time-point t; see Appendix A for more details.

4.2. Learning

Learning the dLCM model involves estimating the parameters of the model, the number of latent variables, and the number of mixture components. With the number of latent variables and the number of mixture components specified, parameter learning becomes simple: since the class variables are always observed during learning, the resulting model is similar to a linear state- space model for which an EM-algorithm [4] can be applied.

Following the approach of [21], we determine the number of latent variables, k, and the number of mixture components, jspðMÞj, by performing a systematic search over a selected subset of candidate structures (note that a model structure is completely specified by its number of latent variables and number of mixture components). For each candidate structure we learn the model parameters by applying the EM-algorithm, and the quality of the resulting model is then evaluated using the wrapper approach [19]. Finally, we select the model with the highest score.

Algorithm 1. Learn a dLCM classifier with the wrapper approach.

Input: A dataset D1:T . Number of wrapper folds to use: γ. Output: A dLCM classifier. 1: Partition the dataset into γ wrapper folds W1; …; Wγ so

that each Ww is a (time-)continuous part of D 1:T .

2: for possible values of k and jspðMÞj do 3: for w ¼ 1; …; γ do 4: Learn a classifier from the D1:T\Ww. 5: Calculate the accuracy on the validation dataset Ww. 6: end for 7: Score the parameter-pair (k, jspðMÞj) by the average

accuracy obtained over the wrapper folds. 8: end for 9: Select the values of k and jspðMÞj with highest accuracy. 10: return classifier learned with these parameters.

In data-rich situations one may alternatively choose to use a validation-set for model selection instead of employing the wrapper-approach. The loop starting on Line 3 in Algorithm 1 is in this case replaced by a single call to the parameter learning routine (using the full training-set D1:T), and the model selection part (Line 7) will be based on the accuracy obtained on the validation-set. This is the approach we will use in Section 5.

The parameters of the dLCM model are learned using the EM algorithm. At each iteration, the model parameters are obtained by maximizing the expected log likelihood function, resulting in the updating rules given in Appendix B.

Fig. 5. The state specific dynamics are encoded at the level of the latent variables with k¼2 and n¼4.

Fig. 6. A mixture variable M is introduced at each time slice to extend of the model, with k¼2, n¼4.

S. Zhong et al. / Reliability Engineering and System Safety 121 (2014) 61–71 65

5. Experiments results

5.1. Setup

In this section we empirically analyze the performance of the proposed dynamic latent classification model. The analysis is based on the oil drilling data described in Section 2. The data consists of sensor readings from 62 sensors captured at a sampling frequency of every 5 s. Out of the 62 sensors, domain experts deem that only 9 of the sensors are important for recognizing the activities listed in Section 2. Consequently, we focus on this reduced set of sensor readings in the experiments. The data was collected during approximately 31 h, yielding time series data containing 220 000 observation vectors, where each observation vector consists of a value for the class variable (encoding the type of activity being performed) and a configuration of the nine sensor variables.

The full dataset covers a total of five oil drilling activities, namely drilling, connection, tripping in, tripping out, and others; others is an abstract activity covering all other activities besides the four mentioned previously, c.f. Section 2. In our experiments, the data was divided into training, validation, and test datasets consisting of 90 000, 80 000, and 50 000 time slices, respectively. The training data was chosen as the initial segment of the time series, whereas the validation and test sets were chosen as the intermediate and end segments, respectively. We note that since the different segments represent different phases of the drilling, the data generation process may be different in the three seg- ments (even for the same activity), and furthermore, the fraction of time spent doing the different activities may also change between the segments. As an example, tripping into the well will take more time when the length of the well increases, and correspondingly, the marginal probability for doing the Tripping In activity changes from 12% during the initial phase (the training set) to 21% in the last phase (the test set); see Table 1 for further details. Obviously, this complicates the classification problem further.

For the classification of the five activities, we have followed the recommendations of the domain experts and used a two-step hierarchical classification process. At the first step we construct two abstract activities by merging together the three activities drilling, connection, and others into one group, and the activities tripping in and tripping out into another. For the experiments, this corresponds to constructing a new data set, where the actual activities have been replaced by the two activities drilling/connection and tripping in/out. After having classified an activity as e.g. tripping in/out, we proceed with a second step and attempt to refine the classification by reclassifying the activity as either tripping in or tripping out. Consequently, we have trained three distinct classifiers (Mtop, Mi=o, and Md=c) corresponding to the two-step classification procedure of the activities. When conducting the experiments we first deployed the Mtop model to classify an activity at time t as either drilling/connection or tripping in/out. If the activity was classified as tripping in/out we used the model Mi=o to refine the classification based on the largest consecutive sequence of data

points yt′:t classified as tripping in/out. The process for reclassifying drilling/connection is analogous, but using Md=c.

5.2. Learning procedure for the dLCM models

The dLCM learning framework consists of two components: learning the parameters for a given model structure and finding an appropriate model structure. In our learning setup these two activities are interleaved (see Section 4.2). For learning the model parameters, we employ the EM algorithm described in Section 4.2, where the termination condition was set to either 50 iterations or if the relative change in log-likelihood over two consecutive iterations falls below ϵ ¼ 10�4. In order to find the structure of the model (i.e., the number of latent variables and the state space of the mixture variable M) we adopt a greedy search strategy. More specifically, the search strategy is characterized by (i) a systematic approach for selecting values for jspðMÞj and jspðZÞj, and, given such a pair of values, (ii) learning the parameters in the model. Each candidate model is then scored by estimating its classification accuracy using a separate validation dataset.

When learning classification models for the second hierarchical step (the Mi=o and Md=c models), we first extract the relevant training data, e.g., those examples that are classified as either tripping in or tripping out are relevant for Mi=o, and subsequently use this extracted data for learning. We thereby obtain a collection of time series, but for the purpose of the experiments reported in the present paper, we have treated these time series as a single time series during learning. For future work, the learning algo- rithm will be adapted to allow for multiple time series.

5.3. Results

To analyze the classification performance of the dLCM classifier, we have compared the following list of classifiers1:

NB: The naïve Bayes classifier discussed in Section 3.1. dNB: The naïve Bayes classifier extended with dynamics on

the class variable (described in Section 3.2). dLCM1: The dLCM classifier without the mixture component, as

described in Section 3.3. The name is chosen to signify that the model is identical to a dLCM with only one state for the mixture variable (i.e., jspðMÞj ¼ 1Þ.

dLCM: The full dLCM classifier, with structure learned as described above.

LGL: The local-global learning extension of the naïve Bayes classifier towards discriminative learning [38].

J48: The J48 decision tree implementation in Weka [14] using standard parameter settings.

Note that the static classifiers (NB, LGL, J48) base their activity- classification at time t only on yt, whereas the dynamic classifiers (dNB, dLCM1, dLCM) use observations up until time t.

The results of the first hierarchical classification step for drilling/connection and tripping in/out are summarized in Table 2. We see that the classification accuracy increases with the expres- siveness of the model: the poorest results are obtained by the naïve Bayes and LGL models. Adding dynamics at the class level (dNB) improves the results, but somewhat surprisingly, only by 0.03%. dLCM1 is clearly better than dNB, hence it seems that representing the dynamics only at the class level (as the dNB does) is not sufficient to recognize the different activities; one must also capture the dynamics among the attributes to faithfully

Table 1 Empirical distribution over the activities for training-set, validation-set and test-set.

Activity Training-set (%) Validation-set (%) Test-set (%)

Drilling 42.7 42.1 42.6 Connection 32.5 21.4 23.3 Tripping in 12.2 15.7 21.3 Tripping out 10.8 20.8 12.8 Others 1.8 0.0 0.0

1 Observe that all the intermediate models considered in the development of the dLCM model are included among the straw-men.

S. Zhong et al. / Reliability Engineering and System Safety 121 (2014) 61–7166

represent the important properties of the data. Introducing the ability to model non-linear systems in the full dLCM model also contributes significantly by reducing the error from 2.86% to 1.28%. Finally, it is interesting to see that J48 also fares very well at this level of classification, being able to separate the combined activ- ities almost at the same level as dLCM1.

The detailed classification results using the dLCM classifier are shown in Fig. 7, with time on the x-axis and the class label on the y-axis. The combination drilling/connection corresponds to y¼1 and tripping in/out corresponds to y¼2. The correct classifications are shown in the topmost plot and the dLCM classifications are given in the lowermost plot. The results appear to be quite good overall, but with some patches of observations erroneously classi- fied as tripping in/out.

Next, we trained separate classifiers for recognizing the five activities drilling, connection, tripping in, tripping out, and other in a hierarchical process as outlined above. The results are given in Table 3. We can see the same tendency as we did for the aggregated activities (Table 2) apart from two important issues: firstly, we see a clearer benefit of the dynamic model at the class level, as the dNB is now much better than the NB. Secondly, J48 does not produce good results for the full hierarchical classifica- tion procedure. Where J48 could separate the aggregated activities using combinations of attribute values inside a single time-step that were impossible (NB, dNB, LGL) or costly (dLCM1, dLCM) to capture for the other models, this is clearly not sufficient for the overall classification task.

For completeness, Table 4 shows the results of each refinement-classifier. In this table, each number gives the accuracy of that sub-classifier given that the aggregated class is correct. For

instance, 60.72 for “NB – tripping in/out” means that out of all examples that are either tripping in or tripping out, NB classifies 60.72% correctly. Apparent discrepancies between Tables 3 and 4 are thus to be understood in light of misclassifications at the top- level of the hierarchical process (see also Table 2). The conclusions we can draw from Table 4 correspond well with those drawn from Table 3: firstly, the experimental results show that the dLCM classifiers (dLCM1, dLCM) achieve significantly better accuracy results than the static classifiers (NB, J48, and LGL). This verifies the need for dynamic classification models as outlined in this paper. Next, the results also show that the full dLCM is more effective than the two intermediate dynamic straw-men (dNB, dLCM1), justifying the added model complexity of the dLCM classifier. Finally, we would like to note that the accuracy results for tripping in/out are consistently lower than those for drilling/ connection. A potential contributing factor to this difference is the change in data characteristics (in particular, the distribution of the tripping in/out activities) that we observe when comparing the training, validation, and test set; see also Table 1.

Table 5 lists the models that were selected by the learning procedure. We chose the number of latent variables in the interval from k¼3 and up to k¼27 while keeping jspðMÞj ¼ 1 for dLCM1; the full dLCM models were chosen with kA½3; 27� and jspðMÞjA½1; 3�.

Finally, the detailed behavior of the dLCM classifier can be seen in Fig. 8, where again time runs along the x-axis and activity is encoded on the y-axis. For the activities, Drilling corresponds to a value of 1, Connection is given value 2, Tripping In and Tripping Out are encoded as y¼3 and y¼4 respectively, and finally Others maps to a value of 5. The top-most plot shows the correct classifications, whereas the bottom-most plot shows the results of the dLCM classifier. It is apparent that many of the mistakes made by the classifier is due to an over-reliance on the Other-activity, which is in fact not present in the test-set at all. A reason for this behavior is that since Other is made up by a combination of several activities, its data pattern is not very well defined, and therefore basically used as a “default” when no activity seems to fit the data well.

6. Conclusions

Systems for analyzing streaming data are of great importance for reliability engineering, where an obvious application area is

Table 2 Summary of classification accuracy for the first hierarchical step.

Model Classification accuracy

NB 84.77 dNB 84.80 dLCM1 97.14 dLCM 98.72 LGL 84.77 J48 96.83

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 1

2

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 1

2

D ril

lin g−

co nn

ec tio

n v

s. tr

ip pi

ng in

/o ut

D ril

lin g−

co nn

ec tio

n v

s. tr

ip pi

ng in

/o ut

Fig. 7. The dLCM classification results for the first hierarchical step: drilling/connection activity (1) vs. tripping in/out activity (2).

S. Zhong et al. / Reliability Engineering and System Safety 121 (2014) 61–71 67

process monitoring; when monitoring a process the typical task is to determine the state of the process based on streaming data consisting of current and past sensor readings. In this paper we have described a new family of models specifically designed for the analysis and classification of such streaming data. We have specified data-driven learning and inference procedures for this model class and exemplified its use by looking at online activity recognition for an oil drilling facility, where we empirically showed that our classification model significantly outperforms other standard candidate straw-men classifiers.

Our dynamic latent classification model (dLCM) generalizes the latent classification model by Langseth and Nielsen [21] to dynamic

domains, and is closely related to switching stat-space models [13]. The model class is sufficiently expressive to capture the underlying dynamics of the oil drilling data, but is unfortunately not amenable to exact inference. Instead we have employed an approximate inference scheme inspired by [2]. We have already initiated an investigation into the appropriateness of the approximate inference scheme by compar- ing our approach to traditional sampling techniques (e.g., [11]), see [37] for details, and we plan to continue this investigation as part of our future work.

The dLCM model is a general purpose classification model for dynamic data, and even though we have exemplified its use for activity detection in the oil production domain, other risk and reliability applications are also relevant. As an example, we plan to use the classification model to do event detection directly, i.e., to foresee – and therefore help the operators prevent – undesired situations. This can be a difficult problem if the events one tries to detect are very rare, as that would lead to unbalanced classification problems [16], and we thus plan to devise a semi-discriminative strategy in the spirit of [38] to examine the appropriateness of maximum-likelihood based learning in this setting. Finally, traditional probabilistic classification based on Eq. (1) requires the specification of a meaningful loss-function, and in this paper we have utilized the 0/1-loss, which is identical to maximizing the classification accuracy. In the dynamic classification setting, one could also want to encode more advanced statements that take the dynamics into consideration, like for instance “Detecting an event before a minute has past is worth $1, but detection after 90 seconds is useless”. Such statements require a richer representation than a (static) loss function, and we are planning to look into formal languages for describing them.

Acknowledgments

We would like to thank Ana M. Martínez for her participation at the offset of this work [39], and Sigve Hovda at Verdande Technology, who helped preparing and understanding the data.

Appendix A. Inference

A.1. Forward recursion: filtering

A simple decomposition of pðzt; mt; ctjy1:tÞ gives us pðzt; mt; ctjy1:tÞ ¼ pðzt; mt; ct; ytjy1:t�1Þ=pðytjy1:t�1Þ

ppðzt; mt; ct; ytjy1:t�1Þ:

Disregarding the normalisation constant pðytjy1:t�1Þ for now, we examine pðzt; mt; ct; ytjy1:t�1Þ further. Using the law of total prob- ability, we get

pðzt; mt; ct; ytjy1:t�1Þ

¼ Z zt�1

∑ mt�1

∑ ct�1

pðzt�1:t; mt�1:t; ct�1:t; ytjy1:t�1Þ dzt�1

¼ Z zt�1

∑ mt�1

∑ ct�1

pðzt�1; mt�1; ct�1jy1:t�1Þ�

pðyt; zt; mt; ctjy1:t�1; zt�1; mt�1; ct�1Þ dzt�1: ðA:1Þ

Eq. (A.1) can be simplified, but first we note that the term pðzt�1; mt�1; ct�1jy1:t�1Þ is already available form the previous time- step, and thus needs no further recalculation. Next,

pðyt; zt; mt; ctjy1:t�1; zt�1; mt�1; ct�1Þ ¼ pðytjzt�1:t; mt�1:t; ct�1:t; y1:t�1Þ � pðztjzt�1; mt�1:t; ct�1:t; y1:t�1Þ�

pðmtjzt�1; mt�1; ct�1:t; y1:t�1Þ � pðctjzt�1; mt�1; ct�1; y1:t�1Þ;

Table 3 Summary of the accuracy results for the full hier- archical classification process.

Model Accuracy

NB 58.53 dNB 61.69 dLCM1 76.29 dLCM 82.09 LGL 60.41 J48 66.75

Table 4 Accuracy results for the second-level classifications.

Model Tripping in/out Drilling/connection

NB 60.72 67.64 dNB 59.96 76.23 dLCM1 75.57 76.66 dLCM 75.57 85.47 LGL 61.07 70.18 J48 57.88 69.56

Table 5 The complexity of the learned LCM models.

Classifier dLCM1 dLCM

Mtop k¼24 k¼24, jspðMÞj ¼ 2 Mi=o k¼21 k¼21, jspðMÞj ¼ 1 Md=c k¼24 k¼21, jspðMÞj ¼ 3

0 1 2 3 4 5 1

2

3

4

5

A ct

iv ity

0 1 2 3 4 5 1

2

3

4

5

A ct

iv ity

Fig. 8. The detailed accuracy results for the full hierarchical classification process.

S. Zhong et al. / Reliability Engineering and System Safety 121 (2014) 61–7168

and now we can utilize the conditional independence statements encoded in the model. We notice (see also Fig. 6) that

Yt⫫fY1:t�1; Ct�1:t; Zt�1; Mt�1gjfZt; Mtg;

thus pðytjzt�1:t; mt�1:t; ct�1:t; y1:t�1Þ ¼ pðytjzt; mtÞ, which is simply a parameter of the model, and therefore requires no further calcula- tion. In a similar way, we can find that, pðztjzt�1; mt�1:t; ct�1:t; y1:t�1Þ ¼ pðztjzt�1; ctÞ, pðmtjzt�1; mt�1; ct�1:t; y1:t�1Þ ¼ pðmtjctÞ and pðctjzt�1; mt�1; ct�1; y1:t�1Þ ¼ pðctjct�1Þ. Eventually, we obtain pðzt; mt; ct; jy1:tÞppðytjzt; mtÞpðmtjctÞ�

∑ ct�1

pðctjct�1Þ Z zt�1

pðztjzt�1; ctÞ ∑ mt�1

pðzt�1; mt�1; ct�1jy1:t�1Þ dzt�1;

ðA:2Þ which shows how pðzt; mt; ct; jy1:tÞ can be calculated recursively in t.

However, by examining Eq. (A.2) we see that exact filtering is intractable (scaling exponentially with t, see also [24,2]) because neither the class variables nor the mixture variables are observed: At time t¼1, pðz1; m1; c1jy1Þ is built up by a single Gaussian. However, at time-step t¼2, due to the summation over the class c1 and mixture variable m1 in Eq. (A.2), pðz2; m2; c2jy1:2Þ will contain a mixture of jspðCÞj � jspðMÞj Gaussians; the model contains a mixture of jspðCÞj2 � jspðMÞj2 Gaussians at t¼3, and so on. To control this explosion in computational complexity, we will resort to Gaussian collapse [3,2]. The Gaussian collapse guarantees that the distribution pðzt; mt; ct; jy1:tÞ is represented by a single Gaussian at any time-step t.

A.2. Backward recursion: smoothing

To compute pðzt; mt; ctjy1:TÞ, we factorize it as

pðzt; mt; ctjy1:TÞ ¼ ∑ mtþ1

∑ ctþ1

Z ztþ1

pðzt:tþ1; mt:tþ1; ct:tþ1jy1:TÞ dztþ1

¼ ∑ mtþ1

∑ ctþ1

pðmtþ1; ctþ1jy1:TÞ � pðmt; ctjmtþ1; ctþ1; y1:TÞ� Z ztþ1

pðztjztþ1; mt:tþ1; ct:tþ1; y1:TÞ � pðztþ1jmt:tþ1; ct:tþ1; y1:TÞ dztþ1:

We note that Ztþ1 fMt; CtgjfY1:TMtþ1; Ctþ1g, but following [2] we assume that the influence of fMt; Ctg on Ztþ1 is “weak” compared to the influence from Y1:T , Mtþ1 and Ctþ1, and we will thus approximate pðztþ1jmt:tþ1; ct:tþ1; y1:TÞ by pðztþ1jmtþ1; ctþ1; y1:TÞ. We are left with the approximation

pðzt; mt; ctjy1:TÞ � ∑ mtþ1

∑ ctþ1

pðmt; ctjmtþ1; ctþ1; y1:TÞ� Z ztþ1

pðztþ1; mtþ1; ctþ1jy1:TÞ � pðztjztþ1; mt:tþ1; ct:tþ1; y1:TÞ dztþ1:

ðA:3Þ

Noticing that Zt⫫fYtþ1:T; Mtþ1gjfY1:t; Mt; Ct:tþ1; Ztþ1g,

pðztjztþ1; mt:tþ1; ct:tþ1; y1:TÞ ¼ pðztjztþ1; mt; ct:tþ1; y1:tÞ ppðzt; ztþ1jmt; ct:tþ1; y1:tÞ ¼ pðztþ1jzt; ctþ1Þ � pðztjmt; ct; y1:tÞ;

where pðztþ1jzt; ctþ1Þ is a parameter in the model, and pðztjmt; ct; y1:tÞ is calculated during the forward phase.

Since pðztþ1; mtþ1; ctþ1jy1:TÞ in Eq. (A.3) is known from the previous step in the backwards recursion, the last piece of the puzzle is to determine how to calculate

pðmt; ctjmtþ1; ctþ1; y1:TÞ

p Z ztþ1

pðmt; ctjztþ1; mtþ1; ctþ1; y1:TÞ � pðztþ1; mtþ1; ctþ1jy1:TÞ dztþ1:

Again, pðztþ1; mtþ1; ctþ1jy1:TÞ is known from the previous step in the backwards recursion, and

pðmt; ctjztþ1; mtþ1; ctþ1; y1:TÞ ¼ pðmt; ctjztþ1; mtþ1; ctþ1; y1:tÞ ppðztþ1jct:tþ1; mt:tþ1; y1:tÞ � pðmt; ctjctþ1; mtþ1; y1:tÞ ¼ Z zt pðztþ1jzt; ct:tþ1; mt:tþ1; y1:tÞpðztjct:tþ1; mt:tþ1; y1:tÞ dzt�

pðmt; ctjctþ1; mtþ1; y1:tÞ ¼ Z zt pðztþ1jzt; ctþ1Þpðztjct; mt; y1:tÞ dzt�

pðmtjct; y1:tÞ � pðctjy1:t; ctþ1Þ: ðA:4Þ Since pðctjy1:t; ctþ1Þppðctþ1jctÞ � pðctjy1:tÞ, it follows that Eq. (A.4) only contains terms that are readily available from the model definition or have already been calculated during the forward phase.

The calculations above show how we achieve the recursion for the backward pass using quantities that can be computed from previous results or from the forward recursion. There are two approximations in the backward recursion: firstly, we approxi- mated pðztþ1jmt:tþ1; ct:tþ1; y1:TÞ by pðztþ1jmtþ1; ctþ1; y1:TÞ. Secondly, we note that also pðzt; mt; ctjy1:TÞ in Eq. (A.3) is a mixture of Gaussians, and this time the number of components increases exponentially in T�t. We employ the same solution strategy as for the forward phase, and approximate pðzt; mt; ctjy1:TÞ by a single Gaussian at each time-point t.

Appendix B. Learning

B.1. The M-step

The parameters of dLCM model are A; L; Σ; Θ; Φ; J; K. At each iteration, these parameters can be obtained by taking the corre- sponding partial derivative of the expected log likelihood. The following are the results:

L, the linear dynamics from the latent space to the attribute spaces: Li;mt denotes the ith row of this matrix when the mixture node is mt:

∂Q ∂Li;mt

¼ Θ�1i;mt ∑ T

t ¼ 1 fPðMt ¼ mtjD1:TÞE½ZtðZtÞTjMt ¼ mt; D1:T�Li;mt

�ðyti �Φi;mt ÞPðMt ¼ mtjD1:TÞE½ZtjMt ¼ mt; D1:T�g

L̂i;mt ¼ ∑ T

t ¼ 1 fPðMt ¼ mtjD1:TÞE½ZtðZtÞTjMt ¼ mt; D1:T�g

� ��1

� ∑ T

t ¼ 1 fðyti �Φi;mt ÞPðMt ¼ mtjD1:TÞE½ZtjMt ¼ mt; D1:T�g

Φ, the offset from the latent space to the attribute spaces: Φi;mt denotes the ith row of this matrix when the mixture node is mt:

∂Q ∂Φi;mt

¼ ∑ T

t ¼ 1 f2PðMt ¼ mtjD1:TÞ yti �Φi;mt

� � �2PðMt ¼ mtjD1:TÞLTi;mt

�E½ZtjMt ¼ mt; D1:T�g

Φ̂i;mt ¼ 1

∑Tt ¼ 1fPðM t ¼ mtjD1:TÞg

∑ T

t ¼ 1 fPðMt ¼ mtjD1:TÞyti g

�

� ∑ T

t ¼ 1 fPðMt ¼ mtjD1:TÞLTi;mt E½ZtjMt ¼ mt; D1:T�g

�

Θ, the covariance matrices of the attribute spaces: Θi;mt denotes the ith element on the diagonal of the matrix when the mixture

S. Zhong et al. / Reliability Engineering and System Safety 121 (2014) 61–71 69

node is mt:

∂Q ∂Θi;mt

¼ �1 2 Θ�1i;mt ∑

T

t ¼ 1 fPðMt ¼ mtjD1:TÞg þ 1

2 Θ�2i;mt ∑

T

t ¼ 1 fðyti �Φi;mt Þ2

�PðMt ¼ mtjD1:TÞ�2ðyti �Φi;mt ÞPðMt ¼ mtjD1:TÞLTi;mt �E½ZtjMt ¼ mt; D1:T� þ PðMt ¼ mtjD1:TÞLTi;mt �E½ZtðZtÞTjMt ¼ mt; D1:T�Li;mt g

Θ̂i;mt ¼ � 1

∑Tt ¼ 1fPðM t ¼ mtjD1:TÞg

∑ T

t ¼ 1 fðyti �Φi;mt Þ2PðMt ¼ mtjD1:TÞ

�2ðyti �Φi;mt ÞPðMt ¼ mtjD1:TÞLTi;mt E½ZtjMt ¼ mt; D1:T� þPðMt ¼ mtjD1:TÞLTi;mt E½ZtðZtÞTjMt ¼ mt; D1:T�Li;mt Þg

Σ, the covariance matrices of the latent spaces: Σc denotes the diagonal matrix when the class variable is c:

∂Q ∂Σc

¼ �αc 2 Σ�1c I�

1 αc

∑ t:ct ¼ c

fE½ZtðZtÞTjD1:T�Σ�1c g

þ2 1 αc

∑ t:ct ¼ c

fAc � E½Zt�1ðZtÞTjD1:T�Σ�1c g

� 1 αc

∑ t:ct ¼ c

fAcE½Zt�1ðZt�1ÞTjD1:T�ATc Σ�1c g !

Σ̂c ¼ 1 αc

∑ t:ct ¼ c

fE½ZtðZtÞTjD1:T��2AcE½Zt�1ðZtÞTjD1:T� þ Act

�E½Zt�1ðZt�1ÞTjD1:T�ATct g

A, the linear dynamics within the latent space from one time slice to next time slice: Ac denotes the matrix when the class variable is c:

∂Q ∂Ac

¼ �2 ∑ t:ct ¼ c

Σ�1c E½ZtðZt�1ÞTjD1:T�g n

þ2 ∑ t:ct ¼ c

fΣ�1c AcE½Zt�1ðZt�1ÞTjD1:T�g

Âc ¼ ∑ t:ct ¼ c

fE½ZtðZt�1ÞTjD1:T�g � ∑ t:ct ¼ c

fE½Zt�1ðZt�1ÞTjD1:T�g !�1

J, the transition matrix of class variable is directly obtained by using frequency estimation on Pðctjct�1Þ:

K, the transition matrix of the mixture node is computed based on Pðmtjct ¼ c; D1:TÞ :

Pðmtjct ¼ c; D1:TÞ ¼ ∑t:ct ¼ cfPðM t ¼ mtjD1:TÞg

∑Tt ¼ 1fPðM t ¼ mtjD1:TÞg

:

B.2. The E-step

To complete the maximization step, the following expected terms need to be calculated given Mt ¼ mt (this also applies when Mt is taking on other values):

� PðMt ¼ mtjD1:TÞ, � E½ZtjMt ¼ mt; D1:T�, � E½ZtðZtÞTjMt ¼ mt; D1:T�, � E½Zt�1ðZt�1ÞTjD1:T�, and � E½ZtðZt�1ÞTjD1:T�.

PðMt ¼ mtjD1:TÞ and E½ZtjMt ¼ mt; D1:T� can be directly obtained from a process similar to Rauch–Tung–Striebel smoother [30]. By a

further decomposition, we can also obtain that

E½ZtðZtÞTjMt ¼ mt; D1:T� ¼ Cov½ZtjMt ¼ mt; D1:T� þE½ZtjMt ¼ mt; D1:T�fE½ZtjMt ¼ mt; D1:T�gT

E½ZtðZt�1ÞTjD1:T� ¼ Cov½ZtðZt�1ÞTjD1:T� þE½ZtjD1:T�fE½Zt�1jD1:T�gT

where Cov½ZtðZt�1ÞTjD1:T� can be calculated with the quantities obtained from Rauch–Tung–Striebel smoother process

Cov½ZtðZt�1ÞTjD1:T� ¼ Cov½ZtðZt�1ÞTjD1:t� þðE½ZtjD1:T��E½ZtjD1:t�Þ=E½ZtjD1:T� � Cov½ZtðZt�1ÞTjD1:t�

After obtaining all the expected term required from the max- imization step, the EM step for the dLCM is the complete.

Appendix C. Fitting a linear dynamical system to model any time series

In the following, we will show that a linear dynamical system can model any real-world time series, given sufficient dimension of latent continuous space.

Assume that we have observed the series v1; …; vT . Our goal is to find a linear dynamical system with transition matrix A and emission matrix B that can fit this given series. We call the latent variables at time t Ht. At each time step we require that BHt ¼ vt, and utilizing the definition of the model, we have that Ht ¼ At�1h1, giving the requirement that BAt�1h1 ¼ vt for t ¼ 1; …; T. These requirements can naively be fulfilled by letting Ht have a number of states equal to the observation sequence (i.e., jspðHtÞj ¼ T), define h1 ¼ ½1; 0; …; 0�T, and let A ¼ L1, the lower shift matrix. In this case, ht ¼ At�1h1 is a vector of only zeros, apart from a single “1” at location t. Defining B to hold the observations, B ¼ ½v1v2…vT�, the constraints above are trivially fulfilled.

We note that the representation above is extremely wasteful, and therefore not useful in practice, but it still proves the point that these models can indeed represent any observation sequence if so desired.

References

[1] Bar-Shalom Y, Li X-R. Estimation and tracking: principles techniques and software. Artech House Publishers; 1993.

[2] Barber D. Expectation correction for smoothed inference in switching linear dynamical systems. Journal of Machine Learning Research 2006;7:2515–40.

[3] Boyen X, Koller D. Approximate learning of dynamic models. In: Advances in neural information processing systems 12, 1999. p. 396–402.

[4] Dempster AP, Laird NM, Rubin DB. Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society, Series B 1977;39 (1):1–38.

[5] Duda RO, Hart PE. Pattern classification and scene analysis.New York: John Wiley & Sons; 1973.

[6] Durand J-B, Gaudoin O. Software reliability modelling and prediction with hidden Markov chains. Statistical Modelling 2005;5(1):75–93.

[7] Everitt BS. An introduction to latent variable models.London: Chapmann & Hall; 1984.

[8] Friedman N, Geiger D, Goldszmidt M. Bayesian network classifiers. Machine Learning 1997;29(2–3):131–63.

[9] Friedman N, Goldszmidt M, Lee TJ. Bayesian network classification with continuous attributes: Getting the best of both discretization and parametric fitting. In: Proceedings of the fifteenth international conference on machine learning. San Francisco, CA: Morgan Kaufmann Publishers; 1998. p. 179–87.

[10] Gama J, A linear-Bayes classifier. In: Proceedings of the 7th Ibero-American conference on AI: advances in artificial intelligence, Lecture notes in computer science, vol. 1952. Berlin, Germany: Springer-Verlag; 2000 p. 269–79.

[11] Geman S, Geman D. Stochastic relaxation, Gibbs distribution and the Bayesian restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence 1984;6:721–41.

S. Zhong et al. / Reliability Engineering and System Safety 121 (2014) 61–7170

[12] Ghahramani Z. Learning dynamic Bayesian networks. In: Adaptive processing of sequences and data structures, international summer school on neural networks, “E.R. Caianiello”-tutorial lectures. London, UK: Springer-Verlag; 1998. p. 168–97.

[13] Ghahramani Z, Hinton GE. Variational learning for switching state-space models. Neural Computation 1998;12:963–96.

[14] Hall M, Frank E, Holmes G, Pfahringer B, Reutemann P, Witten IH. The WEKA data mining software: an update, vol. 11, SIGKDD Explorations, 2009.

[15] He Y, Kundu A. 2-D shape classification using hidden Markov model. IEEE Transactions on Pattern Analysis and Machine Intelligence 1991;13(11):1172–84.

[16] Japkowicz N, Stephen S. The class imbalance problem: a systematic study. Intelligent Data Analysis 2002;6(5):429–49.

[17] Jensen FV, Nielsen TD. Bayesian networks and decision graphs. 2nd ed. Springer-Verlag; 2007.

[18] Khakzad N, Khan F, Amyotte P. Safety analysis in process facilities: comparison of fault tree and Bayesian network approaches. Reliability Engineering and System Safety 2011;96(8):925–32.

[19] Kohavi R, John GH. Wrappers for feature subset selection. Artificial Intelligence 1997;97(1–2):273–324.

[20] Kohda T, Cui W. Risk-based reconfiguration of safety monitoring system using dynamic Bayesian network. Reliability Engineering and System Safety 2007;92(12):1716–23.

[21] Langseth H, Nielsen TD. Latent classification models. Machine Learning 2005;59(3):237–65.

[22] Langseth H, Portinale L. Bayesian networks in reliability. Reliability Engineer- ing and System Safety 2007;92(1):92–108.

[23] Lavee G, Rivlin E, Rudzsky M. Understanding video events: a survey of methods for automatic interpretation of semantic occurrences in video. IEEE Transactions on Systems, Man, and Cybernetics Part C 2009;39(5):489–504.

[24] Lerner U. Hybrid Bayesian networks for reasoning about complex systems. PhD thesis, Department of Computer Science, Stanford University; 2002.

[25] Lerner U, Parr R. Inference in hybrid networks: theoretical limits and practical algorithms. In: Proceedings of the seventeenth conference on uncertainty in artificial intelligence (UAI), 2001. p. 310–8.

[26] Martínez M, Sucar LE. Learning dynamic naive Bayesian classifiers. In: Proceedings of the twentyfirst international Florida artificial intelligence research symposium conference, 2008. p. 655–9.

[27] Montani S, Portinale L, Bobbio A, Raiteri DC. Radyban: A tool for reliability analysis of dynamic fault trees through conversion into dynamic Bayesian networks. Reliability Engineering and System Safety 2008;93(7):922–32.

[28] Pearl J. Probabilistic reasoning in intelligent systems: networks of plausible inference. San Francisco, CA, USA: Morgan Kaufmann Publishers Inc.; 1988.

[29] Rabiner L. A tutorial on hidden Markov models and selected applications in speech recognition. Proceedings of the IEEE 1989;77(2):257–86.

[30] Rauch H, Tung F, Striebel CT. Maximum likelihood estimates of linear dynamic systems. Journal of the American Institute of Aeronautics and Astronautics 1965;3(8):1445–50.

[31] Smyth P. Hidden Markov models for fault detection in dynamic system. Pattern Recognition 1994;27(1):149–64.

[32] Suzuki R, Ogawa M, Otake S, Izutsu T, Tobimatsu Y, Izumi S-I, et al. Analysis of activities of daily living in elderly people living alone: single-subject feasibility study. Telemedicine Journal and E-Health 2004;10(2):260–76.

[33] Torres-Toledano JG, Sucar LE. Bayesian networks for reliability analysis of complex systems. In: Proceedings of the 6th Ibero-American conference on AI: progress in artificial intelligence, IBERAMIA '98, London, UK: Springer-Verlag; 1998. p. 195–206.

[34] Vinnem J-E. Offshore risk assessment – principles, modelling and applications of QRA studies, springer series in reliability engineering. 2nd ed.Springer- Verlag; 2007.

[35] Webb GI, Boughton JR, Wang Z. Not so naïve Bayes: aggregating one- dependence estimators. Machine Learning 2005;58(1):5–24.

[36] Zamalieva D, Yilmaz A, Aldemir T. A probabilistic model for online scenario labeling in dynamic event tree generation. Reliability Engineering and System Safety 2013:1–9.

[37] Zhong S. On approximate inference of dynamic latent classification models for oil drilling monitoring. Presented at the 9th Bayesian modeling applications (in conjunction with UAI-2012), ⟨http://www.abnms.org/uai2012-apps-work shop/papers/Zhong.pdf⟩; 2012.

[38] Zhong S, Langseth H. Local-global-learning of naive Bayesian classifier. In: Proceedings of the 2009 fourth international conference on innovative computing, information and control. Washington, DC, USA: IEEE Computer Society; 2009. p. 278–81.

[39] Zhong S, Martínez AM, Langseth H, Nielsen TD. Towards a more expressive model for dynamic classification. In: Proceedings of the twenty-third inter- national Florida artificial intelligence research society conference, 2010. p. 563–4.

S. Zhong et al. / Reliability Engineering and System Safety 121 (2014) 61–71 71