Dynamic hybrid model for comprehensive risk assessment: a case study of train derailment due to coupler failure

Comprehensive risk assessment plays a significant role in railway rolling stock safety planning to prevent accidents, including rail derailment and collision. Several methods of evaluating individual sources of railway system risk, ranging from human factors to inherent system failure and environmental hazards, exist in the literature. However, the lack of a hybrid technique to integrate these multiple sources of risk holistically, including their interdependent effects, as a single framework for robust, accurate, and comprehensive risk assessment can limit risk perception and risk mitigation actions. This report proposes a dynamic hybrid model (DHM) that incorporates the Bayesian convolutional factorization and elimination method as a compound aggregation of frequency and severity distributions. The DHM validates predicted risk using Bayesian expectation–maximization machine learning with evidenced-based propagation from expert knowledge and learned data. It also incorporates sensitivity analysis to improve the predicted risk further by prioritizing the hazards with the maximum impact on the estimated risk due to organization resource constraints. A railway case study in the UK revealed that risk prediction using the DHM provided a holistic view of the risk. The results showed that the quantitative risk prediction using the DHM was significantly more robust, accurate, and holistic than that of the conventional risk-assessment method based on the inherent failure rate. This research will facilitate the comprehensive development of risk-mitigation strategies, such as improvements in staff training and wiring insulation, to decrease the likelihood of train derailment caused by semi-permanent coupler failure.


