Inference and optimal design for the k-level step-stress accelerated life test based on progressive Type-I interval censored power Rayleigh data

1.   Introduction

To estimate the reliability of highly reliable products under the normally used conditions, the tested units have a long time for failure, and only few failure times are available to make inference about their life model. To overcome this problem, the life testing experiment is conducted using the accelerated life tests (ALTs). This can be done by exposing the tested units to successively higher levels of stress including the effect of temperature, pressure, and voltage. Practically, this procedure will provide the requested lifetime's data within a faster acceptable time and with relatively low costs.

In reliability frame works, two types of ALTs are mainly considered: The simple accelerated life test (SALT) by exposing the tested units to only one single change in stress level or exposing the tested units to more than one higher level of stress with the step-stress accelerated life test (SSALT). Using the ordinary Type-I and Type-II censored data, SALT and SSALT are devoted to estimate the parameters of many life models. Point and interval estimation are obtained for the exponential lifetime model by Balakrishnan and Han ^[1], Balakrishnan and Han ^[2] and Haghighi ^[3]. Lu and Rudy ^[4], McSorley et al. ^[5] and Komori ^[6] have considered inference in the Weibull life model. Chung and Bai ^[7] have explored SSALT for the log-normal distribution. Related researches are conducted also by Ebrahem and Al-Masri ^[8] and Bleed and Hasan ^[9] for the Logistic distribution, Saxena et al. ^[10] for the Rayleigh distribution and Ahmadi et al. ^[11] for the generalized half-normal distribution.

When it is difficult to observe the exact lifetimes by continuous monitoring of the tested units using Type-I and Type-II censoring schemes, one feasible way to conduct the life test is by the intermittent inspection of the tested units in pre-determined time intervals. In this case, the only available data are the numbers of failed units within these intervals. In addition, if for some constraints the life testing experiment requests removing some units during the life test before the final termination point, then we will have what is called the progressively Type-I interval-censored data (PT-ICD). The PT-ICD consists of the numbers of failed and the numbers of the removed units within the inspection intervals. This procedure will decrease the cost and reduce the efforts requested for conducting the life testing experiment.

Based on SSALT conducted under PT-ICD, many studies are focused on either estimation of the model parameters or in determining the optimum the life test plans. Wu et al. ^[12] explored the SSALT for the Exponential life model. Planning of the SSALTs are also considered in ^[13] for the log-normal distribution, in ^[14] for the Weibull distribution and in ^[15] for the Burr type XII distribution. Optimum designs of SSALTs are also implemented in ^[16] for Pareto distribution, in ^[17] for the Weibull Poisson distribution and in ^[18] for the Weibull distribution.

As it has a growing interest in many engineering and medical applications, inference in SSALT under different censoring schemes has become a main topic in research. Examples of such research explorations include Budhiraja and Pradhan ^[19] for estimation of the cost model parameters and Bayesian optimum life testing plans analyzed by Roy and Pradhan ^[20]. Wang ^[21] addressed the group effects in the data analysis of SSALT with random removals. EL-Sagheer and Khder ^[22] have also considered estimation in a k-stage partially SSALT when the life model is the generalized Pareto distribution. Parameter estimation in partially SSALT under different types of censored data was also justified by Kamal ^[23].

Based on a variety of progressive censored data, inferences of SSALT are also conducted in other lifetime models by Zhuang et al. ^[24], Alotaibi et al. ^[25], Bai et al. ^[26], Almarashi ^[27], Nassar and Elshahhat ^[28] and Alam et al. ^[29].

The power Rayleigh distribution (PRD), introduced by Bhat and Ahmad ^[30], is a new extension of the standard Rayleigh distribution. To date, only a few researchers have considered this distribution for modeling survival and reliability data. To utilize this distribution in a broader range of engineering and medical applications as a lifetime model, we assume that the SSALT is conducted under PT-ICD when the life model is the PRD. This paper aims to make inferences about the parameters of the derived model. To increase the efficiency of the estimation procedures, we also approach finding both the optimum inspection times and the optimum stress levels. Unlike most of the previous studies that concentrate in equidistant inspection intervals, this work considers an arbitrary selection of time cut points.

