System Maintainability Estimation with Multi-Peak Time Distribution based on the Bayesian Melding Method

Mochao Pei; Jianping Hao; Cuijuan Gao

doi:10.2478/msr-2025-0014

Full Article

1.

INTRODUCTION

Maintainability by design is an important quality characteristic of a product. A high level of maintainability means that repairs are quick, easy and economical [1]. An in-depth study of maintainability can support the development, maintenance and support of equipment, providing a solid material foundation for combat effectiveness. An important aspect of maintainability studies is maintainability estimation, which in many cases is an important part of the operational suitability assessment in testing and evaluation (T&E) of current equipment [2]. Maintainability (a probability measure), mean time to repair (MTTR), and maximum time to repair [3] are commonly used metrics for maintainability evaluation. According to Military Handbooks 470A (MIL-HDBK-470A), the sample size for the maintainability estimation should be at least 30 to ensure a high degree of reliability of the estimates [4]. However, it is practically impossible to obtain sufficient data for maintainability estimation in certain phases of equipment testing. This is partly due to the considerable cost of testing and partly due to the fact that mapping the performance of the equipment is no longer the primary purpose of testing. This may result in insufficient fault samples being available for maintenance operations.

Thanks to advances in data acquisition and storage technologies, it is possible to collect data from other test phases of the equipment, expert knowledge, and data from similar equipment in addition to the current field test data. If we want to integrate this data, Bayesian theory is a natural choice, and the prior distribution specification is a key issue in the Bayesian framework. Zellner et al. [5] introduced a data quality factor to measure the quality of prior information and field test data, which in turn leads to more accurate prior distributions. Ibrahim et al. [6] investigated the power prior distribution and applied it to regression estimation with good results. Zhou C et al. [7] proposed a demonstration method using mixture prior distributions for the problem of small-sample maintainability demonstration with multiple sources of prior information, using credibility weighting based on quality factors to integrate multiple prior distributions into a single one.

In some phases of T&E, such as operational testing (OT), one faces not only the problem that the sample size of maintenance time is small, but probably also the problem of under-representation of maintenance operations. Due to limitations in test time and conditions, some failure modes that take longer to uncover may not show up, while others may lack excitation conditions, resulting in a lack of the types of operations to fix them. When Bayesian theory is used to integrate prior information, this missing information about maintenance operations can be accounted for at the level of the entire system. Consequently, there are two levels of prior information to consider: historical maintenance data observed from a microscopic viewpoint that is similar to the type of maintenance operations in the field test, and historical maintenance data observed from a macroscopic viewpoint for the entire system that contains this maintenance data. In fact, there are many studies in the field of reliability assessment that have focused on the integration of different levels of prior information. Guo J et al. [8] effectively integrated various expert knowledge and data sources at subsystem and system levels using the Bayesian melding method (BMM) in analyzing system reliability. Yang L et al. [9] integrated multilevel prior information using an improved BMM, which flexibly balanced the contributions of the prior distributions involved in the integration by setting the weighting factors as hyperparameters. From the perspective of system theory, the integration of different levels of information manifests the dialectical relationship between the whole and the local, which helps to obtain a more reasonable prior distribution. However, no corresponding studies have yet been conducted in the field of maintainability assessment.

When determining the prior distribution, the addition of inappropriate prior information can lead to a poor posterior distribution, which can have a negative impact on the maintainability estimation. It is therefore also necessary to check whether there is any significant conflict between the prior information and the field test data. The detection of such a conflict is generally performed by the consistency test in previous studies. Typical methods include the graph comparison method [10] and the Bayesian credible interval method [11], etc. Zhang Z [12], Zhu Z [13], and Wang J [14] all screened prior information from the perspective of data consistency. They used parametric or non-parametric methods to check whether the prior information and the field data were statistically derived from the same population. However, the consistency test is performed at a certain significance level. If this level is exceeded, the prior distribution is discarded, which is a “black or white” approach that results in the loss of useful prior information.

To solve the aforementioned problems, this paper proposes a maintainability estimation method that considers the integration of two different levels of prior information with the field data. The four main contributions of this work can be summarized as follows:

In response to the problems of insufficient maintenance data and underrepresentation of maintenance operation types that existed in some test phases, a BMM-based maintainability estimation method that considers the prior information on different maintenance operations and the one about system level is proposed.
In view of a possible conflict between a prior and the field data, a checking method for a prior-data conflict is given when the maintenance time follows the exponential, normal, and lognormal distribution, respectively, as well as a method for avoiding a prior-data conflict based on mixture priors.
Adaptive sampling importance resampling (ASIR) is flexibly applied to several links of the method proposed in this paper to solve computational difficulties.
An interesting perspective is provided to gain a deeper understanding of the maintenance time distribution for complex equipment.

The rest of this paper is organized as follows: Section 2 presents the motivation for the research in this paper. Section 3 discusses the BMM-based maintainability estimation framework in detail. Sections 4 and 5 present a numerical case and a practical test case to validate and illustrate the benefits, respectively. The last section briefly summarizes this paper and draws three conclusions.

2.

MOTIVATION

The maintenance time is a random variable that is usually considered to follow an exponential, normal or lognormal distribution [15]. However, it follows a mixture distribution when there is more than one time-dependent factor in the maintenance process, and its probability density curve then takes on a multi-peaked shape [16]. Literature [17] and [18] point out that the maintenance time sample set for complex equipment is a combination of samples from different types of maintenance operations, which is a typical heterogeneous dataset. If there are significant differences in the time expectations for these different types of maintenance operations, the sample set histogram will show multiple peaks, the typical sign of a mixture distribution. This representation is not incorrect, but it needs to be explained in more detail.

In general, simple or basic maintenance operations follow a normal distribution, maintenance operations that can be completed after a short adjustment time or a quick replacement of parts follow an exponential distribution, and the lognormal distribution is suitable for describing the maintenance time of all types of complex equipment [19]. The maintenance of complex equipment consists of many different types of maintenance operations, with the types that occur less frequently taking longer, and vice versa. This mechanism leads to a distribution of maintenance time that is skewed to the "left" and has a long "tail" to the right. Ideally, there should be a standard lognormal distribution to model the maintenance time of complex equipment. In practice, however, this is not the case in some test phases. For example, for equipment intended for a specific combat mission, the mission-critical subsystems are only some of its components, which, together with factors such as the limitations of test conditions and test time, means that certain types of maintenance operations are not performed, so some "small populations" are missing in the lognormal distribution "large population" that should be present, thus creating the phenomenon of "multiple peaks". In other words, the lognormal distribution is the "intrinsic nature" of the maintenance time of complex equipment, while the multi-peak phenomenon is the "extrinsic manifestation" under certain conditions.

