Bayesian method to solve the early failure period of numerical control machine tool

Abstract

The hazard rate curve of the numerical control machine tool is a bathtub curve. The change point between the early failure period and the random failure period of the curve is difficult to obtain with a small data sample; thus, a Bayesian method is proposed. A method to build the prior distributions of the Weibull parameters is developed, which integrates the multi-source prior information of the target numerical control machine tool and the reference numerical control machine tool. The Markov chain Monte Carlo method is adopted to calculate the estimators of the Weibull parameters corresponding to each failure, which solves the problem of the absence of an analytical solution. The total working time of the numerical control machine tool when the estimator of the shape parameter is equal to 1 is estimated by taking the estimator of the shape parameter as the function of time. As a result, the change point and the early failure period are obtained. Comparison result shows that the result obtained through an existing change point solving method with a large dataset is close to the result generated through the proposed method with a small dataset. The change point and the early failure period obtained with the proposed method can be used to guide the early failure test and to design a rational maintenance strategy, which are of vital engineering significance.

Introduction

Numerical control machine tool (NCMT) is a complex mechanical–electrical–hydraulic system, and its hazard rate (HR) function presents a bathtub curve as the working time changes.^1–3 According to the failure frequency, the life cycle of the NCMT is divided into three stages, namely, the early failure period, the random failure period, and the wear-out failure period. The early failure period is caused by the defects in the design, poor processing precision of consistency, welding flaws, cracks, and poor assembly process consistency control,^3–10 and the HR in this period is high and rapidly decreases as the working time increases. After reaching the random failure period, causing by random loads, human error, environment,^3,5,7,8 and the HR is approximately constant.^11–17 The time point when NCMT enters into the random failure period from the early failure period is defined as the change point of the bathtub curve. The change point is a significant reference for manufacturing enterprises to expose the potential failures of NCMTs by conducting an early failure test. Ideally, the NCMTs have totally entered the random failure period before the NCMTs are shipped to users or put into field operation. ³ Based on the change point, the manufacturing enterprises can establish rational maintenance strategies for different failure periods, which is significance to productivity improvement and maintenance cost reduction.^17–19

Several studies have solved the change point of the bathtub curve. For example, Clarotti et al. ¹¹ used a piecewise model to approximate the bathtub curve and to find the change point. Ülkü and Yenigün ¹² used a piecewise constant hazard function with a change point; two estimation methods based on the maximum likelihood ideas are considered. Nguyen ¹³ considered the HR function as a two-section piecewise function to determine the change point. Gupta et al. ¹⁴ utilized the extended Weibull distribution to estimate the change point of the Bathtub curve by the maximum likelihood. Bebbington et al. ¹⁵ and Kim ¹⁶ defined the change point of the bathtub curve as the point where the curvature of HR curve has the maximum value. Jiang ¹⁷ proposed a non-parametric approach to estimate the empirical HR function and a cluster analysis method to find the change point. Bebbington et al. ²⁰ used the empirical estimation method to obtain the change point of the HR function of the modified Weibull distribution. Pham and Nguyen ²¹ used a non-regular model and the bootstrapping method to find the change point. Jiang ²² proposed to use spline functions to approximate the degradation process and to find the change point.

The above methods are all based on large data samples. However, NCMTs have recently presented an obvious trend in production as shown by small lot and customization, rapid model upgrading, and improving reliability level. All these features result in the inadequate failure data of NCMTs in a reliability test and in the difficulty in obtaining the large data sample required by the classical statistic methods in a limited test time. Therefore, the existing methods are unsuitable for finding the change point of the bathtub curve with a small data sample. Thus, a Bayesian method that utilizes prior data to estimate the change point of the HR curve for NCMT is proposed. This method comprises four steps: (1) establishing a Bayesian model for solving the change point, (2) solving the prior distributions of the parameters of the reliability model, (3) calculating the posterior distributions of the parameters and the parameter estimators, and (4) solving the change point.

Bayesian model

Building the prior distributions of Weibull parameters

Weibull distribution and its function