I. INTRODUCTION
Despite the technological advancements and developments in autonomous machines, uncertainties caused by internal and external factors still influence the risk profiles of critical industrial physical assets, such as rolling stock (RS). Recently, RS systems and their associated components have undergone incredible transformations to accommodate the increasing need for speed and passenger comfort [1]. Although such transformations have significantly enhanced customer satisfaction, they have also increased the overall system complexity, uncertainty, and risk profiles of most critical components [2]- [5]. One of these critical components is the semi-permanent coupler, which is designed to guarantee permanent mechanical and pneumatic linkage between the individual cars in a typical RS system. Semi-permanent couplers are often fitted using metal components with vulcanized rubber to facilitate relative movements between the connected cars. Further, the coupler functions as a mechanism for resisting vertical, horizontal, and rotational movements. Fig. 1 depicts a typical three-car electric multiple unit (EMU) with an exploded view of a semi-permanent coupler and its associated major components for better visualization. Typical three-car EMU with exploded view of a semi-permanent coupler and its associated major components [6], [7].
The Federal Railroad Administration reportable train accident data for Class I railroads over a period of nine years show that couplers contributed to a staggering 133 derailments involving 771 cars in the United States alone, making it a significant cause of RS accidents [8]. Data on railroad accidents from the UK Railway Archives database also identifies couplers as a major cause of derailment [9]. Owing to these undesirable statistics, adopting a holistic view for the risk management of complex systems is imperative, rather than using the currently dominant individualized approaches [10]- [14]. Traditionally, risk estimation has been used to assist engineers and operators in their decision-making regarding acceptable and unacceptable risks. An acceptable risk is one that can be tolerated by the final user given the constraints related to the organization, limitations of system functions, and governmental factors, such as cost and regulations [15]- [17]. As uncertainty is part of the ecosystem associated with dynamic complex systems (e.g., RS), risk will always exist. Furthermore, significant sources of uncertainty exist in any engineering design, including noise, errors, bias in the sample data, and inaccuracies in the risk estimation model [18], [19]. Overall risk estimation for complex systems is mostly undertaken independently for internal and external uncertainties, without considering the effects of interdependencies that exist among them [15], [16], [20]- [30]. Due to the independent application of these techniques to risk estimation, the combined effects of internal and external uncertainties are often overlooked [31]. For example, despite the low intrinsic failure rates of signaling and RS subsystems, Health and Safety Executive (HSE) rail accident reports in the UK highlight the critical roles of human factors in accidents that have occurred at Clapham Junction [32], Southall [33], and Ladbroke Grove [34]. In addition, analysis of the causes of train derailment in Victoria, Australia, between 2007 and 2013 revealed that infrastructure-related defects accounted for 56% of derailments, followed by operation-related issues at 20% and RS at 12%, with the remaining 12% related to other external causes [35]. Thus, derailment is caused by both internal and external factors. Furthermore, the results indicate the importance of considering risk factors comprehensively to implement prioritized accident mitigation factors sufficiently and robustly.
Given the ubiquitous nature of risk, realistic risk modelling of any engineering system includes estimating system risks in dynamic operating conditions and considering external risk factors, such as human errors and biases in environmental conditions [36]. Over the last decade, human reliability analysis (HRA) [37]- [42], human factors analysis (HFA) [19], [43], [44], and human and organizational factors (HOF) analysis [18], [45] have gained attention and are often performed in isolation from the assessment of technical risks arising from the system design [36]. HRA detects the deviations in computational human-error probabilities using steady-state risk-assessment techniques such as fault trees. However, the absolute risk determination in HRA tends to contradict the actual source of human failure and the available HRA data lack accuracy [36], [46]. Thus far, several technical and human error risk assessments have been conducted independently of each other or in isolation [24], [25], [47]- [50]. Moreover, technical risk analysis dominates most industrial risk assessments and tends to rely mainly on qualitative risk matrix tables that are based on individual operational techniques. Consequently, the analysis sometimes emphasizes the subjective opinions of the analysts, organizational risk tolerance, and operational consequences of failures [17], [51]. Although Kariuki and Löwe [31] proposed a means of integrating human factors into process hazard analysis, the approach focused solely on qualitative techniques without significant emphasis on integrated quantitative assessment of inherent machine errors, human errors, and other external factors during the overall risk estimation.
To enhance the representativeness of overall risk estimation approaches, especially when handling complex systems, researchers have proposed various risk aggregation methods for comprehensive safety management of safety-critical systems [52], [53]. However, owing to random interactions that exist among independent variables, the current risk aggregation models have compatibility issues and inaccuracies considering multiple data types for frequency factors [54]. Therefore, the existing methods of risk aggregation are inadequate for dependence correlation analysis as they struggle to integrate perspectives, learned experiences, and knowledge of subject matter regarding the relationships between causal influences and observed data [55]- [57]. The − and − convolution techniques are the two most popular risk quantification and accumulation techniques [58], [59]. These convolution methods are characterized as the summation of severity and frequency distribution , where for all the contributing factors is stochastically evaluated using the number of fixed −fold convolutions. The existing riskaggregation techniques can derive variables from past data that can evaluate and ; they encounter difficulty in integrating various causal factors to update the predicted risk through evidence propagation with learned information and expert knowledge [59]. Thus, they have difficulty integrating various causal factors to update the predicted risk through evidence propagation with learned information and expert knowledge.
The concept of information or technique hybridization is well-established as a means of combining various techniques to ensure that the weaknesses in some can be compensated by the strengths of others, creating a synergy that enhances the overall robustness of the outcomes. Following this idea, this report proposes a novel dynamic hybrid model (DHM) and demonstrates its proficiency as a comprehensive risk estimation approach for complex systems under multiple dynamic conditions. The proposed method incorporates − and − convolutions in a causal Bayesian dependency model. It uses the advanced features of Bayesian factorization and elimination (BFE) theory to ensure that expert knowledge and information learned from the training data can be validated automatically using the expectation-maximization (EM) technique. Simultaneously, the model enables the quantitative integration of multiple data types, such as discrete and continuous variables, via dynamic Bayesian discretization. The method also incorporates sensitivity analysis features to allow targeted improvements of the overall risk by prioritizing the input frequency variables with significant impacts on the overall estimated risk considering organizational resource constraints.
Similarly, researchers in numerous fields have used multi criteria decision-making (MCDM) frameworks to support risk assessment strategies, and their accuracies have been shown to depend on the quality of data used for the analysis [14], [60]. For instance, the analytical hierarchy process (AHP) has been individually applied as an MCDM tool to analyze complex projects and dynamic issues [61], [62]. One of the widely recognized drawbacks of AHP is the consistency of outcomes across different questions, even when the goal or target remains constant. However, while Saaty used AHP to provide a flexible, systematic, and repeatable evaluation process for selecting optimal alternatives among multiple criteria, the process remains subjective depending on the selected targets, which could be problematic for risk assessment of safety-critical systems such as railway rolling stock systems [63], [64].
The remainder of this paper is organized as follows. Section II describes the proposed DHM technique. Section III presents a comprehensive demonstration of the DHM through a selected case study. Section IV provides an analysis, a discussion, and the implications of the results. Finally, Section V summarizes the conclusion of the study.

