# A Parameter Estimation Method of Shock Model Constructed with Phase-Type Distribution on the Condition of Interval Data

Authors: Long Wang; Yue Li; Yanling Qian; Xu Luo
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

The phase-type distribution (also known as PH distribution) has mathematical properties of denseness and closure in calculation and is, therefore, widely used in shock model constructions describing occurrence time of a shock or its damage. However, in the case of samples with only interval data, modeling with PH distribution will cause decoupling issues in parameter estimation. Aiming at this problem, an approximate parameter estimation method based on building PH distribution with dynamic order is proposed. Firstly, the shock model established by PH distribution and the likelihood function under samples with only interval data are briefly introduced. Then, the principle and steps of the method are introduced in detail, and the derivation processes of some related formulas are also given. Finally, the performance of the algorithm is illustrated by a case with three different types of distributions.


---

## Body

## 1. Introduction

The shock model is usually used in the study of a system’s reliability. The model was first proposed by Esary et al. [1] to describe time intervals between adjacent shocks and its damage by random distributions, as shown in Figure 1. When the cumulative damage (blue line) caused by the shock sequence (green lines) exceeds a certain threshold (red line), the system fails.Figure 1
Traditional shock model.In the shock model, the occurrence time of shocks is usually expressed as a Poisson process [2, 3]; that is, the time interval of adjacent shocks satisfies the exponential distribution, and the one shock damage is also assumed to satisfy the same cluster. Of course, other forms of distributions [4–6], such as gamma distribution, log-normal distribution, and Weibull distribution, also have been used in constructing shock models. While these distribution clusters enhance the model’s adaptability, it also increases the difficulties in estimating their parameters.Observed sample data are required to obtain the model’s parameters. When facing point sample data, their form of log-likelihood function is usually a sum form of its distribution with the logarithmic operator in the front, that is to say, the parameters are easy to be decoupled and methods of MLE (maximum likelihood estimation), including various types of Newton iterative methods [7], can be used for estimating parameters and obtaining a result with a high accuracy. When the sample data are interval data, the likelihood function often contains a large number of convolution operations since the data are incomplete. At this time, exponential distribution is often assumed to simplify the calculation, while it may cause the model to be inconsistent with the actual situation. Of course, some Monte Carlo methods [8] have also been proposed to solve the computational problems caused by convolution, but the estimation accuracy of the methods is usually not controllable.Phase-type distribution (known as PH distribution) is a kind of probability distribution clusters based on continuous-time Markov process [9]. Suppose a continuous-time Markov process Xt,t≥0 with k+1 states, wherein the states 1,...,k are transient states, the state k+1 is an absorbing state, and Xt represents the state of the random process at moment t. Then, the infinitesimal generator matrix of the Markov process can be expressed as Q=D0d100, wherein D0 is a k-dimensional invertible matrix and d1 is a k-dimensional column vector, satisfying d1=−D0e⇀, where e⇀ is a column vector and full elements of 1 with corresponding dimensional.Based on the above Markov process, ak-order continuous PH distribution α,D0 can be defined. A random variable T is distributed as an absorption time taken to transition from any initial state i=1,...,k to the absorbing state k+1, wherein α=α1...αk represents a row vector of probability distribution of each initial phase, D0 is called the phase-type generator, and d1 represents a column vector of absorption rate in which the respective phase transfers to the absorbing state.The cumulative distribution function (CDF) of the PH distribution is(1)Ft=PT≤t=1−αexpD0te⇀,t>0.The probability density function (PDF) of the PH distribution is(2)ft=αexpD0td1,t>0.The PH distribution is dense in the positive abscissa axis; that is, it can fit any type of probability distributions with the random variable located on the interval0,∞. In addition, the PH distribution also has the feature of closure in convolution operations, which greatly simplifies the calculation. Due to the above properties, it is widely used in shock modeling. Neuts and Bhattacharjee [10] first used PH distribution in shock model construction. Ochoa et al. [11, 12] used PH distribution to characterize time intervals of adjacent shocks and one shock damage to analyze reliability of systems. These works are constantly expanding the adaptability of shock models, but it can also be seen that the PH distribution is currently mainly used in shock model constructions under the samples with point data or partial interval data.There are three general parameter estimation methods for PH distribution [13]: maximum likelihood method, Monte Carlo method, and EM method. Among them, the maximum likelihood method, including various types of Newton iterative methods, is often used to solve some specific types of simple PH distributions; the Monte Carlo method is mainly used for specific forms of PH distribution with point data; the expectation maximization method (also known as the EM method) is the most widely used method for parameter estimation; however, it is mainly used in PH distribution with point data [14] or partial interval data [15]. If only interval data can be observed, such as only censored data and interval failure data can be obtained, the model constructed with PH distribution will face the decoupling problem when performing parameter estimation. At this time, the methods described previously cannot be used directly.Aiming at this problem, an approximate method for parameter estimation of the shock model constructed with PH distribution is proposed. By giving an estimated error and constructing PH distribution with dynamic order, the parameters of the shock model with only interval data samples can be estimated. The paper is organized as follows: in Section2, the shock model constructed by PH distribution and the likelihood function with only interval data samples are briefly introduced. Section 3 reconstructs the process statistics based on the framework of the traditional EM method [14], explains the approximation principle in detail, derives the approximate estimated expression of the model parameters, and also gives the total steps of the entire algorithm. Section 4 tests the performance of the algorithm through a case with three different types of distributions; Section 5 summarizes the method proposed in this paper and briefly introduces the follow-up work.


## 2. Shock Model and Likelihood Function

### 2.1. Shock Model Based on PH Distribution

Some assumptions:(1)
The failure threshold of a component isXth, and the censored time is tc.(2)
The shock damageXn>0 caused by the nth n=1,2,... shock is an iid random variable, and the corresponding PDF is fx, x>0, and the CDF is Fx, x>0. The cumulative damage caused by the number of shock n is Xn∗=∑i=1nXi, let X0∗=0, and the corresponding CDF is Fnx,x>0, and the PDF is fnx,x>0.(3)
The occurrence time of thenth n=1,2,... shock is Tn>0, let T0=0, and the corresponding CDF is Gnt, t>0, and the PDF is gnt, t>0. The time interval ΔTn=Tn−Tn−1 is also an iid random variable, and the corresponding PDF is gΔt, t>0, and the CDF is GΔx, t>0.(4)
The occurrence time of a shock is independent of its damage and also independent of the cumulative damage that has been caused by shock sequence.LetGn∗t=PTn−1≤t,Tn>t, t>0, represent the probability that the occurrence time of the nth shock exceed t for the first time. Based on the previous assumptions, the reliability at moment t is(3)Rt=PT>t=P∑n=0NtXn<Xth=∑n=0∞FnXth−Gn+1∗t,where Nt represents the number of shocks that the component receives at the moment t and Xth− represents the left limit of the threshold.The following describes the specific forms of the previous parameters by PH distributions:(1)
Assume that the one shock damageXn satisfies υ order υ=1,2,... PH distribution β,R0, wherein β=β1...βυ represents the probability distribution vector of initial phases, which satisfies ∑i=1υβi=1, and R0=r110...r1υ0.........rυ10...rυυ0 represents the phase-type generator. Let r1=r11...rυ1T=−R0e⇀ represent the column vector of absorption rate. Then, the CDF of Xn is Fx=PX≤x=1−βexpR0xe⇀, x>0. The corresponding PDF is fx=βexpR0xr1, x>0.(2)
The cumulative damageXn∗ caused the number of shock n=1,2,... satisfies the PH distribution γn,Yn0, wherein γn=β0...01,υn represents the distribution probability vector of initial phases and Yn0=R0R1...00R0...0.........R100...R0υn,υn represents the phase-type generator, in which R1=r1β=r111...r1υ1.........rυ11...rυυ1 represents the transfer rate matrix at the occurrence of a shock and the angular scale, respectively, represents the number of rows and columns of a combined vector or matrix. The corresponding column vector of absorption rate is yn1=−Yn0e⇀. Then, the CDF of Xn∗ is Fnx=1−γnexpYn0xe⇀, x>0. The corresponding PDF is fnx=γnexpYn0xyn1, x>0. Let Fn∗x=PXn−1∗<x,Xn∗≥x present the probability that the cumulative damage is no less than x after the nth shock for the first time. Then, Fn∗x=Fn−1x−Fnx=γnexpYn0xe⇀n,Σ, x>0, n=1,2,..., where e⇀n,Σ represents a column vector that the n−1υ+1th to the υmth elements are 1 and others are 0, with corresponding order.The distribution of the adjacent shock time intervalΔTn and shock time Tn is α,D0 (m order) and μn,Mn0. The assumption forms of these two distributions are similar to the shock damage distributions described previously.


### 2.2. Likelihood Function of Interval Data