According to the previous research experience on the reliability of NCMT, the two-parameter Weibull distribution is selected as the reliability model of NCMT.^22–27 TBF _r (time between failures, r = 1, 2, …, n), the time between the rth and r+ 1th failure, is a random variable. TBF _r is independent and comes from the same population TBF, which follows a two-parameter Weibull distribution. The cumulative distribution function (CDF) and the probability density function (PDF) are expressed as follows

$$F(t_r|\alpha ,\beta )=F(t_r|\mathbf{\theta })=1-\exp [-{\left(\frac{t_r}{\alpha }\right)}^{\beta }]$$

(1)

$$f(t_r|\alpha ,\beta )=f(t_r|\mathbf{\theta })=\frac{\beta }{\alpha }{\left(\frac{t_r}{\alpha }\right)}^{\beta -1}\exp [-{\left(\frac{t_r}{\alpha }\right)}^{\beta }]$$

(2)

where t_r is the observed value of TBF _r , α> 0 is the scale parameter, β> 0 is the shape parameter, and θ = (α, β) is the parameter vector.

This section aims to establish the prior distributions of α and β. However, α and β do not have obvious physical meanings and cannot be understood directly by experts; ²⁸ thus, some other quantities that can be easily judged by experts should be given.

MTBF (mean time between failures) and STD² (square of standard deviation) represent the expectation and variance of TBF, respectively, which can be obtained by the following formulas

$$\mathrm{MTBF}=E\left(T\right)=\alpha \mathrm{\Gamma }(1+\frac{1}{\beta })$$

(3)

$$\mathrm{ST}{D}^2={\alpha }^2\{\mathrm{\Gamma }(1+\frac{2}{\beta })-{\mathrm{\Gamma }}^2(1+\frac{1}{\beta })\}$$

(4)

The unit of MTBF and STD is hours. Gutiérrez-Pulido et al. ²⁹ considered that describing TBF with mean value and standard deviation (SD) is easy. Therefore, MTBF and STD are selected as the quantities whose value ranges can be given by experts.

Expert judgment procedure

The expert judgment procedure includes four stages, and each stage includes one or more steps:

Stage 1 is the preparation stage, which includes four steps:

• Step 1.1 selects the reference NCMT according to the target NCMT (NCMT to be tested). The type, model, structure, and function of the reference NCMT should be similar to the target NCMT. The reference NCMT should have a large sample of field reliability test data. The best situation is where the target NCMT is upgraded based on the reference NCMT.

• Step 1.2 collects the related information of the reference and target NCMT.

• Step 1.3 selects the reliability model (two-parameter Weibull distribution) and model parameters (α and β) and determines the judgment quantities (MTBF and STD).

• Step 1.4 defines the meaning of MTBF and STD and determines the multi-source of prior information.

Stage 2 includes two steps, in which the technique experts learn and practice their scoring skill:

• Step 2.1 invites the technique experts to learn and practice their scoring skills repeatedly to reduce the subjective bias. The technique experts are the representatives of experienced NCMT designers, manufacturing staffs, outsourcing staffs, after-sales service staffs, NCMT operators, and NCMT maintenance staffs.

• Step 2.2 asks the technique experts to master the multi-source prior information on the reference and target NCMTs. The multi-source prior information includes the following:

• 1. Basic information of NCMT (type, model, manufactory, manufacture date, and price).

• 2. Comparative analysis of the reliability levels of the sub-systems or components between the reference and target NCMTs.

• 3. Information on the reliability test of the reference NCMT and the value ranges of its MTBF and STD.

• 4. Operation behavior of the target NCMT in the reliability test, which includes stability, failure occurrence, failure cause, and failure mode.

Stage 3 asks the technique experts to score and obtain the judgment result. In this stage, the technique experts formally provide the value ranges of MTBF and STD for the target NCMT based on the practice in Stage 2 in combination with the multi-source prior information.

We assume that the number of the expert is m, and the value ranges of MTBF and STD given by the jth (j = 1, 2, …, m) expert are [M_Lj, M_Uj] and [S_Lj, U_Sj]. A weight W_j% is first assigned to the jth expert to indicate his importance or his judgment accuracy based on his experience and degree of familiarity with judgment skill to integrate the judgment results of all experts into two single intervals, [M_L, M_U] and [S_L, U_S]. The sum of all W_j% should be equal to 1. Thus, the final judgment is given as follows