II. PROPOSED DHM FOR QUANTIFIED RISK ASSESSMENT
The proposed DHM for overall risk assessment with uncertainty consists of eight fundamental steps, as depicted in Fig. 2. The process commences by collecting and defining the input training data, consisting of frequency and severity variables and distributions. The datasets can represent internal (component inherent failure rates) and external error conditions (such as human errors, electromagnetic interference, and weather conditions). Node probability tables (NPTs) are established using the probability distribution and density functions of the input data. Next, the aggregate risk function is established using the dependency and interdependency factors. The structure is defined using a Bayesian network (BN) with discretization to aggregate both discrete and continuous data. The frequency and severity input data are aggregated using BFE with − and − convolution to predict the overall risk, considering the interdependency that exists among the variables and distributions. The predicted risk is validated and updated with expert knowledge and learned data using Bayesian EM by forward and backward evidence-based propagation. To improve the results, sensitivity analysis can be conducted to identify and prioritize the input frequency, severity variables, and distributions that improve the overall risk, given the organizational resource constraints. Finally, the predicted results can be accepted, or further iterations can be implemented to improve the overall predicted risk using additional expert knowledge or input data.

A. DATA INPUT (TRAINING DATA, NEW VARIABLES, AND NPT DATA)
The model begins with data acquisition and definition of the input variables. This step is the most critical component of the DHM technique because the accuracy of outputs is often indicative of the inputs used to generate them. It includes establishing the system failure rates as well as defining the environmental and human error condition data. Data for a particular task (e.g., a particular repair or replacement activity), can be derived using conventional HRA [37]- [39], HFA [19], [43], [44], and HOF methods [18], [45] In this study, secondgeneration HRA techniques were used due to their ability to mitigate against calculation bias and errors. These abilities are due to the facts that these approaches can adequately overcome some of the limitations that plague first-generation tools, especially data scarcity, inconsistencies associated with the treatment of commission errors, insufficient treatment of performance-shaping factors, cognitive abilities, and lack of systematic task analysis structure [65].

B. ESTIMATION OF RISK FACTORS
In this step, the input data for the frequency and severity factors are defined based on the organization operational context, needs, and constraints. The input data for the frequency and severity factors can be a discrete variable or continuous distribution. Some of the external error conditions for the frequency input data that can be defined include weather patterns, human error conditions, earthquakes, and electromagnetic interference. Similarly, the internal error conditions can consider system failure rates, repair rates, degradation times, etc. The severity or consequence factors can consider the mortality rate (e.g., catastrophic, critical, and moderate), legal cost, financial penalty, business reputation, etc.

C. BUILDING OF THE OVERALL RISK AGGREGATION MODEL
Based on the factors identified in Step 2, an estimation of the overall risk model can be defined using the general risk aggregation formula for fixed systems as [58], [59], [66], [67].
Here, is the sum of system failure consequences, and each for = 0, … , represents independent identically distributed (IID) severity distributions. If originates from a common continuous distribution function , then can be considered a return distribution from written as ~. Thus, ~ is referred to as arbitrary constant − convolution. Now, considering that there are systems and if ~, then (1) can be reformulated as − convolution.
is a recursive − convolution on severity . Therefore, (2) can be modelled as a discrete random variable, ( = ) = , for = 0, . . , , where is the range of frequency distribution . Thus, the distribution of the aggregate consequences can be represented as where each is a constant − convolution, ( ) represents discrete random variable, represents coefficient of the random variable ( ) and (3) represents a compound distribution consisting of mutually exclusive parameters.
Step 5 illustrates how the BFE technique can be implemented to estimate the complex risk model in (3) with various dependencies such as common cause, sequential effects, and nonlinear effects. In addition, mixture variables are computed by utilizing the dependent predetermined functions defined in (4).

D. BN STRUCTURE CONSTRUCTION AND APPLICATION OF DYNAMIC DISCRETISATION
The BN structure was constructed for all the nodes that represent the random variables that influence frequency and severity to establish the NPT [68], [69]. Several scholars [57], [59], [70] have already indicated that most of the existing risk aggregation approaches are limited in their ability to recreate causal interactions between hybrid variables accurately. Therefore, the incorporation of the BN representation here enables the DHM to simulate cause-and-effect dependencies between hybrid parameters in the hypothesis to estimate the marginal probability distributions for the parameters under study. Likewise, the Bayes theorem is agnostic to the direction of dependent factors, with the capability to hypothesize and infer the cause and effect regardless of the form, structure, and path existing among the parameters (or convolution to deconvolution). This operation typically reduces the complexity of inference as described by (5) [71]- [73]: Dynamic discretization entails searching for regions of high density during hypothesizing by summing other specific intervals while eliminating unnecessary intervals through merging and deletion, thus overcoming the deficiencies in static discretization [74]. The entropy error-convergence technique repetitively discretized the target parameters [71], [74], [75]. For example, for a continuous node , if the length of is represented as Ω , then the density of random continuous function can be represented as . Discretisation can be achieved by estimating via the following: • Allocating each interval group with a fixed function . Utilizing the maximum bound of the Kullback-Leibler (KL) value that exists between two primary density functions and , the comparative entropy error stimulated by the discretised parameter can be estimated as With the KL parameter, the maximum parameter of a discretized function ̃ can be specified by the average of the parameter contained in the individual intervals of the discretized region. The comparative entropy error obtained from the probability density function ̃ due to discretisation Ψ = { } is evaluated, tested, and checked during the iterative process to ensure that it is below a stated threshold as defined by the convergence stopping rule. The joint marginal distributions for all parameters in the model are evaluated using the inference algorithm upon completing the discretization of each parameter.