The rest of this paper is organized as follows: The model and the test procedures are introduced in the next section. The maximum likelihood estimators (MLEs) of the derived model parameters are obtained in Section 3. The Fisher information matrix is obtained in Section 4. Asymptotic and bootstrap confidence intervals of the model parameters are obtained in Section 5. An extensive simulation study is conducted in Section 6 to evaluate the performance of the obtained estimates. Two optimality criteria are employed in Section 7 to obtain both the optimum inspection times and the optimum stress levels. The estimation procedures are investigated in a real dataset in Section 8. Finally, Section 9 explores the general conclusions of this work.

2.   Model and Test Procedures

2.1.   Power Rayleigh distribution

Using the power transformation method, Bhat and Ahmad ^[30] have derived a new generalization of the one parameter Rayleigh distribution called the PRD with probability density (PDF), cumulative distribution function (CDF) and reliability functions given respectively by

The hazard rate function (HRF) of the PRD is given by

which is constant when the shape parameter is $\beta = 0.5$, increasing if $\beta < 0.5$ and decreasing if $\beta > 0.5$.

This flexibility in the hazard function indicated that PRD can be effectively used in place of the traditional Weibull, Gamma, log-normal, and other distributions to model life data.

Additionally, PRD can be used as a general case of the standard Rayleigh distribution, which is employed in modeling the lifetimes of resisters, transformers, and capacitors (see Ateeq et al. ^[31] and Kundu and Raqab ^[32]). The PDF and the HRF of the PRD with different combinations of ($\theta$, $\beta$) are presented in Figure 1.

Recent inference of the PRD in the analysis of survival data is explored by Tolba et al. ^[33] and Migdadi et al. ^[34].

2.2.   Assumptions

Consider the life-testing experiment in which ${n}$ units are placed in a ${k}$-level SSALT based on progressively Type-I interval censored data. The units are run at stress level ${S_{1}}$ until the time ${\tau_{1}}$, and within the interval $({0}, {\tau_{1}}]$ the number ${n_1}$ of failed units and the number ${r_1}$ of the removed units are recorded. Now, starting at the time ${\tau_1}$, the remaining $({n-n_1-r_1})$ survival units are subject to a higher level of stress ${S_2}$ and run until the time ${\tau_2}$, at which point the number ${n_2}$ of failed units and the number ${r_2}$ of the removed units within the interval $({\tau_1, \tau_2}]$ are recorded and so on. The model assumptions are described in the following constraints:

1. There are successively $k$ levels of higher stresses $S_1 < S_2 < \ldots < S_k$ and $k-1$ predetermined inspection times: $\tau_1 < \tau_2 < \ldots < \tau_{k-1} < \tau_k = \infty$, where the normal used stress is assumed to be $S_0 = 0$.

2. The termination test time is the point $\tau_{k-1}$ at which the remaining tested units are either failed or removed from the test. This implies the inspection intervals are: $\left(0, \tau_1\right], \left(\tau_1, \tau_2\right], \ldots, \left(\tau_{k-1}, \infty\right]$.

3. The probabilities of the removed units $r_j$, $\pi_j$, $j = 1, 2, \ldots, k-1$ are randomly determined by the experimenter with the constraint $\pi_k = 1$.

4. The lifetimes of the n tested units are identically distributed according to the PRD indexed by the scale parameter $\theta_j$, which has at the level stress $S_{j}$, $j = 1, 2, \ldots, k$ the inverse power function given by $\theta_j = c S^p_j, c > 0$, where the parameter $c$ is the constant of proportionality and the parameter $p > 0$ is the power of the applied stress.

5. The shape parameter $\beta$ of the PRD is assumed to be independent of the stress levels.

Figure 2 presents the k-Level SSAT based on PT-ICD.

2.3.   The cumulative exposure model

From the previous assumptions and since failures have occurred according to the cumulative exposure model, this implies, as shown by Nelson ^[35], that the CDF of the lifetimes of the tested units is given by

where $u_i$ is the solution to the equation:

From Eq (2.2), the general solution of $u_j$ is given by:

Now, substituting for $\theta_j = cS_j^p$, $j = 1, 2, \ldots, k$ in Eq (2.5) implies

Therefore, The CDF and the PDF of the lifetimes of the tested units are given, respectively, by:

3.   Maximum likelihood estimation

Based on PT-ICD and following Wu et al. ^[36] and Aly ^[37], the likelihood function for the parameters $(c, p, \beta)$ has the following form:

Substituting for $F(\tau_j), j = 1, 2, \ldots, k$, from equation (2.1), the likelihood function can be written as:

where:

This implies that the log likelihood function is given by:

then, $\ln(L(c, p, \beta))$ can be simplified as:

Taking the first partial derivatives of Eq (3.3) with respect to $c$, $p$ and $\beta$ and then equating each to zero, we have the following normal equations:

where

The system of nonlinear Eqs (3.4)–(3.6) cannot be solved analytically. Therefore, the Newton-Raphson numerical iteration method can be employed to approximate the MLEs of the unknown parameters ($c$, $p$, $\beta$) denoted as ($\hat{c}$, $\hat{p}$, $\hat{\beta}$).

4.   The Fisher information matrix

To quantify the information contained in the data obtained from a life test plan in the long run, the Fisher information matrix is frequently used. The Fisher information matrix is defined by the negative expectation of the second derivatives of the log-likelihood function as:

To find the elements of $\mathcal{F}$, we have to find the expectations of $n_j$, $r_j$, $j = 1, 2, \ldots, k$ and the second derivatives of the log-likelihood function.

Based on the multinomial distribution, the numbers of failures are distributed as:

where $P_j^\ast$ is the conditional probability of failure in the interval ($\tau_{j-1}$, $\tau_j$) given of survival at $\tau_{j-1}$ given by

where $W_j = e^\frac{-(\tau_j-\tau_{j-1+}\sum_{i = 1}^{j-1}{\left(\frac{S_j}{S_i}\right)^\frac{p}{\beta}\tau_i})^{2\beta}}{2c^2S_j^{2p}}, j = 2, \ldots, k$.

This implies that the numbers of removals are distributed as:

Therefore,

By induction,

Since $\ln{\left(L\left(c, p, \beta\right)\right)} = \ln{(L}_1\left(c, p, \beta\right))+ \ln{(L}_2\left(c, p, \beta\right))$, we set

This implies

Hence, $\frac{\partial^2L_1}{\partial p^2}$, $\frac{\partial^2L_1}{\partial\beta^2}$ are obtained by substituting $p$, $\beta$ in place of $c$ in equation (4.6). To find the second derivatives of $U$ with respect to ($c$, $p$, $\beta$), see Appendix A.1, and to find the second derivatives of $\ln{(L}_2\left(c, p, \beta\right))$ with respect to $c$, $p$, $\beta$, set:

This implies

and therefore,

where $\frac{\partial^2L_2}{\partial p^2}$, $\frac{\partial^2L_2}{\partial\beta^2}$ can be obtained by substituting $p$, $\beta$ in place of $c$ in equation (4.7), respectively. Also, $\frac{\partial^2L_2}{\partial c\partial\beta}$, $\frac{\partial^2L_2}{\partial p\partial\beta}$ can also be obtained by replacing $\beta$ in place of $p$ and in place of c in equation (4.8), respectively.

To find the second derivatives of $P_j$, $j = 2, \ldots, k$ and $f_j$, $j = 1, \ldots, k$ with respect to $c$, $p$ and $\beta$, see Appendix A.2.

Now, substituting for $E\left(n_j\right)$ and $E\left(r_j\right)$, $j = 1, \ldots, k$ from the Eqs (4.4)–(4.5) and the corresponding second derivatives of the log-likelihood function from the Eqs (4.6)–(4.8) in Eq (4.1), we will have the elements of the Fisher information matrix ($\mathcal{F}$).

5.   Confidence Intervals

Another desired inference is to obtain interval estimates for the model parameters. In this section, two methods are employed for constructing confidence intervals (CIs) for the parameters $c$, $p$ and $\beta$. The first method is based on the asymptotic normality of the MLEs, while the second method is based on the bootstrap resampling technique.

5.1.   Asymptotic confidence intervals

The estimated variance-covariance matrix $\hat{V}$ for the parameters ($c$, $p$, $\beta$) is given by:

where $\mathcal{F}(\hat{c}, \hat{p}, \hat{\beta})$ is the Fisher information matrix given the estimated values of $(\hat{c}, \hat{p}, \hat{\beta})$. Consequently, based on the large sample theory of the MLEs (see Lawless ^[38]) the approximate (asymptotic) $100(1-\alpha)\%$ confidence intervals (ACIs) for the parameters $(c, p, \beta)$ are obtained respectively as:

where $z_{\frac{\alpha}{2}}$ is the $(1-\alpha)\%$ quantile of the standard normal distribution and $\hat{v}_{11}$, $\hat{v}_{22}$ and $\hat{v}_{33}$ are the diagonal elements of $\hat{V}$.