3.

METHODOLOGIES

The proposed method is described in detail in this section and its flowchart is shown in Fig. 1. First, the natural prior distributions for the system and each type of maintenance operation are determined based on historical data, and the induced system prior is evaluated based on the structural model and the natural priors for all maintenance operations. Then the BMM is used to obtain an updated prior for each maintenance operation, while a prior-data conflict elimination is used to determine an uninformative prior for each maintenance operation. The uninformative prior and the update prior are combined to generate the mixture prior and thus the posterior for each maintenance operation. Finally, the system maintainability is estimated based on the structural model and the posterior of all maintenance operations.

3.1.

Model for system maintainability estimation

The system maintainability M_s(t) at a given time t is a weighted mean of n probabilities [20]: (1) $M_{s} (T) = \sum_{i = 1}^{n} w_{i} M_{i} (T)$ {M_s}(T) = \sum\limits_{i = 1}^n {{w_i}{M_i}(T)} where M_i(t) is the probability that the i^th maintenance operation is completed at time t. The weight $w_{i} = λ_{i} / \sum_{i = 1}^{n} λ_{i}$ {w_i} = {\lambda_i}/\sum\nolimits_{i = 1}^n {{\lambda_i}} , and λ_i denotes the failure rate corresponding to the failure object of the i^th maintenance operation. To illustrate the methodology presented here, we consider three maintenance operations that belong to the same system and whose operation times follow exponential, lognormal or normal distributions. Thus, n=3, and the following equation can easily derived: (2) ${\begin{array}{l} M_{1} (t) = 1 - e^{- θ_{1} t} \\ M_{2} (t) = \frac{1}{θ_{3} \sqrt{2 π}} \int_{0}^{t} \frac{1}{x} exp [- \frac{1}{2} {(\frac{ln x - θ_{2}}{θ_{3}})}^{2}] d x \\ M_{3} (t) = \frac{1}{θ_{5} \sqrt{2 π}} \int_{0}^{t} exp [- \frac{1}{2} {(\frac{x - θ_{4}}{θ_{5}})}^{2}] d x \end{array}$ \left\{{\matrix{{{M_1}(t) = 1 - {e^{- {\theta_1}t}}} \hfill \cr {{M_2}(t) = {1 \over {{\theta_3}\sqrt {2\pi}}}\int_0^t {{1 \over x}\exp \left[ {- {1 \over 2}{{\left({{{\ln x - {\theta_2}} \over {{\theta_3}}}} \right)}^2}} \right]{\rm{d}}x}} \hfill \cr {{M_3}(t) = {1 \over {{\theta_5}\sqrt {2\pi}}}\int_0^t {\exp \left[ {- {1 \over 2}{{\left({{{x - {\theta_4}} \over {{\theta_5}}}} \right)}^2}} \right]{\rm{d}}x}} \hfill \cr}} \right. where θ₁,θ₂,θ₃,θ₄,θ₅ are the parameters. Historical system-level data are less prone to “multi-peak” phenomena, because they are often richer and more patterned than field test data. They can naturally be modeled by a lognormal distribution, and the system maintainability function in this case can be expressed as follows: (3) ${M^{'}}_{s} (t) = \frac{1}{ϕ_{2} \sqrt{2 π}} \int_{0}^{t} \frac{1}{x} exp [- \frac{1}{2} {(\frac{ln x - ϕ_{1}}{ϕ_{2}})}^{2}] d x$ {M'_S}\left(t \right) = {1 \over {{\phi_2}\sqrt {2\pi}}}\int_0^t {{1 \over x}\exp \left[ {- {1 \over 2}{{\left({{{\ln x - {\phi_1}} \over {{\phi_2}}}} \right)}^2}} \right]{\rm{d}}x} where ϕ₁ and ϕ₂ are the parameters.

3.2.

Determination of prior distributions using historical data

The parameters θ₁,θ₂,θ₃,θ₄,θ₅ in (2) are all random variables in the Bayesian context whose prior distributions can be assumed to be in the form of conjugate priors, and probability density functions (p.d.f.s) of the following form can be derived [21]: (4) ${\begin{array}{l} π_{1} (θ_{1}) = \frac{β_{1}^{α_{1}}}{Γ (α^{1})} θ_{1}^{α^{1} - 1} e^{- β^{1} θ^{1}} \\ π_{2} (θ_{2}, θ_{3}^{2}) = \frac{\sqrt{k_{2}}}{2 π} \frac{{(\frac{λ_{2}}{2})}^{\frac{r_{2}}{2}}}{Γ (\frac{r_{2}}{2})} {(θ_{3}^{2})}^{\frac{- r_{2} - 3}{2}} e^{- \frac{λ_{2} + k_{2} {(θ_{2} - μ_{2})}^{2}}{2 θ_{3}^{2}}} \\ π_{3} (θ_{4}, θ_{5}^{2}) = \frac{\sqrt{k_{3}}}{2 π} \frac{{(\frac{λ_{3}}{2})}^{\frac{r_{3}}{2}}}{Γ (\frac{r_{3}}{2})} {(θ_{5}^{2})}^{\frac{- r_{3} - 3}{2}} e^{- \frac{λ_{3} + k_{3} {(θ_{4} - μ_{3})}^{2}}{2 θ_{5}^{2}}} \end{array}$ \left\{{\matrix{{{\pi_1}({\theta_1}) = {{\beta_1^{{\alpha_1}}} \over {\Gamma \left({{\alpha^1}} \right)}}\theta_1^{{\alpha^1} - 1}{e^{- {\beta^1}{\theta^1}}}} \hfill \cr {{\pi_2}\left({{\theta_2},\theta_3^2} \right) = {{\sqrt {{k_2}}} \over {2\pi}}{{{{\left({{{{\lambda_2}} \over 2}} \right)}^{{{{r_2}} \over 2}}}} \over {\Gamma \left({{{{r_2}} \over 2}} \right)}}{{(\theta_3^2)}^{{{- {r_2} - 3} \over 2}}}{e^{- {{{\lambda_2} + {k_2}{{({\theta_2} - {\mu_2})}^2}} \over {2\theta_3^2}}}}} \hfill \cr {{\pi_3}\left({{\theta_4},\theta_5^2} \right) = {{\sqrt {{k_3}}} \over {2\pi}}{{{{\left({{{{\lambda_3}} \over 2}} \right)}^{{{{r_3}} \over 2}}}} \over {\Gamma \left({{{{r_3}} \over 2}} \right)}}{{(\theta_5^2)}^{{{- {r_3} - 3} \over 2}}}{e^{- {{{\lambda_3} + {k_3}{{({\theta_4} - {\mu_3})}^2}} \over {2\theta_5^2}}}}} \hfill \cr}} \right. which contains one Gamma p.d.f. and two Normal Inverse Gamma (N-IG) p.d.f.s. Similarly, the joint prior density form of the parameters ϕ₁ and ϕ₂ in (3) is: (5) $π (ϕ_{1}, ϕ_{2}^{2}) = \frac{\sqrt{k}}{\sqrt{2 π}} \frac{{(\frac{λ}{2})}^{\frac{r}{2}}}{Γ (\frac{r}{2})} {(ϕ_{2}^{2})}^{\frac{- r - 3}{2}} e^{- \frac{λ + k {(ϕ_{1} - μ)}^{2}}{2 ϕ_{2}^{2}}}$ \pi \left({{\phi_1},{\phi_2}^2} \right) = {{\sqrt k} \over {\sqrt {2\pi}}}{{{{\left({{\lambda \over 2}} \right)}^{{r \over 2}}}} \over {\Gamma \left({{r \over 2}} \right)}}{({\phi_2}^2)^{{{- r - 3} \over 2}}}{e^{- {{\lambda + k{{({\phi_1} - \mu)}^2}} \over {2{\phi_2}^2}}}}

The hyperparameters α₁ and β₁ of π₁(θ₁) can be deter-mined using the moment method if historical test data are available [22]. If the maintenance time X follows an exponential distribution, the expected value μ_m and the variance $σ_{m}^{2}$ \sigma_m^2 of its marginal density m(x|α₁, β₁) are each represented as follows: (6) $μ_{m} = \int_{0}^{\infty} μ (θ) π_{1} (θ_{1}) d θ_{1} = \frac{β_{1}}{α_{1} - 1}$ {\mu_m} = \mathop \smallint \nolimits_0^\infty \mu \left(\theta \right){\pi_1}\left({{\theta_1}} \right){\rm d}{\theta_1} = {{{\beta_1}} \over {{\alpha_1} - 1}} (7) $σ_{m}^{2} = \int_{0}^{+ \infty} (σ^{2} (θ) + {[μ (θ) - μ_{m}]}^{2}) π_{1} (θ_{1}) d θ_{1} = {(\frac{β_{1}}{α_{1} - 1})}^{2} \frac{α_{1}}{α_{1} - 2}$ \sigma_m^2 = \mathop \smallint \nolimits_0^{+ \infty} ({\sigma^2}(\theta) + {[\mu (\theta) - {\mu_m}]^2}){\pi_1}({\theta_1}){\rm{d}}{\theta_1} = {\left({{{{\beta_1}} \over {{\alpha_1} - 1}}} \right)^2}{{{\alpha_1}} \over {{\alpha_1} - 2}}

After replacing μ_m and $σ_{m}^{2}$ \sigma_m^2 with the sample mean x̄ and sample variance v² of X, respectively, the estimates of the hyperparameters α₁ and β₁ can be derived by a simple calculation as follows: (8) ${\begin{array}{l} {\hat{α}}_{1} = \frac{2 v^{2}}{v^{2} - {\bar{x}}^{2}} \\ {\hat{β}}_{1} = \frac{(v^{2} + {\bar{x}}^{2}) \bar{x}}{v^{2} - {\bar{x}}^{2}} \end{array}$ \left\{{\matrix{{\hat \alpha_1 = {{2{v^2}} \over {{v^2} - {{\bar x}^2}}}} \hfill \cr {\hat \beta_1 = {{({v^2} + {{\bar x}^2})\bar x} \over {{v^2} - {{\bar x}^2}}}} \hfill \cr}} \right.

Since π₂(θ₂, θ₃²), π₃(θ₄, θ₅²) and π(ϕ₁, ϕ₂²) are all N-IG densities, their hyperparameters are determined by the same approach. Taking π(ϕ₁, ϕ₂²) as an example, the hyperparameters μ, k, r and λ are estimated as follows [23]: (9) ${\begin{array}{l} \hat{μ} = A_{1} \\ \hat{k} = 1 \\ \hat{r} = \frac{- \frac{14}{3} A_{1}^{4} + 4 A_{1}^{2} A_{2} + 2 A_{2}^{2} - \frac{4}{3} A_{4}}{- \frac{2}{3} A_{1}^{4} + A_{2}^{2} - \frac{1}{3} A_{4}} \\ \hat{λ} = \frac{- \frac{\hat{r}}{6} (A_{2} - A_{1}^{2}) (5 A_{1}^{4} - 6 A_{1}^{2} A_{2} + A_{4})}{- \frac{7}{3} A_{1}^{4} + 2 A_{1}^{2} A_{2} + A_{1}^{2} - \frac{2}{3} A_{4}} \end{array}$ \left\{{\matrix{{\hat \mu = {A_1}} \hfill \cr {\hat k = 1} \hfill \cr {\hat r = {{- {{14} \over 3}{A_1}^4 + 4{A_1}^2{A_2} + 2{A_2}^2 - {4 \over 3}{A_4}} \over {- {2 \over 3}{A_1}^4 + {A_2}^2 - {1 \over 3}{A_4}}}} \hfill \cr {\hat \lambda = {{- {{\hat r} \over 6}\left({{A_2} - {A_1}^2} \right)\left({5{A_1}^4 - 6{A_1}^2{A_2} + {A_4}} \right)} \over {- {7 \over 3}{A_1}^4 + 2{A_1}^2{A_2} + {A_1}^2 - {2 \over 3}{A_4}}}} \hfill \cr}} \right. where (10) $A_{k} = \frac{1}{n} \sum_{i = 1}^{n} X_{i}^{k}, k = 1, 2, 3, 4$ {A_k} = {1 \over n}\sum\limits_{i = 1}^n {X_i^k,k = 1,2,3,4} is the k^th moment of the sample of maintenance time X.

3.3.

Updating the prior using BMM

In the Bayesian context, M₁, M₂ and M₃ in (2) can also be regarded as functions of θ₁, (θ₂, θ₃) and (θ₄, θ₅), respectively. Then the mapping relationship from the input vector θ = (θ₁, θ₂, θ₃, θ₄, θ₅) to the output variable M_s can be established, denoted as W: θ → M_s, and the model can be further defined as M_s = W(θ), which is exactly the system maintainability model shown in (1) and (2). Since one M_s can be mapped by more than one θ, the model is irreversible. Take the distribution obtained from expert knowledge or historical experimental data as the natural prior and let q₁(θ) be the natural prior of θ. In (3), M^′_s (t) is strictly monotonically decreasing with respect to ϕ₁. Therefore, the following inverse function can be easily determined: (11) $ϕ_{1} = h (ϕ_{2}, M_{s}^{'}) = ln t - ϕ_{2} \cdot Φ^{- 1} (M_{s}^{'})$ {\phi_1} = h({\phi_2},{M'_s}) = \ln t - {\phi_2} \cdot {\Phi^{- 1}}({M'_s}) where Φ⁻¹(⋅) is the inverse function of the cumulative distribution function (c.d.f.) of the standard normal distribution. The natural prior q₂(M_s ) of the system maintainability can be solved by combining (5) and (11) as follows: (12) $\begin{array}{l} q_{2} (M_{s}) = \\ = \int_{0}^{+ \infty} π (ln t - ϕ_{2} \cdot Φ^{- 1} (M_{s}), ϕ_{2}^{2}) | \frac{\partial [ln t - ϕ_{2} \cdot Φ^{- 1} (M_{s})]}{\partial M_{s}} | d ϕ_{2}^{2} = \\ = \frac{{(\frac{λ}{2})}^{\frac{r}{2}} \sqrt{k}}{Γ (\frac{r}{2})} e^{\frac{{(Φ^{- 1} (M_{s}))}^{2}}{2}} \int_{0}^{+ \infty} {(ϕ_{2}^{2})}^{\frac{- r - 2}{2}} e^{- \frac{λ + k {(ln t - ϕ_{2} \cdot Φ^{- 1} (M_{s}) - μ)}^{2}}{2 ϕ_{2}^{2}}} d ϕ_{2}^{2} \end{array}$ \matrix{{{q_2}({M_s}) =} \hfill \cr {= \mathop \smallint \nolimits_0^{+ \infty} \pi \left({\ln t - {\phi_2} \cdot {\Phi^{- 1}}\left({{M_S}} \right),{\phi_2}^2} \right)\left| {{{\partial \left[ {\ln t - {\phi_2} \cdot {\Phi^{- 1}}\left({{M_S}} \right)} \right]} \over {\partial {M_S}}}} \right|{\rm{d}}{\phi_2}^2 =} \hfill \cr {= {{{{\left({{\lambda \over 2}} \right)}^{{r \over 2}}}\sqrt k} \over {\Gamma \left({{r \over 2}} \right)}}{e^{{{{{({\Phi^{- 1}}({M_S}))}^2}} \over 2}}}\int_0^{+ \infty} {{{\left({{\phi_2}^2} \right)}^{{{- r - 2} \over 2}}}{e^{- {{\lambda + k{{(\ln t - {\phi_2} \cdot {\Phi^{- 1}}({M_S}) - \mu)}^2}} \over {2{\phi_2}^2}}}}{\rm{d}}{\phi_2}^2}} \hfill \cr}

The integral in (12) does not yield a closed expression and can be approximated by numerical integration. The structural model M_s = W(θ) corresponds to a transformation of θ. Consequently, a distribution of M_s, called the induced prior distribution $q_{1}^{*} (M_{s})$ q_1^*\left({{M_s}} \right) of system maintainability, can be determined based on q₁(θ) and M_s = W(θ). Since the model M_s = W(θ) is irreversible, numerical calculations are required to derive $q_{1}^{*} (M_{s})$ q_1^*\left({{M_s}} \right) . To do this, a large number of random samples of θ are first generated from q₁(θ). Then, for each random sample, the model M_s = W(θ) is used to evaluate M_s, resulting in a large number of random samples of M_s. Finally, $q_{1}^{*} (M_{s})$ q_1^*\left({{M_s}} \right) is estimated based on these random samples using kernel density estimation (KDE) [24]. In BMM, $q_{1}^{*} (M_{s})$ q_1^*\left({{M_s}} \right) and q₂(M_s ) are integrated by geometric pooling, which is: (13) $q_{3}^{*} (M_{s}) = ε_{α} q_{1}^{*} {(M_{s})}^{α} q_{2} {(M_{s})}^{1 - α}$ q_3^*\left({{M_s}} \right) = {\varepsilon_\alpha}q_1^*{({M_s})^\alpha}{q_2}{({M_s})^{1 - \alpha}} where $q_{3}^{*} (M_{s})$ q_3^*\left({{M_s}} \right) is the pooled prior of system maintainability, α is the weight, and ɛ_α is a normalization constant. $q_{3}^{*} (M_{s})$ q_3^*\left({{M_s}} \right) can be further transformed into the updated prior q_θ(θ) of θ: (14) $q_{θ} (θ) = q_{3}^{*} (M_{s}) \frac{q_{1} (θ)}{q_{1}^{*} (M_{s})} = ε_{α} q_{1} (θ) {[\frac{q_{2} [W (θ)]}{q_{1}^{*} [W (θ)]}]}^{1 - α}$ {q_\theta}\left({\boldsymbol{\theta}} \right) = q_3^*\left({{M_s}} \right){{{q_1}\left({\boldsymbol{\theta}} \right)} \over {q_1^*\left({{M_s}} \right)}} = {\varepsilon_\alpha}{q_1}\left({\boldsymbol{\theta}} \right){\left[ {{{{q_2}\left[ {W\left({\boldsymbol{\theta}} \right)} \right]} \over {q_1^*\left[ {W\left({\boldsymbol{\theta}} \right)} \right]}}} \right]^{1 - \alpha}}

The foregoing details show that based on the deterministic structural model and natural priors, the BMM integrates the prior information of each maintenance operation to the system level and further transfers the pooled system information to the maintenance operation level. The information integration and updating for different levels is thus realized. q_θ(θ) integrates the prior information of M_s and θ, and also integrates the structural information of the model M_s = W(θ).

3.4.

Checking and avoiding prior-data conflict

There is a significant discrepancy between the prior information at system level and the field test data. When integrated according to (14), there may still be a prior-data conflict, although this discrepancy is mitigated to some extent by the structural model M_s = W(θ). It is necessary to avoid such a conflict while preserving useful information in the prior.

3.4.1.

Conflict check

In the Bayesian framework, a prior-data conflict arises whenever the prior concentrates most of its mass in the low-density region of the likelihood [25], which is then measured by P(s₀) in the following equation: (15) ${\begin{array}{l} P (s_{0}) = F_{s} (m_{s} (s) \leq m_{s} (s_{0})) \\ m_{s} (s) = \int_{Θ} p (s | θ) p (θ) d θ \end{array}$ \left\{{\matrix{{P\left({{s_0}} \right) = {F_s}\left({{m_s}\left(s \right) \le {m_s}\left({{s_0}} \right)} \right)} \hfill \cr {{m_s}\left(s \right) = \int_\Theta {p(s|\theta)p\left(\theta \right){\rm{d}}\theta}} \hfill \cr}} \right. where S denotes a minimum sufficient statistic (MSS) for the parameter θ in the sampling model, which in a sense condenses all the information about θ in the data. s₀ denotes the observation of S, p(θ) is the prior p.d.f. of θ, and p(s|θ) is the marginal p.d.f. of S. m_s (s) is known as the p.d.f. of the prior predictive distribution. The smaller P(s₀) is, the stronger is the tendency that the main mass of the data is concentrated in the tails of m_s(s), i.e., the stronger is a prior-data conflict.

For a sample x = (x₁, …, x_n ) of maintenance time X, S = x̄⁻¹ is an MSS for θ₁ in (2) if X follows an exponential distribution. Since x̄⁻¹ ∼ IG(n, nθ₁) [26], where IG is an abbreviation for Inverse Gamma, so: (16) $p ({\bar{x}}^{- 1} | θ_{1}) = \frac{{(n θ_{1})}^{n}}{Γ (n)} {({\bar{x}}^{- 1})}^{- n - 1} e^{- \frac{n θ_{1}}{{\bar{x}}^{- 1}}}$ p({\bar x^{- 1}}|{\theta_1}) = {{{{(n{\theta_1})}^n}} \over {\Gamma \left(n \right)}}{({\bar x^{- 1}})^{- n - 1}}{e^{- {{n{\theta_1}} \over {{{\bar x}^{- 1}}}}}}

Since the prior distribution of θ₁ is Gamma(α₁, β₁), the posterior distribution is then given by θ₁|x ∼ Gamma(α₁ + n, β₁ + nx̄) [27]. Consequently, the prior predictive density of x̄⁻¹ can be derived according to Bayes' theorem: (17) $m ({\bar{x}}^{- 1}) = \frac{π_{1} (θ_{1}) p ({\bar{x}}^{- 1} | θ_{1})}{π_{1} (θ_{1} | x)} \propto \frac{{({\bar{x}}^{- 1})}^{- n - 1}}{{(β_{1} + n / {\bar{x}}^{- 1})}^{α_{1} + n}}$ m\left({{\bar x^{- 1}}} \right) = {{{\pi_1}\left({{\theta_1}} \right)p({\bar x^{- 1}}|{\theta_1})} \over {{\pi_1}({\theta_1}|x)}} \propto {{{{({\bar x^{- 1}})}^{- n - 1}}} \over {{{({\beta_1} + n/{\bar x^{- 1}})}^{{\alpha_1} + n}}}} where π₁(θ₁|x) is the p.d.f. of Gamma(α₁ + n, β₁ + nx̄). S = (x̄, v²) is an MSS for (θ₄, θ₅²) in (2) if X follows a normal distribution, where v² is the sample variance of X. Then x̄ ∼ N(θ₄, θ₅²/n) is independent of v² ∼ (θ₅²/(n − 1))χ²_(n−1), so that the marginal density of (x̄, v²) is: (18) $p (\bar{x}, v^{2} | θ_{4}, θ_{5}^{2}) = \frac{1}{\sqrt{2 π \frac{θ_{5}^{2}}{n}}} e^{- \frac{{(\bar{x} - θ_{4})}^{2}}{\frac{2 θ_{5}^{2}}{n}}} \frac{n - 1}{θ_{5}^{2}} \frac{e^{- \frac{ν^{2} (n - 1)}{2 θ_{5}^{2}}} {(\frac{v^{2} (n - 1)}{θ_{s}^{2}})}^{\frac{n - 3}{2}}}{2^{\frac{n - 1}{2}} Γ (\frac{n - 1}{2})}$ p(\bar x,{v^2}|{\theta_4},{\theta_5}^2) = {1 \over {\sqrt {2\pi {{{\theta_5}^2} \over n}}}}{e^{- {{{{(\bar x - {\theta_4})}^2}} \over {{{2{\theta_5}^2} \over n}}}}}{{n - 1} \over {{\theta_5}^2}}{{{e^{- {{{\nu^2}\left({n - 1} \right)} \over {2\theta_5^2}}}}{{\left({{{{v^2}\left({n - 1} \right)} \over {\theta_s^2}}} \right)}^{{{n - 3} \over 2}}}} \over {{2^{{{n - 1} \over 2}}}\Gamma \left({{{n - 1} \over 2}} \right)}}

Actually, the following equation is clearly known from the form of π₃(θ₄, θ₅²) in (4): (19) $θ_{4} | θ_{5}^{2} \sim N (μ_{3}, θ_{5}^{2} / k_{3}), θ_{5}^{2} \sim IG (r_{3} / 2, λ_{3} / 2)$ {\theta_4}|{\theta_5}^2 \sim {\rm{N}}({\mu_3},{\theta_5}^2/{k_3}),{\theta_5}^2 \sim {\rm IG} \left({{r_3}/2,{\lambda_3}/2} \right)

The posterior distribution of (θ₄, θ₅²) is then given as: (20) ${\begin{array}{l} θ_{4} | θ_{5}^{2}, x \sim N (μ_{x}, (n + k_{3})^{- 1} θ_{5}^{2}) \\ θ_{5}^{2} | x \sim IG (r_{3} / 2 + n / 2, λ_{x}) \end{array}$ \left\{{\matrix{{{\theta_4}|{\theta_5}^2,x \sim {\rm{N}}({\mu_x},\left({n + {k_3}{)^{- 1}}{\theta_5}^2} \right)} \hfill \cr {{\theta_5}^2|x \,\,\,\,\, \sim {\rm{IG}}\left({{r_3}/2 + n/2,{\lambda_x}} \right)} \hfill \cr}} \right. where μ_x = (n + k₃)⁻¹(μ₃k₃ + nx̄) and λ_x = λ₃/2 + (n − 1)v²/2 + n(x̄ − μ₃)²/2(n/k₃ + 1).

The joint prior predictive density of (x̄, v²) can be derived according to Bayes' theorem as: (21) $m (\bar{x}, v^{2}) = \frac{π_{3} (θ_{4}, θ_{5}^{2}) p (\bar{x}, v^{2} | θ_{4}, θ_{5}^{2})}{π_{3} (θ_{4}, θ_{5}^{2} | x)} \propto \frac{{(v^{2})}^{\frac{n - 3}{2}}}{{(λ_{x})}^{\frac{r_{3} + n}{2}}}$ m\left({\bar x,{v^2}} \right) = {{{\pi_3}\left({{\theta_4},{\theta_5}^2} \right)p\left({\bar x,{v^2}|{\theta_4},{\theta_5}^2} \right)} \over {{\pi_3}\left({{\theta_4},{\theta_5}^2|x} \right)}} \propto {{{{({v^2})}^{{{n - 3} \over 2}}}} \over {{{({\lambda_x})}^{{{{r_3} + n} \over 2}}}}}

In the above description, P(s₀) in (15) does not yield a closed analytical expression, whose Monte Carlo estimate is therefore given here:

➢
Step 1.
For given α₁ and β₁[μ₃, k₃, r₃ and λ₃], generate θ₁[(θ₄, θ₅²)] from π₁(θ₁)[π₃(θ₄, θ₅²)].
➢
Step 2.
For θ₁ [(θ₄, θ₅²)] generated from step 1, generate x̄⁻¹ [(x̄, v²)] from p(x̄⁻¹|θ₁) [p(x̄, v²|θ₄, θ₅²)], and evaluate m(x̄⁻¹)[m(x̄, v²)] according to (17), (18).
➢
Step 3.
Repeat step 1 and step 2 many times and record the proportion of $m ({\bar{x}}^{- 1}) \leq m ({\bar{x}}_{0}^{- 1}) [m (\bar{x}, v^{2}) \leq m ({\bar{x}}_{0}, v_{0}^{2})],$ m({\bar x^{- 1}}) \le m({\bar x_0}^{- 1})[m(\bar x,{v^2}) \le m({\bar x_0},v_0^2)], where ${\bar{x}}_{0}^{- 1} [({\bar{x}}_{0}, v_{0}^{2})]$ {\bar x_0}^{- 1}[({\bar x_0},v_0^2)] is an observation of x̄⁻¹[(x̄, v²)] based on the sample x = (x₁, …, x_n ). This proportion is a Monte Carlo estimate of P(s₀).

The expressions “[⋅]” in the above three steps correspond to the case where X follows a normal distribution. If X follows a lognormal distribution, it is sufficient to set Y = log( X) and then proceed as in the case of a normal distribution.

3.4.2.

Conflict avoidance

To avoid prior-data conflict, a mixture prior distribution π(θ|D_h, ψ) with the following form is used in this paper: (22) $π (θ | D_{h}, ψ) = ψ π (θ) + (1 - ψ) π (θ | D_{h})$ \pi (\theta |{D_h},\psi) = \psi \pi \left(\theta \right) + \left({1 - \psi} \right)\pi (\theta |{D_h}) where π(θ) is the non-informative prior for the parameter θ, π(θ|D_h) is the informative prior constructed from the historical data D_h, and the weight ψ is adjusted to incorporate as much prior information as possible without prior-data conflict, i.e., (23) $ψ = inf {ψ | P^{π (θ | D_{h}, ψ)} (s_{0}) \geq β}$ \psi = \inf \{\psi |{P^{\pi (\theta |{D_h},\psi)}}({s_0}) \ge \beta \} where P^{π(θ|D_h,ψ)}(s₀) denotes the conflict measure between the mixture prior π(θ|D_h, Ψ) and the field test data calculated according to (15), and β is a tuning parameter whose setting is related to the degree of flexibility assumed by the tester.

In this paper, prior-data conflict checks are required for each of the three maintenance operations. The informative prior π(θ|D_h ) of (22) should be sequentially replaced by the corresponding components in the updated prior q_θ(θ), i.e., q(θ₁), q(θ₂, θ₃) and q(θ₄, θ₅). As for the choice of the non-informative prior π(θ), the absence of the possibility of any prior-data conflict should be considered as a necessary characteristic of any non-informative prior for our purposes, as suggested in [25]. The choice of π(θ) is not unique and should be subjected to a process of trial and error, the result of which is that the variance of π(θ) is much larger than that of π(θ|D_h ). This choice should be such that μ_m ≈ x̄ in (6) if the maintenance time X follows an exponential distribution. If X follows a normal distribution, its marginal distribution is a generalized Student's t distribution, i.e., [23]: (24) $X | μ_{3}, k_{3}, r_{3}, λ_{3} \sim t_{r_{3}} (μ_{3}, \frac{λ_{3} (1 + k_{3})}{r_{3} k_{3}})$ X|{\mu_3},{k_3},{r_3},{\lambda_3} \sim {t_{{r_3}}}\left({{\mu_3},{{{\lambda_3}\left({1 + {k_3}} \right)} \over {{r_3}{k_3}}}} \right) where μ₃ is a mean parameter, r₃ is a degree of freedom parameter, and $\sqrt{λ_{3} (1 + k_{3}) / (r_{3} k_{3})}$ \sqrt {{\lambda_3}\left({1 + {k_3}} \right)/\left({{r_3}{k_3}} \right)} is a scale parameter. The choice should result in μ₃ ≈ x̄ and λ₃(1 + k₃)/r₃k₃ ≈ v². The aforementioned approaches make prior-data conflict extremely weak when a non-informative prior is used. If X follows a lognormal distribution, it is sufficient to set Y = log( X) and then proceed as in the case of a normal distribution.

3.5.

Posterior inference of system maintainability

According to Bayes' theorem, the posterior distribution $π_{θ}^{sub}$ \pi_{\boldsymbol{\theta}}^{sub} can be evaluated on the basis of the mixture prior distribution and the field test data as follows: (25) $π_{θ}^{sub} \propto π (θ | D_{h}, ψ) L (θ)$ \pi_{\boldsymbol{\theta}}^{sub} \propto \pi (\theta |{D_h},\psi)L\left(\theta \right) where L(θ) is the likelihood function for the maintenance time observations. The posterior distribution of the parameters θ₁, θ₂, θ₃, θ₄ and θ₅ can be evaluated according to (25). By combining (1) and (2), the distribution and the estimate of the system maintainability can be determined. The results of the evaluation are determined by simulation algorithms, which are described later.

3.6.

The application of adaptive sampling importance resampling

Sampling importance resampling (SIR), originally developed by Rubin [28], was used for posterior inference to overcome the challenges of sampling from complex priors and likelihood functions. The basic principle of SIR is to take a large number of samples from a known distribution, then reweight the samples, and finally resample the samples according to the weights to form the target distribution. However, the occasional presence of a small number of large importance weights may dominate the resampling process during the algorithm run, which may cause the SIR to not perform well. Liang [29] extended the pruned-enriched Rosenbluth method by setting upper bounds on the weights to limit the use of certain samples. The enrichment method explores a broader range of samples by splitting large weights into small ones. In this paper, ASIR is used for three contexts: the first is the updated prior q_θ(θ) in (14), the second is the mixture prior π(θ|D_h, Ψ) in (22), and the third is the posterior $π_{θ}^{sub}$ \pi_{\boldsymbol{\theta}}^{sub} in (25). The ASIR procedure for q_θ(θ) is shown as follows:

➢
Step 1.
Generate initial random samples. That is, generate a set of samples {θ_j, j = 1, …, m} from q₁(θ).
➢
Step 2.
Calculate resampling weights.
- Generate a set of samples {M_{s j}, j = 1, …, m} using the mapping relationship W: θ → M_s and calculate the resampling weight c_j as follows: (26) $c_{j} = {[\frac{q_{2} [W (θ_{j})]}{q_{1}^{*} [W (θ_{j})]}]}^{1 - α}$ {c_j} = {\left[ {{{{q_2}\left[ {W\left({{{\boldsymbol{\theta}}_j}} \right)} \right]} \over {q_1^*\left[ {W\left({{{\boldsymbol{\theta}}_j}} \right)} \right]}}} \right]^{1 - \alpha}} where W(θ_j) = M_{s j}. Based on {M_{s j}, j = 1, …, m}, KDE is used to obtain $q_{1}^{*} [W (θ_{j})]$ q_1^*[W({{\boldsymbol{\theta}}_j})] , and q₂[W(θ_j )] is calculated according to (12).
- Choose d weights randomly from {c_j, j = 1, …, m} and sort them in descending order to obtain ${c_{j}^{'}, j = 1, \dots, d}$ \{c_j^{'},j = 1, \ldots,d\} .
- Define a set {y_k, k = 1, …, d} with $y_{k} = \sum_{j = 1}^{k} c_{j}^{'}$ {y_k} = \sum\nolimits_{j = 1}^k {c_j^{'}} and assign the element y_⌊γ×d⌋ as the threshold value, where γ is an empirical threshold percentage, e.g., γ = 80 %. ⌊γ × d⌋ is the floor value of γ × d.
- Split any weight c_b > y_⌊γ×d⌋ in {c_j, j = 1, …, m} into g_b = ⌊c_b /y_⌊γ×d⌋ + 1⌋ number of weights c_b /g_b used for resampling the corresponding θ_b.
- Repeat (b)-(d) until the predefined condition, such as a certain standard deviation of the sampling weights, is met.
➢
Step 3.
Resample from {θ_j, j = 1, …, m} with the importance weights derived in step 2.

The ASIR procedure for π(θ|D_h, ψ) is similar to the above procedure, where the initial samples are generated from π(θ) and π(θ|D_h) and the initial resampling weights have only two values, with ψ assigned to the samples from π(θ) and (1 − ψ) assigned to the samples from π(θ|D_h). The distribution of the new samples is π(θ|D_h, ψ). Furthermore, these new samples are resampled with their likelihood function values in (25) as initial weights, resulting in the samples of $π_{θ}^{sub}$ \pi_{\boldsymbol{\theta}}^{sub} .

4.

NUMERICAL CASE FOR VALIDATION

In this section, a synthetic dataset is provided for demonstration purposes. Suppose there are 9 different types of maintenance operations for a system whose times follow normal distributions with different parameters. A set of randomly generated time samples from each of the 9 normal distributions is combined into a large sample set. The parameters are set as shown in Table 1. Obviously, the sample size in Table 1 is also the number of failures corresponding to the maintenance operation, and therefore can be used to calculate the weight w_i in (1).

Table 1.

Parameter settings for maintenance operations.

Number of maintenance operation type	Mean parameter	Variance parameter	Sample size
1	33.49	15.07	132
2	75.35	64.59	24
3	17.88	8.33	187
4	44.45	23.13	96
5	97.05	160.39	12
6	25.29	10.81	172
7	129.12	942.91	4
8	11.21	8.98	154
9	58.00	39.06	53

Three datasets are generated for demonstration. Dataset 1 was generated using the settings in Table 1 with 834 maintenance time samples. The data in dataset 1 belonging to maintenance operations 1, 2 and 4 are filtered out and combined into dataset 2. The time samples for maintenance operations 1, 2 and 4 are generated again to form dataset 3, with sample sizes of 20, 4 and 14, respectively. These three data sets are treated as historical data of the system, historical data of maintenance operations 1, 2, and 4 and field test data, whose distribution characteristics are illustrated with histograms, KDE curves and probability plots, as shown in Fig. 2.

As can be seen in Fig. 2, dataset 1 shows a clear positive skewness and fits well with a lognormal distribution. The other datasets are multimodal. Based on datasets 1 and 2, the natural prior distributions of the system and the three maintenance operations with N-IG forms can be determined using (9). All hyperparameter estimates are listed in Table 2.

Table 2.

Hyperparameter estimates for the natural prior distributions.

Hyperparameter source	$\hat{μ}$ \hat \mu	k̂	r̂	$\hat{λ}$ \hat \lambda
System	3.22	1	13.06	1.97
Maintenance operation 1	33.55	1	8.50	49.87
Maintenance operation 2	74.62	1	4.21	75.17
Maintenance operation 4	44.90	1	6.45	45.98

For a given time t = 25, the induced system prior can be derived based on the structural model M_s = W(θ) and the natural prior distributions of maintenance operations 1, 2, and 4. Subsequently, the updated priors of the three operations can be evaluated according to (14), which are also treated as informative priors in Section 3.4.2. In addition, the non-informative priors of the maintenance operations can also be determined according to the method described in Section 3.4.2, whose hyperparameter settings are listed in Table 3. The p.d.f. plots for the three maintenance operations with informative and non-informative priors are shown in Fig. 3, where the informative priors were determined using KDE.

Table 3.

Hyperparameter settings for the non-informative priors.

Number of maintenance operation type	$\hat{μ}$ \hat \mu	k̂	r̂	$\hat{λ}$ \hat \lambda
1	32.32	0.8	1.3	41
2	72.42	0.4	1.2	50
4	43.10	0.2	0.1	8

As can be seen in Fig. 3, the regions of the heat maps in (b), (d) and (f) are more dispersed and flatter, which is a characteristic of non-informative priors. Using the methods described in Section 3.4.2, the mixture posteriors can be evaluated on the basis of non-informative and informative priors. In addition, the system maintainability at t = 25 can be estimated based on the methods described in Section 3.5. Similarly, by changing the value of t, the system maintainability at other times can also be estimated. For comparison purposes, system maintainability is also estimated using the following three methods:

➢
Method 1:
Dataset 1 (which is treated as the historical data of the system) is fitted with a lognormal distribution whose parameters are estimated using maximum likelihood estimation (MLE), and the values of its c.d.f. are treated as the estimates of system maintainability.
➢
Method 2:
The time samples of each maintenance operation in dataset 3 (which is treated as field test data) are fitted with a normal distribution whose parameters are estimated using MLE. According to (1), the weighted mean of all normal distribution c.d.f values is used as an estimate of the system maintainability.
➢
Method 3:
In the proposed method, the system-level prior information is ignored, i.e., dataset 1 is not used.

The three methods and the proposed method for estimating system maintainability are shown in Fig. 4.

As can be seen from Fig. 4, the estimates of Methods 2 and 3 are not far apart, as the same three normal components were used to generate datasets 2 and 3. Compared to Method 2, the estimates of Method 3 are slightly closer to those of Method 1 because Method 3 uses more time samples than Method 2 and therefore should provide more accurate estimates. The estimates of the proposed method are closer to those of Method 1 than to those of Methods 2 and 3 because it takes into account the system-level prior information compared to Method 3. Method 1 uses 834 time samples. Without taking into account the differences in test conditions between current tests and historical tests, it can be assumed that the estimates of Method 1 are closest to the true results. However, there are always differences between the various tests, and some of these are quite significant. In this study, this is thoroughly considered. Therefore, the proposed method provides the most reasonable results.

5.

A PRACTICAL TEST CASE

The data used in this case comes from the performance tests (PT) and OT for a specific type of wheeled armored vehicle (WAV) in 2023. There are several differences between PT and OT, such as those related to vehicle drivers, road conditions and individual task durations. In addition, the OT covered a much shorter period of time than the PT. In addition to several WACs of this type, there was also some accompanying test equipment that was used to form a tactical system to complete the planned combat missions. For the purpose of subsequent demonstrations and to meet the requirements of proprietary information protection, the actual data used in this case was modified accordingly. When analyzing the maintenance operations in the tests, it was found that some maintenance operations are similar and can be grouped into the same categories. Overall, the 28 maintenance operations that occurred during the OT are grouped into three categories and their maintenance times are found to follow exponential, normal, and lognormal distributions, respectively. In addition, prior information on the three maintenance operations was screened in the PT, including 58 maintenance time samples. Similar to the method in Section 4, the characteristics of the data are shown in Fig. 5. The three methods and the proposed method for estimating system maintainability are shown in Fig. 6.

As can be seen from Fig. 5 and Fig. 6, the characteristics of the datasets are similar to those in Section 4, and the result of comparing the maintainability estimates is also similar, although the practical maintenance operation times contain not only a normal distribution but also an exponential distribution and a lognormal distribution.

6.

CONCLUSION

In some test phases of equipment, due to the limited test time and conditions, the test data used for the maintainability evaluation not only faces the problem of small sample size, but also the problem of insufficient representativeness of maintenance operations, which may cause the lognormal data that should have appeared to have a multi-peak state, which in turn poses a challenge to the maintainability estimation based on Bayesian information fusion. To overcome this challenge, two levels of prior information, the system level and the maintenance operation level, are integrated with the field test data via BMM. Mixture priors are used to avoid prior-data conflicts and a Bayesian posterior distribution is used to estimate system maintainability. The main conclusions are as follows:

BMM can consider prior information at different levels, expanding the channels for data sources.
The mixture priors provide a balance between the non-informative prior and the informative prior, avoiding prior-data conflicts in the Bayesian framework.
Compared with other traditional methods, the proposed method can provide more reasonable estimates for maintainability.

In addition to maintainability, MTTR and maximum repair time are also important metrics for maintainability estimation. The next main task will be to investigate the estimation of other metrics.

System Maintainability Estimation with Multi-Peak Time Distribution based on the Bayesian Melding Method

Full Article

Paradigm

My account