E. INITIAL RISK ESTIMATION
The BFE method that performs convolution on the hybrid frameworks requires risk accumulation with or without causal effects that can be sequential, resulting from a common cause, and other nonlinear dependences. To reduce the conditional probability tables, the BFE method consists of three separate steps, each performing a specific optimization task during the convolution process [71], [76]- [79].

1) LOG-BASED AGGREGATION (LBA)
LBA computes the − convolution based on a log-based pattern, thereby improving computational efficiency, rather than aggregation by straight summation in (4). According to (3), , where = 1, … , represents the aggregation of specific parent parameters −1 and , where the summation consists of repetitive sums of the same parameter . The overall aggregate parameter, , is estimated when two intermediate parameters are produced by two parents forming a hierarchy chain during the factorization process. To avoid the computational delay and reduce the cost, LBA first calculates and exploits prior estimated results inductively. Therefore, an individual stage in the process enables the iteration of results from the posterior steps without the creation of a different BN structure.

2) VARIABLE ELIMINATION (VE)
In this stage, parameters are repetitively deleted in the LBA iteration step, significantly enhancing the computational efficiency for evaluating arbitrary fixed − convolutions. The objective of VE is to eliminate nodes in the BN structure through marginalization for mutually exclusive variables outside the boundary of the query sets. However, VE reduces the number of variables in the analysis but integrates additional steps to capitalize on the repeated structures and frameworks within the model. For an − convolution of IID severity parameters, the graph in Fig. 3 = ∑ ( ( 0 ,.., , 1,.., −1) However, by exploiting the conditional independence relations in Fig. 3(a), it can be observed that all the pairs of and +1 are independent of each other. Therefore, it is feasible to marginalize each pair of and +1 individually, as shown in the model. Thus, based on the "query blocks," (7) can be rewritten as With (9), it is feasible to marginalize each pair of parents and +1 iteratively out of the model through elimination or pruning. For example, the elimination order in (9) could be { 0 , 1 }, { 1 , 2 },…,{ −1 , }. The marginal probability distribution of , which refers to the last query set, is then determined at the final deletion stage.
Nevertheless, in the scenario in which BN exhibits common cause effects including dependency features as illustrated by 2 in Fig. 3(b), the set to be directly eliminated during the VE process includes the leaf nodes. The common parent node is the desired node that should be preserved in the query at each step. Marginal distribution 2 can be evaluated and expressed by computing posterior distribution 2 in BN 2 * , as illustrated in Fig. 3(b), as follows: (11) To eliminate 0 and 1 by marginalising them from (11), The posterior marginal of 2 can also be defined in combination with , 1 , and 2 as Next, 2 and 1 can be expressed as The conditional probability distribution for each parameter −1 can be obtained with the following expression: Therefore, the target − convolution with variable −1 aggregate can be achieved for the posterior probability distribution through the marginalization of by specifying the conditional distribution for variable −1 | . The output from (15) can be considered a universal expression for common causes, dependency, and nonlinear effects with multiple influencing factors for risk prediction.

3) COMPOUND DENSITY FACTORISATION (CDF)
CDF is the final process in the BFE method, by which the algorithm factorizes the compound sum of (3) to reduce large NPTs into smaller ones. Consider the compound probability density function, ( ), as expressed in (3) Therefore, the conditionally deterministic expression for variable −1 , which is referred to in BN parlance as the partitioned node, can be expressed as In addition, 0 and 1 are mutually exclusive variables; the marginal probability distribution for variables 0 is This expression is similar to the first two terms for compound density function ( ) in (3). Hence, the marginal probability for assumes Therefore, through the BFE method, the compound density function in (3) can be computed efficiently. Note that the robustness of BN enables any type of labelling node to be considered. Hence, other labelling such as "Yes" and "No" or "On" and "Off," among others, could be used instead of "True" and "False." Appoh and Yunusa-Kaltungo: Dynamic hybrid model for comprehensive risk assessment

F. INCORPORATION OF EXPERT KNOWLEDGE INTO RISK ESTIMATION
Assume that we wish to express a joint probability density function through maximum likelihood estimation (MLE); however, some data or information is missing. In this case, we can use expert knowledge and learned data as observable variables to update the model for greatly improved overall risk prediction. Let be the missing or hidden data, be the dimensionality of the data vector, be the vector of the unknown parameter, and training data Ɗ = { | = 1: }, Ɲ( , ) be a multivariate Gaussian distribution of the hybrid model with as the mean of the multivariate distribution with variance . Then, we can use the EM in two steps to update the model. The goal is to perform MLE for , which becomes ̂ as [80]- [82] ̂= ( , −1 ) .
At the expectation ( ) step, we can evaluate the expected comprehensive data log-likelihood at iteratively as follows: where are the expected sufficient statics, ∑ −1 represents the inverse covariance matrix, ∑ represents the covariance matrix and ( − ) represents transpose of vector ( − ).
For the maximation (M) step, we optimize the auxiliary function ∇ ( , ( −1) ) = 0, where Equations (31) and (32) provide the inference for EM application for a dynamic hybrid model with missing and unknown variables, where we cannot replace variables by their expectations and apply MLE. Instead, we calculate the expectation of sufficient statics and insert that into the MLE in (20).

G. SENSITIVITY ANALYSIS
Upon completing the EM learning process to validate and improve the predicted overall risk, further improvement can be achieved through prioritization and targeting the frequencies with the highest impact on aggregated risk. Let ( ̅ , ) be a discrete frequency node in a compound risk estimation; then, the sensitivity node can be expressed with a generalized approach [71], [83]: where ( = | ) is the current probability value for , given evidence, and ( = | , ̅ = ̅ ) is the new value taken by when values for the set of observable variables, ̅ , are inserted into the BN model. Note that the current instantiated nodes denoted by evidence will be excluded from ̅ . On evaluation, the joint sensitivity of the target, being the overall risk to perturbations in the source nodes, is exponential in time and space. Therefore, pairwise sensitivity can be adopted with a generalized approach to (33): Considering inward analysis whereby all the values are set on the source variable ̅ , assessing the change in and broadcasting involves changing only the target node and evaluating the changes in the source set ̅ [83]. In addition, the broadcasting evaluates changes to the source nodes in parallel and provides time-saving opportunities as follows: In the case of continuous data, sensitivity analysis is derived by broadcasting the changes in value for and storing the full distributions: ( = | , = ), ( = | ), ( = | ), as required in the discrete case. The statistics and percentiles derived at the end will give the variety of results proportional to the corresponding impact on risk, and these can be represented in a tornado graph for different observable variables with respect to the target node (i.e., the overall risk) [71].

H. OVERALL SYSTEM RISK ESTIMATION
Further iterations of the processes in Section II.F can be initiated by iterating and updating the model with new expert knowledge and data learned from the training data as more information becomes available via design, operation, and model maturity. The DHM can be repeated and iterated from the initial step described in Section II.B with a new observable variable for frequency or severity, and from the descriptions in Section II.F with new expert knowledge and learned data. The iterations and computations for the DHM follow the described stable-entropy and low-entropy error-convergence stopping rules [74], [75], [84].
where is the total number of iterations, ( ) = ∑ is the approximate relative entropy error, and is the observable three-sequential iterations to establish whether the entropy has converged to a stable parameter with the restricted region (1 − , 1 + ). The LEE stopping rule determines the absolute entropy error threshold convergence for NPT.