Assume that the samples are independent of each other and the data can only be obtained at the censored timetc. The information of observed sample data is as follows:(1)
The total number of the samples isS>0(2)
The failure mode set of the samples isδ=δ1...δS, wherein δs=0 represents the censored sample s=1,2,...,S and δs=1 represents the failure sample s(3)
The number of shocks received by the samples isN=n1...nSLetθ=θx,θt represent all parameters to be estimated, where θx=β,R0 represents the parameter set of the distribution of the one shock damage and θt=α,D0 represents the parameter set of the distribution of the time interval. Based on the assumption that the occurrence time of the shock is independent of its damage, the likelihood function can be expressed separately as shown in Table 1.Table 1
Likelihood function of samples.
Sample typeLH typeExpressCensored sampleδs=0Damage functionLxsθx;Xth,ns,δs=∫0XthγnsexpYns0xyns1dxTime functionLtsθt;tc,ns,δs=μns+1expMns+10tce⇀ns+1,ΣFailure sampleδs=1Damage functionLxsθx;Xth,ns,δs=γnsexpYns0Xthe⇀ns,ΣTime functionLtsθt;tc,ns,δs=∫0tcμnsexpMns0tmns1dtIn summary, the overall likelihood function based on the sample data is(4)Lθ;Xth,N,δ=∏s=1SLxsθx;Xth,ns,δs∏s=1SLtsθt;tc,ns,δs.It can be seen in Table1 that parameters in equations with the integral cannot be decoupled easily with the general methods proposed in Section 1. Combined with the assumed model above, observation data, and likelihood function proposed in this section, a specific method on estimating parameters θ is introduced as follows.


## 2.1. Shock Model Based on PH Distribution

Some assumptions:(1)
The failure threshold of a component isXth, and the censored time is tc.(2)
The shock damageXn>0 caused by the nth n=1,2,... shock is an iid random variable, and the corresponding PDF is fx, x>0, and the CDF is Fx, x>0. The cumulative damage caused by the number of shock n is Xn∗=∑i=1nXi, let X0∗=0, and the corresponding CDF is Fnx,x>0, and the PDF is fnx,x>0.(3)
The occurrence time of thenth n=1,2,... shock is Tn>0, let T0=0, and the corresponding CDF is Gnt, t>0, and the PDF is gnt, t>0. The time interval ΔTn=Tn−Tn−1 is also an iid random variable, and the corresponding PDF is gΔt, t>0, and the CDF is GΔx, t>0.(4)
The occurrence time of a shock is independent of its damage and also independent of the cumulative damage that has been caused by shock sequence.LetGn∗t=PTn−1≤t,Tn>t, t>0, represent the probability that the occurrence time of the nth shock exceed t for the first time. Based on the previous assumptions, the reliability at moment t is(3)Rt=PT>t=P∑n=0NtXn<Xth=∑n=0∞FnXth−Gn+1∗t,where Nt represents the number of shocks that the component receives at the moment t and Xth− represents the left limit of the threshold.The following describes the specific forms of the previous parameters by PH distributions:(1)
Assume that the one shock damageXn satisfies υ order υ=1,2,... PH distribution β,R0, wherein β=β1...βυ represents the probability distribution vector of initial phases, which satisfies ∑i=1υβi=1, and R0=r110...r1υ0.........rυ10...rυυ0 represents the phase-type generator. Let r1=r11...rυ1T=−R0e⇀ represent the column vector of absorption rate. Then, the CDF of Xn is Fx=PX≤x=1−βexpR0xe⇀, x>0. The corresponding PDF is fx=βexpR0xr1, x>0.(2)
The cumulative damageXn∗ caused the number of shock n=1,2,... satisfies the PH distribution γn,Yn0, wherein γn=β0...01,υn represents the distribution probability vector of initial phases and Yn0=R0R1...00R0...0.........R100...R0υn,υn represents the phase-type generator, in which R1=r1β=r111...r1υ1.........rυ11...rυυ1 represents the transfer rate matrix at the occurrence of a shock and the angular scale, respectively, represents the number of rows and columns of a combined vector or matrix. The corresponding column vector of absorption rate is yn1=−Yn0e⇀. Then, the CDF of Xn∗ is Fnx=1−γnexpYn0xe⇀, x>0. The corresponding PDF is fnx=γnexpYn0xyn1, x>0. Let Fn∗x=PXn−1∗<x,Xn∗≥x present the probability that the cumulative damage is no less than x after the nth shock for the first time. Then, Fn∗x=Fn−1x−Fnx=γnexpYn0xe⇀n,Σ, x>0, n=1,2,..., where e⇀n,Σ represents a column vector that the n−1υ+1th to the υmth elements are 1 and others are 0, with corresponding order.The distribution of the adjacent shock time intervalΔTn and shock time Tn is α,D0 (m order) and μn,Mn0. The assumption forms of these two distributions are similar to the shock damage distributions described previously.


## 2.2. Likelihood Function of Interval Data

Assume that the samples are independent of each other and the data can only be obtained at the censored timetc. The information of observed sample data is as follows:(1)
The total number of the samples isS>0(2)
The failure mode set of the samples isδ=δ1...δS, wherein δs=0 represents the censored sample s=1,2,...,S and δs=1 represents the failure sample s(3)
The number of shocks received by the samples isN=n1...nSLetθ=θx,θt represent all parameters to be estimated, where θx=β,R0 represents the parameter set of the distribution of the one shock damage and θt=α,D0 represents the parameter set of the distribution of the time interval. Based on the assumption that the occurrence time of the shock is independent of its damage, the likelihood function can be expressed separately as shown in Table 1.Table 1
Likelihood function of samples.
Sample typeLH typeExpressCensored sampleδs=0Damage functionLxsθx;Xth,ns,δs=∫0XthγnsexpYns0xyns1dxTime functionLtsθt;tc,ns,δs=μns+1expMns+10tce⇀ns+1,ΣFailure sampleδs=1Damage functionLxsθx;Xth,ns,δs=γnsexpYns0Xthe⇀ns,ΣTime functionLtsθt;tc,ns,δs=∫0tcμnsexpMns0tmns1dtIn summary, the overall likelihood function based on the sample data is(4)Lθ;Xth,N,δ=∏s=1SLxsθx;Xth,ns,δs∏s=1SLtsθt;tc,ns,δs.It can be seen in Table1 that parameters in equations with the integral cannot be decoupled easily with the general methods proposed in Section 1. Combined with the assumed model above, observation data, and likelihood function proposed in this section, a specific method on estimating parameters θ is introduced as follows.


## 3. Principles, Steps, and Derivation of the Method

It is not difficult to see from Table1 that forms of the likelihood functions are essentially the same, so only the estimating process of parameters θx is described in detail here, which can also be used for estimating θt. According to the main steps of the EM method [14], it is necessary to construct a complete data likelihood function.


### 3.1. Process Statistics Construction

Four variables and an indicator function are defined in the following:(1)
dn,ks represents the phase of the sample s=1,2,...,S after the kth phase transition in the nth n=1,2,... shock damage period;(2)
bn,ks represents the damage caused by phase dn,ks;(3)
nns represents the number of phase transitions of the sample s=1,2,…,S during the nth n=1,2,… shock damage period, including the times when the damage period changes;(4)
Nsh represents the number of additional shocks that may cause failure if the sample sδs=0 is censored with shock number ns. It can be seen that Nsh>0 is a random variable and nsh represents its specific sample value;(5)
Indicator functionδn,isk=1dn,ks=i0dn,ks≠i represents whether the phase of the sample s=1,2,…,S is i=1,2,…,υ after the kth phase transition in the nth n=1,2,… shock damage period.With the above variables and function, four process statistics can be constructed as follows:(1)
The set of statisticsBn,is represents whether the initial phase of the sample s=1,2,…,S is i=1,2,…,υ during the nth n=1,2,… shock damage period. The values of Bn,is are shown in Table 2.Table 2
Value of process statisticsBn,is.
Sample typeShock numberProcess statistics valueCensored sampleδs=01≤n≤ns+nshBn,is=δn,is0, satisfies ∑i=1υBn,is=1n>ns+nshBn,is=0Failure sampleδs=11≤n≤nsBn,is=δn,is0, satisfies ∑i=1υBn,is=1n>nsBn,is=0(2)
The set of statisticsZn,is represents the damage caused by the sample s=1,2,…,S in phase i=1,2,…,υ during the nth n=1,2,… shock damage period. The values of Zn,is are shown in Table 3.Table 3
Value of process statisticsZn,is.
Sample typeShock numberProcess statistics valueCensored sampleδs=01≤n≤ns+nshZn,is=∑k=0nns−1δn,iskbn,ks, satisfies Xth=∑n=1ns+nsh∑i=1υZn,isn>ns+nshZn,is=0Failure sampleδs=11≤n≤nsZn,is=∑k=0nns−1δn,iskn,ks, satisfies Xth=∑n=1ns∑i=1υZn,isn>nsZn,is=0(3)
The set of statisticsMn,ijs represents the total number of phase transition times of the sample s=1,2,…,S when the phase is transferred from i=1,2,…,υ to j=1,2,…,υi≠j during the nth n=1,2,… shock damage period, which does not include the times of the phase transition when the damage period changes. The values of Mn,ijs are shown in Table 4.Table 4
Value of process statisticsMn,ijs.
Sample typeShock numberProcess statistics valueCensored sampleδs=01≤n≤ns+nshMn,ijs=∑k=0nns−2δn,iskδn,jsk+1n>ns+nshMn,ijs=0Failure sampleδs=11≤n≤nsMn,ijs=∑k=0nns−2δn,iskδn,jsk+1n>nsMn,ijs=0(4)
The set of statisticsYn,is represents whether the phase of the sample s=1,2,…,S is transferred from i=1,2,…,υ to the absorbing state during the nth n=1,2,… shock damage period. The values of Yn,is are shown in Table 5.Table 5
Value of process statisticsYn,is.
Sample typeShock numberProcess statistics valueCensored sampleδs=01≤n≤ns+nshYn,is=δn,isnns−1, satisfies∑n=1ns+nsh∑i=1υYn,is=ns+nsh−1n>ns+nshYn,is=0Failure sampleδs=11≤n≤nsYn,is=δn,isnns−1, satisfies∑n=1ns∑i=1υYn,is=ns−1n>nsYn,is=0LetH=H1,…,HS represent a set of process statistics of all samples, where Hs=Bn,is,Zn,is,Mn,ijs,Yn,is represents the set of process statistics of the sample s=1,2,…,S.