5.2.   Bootstrap confidence interval

The bootstrap resampling technique is a common method to construct confidence intervals for the unknown parameters and provides more efficient estimation than conventional asymptotic intervals for relatively small sample sizes. For more details, refer to Hall ^[39] and Efron and Tibshirani ^[40]. Given initial values of the parameters $(c, p, \beta)$, initial stress levels, initial inspection times and a progressive censoring scheme (PCs), the procedure to obtain interval estimates for the parameters $(c, p, \beta)$ by this method is described in the following steps:

Step 1: Compute $(\hat{c}, \hat{p}, \hat{\beta})$

Step 2 Using $(\hat{c}, \hat{p}, \hat{\beta})$, generate a bootstrap resample from equation (2.6).

Step 3: Calculate the bootstrap estimates $(\hat{c}_1, \hat{p}_1, \hat{\beta}_1)$.

Step 4: Repeat steps two to three up to $k$ times to obtain: $\hat{c}_1^B, \hat{c}_2^B, \ldots, \hat{c}_k^B, \hat{p}_1^B, \hat{p}_2^B, \ldots, \hat{p}_k^B, \hat{\beta}_1^B, \hat{\beta}_2^B, \ldots, \hat{\beta}_k^B$.

Step 5: Rearrange the bootstrap estimates in ascending order to have $100(1-\alpha)\%$ percentile confidence intervals for the parameters $(c, p, \beta)$, respectively, as:

6.   Simulation Study

This study aims to evaluate the performance of the obtained point and interval estimates of the model parameters. For this purpose, a Monte Carlo simulation of 1000 repetitions is conducted for the SSALT with $k = 4$ levels under PT-ICD assuming PRD is the underlying life model of the tested units. The performance of the ML point estimates is measured in terms of their average absolute biases (ABs) and their mean squared errors (MSEs), while the performance of interval estimates is measured in terms of their average lengths (ALs) and their coverage probabilities (CPs). Different settings in the simulation processes are considered as listed in Table 1. The simulation process is performed using MATLAB9.7R2019 software code through the following steps:

Step 1: Generate a random sample of size $n$ from the uniform distribution $U(0, 1)$.

Step 2: Using the inverse-CDF method and the random sample generated in step one, generate a random sample of size $n$ from model (2.6) based on the true parameter values, initial stress values $S_j$, $\tau_{j}$, $j = 1, 2, \ldots, 4$ and initial inspection times $\tau_j$, $j = 1, 2, \ldots, 4$.

Step 3: Based on the indicated progressively censored scheme (PSC), compute the number of failures $n_j$ and the number of removals $r_j$, $j = 1, 2, \ldots, 4$.

Step 4: Obtain the MLEs of the parameters $(c, p, \beta)$.

Step 5: Based on Step 4, obtain the Fisher information matrix.

Step 6: Obtain the 95% asymptotic confidence intervals.

Step 7: Repeat Steps (1–4) 1000 times to obtain the 95% bootstrap confidence intervals.

Step 8: Repeat Steps (1–7) 1000 times and compute ABs and MSEs of the MLEs for the parameters $(c, p, \beta)$, as well as ALs and CPs of the interval estimates.

The ABs and the MSEs of the MLE point estimates are listed in Table 2 and the ALs and the CPs of the 95% asymptotic and bootstrap interval estimates are listed in Tables 3–4.

From the simulation results listed in Tables 2–4, the following concluding remarks are clearly observed:

1. The obtained point estimates have satisfactory performance since they have relatively small ABs and MSEs under the three PCs.

2. The ABs and the MSEs of the point estimates are uniformly decreasing with an increase in the sample size and increasing with an increase in their indicated initial values.

3. The ABs and the MSEs of the point estimates are increasing as the number of removals is increasing. This is expected because point estimates obtained under the conventional Type-II (PC1) censored scheme have more information about the failure process than point estimates obtained under the PC2 and PC3 (progressively censored schemes).

4. The ALs of the asymptotic and bootstrap 95% confidence intervals are also decreasing with an increase in the sample size and increasing with an increase in their indicated initial values.

5. The bootstrap 95% confidence intervals have smaller ALs than their corresponding 95% asymptotic confidence intervals. (The intuitive reason for this result is that the bootstrap technique can be considered a modification of the ML method).