$$\left[M_L,M_U\right]=\left[{\displaystyle \sum _{j=1}^m(W_j\%\times M_{Lj})},{\displaystyle \sum _{j=1}^m(Wj\%\times M_{Uj})}\right]$$

(5)

$$\left[S_L,S_U\right]=\left[{\displaystyle \sum _{j=1}^m(W_j\%\times S_{Lj}),}{\displaystyle \sum _{j=1}^m(W_j\%\times S_{Uj})}\right]$$

(6)

Stage 4 is the feedback stage of the reliability test of the target NCMT. According to the reliability test status, testers determine whether the test should be continued. If the test should be continued, testers feed the failure and maintenance data in the test back to the experts; if the test should be ended, the prior distributions should be obtained, and the final report should be made.

MTBF and STD

MTBF is a direct indicator of the reliability of NCMT. If experts assert that NCMT A is more reliable than NCMT B, then the MTBF of A should be higher than the MTBF of B. STD is the indicator of the difference between t_r and MTBF, as well as the indicator of stability or uncertainty of the NCMT. For two products A and B, experts can compare the reliability levels and the degree of difference between t_r and MTBF based on experience and prior information.

From the point of HR, when β<1, NCMT is in the early failure period; when β = 1, NCMT is in the random failure period; when β> 1, NCMT is in the wear-out failure period. ⁹ The ratio between MTBF and STD depends only on β, which is demonstrated by the following formula

$$\frac{\mathrm{STD}}{\mathrm{MTBF}}=\frac{\alpha \sqrt{\mathrm{\Gamma }(1+\frac{2}{\beta })-{\mathrm{\Gamma }}^2(1+\frac{1}{\beta })}}{\alpha \mathrm{\Gamma }(1+\frac{1}{\beta })}=\frac{\sqrt{\mathrm{\Gamma }(1+\frac{2}{\beta })-{\mathrm{\Gamma }}^2(1+\frac{1}{\beta })}}{\mathrm{\Gamma }(1+\frac{1}{\beta })}$$

(7)

In (0, ∞), STD/MTBF is a monotonically decreasing function of β and approximates to 0 infinitely as β increases. Therefore, when STD/MTBF> 1, NCMT is in the early failure period; when STD/MTBF=1, NCMT is in the random failure period; when STD/MTBF<1, NCMT is in the wear-out failure period.

Transforming expert judgment into prior distributions

Building the prior distributions of Weibull parameters consists of three steps. First, a pair of point estimators of MTBF and STD is transformed into a pair of point values of α and β. Second, the estimation intervals of MTBF and STD are transformed into the intervals of α and β. Third, the prior distributions of α and β are built based on their intervals.

Transformation of point values

M and S are assumed as the two point estimators of MTBF and STD and (M, S) can be transformed to a pair of point estimators (α*, β*). From equations (3) and (4), (α*, β*) should meet the following two equations

$$M={\alpha }^{*}\mathrm{\Gamma }(1+\frac{1}{{\beta }^{*}})$$

(8)

$$S^2=({\alpha }^{*}{)}^2[\mathrm{\Gamma }(1+\frac{2}{{\beta }^{*}})-{\mathrm{\Gamma }}^2(1+\frac{1}{{\beta }^{*}})]$$

(9)

α is removed from equations (8) and (9)

$$\frac{S^2+M^2}{M^2}=\frac{\mathrm{\Gamma }(1+\frac{2}{{\beta }^{*}})}{{\mathrm{\Gamma }}^2(1+\frac{1}{{\beta }^{*}})}$$

(10)

To obtain β*, we define the following function of β

$$f\left(\beta \right)=\frac{\mathrm{\Gamma }(1+\frac{2}{\beta })}{{\mathrm{\Gamma }}^2(1+\frac{1}{\beta })}$$

(11)

The solution of the following equation will be the value of β*

$$f\left(\beta \right)-\frac{S^2+M^2}{M^2}=0$$