I. DHM ALGORITHM FOR COMPREHENSIVE RISK ASSESSMENT
Given input training data with frequency distribution and variables and severity distributions and variables, the DHM algorithm is as follows: 1. Initialize the function ( ) in (3). 2. Construct and compute the prior probability variables and distributions for the input NPTs based on the input data. 3. Compute the posterior probability for the NPTs using (5) and (6). 4. Compute the marginal distribution based on the BN structure and interdependency factors in Step 2 using (7)-(15). 5. Compute the initial overall risk ( ) using the convergence stopping rules in (36) and (37). 6. Validate the predicted risk via Bayesian EM in (20) with learned data and expert knowledge using the convergence stopping rules in (36) and (37). 7. Compare and contrast the predicted initial risk in Step 4 and validated risk in Step 5. 8. Accept the predicted overall risk or continue iteration with new expert knowledge and learned data. 9. Identify and prioritize the highest input frequency variable and distributions for further improvement in the predicted risk via sensitivity analysis using (33)- (35). 10. After completing Step 9, iterate Steps 4-7. 11. With new input data for severity and frequency, iterate Steps 1-9. 12. End.

III. CASE STUDY: TRAIN COUPLER FAILURE
Section II provides a full description of the theoretical framework of the proposed approach. For a practical demonstration of its proficiency, a semi-permanent coupler of a three-car EMU train operating in the UK was used to demonstrate the overall risk assessment. This approach enables the implementation of robust mitigation factors and facilitates the adoption of safety management plans. The coupler incorporates a center pivot, which consists of male and female mechanical couplers with two-socket joints that are placed over the flange of the coupler halves. The coupler joins two train carriages together and has a deformation tube that absorbs excess forces during collision. The core is pressed into the deformation tube, which in turn expands, thereby enabling the coupler to be compressed and for energy to be absorbed. A rubber doughnut anchor absorbs energy, which enables pivoting and transfer of forces between the center section and vehicle underframe. Pneumatic jumper hoses are provided for the main reservoir pipe and brake pipe between intermediate ends. Jumper connections for train carriages are permanently mounted between the intermediate ends of each car to enable electric power and pneumatic hoses to be connected. Failure of the semi-permanent coupler can be catastrophic, especially for high-speed trains, and could lead to derailment, collision, and fires with multiple casualties. The objective of this case study was to assess the overall risk for a three-car EMU fleet based in the UK similar to Fig. 1 considering coupler failure rate, human error conditions during maintenance, and electromagnetic interference (EMI) from the rail infrastructure. Therefore, the sources of risks considered for this case study include human or maintenance errors, electromagnetic interference (EMI) from railway infrastructure, and couplers inherent failure rates.
The analysis began with data curation and NPT definitions for semi-permanent coupler failure frequency and severity variables. The frequency data have four subcategories: coupler 1 failure rates, coupler 2 failure rates, human error probabilities (HEP) due to coupler maintenance activities, and EMI. For semi-permanent couplers, EMI is a product of rail infrastructure and other trains, which cause intermittent loss of electrical and signal transmissions along the jumper cables, affecting critical subsystem functionalities such as door control systems, braking systems, and lighting systems. Three severity variables were identified: severity related to fatalities of up to five people, severities related to more than five fatalities, and severities related to financial loss to the operator due to settlements, penalties, and rail line closures. The HEP for the semi-permanent coupler was estimated for experienced and inexperienced technicians using the second-generation cognitive reliability and error analysis method (CREAM). The semi-permanent coupler maintenance activity encompassed the following tasks: • Task 1: Visual inspection of the coupler and associated electrical devices for damage. • Task 2: Inspection of the coupler with NDT for cracks and damage.
• Task 3: Replacement of any damaged jumper cables including wiring. • Task 4: Ensuring that all jumpers and locks are securely fitted and test the system for correct functionality. Tables I and II show the common performance conditions (CPCs) that can be defined for Tasks 1-4, which enable the estimation of the context influence index (CII) [27], [85]:  The overall HEPs for inexperienced and experienced technicians were estimated to be 0.820 and 0.182, as indicated in Tables III and IV. Tasks 2 and 3 demand more significant cognitive aspects; although there are cognitive aspects in Tasks 1 and 4, they are less significant. In addition, the probabilities of failure and success of the couplers were predicted at a train operating time t of 2000 h. All models, including the NPTs, were run for 30 iterations with a convergence error of 0.01, as expressed in (36) and (37). The remainder of the discrete and continuous NPTs were defined for frequency and severity, respectively, as shown in Table V. The overall severity was divided by a factor of 10 to normalize the overall risk scale from 1 (very low) to 6 (very high) to enable comparison with the risk matrix table of the UK-based operating company. The analysis was conducted using AgenaRisk tool [86].

IV. ANALYSIS AND DISCUSSION OF OUTCOMES
The overall BN structure for the complete model, including the frequency, severity, and synthetic nodes for NPTs of various data types is shown in Fig. 5. The synthetic node coupler status is introduced into the BFE computations. A hidden/synthetic link connects the overall risk node to the overall severity and frequency for the BN structure. Firstly, we estimated the overall risk after a train time t of 2000 h. In scenario 1, we assumed normal operating conditions with severities as IID and with the two couplers working as required, maintained by an experienced technician in the presence of low EMI. The overall risk was estimated using (3)- (19). Further, (5) and (6) enabled the mixture of discrete and continuous nodes to be computed conveniently using dynamic discretization. The corresponding marginal distribution for the query, including the overall risk in scenario one, is shown in Fig. 6. The resulting overall risk for train derailment caused by coupler failure was estimated to have a mean of 3.04 with a 99th percentile interval of 0. 33-7.53. Compared to the quantitative risk assessment conducted based on the coupler failure rate only, the overall risk of the derailment was estimated to be 1, which represents low risk. This finding implies that coupler derailment could occur, but it would be an exceptional event and significantly unlikely. However, when additional frequency factors were considered in combination with the inherent failure rate, the overall risk indicated was moderate risk of 3 (3.04 ≈ 3), indicating the event was a distinct possibility but rare.
In scenario 2, assuming IID for the severities, when the human error NPT was set to an inexperienced technician and with EMI set to high as shown in Fig. 7, the worst-case scenario overall risk of coupler derailment was estimated at 6 (6.02 ≈ 6) very high risk, even though the couplers were considered to be fully operational. The risk rating of 6 indicated that the probability of an accident due to coupler derailment was frequent under scenario 2.  The predicted risk ratings from scenarios 1 and 2 were validated through further comparisons and incorporation of new test data from a similar three-car EMU with different EMI distributions and coupler maintenance error probabilities with expert knowledge of the fleet performance. The ratio of expert knowledge to new test data was considered to be 50%, meaning that there was equal confidence with the new test data and with the current NPTs in scenarios 1 and 2. Hence, learned NPT data from scenarios 1 and 2 were used for missing NPT data. The Bayesian EM learning algorithm in (20)−(32) was utilized to conduct the learning and validation process. This step began with the validation of overall risk from scenario 2 in scenario 3 with new NPTs, as shown in Fig. 8. The high EMI was considered to be 59.98%, contrasting with the 100% in scenario 2. Furthermore, based on expert knowledge of test data, HEP for the inexperienced technician was considered as 40%, with couplers 1 and 2 having success probabilities of 79.073% and 60%, respectively, contrasting with the value of 100% in scenario 2. After a maximum of 30 iterations with a convergence threshold of 0.01, the new overall risk was estimated to be 5.72 (mid-range between high and very high risk), as shown in scenario 3 in Fig. 8. Comparing the new overall risk in scenario 3 to scenario 2 indicates a minor difference between the very high risk predicted using original data (6.02) and that predicted based on test data with expert knowledge (5.72). The confidence level achieved for the predicted overall risk is greater in scenario 2. The node with "M" indicates that the NPTs have been learned as a mixture, consisting partially of learned data from the previous NPT and from expert knowledge about the test data. It also means that the expert knowledge to data ratio is higher than 0% and lower than 100% (50% as indicated). "K" shows that the NPTs were solely based on learned data from existing NPTs.
Similarly, the overall risk for scenario 1 (i.e. normal operation) as shown in Fig. 6 was validated using scenario 4 as shown in Fig. 9. The figure shows the same NPTs from expert knowledge on test data and learned data as those discussed in scenario 3. Thus, the conditions in scenarios 1 and 4 were set identically, using an experienced technician with fully functioning couplers and low EMI. The overall risk for coupler derailment in scenario 4 was estimated to be 3 (2.93 ≈ 3). The original risk predicted in scenario 1 (3.00) under normal conditions is similar to that predicted in scenario 4 (2.93), with learned data as well as expert knowledge.  Finally, sensitivity analysis was conducted based on (33)- (35) to prioritize improvement actions for the highest frequency factor with a significant impact on the overall risk to reduce the likelihood of derailment. The overall risk was selected as the targeted node, and the discrete frequency factors (HEP, EMI, couplers 1 and 2 failure rates) and continuous severities were selected as sensitive nodes. As shown in Fig. 10, the overall coupler frequency (i.e., the synthetic node for all individual frequency factors) has the greatest effect on the predicted overall risk, reaching 25.71 on the Tornado scale for train derailment. This value contrasts with the overall severity of 2.929. Further analysis indicated that the EMI had the greatest effect on the overall frequency, the overall risk with a factor of 7.32, followed by human error factors at 4.882, and the inherent factors of couplers 2 and 1 at 4.881 and 3.704, respectively. The results highlight a significant need for EMI improvements, such as the proper insulation of cables and use of twisted pair cables. These measures can help reduce the loop areas and induced voltage, minimizing EMI and improving overall risk for coupler-against-rail derailment. The results indicate that technicians require further training to improve their competence levels and to reduce human errors during coupler maintenance.

V. CONCLUSION
Herein, a hybrid model was proposed for the comprehensive quantitative risk assessment of multiple hazard sources. Unlike other quantitative risk assessment methods, to predict the overall risk from multiple sources, the proposed approach harnesses and aggregates discrete variable and continuous distribution data dynamically. Some of the theoretical and experimental benefits of DHM include the following: i. It can aggregate various risk sources from previous data and integrate the causative elements among the risk sources to aid overall risk prediction. ii.
It uses Bayesian factorization and elimination technique to incorporate expert knowledge to update predicted risk, depending on changing input data sources stochastically. iii.
It enables mathematical combination of several types of data (i.e. discrete and continuous data) to facilitate overall risk estimate, utilizing the Kullback-Leibler dynamic discretization method. iv.
The DHM method's sensitivity analysis function allows for targeted improvements in overall risk by prioritizing the highest input risk, which substantially impacts the overall predicted risk in light of organizational resource constraints.
The BFE algorithm embodied in the DHM technique enables the utilization of interdependent factors between multiple sources of risk, which complements in the overall risk estimation. The dynamic updating feature of the DHM via the Bayesian expectation-maximization technique achieves accurate risk prediction in dynamic operating contexts by enabling the model to be updated using forward and backward propagation with additional data and new expert knowledge. The addition of sensitivity analysis to the DHM technique significantly enhances holistic risk assessment by prioritizing the input data (hazards) with the greatest effects on the predicted risk for further improvements in the context of organizational resource constraints. As a result, when compared to analyzing individual sources of risk independently of one another, the ability to amalgamate the sources of risk using the DHM technique aids in showing the broader degree of risk and their implications. As demonstrated in the case study and various scenarios, the results indicate that the risk predicted using the DHM is far more robust and presents a more comprehensive view of the risk of semi-permanent rail derailment than conventional risk estimation based only on inherent failure rate. It was also observed that the DHM can quantitatively and dynamically aggregate various sources of risk into a single holistic outcome for effective risk mitigation. However, the availability of quality data may be a shortcoming in the implementation of the DHM. This issue can be addressed through more effective data collection methods. While the DHM technique is versatile and robust, it should be tailored to the specific system under study, available data, and verified assumptions to support the analysis. Furthermore, while the DHM has been demonstrated in its application using a UK rolling stock electrical multiple units as a case study, it is worth noting that the application of a different rolling stock configuration in a different country should take into account the rules and regulations governing risk assessments for fleet-wide acceptance into operation in that country. Future studies can consider software and online development of the DHM algorithm, including further validation of the proposed technique via new case studies.

ACKNOWLEDGMENT
Frederick Appoh is grateful to RAMS Engineering and Asset Management Consultancy Limited for sponsoring his ongoing doctoral study. Frederick Appoh would also like to thank the UK-based train operator for providing the data sets used to validate the proposed framework.
His EPSRC-funded research project focuses on the development of safer drones-based inspection and stock level estimation mechanisms for confined spaces to minimise human exposure to inherent hazards. He has published over 45 technical articles (peer-reviewed top quartile journals and conference articles) with internationally reputable publishers. He has also published a book on the condition monitoring of industrial assets.
Dr. Yunusa-Kaltungo is currently a member of several industrial and academic committees and working groups including Institution of Mechanical Engineers (IMechE) Safety & Reliability Working Group (SRWG) and British Standards Institute (BSI). He has delivered several national and international keynotes and invited talks. He has reviewed several internationally reputable journal articles and engineering standards (e.g., International standards Organisation and International Electro-technical Commission). He also has extensive industrial experience with the world's largest manufacturer of building materials -LafargeHolcim PLC in diverse roles, including Health & Safety (H&S) Manager, Maintenance Manager, Training & Learning (T&L) Manager, Reliability Engineer, Mechanical Execution Engineer, Plant Operations Champion and Project Core Team Lead for a multi-million GBP coal plant project. On this project, he will primarily contribute his immense expertise in H&S and T&L towards the successful realisation of the project milestones.