6. Generally, the CPs of the obtained interval estimates exceed the nominal 95% significance level, but for interval estimates with extremely small ALs, their CPs are strictly less than the 95% significance level because, mathematically, the conditional CPs of the interval estimates are decreasing with a decrease in their ALs.

7.   Optimum Inspection Times and Stress Levels

This section explores two criteria to determine the optimum inspection times and the optimum stress levels that lead to the most precision in the estimation of the model parameters. Both criteria are based on the Fisher information matrix since the determinant of the Fisher information matrix is proportional to the reciprocal of the volume of the joint confidence region of the parameters $\left(c, p, \beta\right)$. Therefore, one way to minimize the asymptotic variances of $\left(\hat{c}, \hat{p}, \ \hat{\beta}\right)$ is to maximize the determinant of the Fisher information matrix with respect to $\tau = (\tau_1, \tau_2, \ldots, \tau_k)^t$ and $= (S_1, S_2, \ldots, S_k)^t$. In reliability contexts, such a criterion is called the D-optimality criterion.

Given the proposed values of the progressively censored scheme $\pi = (\pi_1, \pi_2, \ldots, \pi_k)$, and the initial values of stress levels, applying this criterion yields the optimal values of $\tau = (\tau_1, \tau_2, \ldots, \tau_k)^t$ as solutions to:

where $|\hat{\mathcal{F}}|$ is the determinant of the Fisher information matrix $\mathcal{F}$ computed at $\left(\hat{c}, \hat{p}, \ \hat{\beta}\right)$.

Consequently, after determining the optimum inspection times $\tau^\ast = ({\tau_1}^\ast, {\tau_2}^\ast, \ldots, {\tau_k}^\ast)^t$, the optimum stress levels: $S^\ast = ({S_1}^\ast, {S_2}^\ast, \ldots, {S_k}^\ast)^t$ are obtained simultaneously as the solutions of the equation:

Another criterion, which is called the A-optimality criterion, is based on the trace of the Fisher information matrix. Applying this criterion, the optimum inspection times are determined by maximizing the sum of the diagonal elements of the Fisher information matrix with respect to $\tau = (\tau_1, \tau_2, \ldots, \tau_k)^t$ and $S = (S_1, S_2, \ldots, S_k)^t$ through solving the following two nonlinear equations:

More details about the D-optimality and the A-optimality criteria are explored in Gouno et al. ^[41]. Figure 3 presents the optimum inspection times and optimum stress levels when the sample size $n = 50$ with different true values of the model parameters. Table 5 lists the ABs and the MSEs of the point estimates when the true parameter values are $(c = 2, p = 0.3, \beta = 1.2)$, obtained using the optimum inspection times and optimum stress levels by applying the simulation process in the previous section when the sample size $n = 50$.

From Figure 3 and the results in Table 5, it is clear that:

1. The optimum inspection times and optimum stress levels increase as the true values of the model parameters increase.

2. The optimum inspection times obtained by PC1 are strictly greater than those obtained by PC2 and PC3. This implies that removing some units during the life test will shorten the termination time of the life testing experiment.

3. The optimum inspection times and optimum stress levels obtained by the D-optimality criterion are greater than their corresponding values obtained by the A-optimality criterion.

4. It is noted that the ABs and the MSEs of the point estimates obtained using the optimum inspection times and optimum stress levels by the D-optimality criterion are strictly lower than their corresponding values obtained by the A-optimality criterion.

5. Additionally, the ABs and the MSEs of the point estimates obtained using the optimum inspection times and optimum stress levels by both criteria are lower than their corresponding values in the simulation study. This confirms that the optimum inspection times and optimum stress levels significantly improve the point estimates of the model parameters.

8.   Application to Real Life Data

In this section, the developed estimation procedures are applied to a real-life data. This dataset consists of 66 observations representing the breaking stress of carbon fibers of 50 mm length (measured in GPa), as reported in Al-Aqtash et al. ^[42] and reanalyzed by Bhat and Ahmad ^[30]. Assuming the PRD modeling for this data, the MLEs of the parameters $\theta$ and $\beta$ are $\hat{\theta} = 4.850$, $\hat{\beta} = 1.721$ and the asymptotic 95% confidence intervals for $\theta$ and $\beta$ are $(2.818, 6.883)$ and $(1.396, 2.045)$, respectively. The estimated value of Kolmogorov-Smirnov test statistics is (K-S) = 0.08 with a p-value = 0.768. This confirms validity of the PRD to this dataset.