(12)

According to the property of the gamma function, $\mathrm{\Gamma }(x+1)=x\mathrm{\Gamma }(x)$ , we can have

$$f\left(\beta \right)=\frac{\mathrm{\Gamma }(1+\frac{2}{\beta })}{{\mathrm{\Gamma }}^2(1+\frac{1}{\beta })}=\frac{\frac{2}{\beta }\mathrm{\Gamma }\left(\frac{2}{\beta }\right)}{\frac{1}{{\beta }^2}{\mathrm{\Gamma }}^2\left(\frac{1}{\beta }\right)}=2\beta \frac{\mathrm{\Gamma }(\frac{1}{\beta }+\frac{1}{\beta })}{\mathrm{\Gamma }\left(\frac{1}{\beta }\right)\mathrm{\Gamma }\left(\frac{1}{\beta }\right)}$$

(13)

According to the relation between the beta function and gamma functions

$$B(x,y)=\frac{\mathrm{\Gamma }\left(x\right)\mathrm{\Gamma }\left(y\right)}{\mathrm{\Gamma }(x+y)}$$

(14)

We can have

$$f\left(\beta \right)=\frac{2\beta }{B(\frac{1}{\beta },\frac{1}{\beta })}$$

(15)

According to equation (15), the curve of f(β) is plotted; f(β) is a monotonically decreasing function of β in (0, ∞) and approximates to 0 infinitely as β increases. Therefore, the abscissa of the intersection point of the curve of f(β) and the horizon line y = [(S² + M²)/M²] is the value of β*. The solution of β* can be obtained with a numerical method. α* can be determined by equation (8).

Transformation of intervals

As shown in Figure 1, according to the curve of f(β) and y = [(S² + M²)/M²], the surface of function β(M, S) can be drawn with intervals [M_L, M_U] and [S_L, U_S]. According to equation (8) and by taking (M, β) as independent variables, the surface of function α (M, β) can be drawn, which is shown in Figure 2.

Figure 1.

Surface of function β(M, S).

Figure 2.

Surface of function α(M, β).

Figures 1 and 2 show that the extreme values of M and S determine the extreme values of β; the extreme values of M and β determine the extreme values of α. Therefore, according to the values (M_L, S_L), (M_U, S_L,), (M_L, S_U), and (M_U, S_U), the extreme values of α and β are determined and their intervals [α_L, α_U] and [β_L, β_U] can be obtained.

Prior distributions of Weibull parameters

Although [α_L, α_U] and [β_L, β_U] are obtained, the most possible value and the interaction between α and β have not been determined. According to Berger, ³⁰ both α and β are uniformly and independently distributed. Therefore, the PDFs of α, β, and (α, β) are regarded as prior distributions

$$\pi (\mathbf{\theta })=\pi (\alpha ,\beta )=\frac{1}{({\alpha }_U-{\alpha }_L)}\times \frac{1}{({\beta }_U-{\beta }_L)}$$

(16)

where α∈(α_L, α_U) and β∈(β_L, β_U).

Posterior distributions of Weibull parameters

The posterior distributions of Weibull parameters can be obtained through Bayes theorem. The prior distribution of θ is π( θ ) and the observed data sample is t = {t1, t2, …, tn} = {tr}, r = 1, 2, …, n. The likelihood function of t can be expressed as

$$L(\mathbf{t}|\mathbf{\theta })=\underset{r=1}{\overset{n}{\Pi }}p(t_r|\mathbf{\theta })=\underset{r=1}{\overset{n}{\Pi }}\left\{\frac{\beta }{\alpha }{\left(\frac{t_r}{\alpha }\right)}^{\beta -1}\exp [-{\left(\frac{t_r}{\alpha }\right)}^{\beta }]\right\}$$

(17)

The marginal distribution of t is

$$\begin{array}{c}\multicolumn{1}{c}{m(\mathbf{t})=\int \pi (\mathbf{\theta })L(\mathbf{t}|\mathbf{\theta })d\mathbf{\theta }}\\ \multicolumn{1}{c}{=\int \int \pi (\alpha ,\beta )\coprod _{r=1}^n\left\{\frac{\beta }{\alpha }{\left(\frac{t_r}{\alpha }\right)}^{\beta -1}\exp [-{\left(\frac{t_r}{\alpha }\right)}^{\beta }]\right\}d\alpha d\beta }\end{array}$$

(18)

π(θ|t) is defined as the posterior distribution of θ . According to Bayes theorem

$$\pi (\mathbf{\theta }|\mathbf{t})=\frac{\pi (\mathbf{\theta })f(\mathbf{t}|\mathbf{\theta })}{f(\mathbf{t})}$$

(19)

By taking the prior distribution, the likelihood function, and the marginal distribution of t into Bayes theorem, the expression of the posterior distribution of (α, β) can be obtained as follows

$$\pi (\alpha ,\beta |t)=\frac{\pi (\alpha ,\beta )\underset{i=1}{\overset{n}{\Pi }}\left\{\frac{\beta }{\alpha }{\left(\frac{t_i}{\alpha }\right)}^{\beta -1}\exp [-{\left(\frac{t_i}{\alpha }\right)}^{\beta }]\right\}}{\int \int \pi (\alpha ,\beta )\coprod _{i=1}^n\left\{\frac{\beta }{\alpha }{\left(\frac{t_i}{\alpha }\right)}^{\beta -1}\exp [-{\left(\frac{t_i}{\alpha }\right)}^{\beta }]\right\}d\alpha d\beta }$$

(20)

Equation (20) is complex and has no analytic solution. The Weibull parameters also have these characteristics. To solve this problem, the Markov chain Monte Carlo (MCMC) method is adopted to simulate the theoretical posterior distributions of α and β.^31,32

Solving the posterior distributions of parameters

For θ = (α, β), the task of the MCMC method is to conduct N times of iterations (e.g. N = 10,000). A new simulation value θ⁽ⁱ⁾ = (α⁽ⁱ⁾, β⁽ⁱ⁾) (i = 1, 2, …) will be generated each time. Subsequently {θ⁽ⁱ⁾ = (α⁽ⁱ⁾, β⁽ⁱ⁾)} can constitute a Markov chain. Assume that the Markov chain reaches the steady-state distribution after B times of iterations. The first B simulation value is considered as the burn-in period. The simulation values after the burn-in period are {θ^(B+ k) =(α^(B+ k), β^(B+ k))}(k = 1, 2, …, N-B), and these can be used to approximate to the target distribution (equation (20)). The mean values of {α^(B+ k)} and {β^(B+ k)} are taken as the point estimators, which are calculated by the following equations

$$\hat{\alpha }=\frac{1}{N-B}\sum _{k=B+1}^N{\alpha }^{(k)}$$

(21)

$$\hat{\beta }=\frac{1}{N-B}\sum _{k=B+1}^N{\beta }^{(k)}$$

(22)

Given that equation (20) is too complex, the software WinBUGS is used to develop and conduct the MCMC method. By conducting the MCMC method, WinBUGS generates considerable simulated parameter values and provides a large number of statistics of the parameters automatically, which include the posterior mean values of α and β. To conduct the MCMC method in WinBUGS software, equation (20) should be rewritten as follows

$$\pi (\alpha ,\beta |\mathbf{t})\propto \pi (\alpha ,\beta )\underset{r=1}{\overset{n}{\Pi }}\left\{\frac{\beta }{\alpha }{\left(\frac{t_r}{\alpha }\right)}^{\beta -1}\exp [-{\left(\frac{t_r}{\alpha }\right)}^{\beta }]\right\}$$

(23)

The BUGS model of the target distribution is based on equation (23). An intact BUGS model includes three parts, namely, the prior distributions, the likelihood function, and the data. Each part needs to be described in BUGS language.

Comparison and verification

Two reliability tests are conducted on two batches of machining centers. In the reliability test on the first batch, a large data sample is collected and used to compare the traditional change point estimation method with the Bayes change point estimation method proposed in this study. When the Bayes method is verified, it is applied to the second batch of machining centers with a small data sample.

Obtaining failure data

The first batch of machining centers includes seven products of the same model M1, which are denoted by V01, V02, …, V07. A total of 60 failure data points are obtained, as shown in Table 1.

Table 1.

Reliability test data of seven machining centers.

Type	Number	Observed TBF (h)
M1	V01	20	112	220	285	307	406	897	1005	1321
	V02	53	97	190	427	503	590	724	838	906
	V03	49	251	314	329	702	799	535	554	817
	V04	41	294	305	310	314	1114	1861
	V05	58	154	320	400	476	769	794	1112
	V06	75	264	93	198	627	801	1046	1232
	V07	45	108	102	175	372	260	872	812	768	1105

TBF: time between failures.

Solving the change point with the Bayesian method

According to the expert judgment procedure in section “Building the prior distributions of Weibull parameters,” M1 is taken as the reference NCMT, and V01, V02, …, V07 are taken as the target NCMTs in turn. The experts determine the value ranges of MTBF and STD based on the multi-source prior information, and their judgments are integrated into the final result as follows: [M_L, M_U] = [375, 702], [S_L, U_S] = [549, 732]. Through equation (8), equation (15), and the transformation procedure in “Transforming expert judgment into prior distributions,” the extreme values of α and β are obtained as follows: α_L = 221.99; α_U = 758.75; β_L = 0.55; β_U = 1.29. The prior distributions of α and β are π(θ) = π(α)π(β) = 1/(α_U − α_L)×1/(β_U − β_L) = 1/397.20.

V01 is taken as the target NCMT, and the data samples of V01 corresponding to each failure, t₁ = {20}, t₂ = {20, 112}, t₃ = {20, 112, 220} are taken into WinBUGS model. The corresponding βs are obtained by the MCMC method. As shown in Table 2, the sum of the observed TBFs in each data sample is the total working time of the NCMT, and β changes as the total working time increases. Therefore, β is taken as the function of the total working time, β(t). Table 2 shows that the total working time t_c that satisfies β(t) = 1 can be worked out with the linear interpolation method, and t_c is defined as the length of the early failure period of V01. The early failure period length of V02, …, V07 can be worked out in the same way, which is shown in Table 2.

Table 2.

Estimators of the beta of the first batch of machining centers with different data samples.

Number	β						β = 1
	First 1 TBF	First 2 TBF	First 3 TBF	First 4 TBF	First 5 TBF	First 6 TBF	tc (h)
V01	0.8466	0.8627	0.9041	0.9472	0.9905	1.0230	1062.68
V02	0.9064	0.9107	0.9474	0.9841	1.0122	1.0331	1051.60
V03	0.8867	0.9281	0.9712	1.0105	1.0224	1.0358	855.10
V04	0.8782	0.9220	0.9650	1.0048	1.0401	1.0068	912.61
V05	0.8953	0.9231	0.9669	1.0060	1.0360	1.0430	879.51
V06	0.9070	0.9487	0.9581	0.9952	1.0075	1.0100	873.68
V07	0.8836	0.8980	0.9168	0.9557	0.9935	1.0300	843.67

TBF: time between failures.

Cluster analysis to solve the change point

An earlier study ¹⁷ drew the bathtub curve with a non-parametric method based on a large data sample. The change point of the bathtub curve can then be obtained through the clustering analysis. ³³ The detailed steps are as follows:

1. The observed y_i (i = 1, 2, …, n) are rearranged in an ascending order. The rearranged data are divided into m intervals according to their values. The number k_j (j = 1, 2, …, m) of the data points in each interval is counted.

2. t₁ is used to represent the interval $(0,(t_1+t_2)/2)$ , t_i (1<i<n) to represent the interval $((t_{i-1}+t_i)/2,(t_i+t_{i+1})/2)$ , and $\mathrm{\Delta }{t}_i$ to represent the length of each interval. The HR of each interval can be expressed as

$$r_i\left(t_i\right)=\frac{k_i}{N\mathrm{\Delta }{t}_i}$$

(24)

• 3. The bathtub curve is drawn based on equation (24), as shown in Figure 3.

• 4. A partition point i = i_c is needed to divide the data sample into two groups using the cluster method. The data in one group are close to another and are considerably different from those in the other group. The first group average µ₁(i) and the second group average µ₂(i) are given as follows

$${\mu }_1\left(i\right)=\frac{1}{i}\sum _{j=1}^i{y}_j,{\mu }_2\left(i\right)=\frac{1}{n-i+i}\sum _{j=1}^n{y}_j$$

(25)

Figure 3.

Empirical intensity function.

• 5. The “distance” between y_i and the first (second) group is measured by

  $$d_1\left(i\right)=|{\mu }_1\left(i\right)-y_i|,d_2\left(i\right)=|{\mu }_2\left(i\right)-y_i|$$

  (26)

• The “similarity” between y_i and the first (second) group is measured by

  $$S_1\left(i\right)=1-d_1\left(i\right),S_2\left(i\right)=1-d_2\left(i\right)$$

  (27)

The curves of S₁ and S₂ are drawn, as shown in Figure 4.

Figure 4.

Plots of similarities for the change point.

• 7. Using linear interpolation, the estimate of t_c is derived by

$$t_c=\frac{S_1\left(t_{i_c}\right)-S_2\left(t_{i_c}\right)}{S_2\prime -S_1\prime }+t_{i_c}$$

(28)

where $S_1\prime =\frac{S_1\left(t_{i_c+1}\right)-S_1\left(t_{i_c}\right)}{t_{i_c+1}-t_{i_c}},S_2\prime =\frac{S_2\left(t_{i_c+1}\right)-S_2\left(t_{i_c}\right)}{t_{i_c+1}-t_{i_c}}$ .

Through the six steps above, the change point of the first batch of machining centers can be solved; its value t_c is 974.22 h.

The eighth column of Table 2 presents the estimated values of t_c derived by the Bayesian method proposed by this study, in which the highest and lowest values of t_c are 1062.68 and 843.67 h, respectively. The mean value of t_c in the eighth column is 925.55 h, which is close to 974.22 h. The value is obtained through the clustering analysis and based on the large data sample. The distance between 925.55 and 974.22 h is 2 or 3 days. From the engineering perspective, this distance nearly has no influence on the economic interest of enterprises. The preceding comparison and discussion verify that the Bayes method based on small data sample proposed by this study can be the best way to estimate the early failure period when the traditional method is invalidated by the lack of data. The proposed Bayesian method is applied to a small data sample in the following section.

Example

Target machine tool

The reliability test is conducted on the second batch of machining centers from 26 October 2014 to 28 May 2015. This batch has only one machining center, and its model is represented by M2. During the test, five failures occur, as shown in Table 3.

Table 3.

Small-sample data of the target machining center.

Type	Number	Observed TBF (h)
		First	Second	Third	Fourth	Fifth
M2	TBF (h)	82	582	1092	356	568

TBF: time between failures.

Prior distributions of Weibull parameters

According to the expert judgment procedure in section “Building the prior distributions of Weibull parameters,” an expert group composed of five experts is organized. Given that M2 is upgraded based on M1, we take M1 as the reference machining center. Then, the value ranges of MTBF and STD are [415.64, 628.08] and [359.61, 604.86], respectively, with the significance level 0.05.

According to the data above and the multi-source prior information of M1, the five experts give the value ranges of MTBF and STD about M2, as shown in Table 4.

Table 4.

Results of expert judgment.

Expert index	Expert’s weight	Estimated interval [MLj, MUj] (j = 1, 2, …, 5)	Estimated interval [SLj, SUj] (j = 1, 2, …, 5)
1	W1% = 25%	[ML1, MU1] = [580, 790]	[SL1, SU1] = [540, 840]
2	W2% = 20%	[ML2, MU2] = [595, 800]	[SL2, SU2] = [665, 890]
3	W3% = 20%	[ML3, MU3] = [575, 805]	[SL3, SU3] = [590, 895]
4	W4% = 25%	[ML4, MU4] = [540, 780]	[SL4, SU4] = [640, 875]
5	W5% = 10%	[ML5, MU5] = [600, 810]	[SL5, SU5] = [660, 930]
Combination	100%	[ML, MU] = [578, 797]	[SL, SU] = [619, 886]

According to the transformation method in section “Transforming expert judgment into prior distributions,” (M_L, M_U) and (S_L, S_L) are transformed into intervals (α_L, α_U) = (447.3984, 862.7213), (β_L, β_U) = (0.6855, 1.2983). Assume that α and β are independently and uniformly distributed in the corresponding intervals; then, the prior distribution of the parameter vector is π(θ) = π(α)π(β) = 1/(α_U − α_L)×1/(β_U − β_L) = 1/254.5099.

Posterior distributions of Weibull parameters

Based on equation (23), the MCMC method is conducted in WinBUGS to simulate the posterior distributions derived using equation (20). To make the estimation accurate, three Markov chains are generated in WinBUGS, and the iteration number of each parameter is set as 10,000.

First, the WinBUGS model is established with the first data sample t₁ = 82 h. The statistics of α and β given by the MCMC method are shown in Table 5.

Table 5.

Posterior statistics obtained by MCMC simulation in WinBUGS.

Parameter	Posterior mean	MC	SD
α	637.1	0.7866	119.1
β	0.9694	0.001172	0.1743

MCMC: Markov chain Monte Carlo; SD: standard deviation.

From Table 5, the following equations can be obtained: MC_α = 0.7866<SD_α×5% = 119.1×5% = 5.955; MC_β=0.001172<SD_β×5%=0.1743×5%=0.008715.

The MC values of α and β are all less than 5% of SD. Therefore, the calculation results are accepted. In other words, the estimations of the posterior distributions of α and β are good.

The iteration processes of the dynamic simulation are shown in Figures 5 and 6.

Figure 5.

Iteration processes of α.

Figure 6.

Iteration processes of β.

Figures 7 and 8 present the nuclear density curves of α and β, where the first 500 simulation values are abandoned as burning-in data. The remaining 9500 simulation values are used to simulate the posterior PDFs of α and β.

Figure 7.

Nuclear density curves of α.

Figure 8.

Nuclear density curves of β.

β₂ = 0.9988 is obtained with the MCMC method based on the WinBUGS model, which is built by the first two data (t₁ = 82 h, t₂ = 582 h). The other estimation values of β can be worked out in the same way, based on the first three, first four, and first five data, which are shown in Table 6.

Table 6.

The total working time of NCMT and the corresponding β

Time of failure (h)	82	664	1765	2121	2689
β	0.9694	0.9988	1.0150	1.037	1.063

Solving the change point

According to Table 6 and by taking β as the function of the total working time t, the total working time that satisfies β(t)=1 can be derived easily with the linear interpolation method; the result is t_c=715.56 h. t_c = 715.56 h is regarded as the early failure period of M2. Therefore, the Bayesian method proposed in this study can estimate the early failure period of NCMT with a few machines and less time.

Conclusion

This study proposes a Bayesian method to estimate the early failure period of NCMT with a small data sample, based on the reliability model of NCMT conforms to the two-parameter Weibull distribution. By taking the Weibull shape parameter β as the function of time t, we obtain the total working time that satisfies β(t)=1 as the change point of the bathtub curve. The conclusions are as follows:

1. An expert judgment procedure, which integrates the multi-source prior information of NCMT and the reliability test data of the reference NCMT, is designed and applied. The judgment result is transformed into the prior information of the Weibull parameters. The expert judgment procedure can avoid the subjective bias and significantly promote the accuracy of the modeling.

2. The cluster analysis method is applied to a large data sample of machining center with type M1, working out the early failure period as 974.22 h. The method proposed by this study is applied to the small data sample of M1, working out the mean value of the early failure period as 925.55 h. Comparison shows that the estimation result of the proposed Bayesian method with a small data sample is close to that of the classic method with a large data sample, and the error can be omitted from the engineering perspective.

3. The early failure period of M2, which is upgraded from M1, is obtained by the Bayesian method and the result is 715.56 h. The result of the comparison between M1 and M2 shows that the early failure period decreases with the improvement of the NCMT. This method can be used to guide the early failure test and to design rational maintenance strategies in different failure periods, which are of vital significance.