### 3.2. Parameter Estimation Expression

With the above process statistics, the likelihood function under the complete data can be expressed as follows:(5)Lxθx;Xth,N,H,δ=Lxθx;H=∏s=1SLxsθx;Xth,ns,Hs,δs=∏s=1SLxsθx;Hs,wherein, the complete data likelihood function of the censored sample δs=0 is(6)Lxsθx;Hs=∏n=1ns∏i=1υβiBn,is∏n=1ns∏i=1υerii0Zn,is∏n=1ns∏i=1υ∏j=1,i≠jυrij0Mn,ijs∏n=1ns∏i=1υri1Yn,is∏n=ns+1∞∏i=1υβiBn,is∏n=ns+1∞∏i=1υerii0Zn,is∏n=ns+1∞∏i=1υ∏j=1,i≠jυrij0Mn,ijs∏n=ns+1∞∏i=1υri1Yn,is.The complete data likelihood function of the failure sampleδs=1 is(7)Lxsθx;Hs=∏n=1ns∏i=1υβiBn,is∏n=1ns∏i=1υerii0Zn,is∏n=1ns∏i=1υ∏j=1,i≠jυrij0Mn,ijs∏n=1ns−1∏i=1υri1Yn,is.A log-likelihood function can be obtained from equation (5):(8)lnLxθx;H=∑s=1SlnLxsθx;Hs.According to the EM method, the posterior distribution of process statisticsH is(9)fH|θx,Xth,N=Lxθx;H∫HLxθx;HdH.Assuming that the parameter valuesθ0 have been obtained by previous iteration, the expectation of posterior distribution of process statistics H in equation (8) by equation (9) is(10)Eθ0lnLxθx;HsXth,N=∑s=1SEθ0lnLxsθx;HsXth,ns,where the expectation of the log-likelihood function with the censored sample δs=0 data is(11)Eθ0lnLxsθx;HsXth,ns=∑n=1ns∑i=1υEθ0Bn,isXth,nslnβi+∑n=1nsc∑i=1υEθ0Zn,isXth,nsrii0+∑n=1ns∑i=1υ∑j=1,i≠jυEθ0Mn,ijsXth,nslnrij0+∑n=1ns∑i=1υEθ0Yn,isXth,nslnri1+∑n=ns+1∞∑i=1υEθ0Bn,isXth,nslnβi+∑n=ns+1∞∑i=1υEθ0Zn,isXth,nsrii0+∑n=ns+1∞∑i=1υ∑j=1,i≠jυEθ0Mn,ijsXth,nslnrij0+∑n=ns+1∞∑i=1υEθ0Yn,isXth,nslnri1.And, the expectation of the log-likelihood function with the failure sampleδs=1 data is(12)Eθ0lnLxsθx;HsXth,ns0,1=∑n=1ns∑i=1υEθ0Bn,isXth,nslnβi+∑n=1ns∑i=1υEθ0Zn,isXth,nsrii0+∑n=1ns∑i=1υ∑j=1,i≠jυEθ0Mn,ijsXth,nslnrij0+∑n=1ns−1∑i=1υEθ0Yn,isXth,nslnri1.The parametersθx to be estimated satisfy the following constraints:(1)
The initial phase probability satisfies∑i=1υβi=1(2)
For any phasei=1,2,…,υ, ∑j=1υrij0+ri1=0With the MLE method, the parametersθx can be expressed in the following equations:(13)βi=∑s=1S∑n=1nsEθ0Bn,isXth,ns+∑s=1,δs=0S∑n=ns+1∞Eθ0Bn,isXth,ns∑s=1Sns+∑s=1,δs=0S∑n=ns+1∞∑j=1υEθ0Bn,jsXth,ns,i=1,2,…,υ,rij0=∑s=1S∑n=1nsEθ0Mn,ijsXth,ns+∑s=1,δs=0S∑n=ns+1∞Eθ0Mn,ijsXth,ns∑s=1S∑n=1nsEθ0Zn,is|Xth,ns+∑s=1,δs=0S∑n=ns+1∞Zn,is|Xth,ns,i=1,2,…,υ;j=1,2,…,υ;i≠j,ri1=∑s=1S∑n=1ns−1Eθ0Yn,isXth,ns+∑s=1,δs=0S∑n=ns∞Eθ0Yn,isXth,ns∑s=1S∑n=1nsEθ0Zn,isXth,ns+∑s=1,δs=0S∑n=ns+1∞Eθ0Zn,isXth,ns,i=1,2,…,υ,rii0=−∑j=1,i≠jυrij0−ri1,i=1,2,…,υ.The general steps of the EM method areStep E: solving the expectation sets of posterior distribution of the process statistics proposed in Section3.1:(14)Eθ0Bn,isXth,ns,Eθ0Zn,isXth,ns,Eθ0Mn,ijsXth,ns,Eθ0Yn,isXth,ns.Step M: bringing the result obtained in step E into equation (13) to obtain an iteration value of parameters θx.It is not difficult to see that the set of process statistics constructed in Section3.2 is infinite, and expectations of the posterior distribution of the process statistics cannot be obtained at this time, so an approximate method is introduced in the following sections.


### 3.3. Approximate Principle

Assume that the samples is censored after the nsth ns=1,2,… shock, and the probability of failure after another nshth nsh=1,2,… shock is PXns+nsh−1∗<Xth,Xns+nsh∗≥Xth.Then, the survival probability of the sample after thensth shock can be expressed as(15)PXns∗<Xth=∑nsh=1∞PXns+nsh−1∗<Xth,Xns+nsh∗≥Xth=limn⟶∞∑nsh=1nPXns+nsh−1∗<Xth,Xns+nsh∗≥Xth.LetΔn=PXns∗<Xth−∑nsh=1nPXns+nsh−1∗<Xth,Xns+nsh∗≥Xth, then it is not difficult to see that for any determined PH distribution β,R0 and failure thresholdXth, limn⟶∞Δn=0 is constantly satisfied.That is, given any positive real numberR>0, there must exist at least one nsΔ=1,2,… such that for any n>nsΔ, Δn≤R is constantly satisfied.Now, assume thatR is sufficiently small and nsΔ satisfies ΔnsΔ≤R and ΔnsΔ−1>R, then equation (15) can be approximated by a function of finite summation:(16)PXns∗<Xth≈γns+nsΔexpYns+nsΔ0Xthe⇀ns+1,ns+nsΔΣ,where e⇀m,nΣ represents a column vector that the element from m−1υ+1 to nυ is 1 and the others are 0, with corresponding order.Based on the above description, the parametersθx can be approximated as(17)βi≈∑s=1S∑n=1nsEθ0Bn,isXth,ns+∑s=1,δs=0S∑n=ns+1ns+nsΔEθ0Bn,isXth,ns∑s=1Sns+∑s=1,δs=0S∑n=ns+1ns+nsΔ∑j=1υEθ0Bn,jsXth,ns,i=1,2,…,υ,rij0≈∑s=1S∑n=1nsEθ0Mn,ijsXth,ns+∑s=1,δs=0S∑n=ns+1ns+nsΔEθ0Mn,ijsXth,ns∑s=1S∑n=1nsEθ0Zn,isXth,ns+∑s=1,δs=0S∑n=ns+1ns+nsΔEθ0Zn,isXth,ns,i=1,2,…,υ;j=1,2,…,υ;i≠j,ri1≈∑s=1S∑n=1ns−1Eθ0Yn,isXth,ns+∑s=1,δs=0S∑n=nsns+nsΔEθ0Yn,isXth,ns∑s=1S∑n=1nsEθ0Zn,isXth,ns+∑s=1,δs=0S∑n=ns+1ns+nsΔEθ0Zn,isXth,ns,i=1,2,…,υ,rii0≈−∑j=1,i≠jυrij0−ri1,i=1,2,…,υ.


### 3.4. Expectations of Process Statistics with Posterior Distribution

After approximation, the expectations of process statistics with posterior distribution can be obtained referring to the method in [14, 15] by making some improvements. The results are given directly as follows:(1)
Censored sampleδs=0:(1)
The expectations ofBn,is with posterior distribution:① Whenn=1,(18)Eθ0Bn,isXth,ns≈βiθ0e⇀iTexpYnsT0,θ0Xthe⇀ns+1,nsTΣγnsTθ0expYnsT0,θ0Xthe⇀ns+1,nsTΣ,i=1,2,…,υ,wherensT=ns+nsΔ.② Whenn=2,3,…,ns+1,(19)Eθ0Bn,isXth,ns≈∫0Xthγn−1θ0expYn−10,θ0xyn−11,θ0βiθ0e⇀iTexpYnsT−n+10,θ0Xth−xe⇀ns−n+2,nsT−n+1ΣdxγnsTθ0expYnsT0,θ0Xthe⇀ns+1,nsTΣ,i=1,2,…,υ.① Whenn=1,(20)Eθ0Bn,isXth,ns≈∫0Xthγn−1θ0expYn−10,θ0xyn−11,θ0βiθ0e⇀iTexpYnsT−n+10,θ0Xth−xe⇀dxγnsTθ0expYnsT0,θ0Xthe⇀ns+1,nsTΣ,i=1,2,…,υ.(2)
The expectations ofZn,is with posterior distribution:① Whenn=1,2,…,ns,(21)Eθ0Zn,isXth,ns≈∫0XthγnsTθ0expYnsT0,θ0xe⇀n,ie⇀n,iTexpYnsT0,θ0Xth−xe⇀ns+1,nsTΣdxγnsTθ0expYnsT0,θ0Xthe⇀ns+1,nsTΣ,i=1,2,…,υ.③ Whenn=ns+2,ns+3,…,nsT,(22)Eθ0Zn,isXth,ns≈∫0XthγnsTθ0expYnsT0,θ0xe⇀n,ie⇀n,iTexpYnsT0,θ0Xth−xe⇀n,nsTΣdxγnsTθ0expYnsT0,θ0Xthe⇀ns+1,nsTΣ,i=1,2,…,υ.(3)
The expectations ofMn,ijs with posterior distribution:① Whenn=1,2,…,ns,(23)Eθ0Mn,ijsXth,ns≈∫0XthγnsTθ0expYnsT0,θ0xe⇀n,idij0,θ0e⇀n,jTexpYnsT0,θ0Xth−xe⇀ns+1,nsTΣdxγnsTθ0expYnsT0,θ0Xthe⇀ns+1,nsTΣ,i=1,2,…,υ;j=1,2,…,υ;i≠j.② Whenn=ns+1,ns+3,…,nsT,(24)Eθ0Mn,ijsXth,ns≈∫0XthγnsTθ0expYnsT0,θ0xe⇀n,idij0,θ0e⇀n,jTexpYnsT0,θ0Xth−xe⇀n,nsTΣdxγnsTθ0expYnsT0,θ0Xthe⇀ns+1,nsTΣ,i=1,2,…,υ;j=1,2,…,υ;i≠j.(4)
The expectations ofYn,is with posterior distribution:① Whenn=1,2,...,ns,(25)Eθ0Yn,isXth,Nc≈∫0Xthγnθ0expYn0,θ0xe⇀iyi1,θ0γnsT−nθ0expYnsT−n0,θ0Xth−xe⇀ns−n+1,nsT−nΣdxγnsTθ0expYnsT0,θ0Xthe⇀ns+1,nsTΣ,i=1,2,…,υ.② Whenn=ns+1,ns+2,…,nsT,(26)Eθ0Yn,isXth,ns≈∫0Xthγnθ0expYn0,θ0xe⇀iyi1,θ0γnsT−nθ0expYnsT−n0,θ0Xth−xe⇀dxγnsTθ0expYnsT0,θ0Xthe⇀ns+1,nsTΣ,i=1,2,…,υ.(2)
Failure sampleδs=1:(1)
The expectations ofBn,is with posterior distribution:(27)Eθ0Bn,isXth,ns=βiθ0e⇀iTexpYns0,θ0Xthe⇀ns,Σγnsθ0expYns0,θ0Xthe⇀ns,Σ,i=1,2,…,υ.② Whenn=ns+1,ns+2,…,nsT−1,(28)Eθ0Bn,isXth,ns=∫0Xthγn−1θ0expYn−10,θ0xyn−11,θ0βiθ0e⇀iTexpYns−n+10,θ0Xth−xe⇀ns−n+1,Σdxγnsθ0expYns0,θ0Xthe⇀ns,Σ,i=1,2,…,υ.(2)
The expectations ofZn,is with posterior distribution:(29)Eθ0Zn,isXth,ns=∫0Xthγnsθ0expYns0,θ0xe⇀n,ie⇀n,iTexpYns0,θ0Xth−xe⇀ns,Σdxγnsθ0expYns0,θ0Xthe⇀ns,Σ,n=1,2,…,ns,i=1,2,…,υ.(3)
The expectations ofMn,ijs with posterior distribution:(30)Eθ0Mn,ijsXth,ns=∫0Xthγnsθ0expYns0,θ0xe⇀n,idij0,θ0e⇀n,jTexpYns0,θ0Xth−xe⇀ns,Σdxγnsθ0expYns0,θ0Xthe⇀ns,Σ,n=1,2,…,ns;i=1,2,…,υ;j=1,2,…,υj=1,2,…,υ;j=1,2,…,υ;i≠j.(4)
The expectations ofYn,is with posterior distribution:(31)Eθ0Yn,isXth,ns=∫0Xthγnθ0expYn0,θ0xe⇀n,idi1,θ0γns−nθ0expYns−n0,θ0Xth−xe⇀ns−n,Σdxγnsθ0expYns0,θ0Xthe⇀ns,Σ,n=1,2,…,ns−1;i=1,2,…,υ,wherein βθ0=β1θ0...βυθ0, γnθ0=βθ00...01,υn, n=1,2,...,(32)R0θ0=r110,θ0…r1υ0,θ0………rυ10,θ0…rυυ0,θ0,R1θ0=r111,θ0…r1υ1,θ0………rυ11,θ0…rυυ1,θ0,Yn0,θ0=R0θ0R1θ0…00R0θ0…0………R1θ000…R0θ0υn,υn,n=1,2,…,r1θ0=r11,θ0…rυ1,θ0T=−R0θ0e⇀,yn1,θ0=−Yn0,θ0e⇀,n=1,2,…,and e⇀n,i represents a column vector that the n−1υ+ith element is 1 and others are 0, with corresponding dimensions.② Whenn=2,3,...,ns,


### 3.5. Steps of the Method

It can be seen from the previous description that it is necessary to determine the value of the parameterR first. After that, the corresponding number of shocks nT according to the parameter values θ0 obtained by the current iteration can be determined, so that the survival probability that the number of additional shocks exceeds nsΔ=nT−ns can be suffered for most censored samples δs=0 is less than R.If the number of shocks of a censored sample satisfiesns≥nT, it is not difficult to see that the survival probability of such samples satisfies(33)PXnT∗<Xth=∑n=nT+1nsγnexpYn0Xthe⇀n,Σ+PXns∗<Xth≤R.Therefore, it can be assumed that if the sample suffers one more shock, then it will inevitably fail, which meansnsT=ns+1. However, for the other censored samples, the original value nT remains unchanged.Taking the process of parameter estimation ofθx in damage distribution β,R0 as an example, the specific steps of the method based on the shock model and interval sample data described in Section 2 are as follows:Step 1: given the initial valueθ0 to the distribution β,R0 and ε>0 as the end condition of iteration, determine the approximate threshold R>0Step 2: calculate the correspondingnT by θ0 and RStep 3: calculate the value of the log-likelihood functionlnLxθ0;Xth,N,δ according to θ0Step 4: use the E-step according to the distribution parameterθn−1 obtained in the n−1th iteration and the formulas in Section 3.4 to obtain expectation Eθn−1HXth,N of the process statistic set H with posterior distributionNote that during the processes, for any censored sampleδs=0,①Whenns<nT, nT does not change②
Whenns≥nT, there is nsT=ns+1Step 5: bringEθn−1HXth,N into equation (17) to obtain a new iteration value θnStep 6: calculate the correspondingnT by θn and RStep 7: calculate the value of the log-likelihood functionlnLxθn;Xth,N,δ of the sample according to θnStep 8: compare the value ofΔ=lnLxθn;Xth,N,δ−lnLxθn−1;Xth,N,δ with ε:① WhenΔ>ε, let n=n+1 and repeat steps 4–7② WhenΔ≤ε, end the iteration and return the estimated value θ⌢x=θnThe parameterθt of time interval distribution α,D0 can also be obtained based on the above steps. For censored samples δs=0, note that since the last shock time Tns satisfies Tns<tc, while for failure samples δs=1, the cumulative damage Xns∗ caused by shocks at the time of failure satisfies Xns∗≥Xth, it is necessary to consider that the number of shocks should be ns+1 when estimating θt, which is similar to ns when estimating θx.


## 3.1. Process Statistics Construction

Four variables and an indicator function are defined in the following:(1)
dn,ks represents the phase of the sample s=1,2,...,S after the kth phase transition in the nth n=1,2,... shock damage period;(2)
bn,ks represents the damage caused by phase dn,ks;(3)
nns represents the number of phase transitions of the sample s=1,2,…,S during the nth n=1,2,… shock damage period, including the times when the damage period changes;(4)
Nsh represents the number of additional shocks that may cause failure if the sample sδs=0 is censored with shock number ns. It can be seen that Nsh>0 is a random variable and nsh represents its specific sample value;(5)
Indicator functionδn,isk=1dn,ks=i0dn,ks≠i represents whether the phase of the sample s=1,2,…,S is i=1,2,…,υ after the kth phase transition in the nth n=1,2,… shock damage period.With the above variables and function, four process statistics can be constructed as follows:(1)
The set of statisticsBn,is represents whether the initial phase of the sample s=1,2,…,S is i=1,2,…,υ during the nth n=1,2,… shock damage period. The values of Bn,is are shown in Table 2.Table 2
Value of process statisticsBn,is.
Sample typeShock numberProcess statistics valueCensored sampleδs=01≤n≤ns+nshBn,is=δn,is0, satisfies ∑i=1υBn,is=1n>ns+nshBn,is=0Failure sampleδs=11≤n≤nsBn,is=δn,is0, satisfies ∑i=1υBn,is=1n>nsBn,is=0(2)
The set of statisticsZn,is represents the damage caused by the sample s=1,2,…,S in phase i=1,2,…,υ during the nth n=1,2,… shock damage period. The values of Zn,is are shown in Table 3.Table 3
Value of process statisticsZn,is.
Sample typeShock numberProcess statistics valueCensored sampleδs=01≤n≤ns+nshZn,is=∑k=0nns−1δn,iskbn,ks, satisfies Xth=∑n=1ns+nsh∑i=1υZn,isn>ns+nshZn,is=0Failure sampleδs=11≤n≤nsZn,is=∑k=0nns−1δn,iskn,ks, satisfies Xth=∑n=1ns∑i=1υZn,isn>nsZn,is=0(3)
The set of statisticsMn,ijs represents the total number of phase transition times of the sample s=1,2,…,S when the phase is transferred from i=1,2,…,υ to j=1,2,…,υi≠j during the nth n=1,2,… shock damage period, which does not include the times of the phase transition when the damage period changes. The values of Mn,ijs are shown in Table 4.Table 4
Value of process statisticsMn,ijs.
Sample typeShock numberProcess statistics valueCensored sampleδs=01≤n≤ns+nshMn,ijs=∑k=0nns−2δn,iskδn,jsk+1n>ns+nshMn,ijs=0Failure sampleδs=11≤n≤nsMn,ijs=∑k=0nns−2δn,iskδn,jsk+1n>nsMn,ijs=0(4)
The set of statisticsYn,is represents whether the phase of the sample s=1,2,…,S is transferred from i=1,2,…,υ to the absorbing state during the nth n=1,2,… shock damage period. The values of Yn,is are shown in Table 5.Table 5
Value of process statisticsYn,is.
Sample typeShock numberProcess statistics valueCensored sampleδs=01≤n≤ns+nshYn,is=δn,isnns−1, satisfies∑n=1ns+nsh∑i=1υYn,is=ns+nsh−1n>ns+nshYn,is=0Failure sampleδs=11≤n≤nsYn,is=δn,isnns−1, satisfies∑n=1ns∑i=1υYn,is=ns−1n>nsYn,is=0LetH=H1,…,HS represent a set of process statistics of all samples, where Hs=Bn,is,Zn,is,Mn,ijs,Yn,is represents the set of process statistics of the sample s=1,2,…,S.


## 3.2. Parameter Estimation Expression

With the above process statistics, the likelihood function under the complete data can be expressed as follows:(5)Lxθx;Xth,N,H,δ=Lxθx;H=∏s=1SLxsθx;Xth,ns,Hs,δs=∏s=1SLxsθx;Hs,wherein, the complete data likelihood function of the censored sample δs=0 is(6)Lxsθx;Hs=∏n=1ns∏i=1υβiBn,is∏n=1ns∏i=1υerii0Zn,is∏n=1ns∏i=1υ∏j=1,i≠jυrij0Mn,ijs∏n=1ns∏i=1υri1Yn,is∏n=ns+1∞∏i=1υβiBn,is∏n=ns+1∞∏i=1υerii0Zn,is∏n=ns+1∞∏i=1υ∏j=1,i≠jυrij0Mn,ijs∏n=ns+1∞∏i=1υri1Yn,is.The complete data likelihood function of the failure sampleδs=1 is(7)Lxsθx;Hs=∏n=1ns∏i=1υβiBn,is∏n=1ns∏i=1υerii0Zn,is∏n=1ns∏i=1υ∏j=1,i≠jυrij0Mn,ijs∏n=1ns−1∏i=1υri1Yn,is.A log-likelihood function can be obtained from equation (5):(8)lnLxθx;H=∑s=1SlnLxsθx;Hs.According to the EM method, the posterior distribution of process statisticsH is(9)fH|θx,Xth,N=Lxθx;H∫HLxθx;HdH.Assuming that the parameter valuesθ0 have been obtained by previous iteration, the expectation of posterior distribution of process statistics H in equation (8) by equation (9) is(10)Eθ0lnLxθx;HsXth,N=∑s=1SEθ0lnLxsθx;HsXth,ns,where the expectation of the log-likelihood function with the censored sample δs=0 data is(11)Eθ0lnLxsθx;HsXth,ns=∑n=1ns∑i=1υEθ0Bn,isXth,nslnβi+∑n=1nsc∑i=1υEθ0Zn,isXth,nsrii0+∑n=1ns∑i=1υ∑j=1,i≠jυEθ0Mn,ijsXth,nslnrij0+∑n=1ns∑i=1υEθ0Yn,isXth,nslnri1+∑n=ns+1∞∑i=1υEθ0Bn,isXth,nslnβi+∑n=ns+1∞∑i=1υEθ0Zn,isXth,nsrii0+∑n=ns+1∞∑i=1υ∑j=1,i≠jυEθ0Mn,ijsXth,nslnrij0+∑n=ns+1∞∑i=1υEθ0Yn,isXth,nslnri1.And, the expectation of the log-likelihood function with the failure sampleδs=1 data is(12)Eθ0lnLxsθx;HsXth,ns0,1=∑n=1ns∑i=1υEθ0Bn,isXth,nslnβi+∑n=1ns∑i=1υEθ0Zn,isXth,nsrii0+∑n=1ns∑i=1υ∑j=1,i≠jυEθ0Mn,ijsXth,nslnrij0+∑n=1ns−1∑i=1υEθ0Yn,isXth,nslnri1.The parametersθx to be estimated satisfy the following constraints:(1)
The initial phase probability satisfies∑i=1υβi=1(2)
For any phasei=1,2,…,υ, ∑j=1υrij0+ri1=0With the MLE method, the parametersθx can be expressed in the following equations:(13)βi=∑s=1S∑n=1nsEθ0Bn,isXth,ns+∑s=1,δs=0S∑n=ns+1∞Eθ0Bn,isXth,ns∑s=1Sns+∑s=1,δs=0S∑n=ns+1∞∑j=1υEθ0Bn,jsXth,ns,i=1,2,…,υ,rij0=∑s=1S∑n=1nsEθ0Mn,ijsXth,ns+∑s=1,δs=0S∑n=ns+1∞Eθ0Mn,ijsXth,ns∑s=1S∑n=1nsEθ0Zn,is|Xth,ns+∑s=1,δs=0S∑n=ns+1∞Zn,is|Xth,ns,i=1,2,…,υ;j=1,2,…,υ;i≠j,ri1=∑s=1S∑n=1ns−1Eθ0Yn,isXth,ns+∑s=1,δs=0S∑n=ns∞Eθ0Yn,isXth,ns∑s=1S∑n=1nsEθ0Zn,isXth,ns+∑s=1,δs=0S∑n=ns+1∞Eθ0Zn,isXth,ns,i=1,2,…,υ,rii0=−∑j=1,i≠jυrij0−ri1,i=1,2,…,υ.The general steps of the EM method areStep E: solving the expectation sets of posterior distribution of the process statistics proposed in Section3.1:(14)Eθ0Bn,isXth,ns,Eθ0Zn,isXth,ns,Eθ0Mn,ijsXth,ns,Eθ0Yn,isXth,ns.Step M: bringing the result obtained in step E into equation (13) to obtain an iteration value of parameters θx.It is not difficult to see that the set of process statistics constructed in Section3.2 is infinite, and expectations of the posterior distribution of the process statistics cannot be obtained at this time, so an approximate method is introduced in the following sections.


## 3.3. Approximate Principle

Assume that the samples is censored after the nsth ns=1,2,… shock, and the probability of failure after another nshth nsh=1,2,… shock is PXns+nsh−1∗<Xth,Xns+nsh∗≥Xth.Then, the survival probability of the sample after thensth shock can be expressed as(15)PXns∗<Xth=∑nsh=1∞PXns+nsh−1∗<Xth,Xns+nsh∗≥Xth=limn⟶∞∑nsh=1nPXns+nsh−1∗<Xth,Xns+nsh∗≥Xth.LetΔn=PXns∗<Xth−∑nsh=1nPXns+nsh−1∗<Xth,Xns+nsh∗≥Xth, then it is not difficult to see that for any determined PH distribution β,R0 and failure thresholdXth, limn⟶∞Δn=0 is constantly satisfied.That is, given any positive real numberR>0, there must exist at least one nsΔ=1,2,… such that for any n>nsΔ, Δn≤R is constantly satisfied.Now, assume thatR is sufficiently small and nsΔ satisfies ΔnsΔ≤R and ΔnsΔ−1>R, then equation (15) can be approximated by a function of finite summation:(16)PXns∗<Xth≈γns+nsΔexpYns+nsΔ0Xthe⇀ns+1,ns+nsΔΣ,where e⇀m,nΣ represents a column vector that the element from m−1υ+1 to nυ is 1 and the others are 0, with corresponding order.Based on the above description, the parametersθx can be approximated as(17)βi≈∑s=1S∑n=1nsEθ0Bn,isXth,ns+∑s=1,δs=0S∑n=ns+1ns+nsΔEθ0Bn,isXth,ns∑s=1Sns+∑s=1,δs=0S∑n=ns+1ns+nsΔ∑j=1υEθ0Bn,jsXth,ns,i=1,2,…,υ,rij0≈∑s=1S∑n=1nsEθ0Mn,ijsXth,ns+∑s=1,δs=0S∑n=ns+1ns+nsΔEθ0Mn,ijsXth,ns∑s=1S∑n=1nsEθ0Zn,isXth,ns+∑s=1,δs=0S∑n=ns+1ns+nsΔEθ0Zn,isXth,ns,i=1,2,…,υ;j=1,2,…,υ;i≠j,ri1≈∑s=1S∑n=1ns−1Eθ0Yn,isXth,ns+∑s=1,δs=0S∑n=nsns+nsΔEθ0Yn,isXth,ns∑s=1S∑n=1nsEθ0Zn,isXth,ns+∑s=1,δs=0S∑n=ns+1ns+nsΔEθ0Zn,isXth,ns,i=1,2,…,υ,rii0≈−∑j=1,i≠jυrij0−ri1,i=1,2,…,υ.


## 3.4. Expectations of Process Statistics with Posterior Distribution

After approximation, the expectations of process statistics with posterior distribution can be obtained referring to the method in [14, 15] by making some improvements. The results are given directly as follows:(1)
Censored sampleδs=0:(1)
The expectations ofBn,is with posterior distribution:① Whenn=1,(18)Eθ0Bn,isXth,ns≈βiθ0e⇀iTexpYnsT0,θ0Xthe⇀ns+1,nsTΣγnsTθ0expYnsT0,θ0Xthe⇀ns+1,nsTΣ,i=1,2,…,υ,wherensT=ns+nsΔ.② Whenn=2,3,…,ns+1,(19)Eθ0Bn,isXth,ns≈∫0Xthγn−1θ0expYn−10,θ0xyn−11,θ0βiθ0e⇀iTexpYnsT−n+10,θ0Xth−xe⇀ns−n+2,nsT−n+1ΣdxγnsTθ0expYnsT0,θ0Xthe⇀ns+1,nsTΣ,i=1,2,…,υ.① Whenn=1,(20)Eθ0Bn,isXth,ns≈∫0Xthγn−1θ0expYn−10,θ0xyn−11,θ0βiθ0e⇀iTexpYnsT−n+10,θ0Xth−xe⇀dxγnsTθ0expYnsT0,θ0Xthe⇀ns+1,nsTΣ,i=1,2,…,υ.(2)
The expectations ofZn,is with posterior distribution:① Whenn=1,2,…,ns,(21)Eθ0Zn,isXth,ns≈∫0XthγnsTθ0expYnsT0,θ0xe⇀n,ie⇀n,iTexpYnsT0,θ0Xth−xe⇀ns+1,nsTΣdxγnsTθ0expYnsT0,θ0Xthe⇀ns+1,nsTΣ,i=1,2,…,υ.③ Whenn=ns+2,ns+3,…,nsT,(22)Eθ0Zn,isXth,ns≈∫0XthγnsTθ0expYnsT0,θ0xe⇀n,ie⇀n,iTexpYnsT0,θ0Xth−xe⇀n,nsTΣdxγnsTθ0expYnsT0,θ0Xthe⇀ns+1,nsTΣ,i=1,2,…,υ.(3)
The expectations ofMn,ijs with posterior distribution:① Whenn=1,2,…,ns,(23)Eθ0Mn,ijsXth,ns≈∫0XthγnsTθ0expYnsT0,θ0xe⇀n,idij0,θ0e⇀n,jTexpYnsT0,θ0Xth−xe⇀ns+1,nsTΣdxγnsTθ0expYnsT0,θ0Xthe⇀ns+1,nsTΣ,i=1,2,…,υ;j=1,2,…,υ;i≠j.② Whenn=ns+1,ns+3,…,nsT,(24)Eθ0Mn,ijsXth,ns≈∫0XthγnsTθ0expYnsT0,θ0xe⇀n,idij0,θ0e⇀n,jTexpYnsT0,θ0Xth−xe⇀n,nsTΣdxγnsTθ0expYnsT0,θ0Xthe⇀ns+1,nsTΣ,i=1,2,…,υ;j=1,2,…,υ;i≠j.(4)
The expectations ofYn,is with posterior distribution:① Whenn=1,2,...,ns,(25)Eθ0Yn,isXth,Nc≈∫0Xthγnθ0expYn0,θ0xe⇀iyi1,θ0γnsT−nθ0expYnsT−n0,θ0Xth−xe⇀ns−n+1,nsT−nΣdxγnsTθ0expYnsT0,θ0Xthe⇀ns+1,nsTΣ,i=1,2,…,υ.② Whenn=ns+1,ns+2,…,nsT,(26)Eθ0Yn,isXth,ns≈∫0Xthγnθ0expYn0,θ0xe⇀iyi1,θ0γnsT−nθ0expYnsT−n0,θ0Xth−xe⇀dxγnsTθ0expYnsT0,θ0Xthe⇀ns+1,nsTΣ,i=1,2,…,υ.(2)
Failure sampleδs=1:(1)
The expectations ofBn,is with posterior distribution:(27)Eθ0Bn,isXth,ns=βiθ0e⇀iTexpYns0,θ0Xthe⇀ns,Σγnsθ0expYns0,θ0Xthe⇀ns,Σ,i=1,2,…,υ.② Whenn=ns+1,ns+2,…,nsT−1,(28)Eθ0Bn,isXth,ns=∫0Xthγn−1θ0expYn−10,θ0xyn−11,θ0βiθ0e⇀iTexpYns−n+10,θ0Xth−xe⇀ns−n+1,Σdxγnsθ0expYns0,θ0Xthe⇀ns,Σ,i=1,2,…,υ.(2)
The expectations ofZn,is with posterior distribution:(29)Eθ0Zn,isXth,ns=∫0Xthγnsθ0expYns0,θ0xe⇀n,ie⇀n,iTexpYns0,θ0Xth−xe⇀ns,Σdxγnsθ0expYns0,θ0Xthe⇀ns,Σ,n=1,2,…,ns,i=1,2,…,υ.(3)
The expectations ofMn,ijs with posterior distribution:(30)Eθ0Mn,ijsXth,ns=∫0Xthγnsθ0expYns0,θ0xe⇀n,idij0,θ0e⇀n,jTexpYns0,θ0Xth−xe⇀ns,Σdxγnsθ0expYns0,θ0Xthe⇀ns,Σ,n=1,2,…,ns;i=1,2,…,υ;j=1,2,…,υj=1,2,…,υ;j=1,2,…,υ;i≠j.(4)
The expectations ofYn,is with posterior distribution:(31)Eθ0Yn,isXth,ns=∫0Xthγnθ0expYn0,θ0xe⇀n,idi1,θ0γns−nθ0expYns−n0,θ0Xth−xe⇀ns−n,Σdxγnsθ0expYns0,θ0Xthe⇀ns,Σ,n=1,2,…,ns−1;i=1,2,…,υ,wherein βθ0=β1θ0...βυθ0, γnθ0=βθ00...01,υn, n=1,2,...,(32)R0θ0=r110,θ0…r1υ0,θ0………rυ10,θ0…rυυ0,θ0,R1θ0=r111,θ0…r1υ1,θ0………rυ11,θ0…rυυ1,θ0,Yn0,θ0=R0θ0R1θ0…00R0θ0…0………R1θ000…R0θ0υn,υn,n=1,2,…,r1θ0=r11,θ0…rυ1,θ0T=−R0θ0e⇀,yn1,θ0=−Yn0,θ0e⇀,n=1,2,…,and e⇀n,i represents a column vector that the n−1υ+ith element is 1 and others are 0, with corresponding dimensions.② Whenn=2,3,...,ns,


## 3.5. Steps of the Method

It can be seen from the previous description that it is necessary to determine the value of the parameterR first. After that, the corresponding number of shocks nT according to the parameter values θ0 obtained by the current iteration can be determined, so that the survival probability that the number of additional shocks exceeds nsΔ=nT−ns can be suffered for most censored samples δs=0 is less than R.If the number of shocks of a censored sample satisfiesns≥nT, it is not difficult to see that the survival probability of such samples satisfies(33)PXnT∗<Xth=∑n=nT+1nsγnexpYn0Xthe⇀n,Σ+PXns∗<Xth≤R.Therefore, it can be assumed that if the sample suffers one more shock, then it will inevitably fail, which meansnsT=ns+1. However, for the other censored samples, the original value nT remains unchanged.Taking the process of parameter estimation ofθx in damage distribution β,R0 as an example, the specific steps of the method based on the shock model and interval sample data described in Section 2 are as follows:Step 1: given the initial valueθ0 to the distribution β,R0 and ε>0 as the end condition of iteration, determine the approximate threshold R>0Step 2: calculate the correspondingnT by θ0 and RStep 3: calculate the value of the log-likelihood functionlnLxθ0;Xth,N,δ according to θ0Step 4: use the E-step according to the distribution parameterθn−1 obtained in the n−1th iteration and the formulas in Section 3.4 to obtain expectation Eθn−1HXth,N of the process statistic set H with posterior distributionNote that during the processes, for any censored sampleδs=0,①Whenns<nT, nT does not change②
Whenns≥nT, there is nsT=ns+1Step 5: bringEθn−1HXth,N into equation (17) to obtain a new iteration value θnStep 6: calculate the correspondingnT by θn and RStep 7: calculate the value of the log-likelihood functionlnLxθn;Xth,N,δ of the sample according to θnStep 8: compare the value ofΔ=lnLxθn;Xth,N,δ−lnLxθn−1;Xth,N,δ with ε:① WhenΔ>ε, let n=n+1 and repeat steps 4–7② WhenΔ≤ε, end the iteration and return the estimated value θ⌢x=θnThe parameterθt of time interval distribution α,D0 can also be obtained based on the above steps. For censored samples δs=0, note that since the last shock time Tns satisfies Tns<tc, while for failure samples δs=1, the cumulative damage Xns∗ caused by shocks at the time of failure satisfies Xns∗≥Xth, it is necessary to consider that the number of shocks should be ns+1 when estimating θt, which is similar to ns when estimating θx.


## 4. A Case

In this section, the performance of the parameter estimation method is tested by a case with only interval data samples. The case includes three types of nonexponential distributions (gamma distribution, log-normal distribution, and Weibull distribution) which are commonly used in shock model constructions to generate sample data. The proposed method is used to estimate the parameters of the corresponding PH distributions, and the results are compared with the above three distributions to analyze its performance.


### 4.1. Model Assumption

(1)
Assume that the sample size isS=100, the failure threshold is Xth=30, and the censored time is tc=70. Sample data can only be obtained at the censored time.(2)
Assume that the end condition of iteration isε=0.0001 and the approximate threshold is R=0.0001.(3)
There are two shock models:(1)
Model 1:① The interval of shock time satisfies gamma distribution with the shape parameter of 3 and scale parameter of 2② The one shock damage satisfies log-normal distribution with a log-mean value of 0.5 and log-variance value of 1(3)
Model 2:① The interval of shock time satisfies gamma distribution with the shape parameter of 3 and scale parameter of 2② The one shock damage satisfies Weibull distribution with the shape parameter of 3 and scale parameter of 2


### 4.2. Parameter Estimation and Result Analysis

Based on the method, the initial value, including type and order of PH distribution, should be given first. How to select type and order of a PH distribution can be referred to [16], and each type of PH distribution is described in detail in [17]. The estimation results of parameters of the time interval distributions are only shown in the first model since the distributions in both models are the same.


#### 4.2.1. Parameter Estimation of Time Interval Distribution

The distribution of the time interval is fitted by Erlang-3 distribution with 3-order. The initial value of the PH distribution is as follows:(34)α=100,D0=−3300−3300−3,d1=003.The simulation time is about 5 h, and the parameter estimation results based on the simulation data are as follows:(35)α=100,D0=−0.51030.510300−0.50920.509200−0.5087,d1=000.5087.The estimated PDF of the PH distribution and the original gamma distribution are both shown in Figure2, and the related statistics are shown in Table 6.Figure 2
PDF of the time interval (gamma).Table 6
Statistics of the distribution of the time interval (gamma).
StatisticsSampleEstimationDeviationPercentageMean65.8893−0.1107−1.8451Variance1211.5613−0.4387−3.6561It can be seen from the results shown in Figure2 and Table 6 that the estimated PH distribution fits well with the original gamma distribution based on the interval sample data. When the shape parameter of gamma distribution is an integer, the random variable can be regarded as the sum of the same number of iid random variables from exponential distribution with the parameter equal to the scale parameter of the gamma distribution, so the initial distribution form can be easily determined as Erlang-3 distribution in this case. At this time, the parameter estimation of the PH distribution is an unbiased estimation, and the estimated errors are entirely derived from the randomness of the samples.


#### 4.2.2. Parameter Estimation of Damage Distribution (Log-Normal Damage)

The distribution of the shock damage is fitted by Coxian distribution with 3-order. The initial value of the PH distribution is given as follows:(36)β=100,R0=−3200−2100−1,r1=111.The simulation time is about 14 h, and the parameter estimation results based on the simulation data are as follows:(37)β=100,R0=−3.00173.001700−0.57680.031200−0.1471,r1=00.54560.1471.The estimated PDF of the PH distribution and the original log-normal distribution are both shown in Figure3, and the related statistics are shown in Table 7.Figure 3
PDF of the one shock damage (log-normal).Table 7
Statistics of the distribution of the shock damage (log-normal).
StatisticsSampleEstimationDeviationPercentageMean2.71832.43450.283810.4403Variance12.696510.98041.716113.5163It can be seen from Figure3 that the fitting result of the estimated PH distribution is not much different from the original log-normal distribution, but it also can be seen from Table 7 that the deviation of the statistics is larger than the deviation in Table 6. This is because the parameter estimation result of the PH distribution is affected by the randomness of the sample, the initial type of the PH distribution, and the initial value of the parameter.


#### 4.2.3. Parameter Estimation of Damage Distribution (Weibull Damage)

The distribution of the shock damage is fitted by Coxian distribution with 4-order, and the initial value of the PH distribution is given as follows:(38)α=1000,D0=−44000−42000−42000−4,d1=0224.The simulation time is about 17 h, and the parameter estimation results based on the simulation data are as follows:(39)α=1000,D0=−1.41211.4121000−1.41201.2298000−1.30861.3086000−1.3072,d1=00.182201.3072.The estimated PDF of the PH distribution and the original Weibull distribution are both shown in Figure4, and the relevant statistics are shown in Table 8.Figure 4
PDF of the one shock damage (Weibull).Table 8
Statistics of the distribution of shock damage (Weibull).
StatisticsSampleEstimationDeviationPercentageMean2.65872.74820.08953.3682Variance1.93142.28420.352818.2644It can be found that if the initial value of the first phase in the absorption rate vector is not set to 0, the PDF of the PH distribution obtained by parameter estimation at point 0 in abscissa is not 0. The reason for this is that the EM method always obtains the iterative parameter value which can maximize the likelihood function that means the PH distribution obtained under the nonzero absorption rate of the first phase has a larger likelihood value, but the result does not match the actual situation. Therefore, the initial value of the absorption rate of the first phase should be set to 0.It can be seen from the case that under the interval data, the PH distributions based on the estimation results can fit well with the original distributions mentioned above. Therefore, through the approximate method proposed in this paper, the PH distribution can be used for shock model construction under only interval sample data, which not only improves the integral problem cause by the original distributions, but also can obtain the distribution form with a higher degree of fitting.During the process of the above case, it is found that, in the early stage of the iteration, the value of the log-likelihood function(40)∑s=1,δs=0Sln∫0XthγnsexpYns0xyns1dx,∑s=1,δs=1Sln∫0tcμnsexpMns0tm1ndt,which is calculated with the relative parameter iterative value θt obtained by the interval failure data and θx obtained by the censored sample data, is very small; that is, the difference between the two partial log-likelihood functions through two adjacent iterative loops does not affect the end condition of the iteration. However, it can be seen from equation (17) that performing parameter estimation with these two types of samples needs to calculate a large number of posterior distribution expectation of process statistics. In order to improve the efficiency of the proposed method, some improvements are given here with the estimation of parameters θx as an example.When performing the estimation of parametersθx, a threshold εx>0 is set for the censored samples which satisfies εx≈ε.When∑s=1,δs=0Sln∫0XthγnsexpYns0xyns1dx≥εx, the censored sample data are considered in the parameter estimation.When∑s=1,δs=0Sln∫0XthγnsexpYns0xyns1dx<εx, the parameters are estimated only by the interval failure data.The estimation of parametersθt in distribution α,D0 is similar when dealing with interval failure samples. This improvement is only applicable to the case where the sample set contains both censored samples δs=0 and interval failure samples δs=1.


## 4.1. Model Assumption

(1)
Assume that the sample size isS=100, the failure threshold is Xth=30, and the censored time is tc=70. Sample data can only be obtained at the censored time.(2)
Assume that the end condition of iteration isε=0.0001 and the approximate threshold is R=0.0001.(3)
There are two shock models:(1)
Model 1:① The interval of shock time satisfies gamma distribution with the shape parameter of 3 and scale parameter of 2② The one shock damage satisfies log-normal distribution with a log-mean value of 0.5 and log-variance value of 1(3)
Model 2:① The interval of shock time satisfies gamma distribution with the shape parameter of 3 and scale parameter of 2② The one shock damage satisfies Weibull distribution with the shape parameter of 3 and scale parameter of 2


## 4.2. Parameter Estimation and Result Analysis

Based on the method, the initial value, including type and order of PH distribution, should be given first. How to select type and order of a PH distribution can be referred to [16], and each type of PH distribution is described in detail in [17]. The estimation results of parameters of the time interval distributions are only shown in the first model since the distributions in both models are the same.


### 4.2.1. Parameter Estimation of Time Interval Distribution

The distribution of the time interval is fitted by Erlang-3 distribution with 3-order. The initial value of the PH distribution is as follows:(34)α=100,D0=−3300−3300−3,d1=003.The simulation time is about 5 h, and the parameter estimation results based on the simulation data are as follows:(35)α=100,D0=−0.51030.510300−0.50920.509200−0.5087,d1=000.5087.The estimated PDF of the PH distribution and the original gamma distribution are both shown in Figure2, and the related statistics are shown in Table 6.Figure 2
PDF of the time interval (gamma).Table 6
Statistics of the distribution of the time interval (gamma).
StatisticsSampleEstimationDeviationPercentageMean65.8893−0.1107−1.8451Variance1211.5613−0.4387−3.6561It can be seen from the results shown in Figure2 and Table 6 that the estimated PH distribution fits well with the original gamma distribution based on the interval sample data. When the shape parameter of gamma distribution is an integer, the random variable can be regarded as the sum of the same number of iid random variables from exponential distribution with the parameter equal to the scale parameter of the gamma distribution, so the initial distribution form can be easily determined as Erlang-3 distribution in this case. At this time, the parameter estimation of the PH distribution is an unbiased estimation, and the estimated errors are entirely derived from the randomness of the samples.


### 4.2.2. Parameter Estimation of Damage Distribution (Log-Normal Damage)

The distribution of the shock damage is fitted by Coxian distribution with 3-order. The initial value of the PH distribution is given as follows:(36)β=100,R0=−3200−2100−1,r1=111.The simulation time is about 14 h, and the parameter estimation results based on the simulation data are as follows:(37)β=100,R0=−3.00173.001700−0.57680.031200−0.1471,r1=00.54560.1471.The estimated PDF of the PH distribution and the original log-normal distribution are both shown in Figure3, and the related statistics are shown in Table 7.Figure 3
PDF of the one shock damage (log-normal).Table 7
Statistics of the distribution of the shock damage (log-normal).
StatisticsSampleEstimationDeviationPercentageMean2.71832.43450.283810.4403Variance12.696510.98041.716113.5163It can be seen from Figure3 that the fitting result of the estimated PH distribution is not much different from the original log-normal distribution, but it also can be seen from Table 7 that the deviation of the statistics is larger than the deviation in Table 6. This is because the parameter estimation result of the PH distribution is affected by the randomness of the sample, the initial type of the PH distribution, and the initial value of the parameter.


### 4.2.3. Parameter Estimation of Damage Distribution (Weibull Damage)

The distribution of the shock damage is fitted by Coxian distribution with 4-order, and the initial value of the PH distribution is given as follows:(38)α=1000,D0=−44000−42000−42000−4,d1=0224.The simulation time is about 17 h, and the parameter estimation results based on the simulation data are as follows:(39)α=1000,D0=−1.41211.4121000−1.41201.2298000−1.30861.3086000−1.3072,d1=00.182201.3072.The estimated PDF of the PH distribution and the original Weibull distribution are both shown in Figure4, and the relevant statistics are shown in Table 8.Figure 4
PDF of the one shock damage (Weibull).Table 8
Statistics of the distribution of shock damage (Weibull).
StatisticsSampleEstimationDeviationPercentageMean2.65872.74820.08953.3682Variance1.93142.28420.352818.2644It can be found that if the initial value of the first phase in the absorption rate vector is not set to 0, the PDF of the PH distribution obtained by parameter estimation at point 0 in abscissa is not 0. The reason for this is that the EM method always obtains the iterative parameter value which can maximize the likelihood function that means the PH distribution obtained under the nonzero absorption rate of the first phase has a larger likelihood value, but the result does not match the actual situation. Therefore, the initial value of the absorption rate of the first phase should be set to 0.It can be seen from the case that under the interval data, the PH distributions based on the estimation results can fit well with the original distributions mentioned above. Therefore, through the approximate method proposed in this paper, the PH distribution can be used for shock model construction under only interval sample data, which not only improves the integral problem cause by the original distributions, but also can obtain the distribution form with a higher degree of fitting.During the process of the above case, it is found that, in the early stage of the iteration, the value of the log-likelihood function(40)∑s=1,δs=0Sln∫0XthγnsexpYns0xyns1dx,∑s=1,δs=1Sln∫0tcμnsexpMns0tm1ndt,which is calculated with the relative parameter iterative value θt obtained by the interval failure data and θx obtained by the censored sample data, is very small; that is, the difference between the two partial log-likelihood functions through two adjacent iterative loops does not affect the end condition of the iteration. However, it can be seen from equation (17) that performing parameter estimation with these two types of samples needs to calculate a large number of posterior distribution expectation of process statistics. In order to improve the efficiency of the proposed method, some improvements are given here with the estimation of parameters θx as an example.When performing the estimation of parametersθx, a threshold εx>0 is set for the censored samples which satisfies εx≈ε.When∑s=1,δs=0Sln∫0XthγnsexpYns0xyns1dx≥εx, the censored sample data are considered in the parameter estimation.When∑s=1,δs=0Sln∫0XthγnsexpYns0xyns1dx<εx, the parameters are estimated only by the interval failure data.The estimation of parametersθt in distribution α,D0 is similar when dealing with interval failure samples. This improvement is only applicable to the case where the sample set contains both censored samples δs=0 and interval failure samples δs=1.


## 4.2.1. Parameter Estimation of Time Interval Distribution

The distribution of the time interval is fitted by Erlang-3 distribution with 3-order. The initial value of the PH distribution is as follows:(34)α=100,D0=−3300−3300−3,d1=003.The simulation time is about 5 h, and the parameter estimation results based on the simulation data are as follows:(35)α=100,D0=−0.51030.510300−0.50920.509200−0.5087,d1=000.5087.The estimated PDF of the PH distribution and the original gamma distribution are both shown in Figure2, and the related statistics are shown in Table 6.Figure 2
PDF of the time interval (gamma).Table 6
Statistics of the distribution of the time interval (gamma).
StatisticsSampleEstimationDeviationPercentageMean65.8893−0.1107−1.8451Variance1211.5613−0.4387−3.6561It can be seen from the results shown in Figure2 and Table 6 that the estimated PH distribution fits well with the original gamma distribution based on the interval sample data. When the shape parameter of gamma distribution is an integer, the random variable can be regarded as the sum of the same number of iid random variables from exponential distribution with the parameter equal to the scale parameter of the gamma distribution, so the initial distribution form can be easily determined as Erlang-3 distribution in this case. At this time, the parameter estimation of the PH distribution is an unbiased estimation, and the estimated errors are entirely derived from the randomness of the samples.


## 4.2.2. Parameter Estimation of Damage Distribution (Log-Normal Damage)

The distribution of the shock damage is fitted by Coxian distribution with 3-order. The initial value of the PH distribution is given as follows:(36)β=100,R0=−3200−2100−1,r1=111.The simulation time is about 14 h, and the parameter estimation results based on the simulation data are as follows:(37)β=100,R0=−3.00173.001700−0.57680.031200−0.1471,r1=00.54560.1471.The estimated PDF of the PH distribution and the original log-normal distribution are both shown in Figure3, and the related statistics are shown in Table 7.Figure 3
PDF of the one shock damage (log-normal).Table 7
Statistics of the distribution of the shock damage (log-normal).
StatisticsSampleEstimationDeviationPercentageMean2.71832.43450.283810.4403Variance12.696510.98041.716113.5163It can be seen from Figure3 that the fitting result of the estimated PH distribution is not much different from the original log-normal distribution, but it also can be seen from Table 7 that the deviation of the statistics is larger than the deviation in Table 6. This is because the parameter estimation result of the PH distribution is affected by the randomness of the sample, the initial type of the PH distribution, and the initial value of the parameter.


## 4.2.3. Parameter Estimation of Damage Distribution (Weibull Damage)

The distribution of the shock damage is fitted by Coxian distribution with 4-order, and the initial value of the PH distribution is given as follows:(38)α=1000,D0=−44000−42000−42000−4,d1=0224.The simulation time is about 17 h, and the parameter estimation results based on the simulation data are as follows:(39)α=1000,D0=−1.41211.4121000−1.41201.2298000−1.30861.3086000−1.3072,d1=00.182201.3072.The estimated PDF of the PH distribution and the original Weibull distribution are both shown in Figure4, and the relevant statistics are shown in Table 8.Figure 4
PDF of the one shock damage (Weibull).Table 8
Statistics of the distribution of shock damage (Weibull).
StatisticsSampleEstimationDeviationPercentageMean2.65872.74820.08953.3682Variance1.93142.28420.352818.2644It can be found that if the initial value of the first phase in the absorption rate vector is not set to 0, the PDF of the PH distribution obtained by parameter estimation at point 0 in abscissa is not 0. The reason for this is that the EM method always obtains the iterative parameter value which can maximize the likelihood function that means the PH distribution obtained under the nonzero absorption rate of the first phase has a larger likelihood value, but the result does not match the actual situation. Therefore, the initial value of the absorption rate of the first phase should be set to 0.It can be seen from the case that under the interval data, the PH distributions based on the estimation results can fit well with the original distributions mentioned above. Therefore, through the approximate method proposed in this paper, the PH distribution can be used for shock model construction under only interval sample data, which not only improves the integral problem cause by the original distributions, but also can obtain the distribution form with a higher degree of fitting.During the process of the above case, it is found that, in the early stage of the iteration, the value of the log-likelihood function(40)∑s=1,δs=0Sln∫0XthγnsexpYns0xyns1dx,∑s=1,δs=1Sln∫0tcμnsexpMns0tm1ndt,which is calculated with the relative parameter iterative value θt obtained by the interval failure data and θx obtained by the censored sample data, is very small; that is, the difference between the two partial log-likelihood functions through two adjacent iterative loops does not affect the end condition of the iteration. However, it can be seen from equation (17) that performing parameter estimation with these two types of samples needs to calculate a large number of posterior distribution expectation of process statistics. In order to improve the efficiency of the proposed method, some improvements are given here with the estimation of parameters θx as an example.When performing the estimation of parametersθx, a threshold εx>0 is set for the censored samples which satisfies εx≈ε.When∑s=1,δs=0Sln∫0XthγnsexpYns0xyns1dx≥εx, the censored sample data are considered in the parameter estimation.When∑s=1,δs=0Sln∫0XthγnsexpYns0xyns1dx<εx, the parameters are estimated only by the interval failure data.The estimation of parametersθt in distribution α,D0 is similar when dealing with interval failure samples. This improvement is only applicable to the case where the sample set contains both censored samples δs=0 and interval failure samples δs=1.


## 5. Conclusion

In this paper, an approximate parameter estimation method based on constructing PH distribution with dynamic orders is proposed to solve the decoupling problems in parameter estimation of the shock model when facing the samples with only interval data. Its performance is tested through a simulation case from which can be seen that the method proposed in this paper can obtain well-fitted parameter estimation results using only interval sample data.There are also limitations in the application of this method. For example, it is only applicable to the shock problem under a single stress source which means the shock time interval and the shock damage need to be subjected to the condition of iid. If there exists multistress sources which means the shock interval and shock damage may satisfy different types of distributions, the methods presented in this paper will face some identifiability problems. Later work will focus on solving these.



---

*Source: 1424105-2020-06-30.xml*