Consider a SSALT with $k = 5$ levels of stress. This data is divided into 5 groups, each consisting of 10 observations by the initial inspection times $\tau_1 = 1.81, \tau_2 = 2.46, \tau_3 = 2.75, \tau_4 = 3, \tau_5 = 3.3$, as it observed in Table 6.

Therefore, the termination time of the life test is $\tau_5 = 3.3$, thus we will have 16 units still working and removed from the step stress test after the fifth stage.

Suppose the scale parameter of the PRD, $\theta_j, j = 1, 2, …, 5$ at each stress level has the inverse power low given by $\theta_j = \theta_j = cS_j^p, c > 0$. Setting the initial value of the first stress level ($S_1 = 1$), initial values of $S_j, j = 2, 3, ..., 5$ are indicated such that:

This gives: $S_2 = 1.9197, S_3 = 2.6985, S_4 = 3.2406, S_5 = 5.9292$.

Using the progressively censored schemes:

Employing the D-Optimality and the A-Optimality criteria, MLEs of the model parameters and 95% asymptotic and bootstrap CIs with optimum inspection times and optimum stress levels are listed in Table 7.

From the above table, it is clearly observed that the obtained estimates are closest to their proposed initial values, agreed with the simulation results, and gave efficient estimates for the model parameters.

9.   Conclusions

This paper dealt with the k-level SSALT using PT-ICD as an extension of the traditional Rayleigh distribution. The PRD was considered to model the lifetimes of the tested units. The cumulative exposure distribution is considered to construct the lifetime model, assuming that the scale parameter of the PRD has the inverse power function of the stress level. Point estimates of the formulated model were obtained via the MLE approach. Additionally, interval estimates of the model parameters were obtained using the approximated variance-covariance matrix and the bootstrap resampling method.

To evaluate the performance of the obtained estimates, a simulation study of different combinations of life testing experiments was conducted under three PCs. The results of this study indicate that all obtained estimates have good characteristics and are efficient for estimating the model parameters.

In order to increase the accuracy of the obtained estimates, the test procedures were improved by obtaining the optimum inspection times and the optimum stress levels using the derived Fisher information matrix with two optimality criteria: The D-optimality criterion and the A-optimality criterion. The concluding results confirmed that the optimum inspection times and the optimum stress levels increase the efficiency of the obtained estimates with priority to the D-optimality criterion.

The developed estimation procedures also investigated in a real stress data. In this investigation, concluding results confirmed that the obtained estimates are closest to their initial practical setting and gave an efficient estimation for the model parameters.

In this work, the obtained results also indicate that both point and interval estimates obtained under the ordinary Type-II PC1 censoring scheme are more accurate than their corresponding values obtained under the PC2 and PC3 progressive censoring schemes. One explanation for this result is that, under PC1, there will be no removed survival units in the inspected intervals except for the last one. However, when employing PC2 and PC3 progressive censoring schemes, the life testing experiment will terminate in a reasonable time. Therefore, the experimenters have to quantify the test plans based on the intended objectives and available conditions of the life testing experiment.

As a new extension of the standard Rayleigh distribution, concluding results confirmed that PRD can be effectively used in the step stress-accelerated life tests for more analysis of reliability and survival data. Future research may extend this work, taking into account other estimating methods and acceptance sampling plans.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors would like to express their sincere gratitude to the anonymous referees for their insightful recommendations and remarks.

Conflict of interest

The authors declare there is no conflict of interest.

Appendix

1.   Appendix 1

To find the second derivatives of $L_1(c, p, \beta)$ with respect to the parameters $c$, $p$, and $\beta$, set:

then

2.   Appendix 2

To find the second derivatives of $L_2(c, p, \beta)$ with respect to the parameters $c$, $p$, and $\beta$, set

This implies

Similarly, $\left(\frac{\partial^2P_j}{\partial p^2}, \frac{\partial^2P_j}{\partial\beta^2} \right)$ are obtained by substituting $p$, $\beta$ in place of $c$ in Eq (A1).

$\frac{\partial^2P_j}{\partial c\partial\beta}, \frac{\partial^2P_j}{\partial\beta\partial p}$ can be also obtained by substituting $\beta$ in place of $p$ in Eq (A2) and $\beta$ in place of $c$ in Eq (A2), respectively. So, it is enough to find the second derivatives of $f_j$, $j = 1, \ldots, k$ with respect to $c$, $p$, $\beta$ as:

Where: