A parametric bootstrap control chart for Lindley Geometric percentiles

TL;DR

A new control chart using parametric bootstrap improves quality control for skewed Lindley geometric distributions in industrial processes.

Contribution

A novel parametric bootstrap control chart for Lindley geometric percentiles is introduced for better process monitoring.

Findings

The new control chart is sensitive to changes in Lindley geometric distribution parameters.

Subgroup size, percentiles, and significance levels significantly affect control limits.

The chart performs consistently across different percentiles and distribution parameters.

Abstract

Control charts are vital for quality control and process monitoring, helping businesses identify variations in production. Traditional control charts, like Shewhart charts, may not work well for skewed distributions, such as the Lindley geometric distribution (LG). This study introduces a new control chart that uses parametric bootstrap techniques to monitor percentiles of the LG distribution, providing a more effective quality control method. The LG distribution is useful for modeling material strength and failures, especially in structural design, where lower percentiles indicate reduced tensile strength. We conducted extensive simulations to assess the proposed control chart’s effectiveness, considering various distribution parameters, percentile values, Type I error rates, and sample sizes. Our findings highlight how subgroup size, percentiles, and significance levels affect control limits, stressing the need for careful parameter selection in monitoring processes. The results show that the new control chart is highly sensitive to changes in LG distribution parameters and performs consistently across different percentiles. This suggests its practical relevance and robustness for industrial applications in quality control. Future research should explore its performance in real-world production settings to confirm its efficiency and reliability.

Figures

pdf of LG distribution for different values of θ and p.

ARL for different values of p based on the LGD with (θ = 0.5, p = 0.05).

ARL for different values of θ based on the LGD with (θ = 0.5, p = 0.05).

ARL for different values of u.

Bootstrap percentile control chart with u = 0.5 for randomly generated subgroups of size 5 from the PLD.

Bootstrap percentile control chart with u = 0.1 for randomly generated subgroups of size 5 from the PLD.

The X¯ Shewhart control chart for randomly generated subgroups of size 5 from the PLD.

Histogram (top left), Q-Q plot (top right), CDF plot (bottom left), and P-P plot (bottom right) of survival time dataset along with the LG distribution fitted to it.

Bootstrap percentile control chart with u = 0.05 for subgroups of size 5 of the survival times data set.

Bootstrap percentile control chart with u = 0.05 for phases I and II data sets.

The X¯ Shewhart control chart for phases I and II data sets.

Full text

1. Introduction

The standard procedures of statistical quality control (SQC) often rely on control charts and acceptance sampling plans, which are typically based on the assumption of normal data. However, in practice, this assumption rarely holds true. Upon analyzing various data sets from different applications, such as statistical process control (SPC), we have observed that this type of data often displays moderate to strong asymmetry and light to heavy tails. Therefore, fitting a normal distribution to the data is generally not the most suitable option. Additionally, modeling real data sets, even when potential symmetric models for the underlying data distribution exist, is always challenging due to uncontrollable perturbation factors.

To address this issue, it is necessary to use distributions with more flexible characteristics as the quality characteristic distribution. The introduction of new distributions is driven by the goal of transitioning from former distributions to more flexible structures. The number of parameters, as well as the method of defining the distribution and its structure, are crucial for fitting and flexibility. An ideal distribution should exhibit better fitting, be more flexible in shape, and have easier formulas for implementation. The Lindley distribution [1] is one of the distributions particularly suitable for studying the reliability modeling of stress-strength characteristics due to its excellent properties. Hence, many researchers have proposed new classes of distributions based on extensions and modifications of the Lindley distribution. Table 1 shows a literature review of the distributions introduced based on the Lindley distribution, along with some of their most important applications, including reliability and SPC fields.

Table 1: Review of literature about the distributions introduced based on the Lindley distribution.

Recently, by compounding Lindley and geometric distributions, [13] introduced the Lindley-geometric (LG) distribution. Advantages of using the Lindley-Geometric distribution include:

Flexibility: The LG distribution offers flexibility in modeling continuous data with additional parameters, allowing for a better fit to the data and capturing various patterns and characteristics;Overdispersion and under dispersion: The LG distribution can effectively model data exhibiting overdispersion or under dispersion, providing a more accurate representation of the variance in the data compared to other distributions;Versatility: The LG distribution can be applied in various fields such as finance, economics, and environmental studies, making it a versatile tool for modeling continuous data;Statistical inference: The LG distribution allows for improved statistical inference and analysis of continuous data, enabling researchers to make more precise predictions and draw meaningful conclusions from the data.

Therefore, The Lindley geometric distribution can be utilized in statistical process control methodologies to monitor and improve quality control processes by providing a robust framework for analyzing and interpreting quality control data. So, in this article, we will present control charts for the quality characteristics that follow the LG distribution.

Various researchers have explored the use of bootstrap techniques in statistical process control charts. The non-parametric bootstrap approach can be utilized in control charts, removing the need for conventional parametric assumptions. The bootstrap techniques can also be applied in cases where the distribution of the statistic used for process monitoring is unknown. [37] introduced a bootstrap control chart for monitoring the process mean, offering an alternative to Shewhart’s $[eqn]$ chart, particularly useful for non-normal process data. Researchers such as [38–40] established the groundwork by introducing bootstrap control charts for quality engineering, process control, and reliability analysis. These investigations underscored the advantages of employing bootstrap methods to develop customized control charts tailored to distinct distributions and process parameters. Subsequent research led by [41–50] further advanced the use of bootstrap methods in developing control charts for various distributions and applications in quality control. These studies provided valuable insights into the reliability, efficiency, and novel approaches to control charting using bootstrap resampling techniques.

Innovative strategies introduced by researchers like [51–56] aimed to enhance the performance and accuracy of control charts. These approaches included Hotelling’s T^2^ control charts, model selection methods, and adjustments to control limits, with a focus on improving the detection of out-of-control conditions and upholding quality standards in industrial processes.

Recent studies in the 2020s conducted by researchers such as [36, 57–63] have delved into new frontiers in bootstrap control charting. These studies explore areas such as monitoring non-normal processes, analyzing the availability index, handling proportion data, and examining specific distribution percentiles. They provide valuable insights into overcoming challenges in quality control, manufacturing systems, and reliability data analysis through the application of bootstrap techniques.

The focus of this study lies in the LG distribution, which has been applied in modeling material strength and various types of failures, as demonstrated in studies such as [13, 34].

There is no explicit formula or known sampling distribution of a statistic for the LGD percentile in cases of small sample sizes. Therefore, in this research, we introduce a new control chart using parametric bootstrap for LGD percentiles.

This paper is organized as follows: We provide some of the characteristics of the LGD in Section 2. The procedure of a parametric bootstrap control chart for LGD percentiles is presented in Section 3. The performance of the bootstrap Lindley-geometric percentile control chart is examined in Section 4. Illustrative examples in a simulated data set and in a real-world data set are included in Section 5. The paper ends with discussions and concluding remarks in Section 6.

2. Some of the characteristics of the LGD

In this section, we will first explain some characteristics of the LG distribution, such as the graph of the probability distribution function(pdf), in addition to the mean, variance, and quantile function. Then, we will talk about estimating the parameters, and generating random numbers of the LGD.

[13] introduced the LG distribution with pdf and cumulative distribution function (cdf) given by

$$[eqn]$$ $$[eqn]$$

Several distinct sub-models of the LG distribution (1) can be derived. For the Lindley distribution, the case where p = 0 is considered. As p approaches from left to 1, the LG distribution converges to a degenerate distribution at zero. Consequently, the parameter p is regarded as a concentration parameter. The density function of the LG distribution exhibits the following characteristics: (i) it decreases for all combinations of p and θ where p is greater than (1-θ^2^)/(1 + θ^2^), and (ii) it is unimodal for all combinations of p and θ where p is less than or equal to (1-θ^2^)/(1 + θ^2^). Fig 1 illustrates the pdf for the LGD in Eq (1), showing how it changes at different values of θ and p.

pdf of LG distribution for different values of θ and p.

2.1 Moments of the LGD

Assuming Y follows a LG(p, θ) distribution, utilizing Eq (1) and applying the binomial expression for $[eqn]$ , the r—th moment of Y can be expressed as

$$[eqn]$$

Replacing the r = 1 in the above equation, yields the mean of the LGD as

$$[eqn]$$

Also, the variance of LGD is obtained by substituting r = 2 in the moment expression and using the following relationship

$$[eqn]$$

2.2 Quantile function of the LGD

Let Y be an arbitrary random variable with cdf F(y) = Pr(Y ≤ y), where x ∈ R. For any u ∈ (0,1), the u—th quantile function, Q(u) of Y is the solution of

F (Q(u)) = u,

for Q(u) > 0.

For any fixed θ > 0, from Eq (2), we obtain

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

In the above equation, we note that (θQ(u)+θ + 1) is the Lambert W function of the real argument $[eqn]$ The Lambert W function is defined by

$$[eqn]$$

The Lambert function has two real branches with a branching point located at (-e^-1^,1). The lower branch, W-1(y), is defined in the interval [-e^-1^,1] and has a negative singularity for x → 0^-^. The upper branch, W0 (y), is defined for y ∈ [-e^-1^,∞].

Then, we have

$$[eqn]$$

Clearly, for any θ > 0, 0 < p < 1 and u ∈(0,1), we have (θQ(u)+θ+1) > 1 and then $[eqn]$ . Therefore, considering the lower branch of the Lambert W function, we can write Eq (3) as

$$[eqn]$$

Hence, the quantile function of Y is given by

$$[eqn]$$

2.3 Estimation of the parameters of the LGD

It is necessary to estimate the parameters as accurately as possible to obtain efficient information about the distribution. The maximum likelihood method is generally preferred for parameter estimation.

We consider the maximum likelihood estimation (MLE) about the parameters (θ,p) of the LGD. Suppose yobs = {x1,x2,…,xn} is a random sample of size n from the LGD. Then the log-likelihood function is given by

$$[eqn]$$

The MLEs of p and θ are obtained by solving the following equations:

$$[eqn]$$ $$[eqn]$$

However, they do not lead to explicit analytical solutions for the parameters. Thus, the estimates can be obtained by means of numerical procedures such as Newton-Raphson method. The R software provides the nonlinear optimization routine for solving such problems. Also, to obtain the MLE of the LGD, [13] introduced an EM Algorithm.

The other methods of estimation methods of the LGD parameters such as least-squares, weighted least-squares, Anderson-Darling, and Crámer–von-Mises are presented by [34].

2.4 Generating random numbers from LGD

In this section, we discuss how to generate random numbers from the LGD. One reliable method for generating random numbers for the LGD is using the quantile function for the distribution. Instead of determining the cumulative probabilities for a set of values, the quantile function determines the values for a set of cumulative probabilities. By utilizing the quantile function, we can generate random numbers from a uniform distribution and then transform them into the distribution of interest using its quantile function.

So, we can generate a random number x from the LGD using the following algorithm.

Algorithm 1: Generating a random number from the LGD

Step 1: Generate u from a Uniform (0,1) distribution.Step 2: calculate x = Q(u) where Q(u) is calculated from Eq (4).

3. Bootstrap control charts for percentiles of the LGD

Understanding the strength distribution is crucial for structural design purposes, especially when considering the lower percentiles that signify a decrease in material tensile strength. Given the asymmetric nature of the LG distribution, traditional $[eqn]$ and R Shewhart control charts may not effectively detect shifts in specific lower percentiles, particularly for quality characteristics like the breaking strength of brittle materials. Our focus is on utilizing the parametric bootstrap (percentile) method to develop control charts for LG distribution.

The following algorithm, which is similar to those proposed by [41] can be used to construct the bootstrap Lindley Geometric percentile control chart.

Algorithm 2: Constructing the bootstrap percentile control chart for the LGD

Gather n × m observations from an in-control, stable process assuming they follow a LGD with unknown parameters p and θ. The observations, denoted as xij with i = 1, …,n and j = 1, …,m, are drawn from m independent subgroups of size n.Determine the maximum likelihood estimators (MLEs), $[eqn]$ and $[eqn]$ , for the unknown parameters using all n × m observations. These estimators are obtained by solving the Eqs (6) and (7).Create a parametric bootstrap subgroup sample of size $[eqn]$ , from the LGD using the maximum likelihood estimators, $[eqn]$ and $[eqn]$ which are calculated from step 2, as the estimated parameters.Determine the maximum likelihood estimators from the bootstrap subgroup sample and represent them as $[eqn]$ and $[eqn]$ .For the bootstrap subgroup sample, calculate $[eqn]$ , which is the estimated 100u-th percentile, Q(u).Iterate through steps 3–5 numerous times, denoted as B, to acquire B bootstrap estimates of Q(u), labeled as $[eqn]$ .Sort the B bootstrap estimates $[eqn]$ in ascending order. The LCL is the value of the smallest ordered $[eqn]$ that has (α/2)B values below it. In this context, α represents the probability of considering an observation as out of control when the process is actually in control (commonly set at 0.0027 for Shewhart-type charts). The UCL is the value of the smallest ordered $[eqn]$ that has (α/2)B values above it.

After calculating the control limits, for implementation the phase II of the bootstrap control chart, subsequent subgroup samples of size n are collected from the process at regular intervals, and Q(u) is estimated for each new subgroup using the MLE as described in step 5. If the estimate, $[eqn]$ , falls within the UCL and LCL determined in step 7, the process is considered to be in control. Any $[eqn]$ values below the LCL or above the UCL suggest process may be out of control. Therefore, following the determination of bootstrap control limits, the process is monitored using the statistic $[eqn]$ in the conventional manner.

The subsequent section will evaluate the effectiveness of the suggested bootstrap Lindley-geometric percentile control chart through computer simulations. Specifically, an examination of the statistical properties of the bootstrap Lindley-geometric percentile control limits will be conducted. Additional simulations will assess the performance of the bootstrap percentile control chart in terms of Average Run Lengths (ARLs), and their standard deviations (SDARLs).

4. Performance of the bootstrap Lindley-geometric percentile control chart

This section includes computer simulations that assess the effectiveness of the bootstrap Lindley-geometric percentile control chart. The analysis involves investigating the behavior of the bootstrap control limits by determining the average UCL and LCL (MLCL, and MUCL) along with their respective standard deviations (SDLCL, and SDUCL) from the simulations. Additionally, further simulations are conducted to evaluate the ARL when the process is in control, along with the associated standard deviation for each ARL (i.e., SDRL). These simulations encompass diverse sample sizes, various percentiles of interest, and different levels of α. The ARL is also computed for when the process is out of control, accompanied by its variance. Similar simulations are carried out with varying sample sizes, different percentiles of interest, and various α values.

The mean UCL and LCL along with their respective standard deviations were calculated as follows: A total of n × 25 observations were generated from a LG distribution with parameters p and θ, following the Algorithm 1, and the procedure outlined in step 1 of Algorithm 2. Subsequently, steps 2 to 7 of Algorithm 2 were executed. The value of B was set to 10,000 in this scenario. This entire sequence of steps (1–7) was iterated k = 100 times, and the mean LCL and mean UCL were determined from the 100 data sets created using the Monte Carlo method. Additionally, the standard deviations of the control limits were calculated based on these 100 values. The simulation process was conducted across various sample sizes (n = 4,5,6), different percentiles (u = 0.05,0.10), different parameters of the LG distribution (θ and p), and diverse α values (α = 0.0027,0.002,0.01). The simulations were implemented using the R software, and Tables 2–4 present some of the outcomes. As anticipated, an increase in subgroup size from four to five to six leads to the convergence of control limits, while elevating the percentile from 0.05 to 0.10 results in the widening of the limits, and lower α (significance level) leads to wider control limits. Additionally, an increase in any of the θ or p parameters of the LG distribution narrows the control limits.

Table 2: Average control limits along with their respective standard deviations for α = 0.0027.

Table 3: Average control limits along with their respective standard deviations for α = 0.002.

Table 4: Average control limits along with their respective standard deviations for α = 0.01.

The in-control ARL was determined by creating m = 25 subgroups, each comprising n observations, and calculating the corresponding bootstrap UCL and LCL using steps 2–7 detailed in Algorithm 2. After establishing the control limits, forthcoming subgroup samples were generated from a LG distribution with identical parameters. For every new subgroup, the MLE was used to estimate Q(u). The count of instances where the estimate, $[eqn]$ , fell within the limits was tracked until an estimate exceeded the control limits. The run length was computed as the total instances the estimate remained within the limits plus the initial subgroup indicating an out-of-control situation. The process of determining the run length was repeated 1000 times, and the average of these 1000 run lengths (Average Run Length, ARL) along with their standard errors (SDRL) were computed. These simulations for the in-control ARL were conducted across various values of the parameters for decreasing and unimodal LG densities. The simulations were also carried out for different percentiles, with different values of α across all runs. To determine the out-of-control ARLs and related SDRLs, 25 sets of samples with a size of n were generated from a LG distribution with parameters θ and p. The UCL and LCL were calculated in a similar manner as for the in-control ARLs. Subsequent subgroup samples, each consisting of n observations, were then simulated from a LG distribution with distinct values of θ and p. The alteration in the parameters mimics a shift in the process percentile, signaling an ’out-of-control’ state to be identified. The percentile estimate, $[eqn]$ , was calculated for every subgroup of size n. The count of subgroups where $[eqn]$ lies within the control limits was tallied until an estimate exceeded the limits. The run length was determined by counting the instances where the estimates stayed within the limits, along with the first subgroup indicating an out-of-control situation. This process was iterated 1000 times, and the ARLs and their corresponding standard errors were derived from the 1000 generated run lengths.

Selected results for the LGD parameters of p = (0.25,0.75) and θ = (0.25,0.75), α values of 0.0027, 0.002, and 0.01, and u values of 0.05 and 0.10 for the ARL simulations are presented in Tables 5–7. According to established theory, the reciprocal of the α value is anticipated to represent the theoretical in-control ARL. Given α = 0.0027, the in-control ARLs for these simulations is expected to be approximately 1/α, which equals 370. In Table 5, it can be seen that when the values of θ and p in the LG distribution do not change (i.e. in-control ARLs is calculated), these values are around 370. Furthermore, Tables 6 and 7 show that the in-control ARL values for α values of 0.002 and 0.01 are around 500 and 100 respectively. When the process is in-control, a smaller ARL suggests that the calculated control limits may be overly restrictive, while ARLs exceeding 1/α indicate that the control limits may be too lenient or that the bootstrap control charts result in fewer false alarms. For the out-of-control ARL simulations, the bootstrap LGD percentile control chart detects the shift towards an out-of-control state. In this case, the smaller the calculated value of the ARL, the quicker the control chart detects changes in the process. For example, examining the first row of Table 5 reveals that when α = 0.0027 and n = 4, for u = 0.05, as the parameters of the LG distribution shift from (θ = 0.25,p = 0.25) to (θ = 0.25,p = 0.75), the bootstrap percentile control chart detects this change on average after 8.393 sampling times (with a standard deviation of 7.764). For u = 0.1, this change is detected on average after 7.603 sampling times with an approximate standard deviation of 6.876.

Table 5: ARLs along with their respective standard deviations for α = 0.0027.

Table 6: ARLs along with their respective standard deviations for α = 0.002.

Table 7: ARLs along with their respective standard deviations for α = 0.01.

Generally, based on Tables 5–7, changes in the parameter p are detected faster compared to the same amount of change in the θ parameter. Moreover, if both θ and p parameters change, this change is detected much faster. For example, in the first row of Table 5, it can be observed that if the parameters of the LG distribution change from (θ = 0.25,p = 0.25) to (θ = 0.75,p = 0.75), this change for u = 0.05 is detected after approximately 1.91 sampling iterations with an estimated standard deviation of 1.291, and for u = 0.1 after approximately 1.884 sampling iterations with an estimated standard deviation of 1.192, using the proposed bootstrap control chart. These demonstrate the very good performance of the suggested bootstrap control charts. Increasing the sample size (n) also leads to a reduction in the ARL. This means that with an increase in sample size, the bootstrap percentile control chart detects changes in the quality characteristic distribution parameters more quickly. Finally, the choice of percentile u (whether 0.05 or 0.1) does not have a significant impact on the performance of the control charts introduced.

In the topic of reliability, there is a well-known proverb that states: "A chain is only as strong as its weakest link." Therefore, in this context, lower percentiles are very important. The upper percentiles of the qualitative distribution are also noteworthy in the quality control field. Consequently, in this section, we have examined our control chart concerning the upper percentages of the LG distribution. For this purpose, we considered the two upper percentiles, u = 0.75 and u = 0.9. Then using the same parameters p, θ, and n that we applied for the lower percentiles, we calculated the values of MCL,MUCL,SDLCL,SDUCL,ARL, and SDRL. Since the usual level of significance in the field of quality control is α = 0.0027, and including other α levels would increase the volume of the article without affecting the results, we summarize the calculations in this section in this level only. The results are presented in Tables 8 and 9.

Table 8: Average control limits along with respective standard deviations for α = 0.0027 with upper percentiles.

Table 9: ARLs along with their respective standard deviations for α = 0.0027 with upper percentiles.

The results are quite similar to those obtained for the lower percentiles. Again, the in-control ARLs are around 370. Changes in the parameter p are detected more quickly compared to the same magnitude of change in the θ parameter. When both θ and p parameters change, this change is detected even faster. Additionally, increasing the sample size (n) leads to a reduction in the ARL.

5. Sensitivity analysis

The accurate estimation of LG distribution parameters play a fundamental role in the efficiency and implementation of the control charts presented in this article. Incorrect estimation of these parameters may reduce the ARL of the control chart, resulting in a delayed discovery of changes in the distribution of the qualitative characteristic. Therefore, in this section, we analyze the sensitivity of these control charts to the LG distribution parameters. We also examine the effect of selecting percentiles on the implementation of the control chart. For this purpose, while keeping other parameters constant, we assess the impact of changes in the parameter for which we want to perform sensitivity analysis on the ARL of the control chart.

For the calculations in this section, we considered the values of n = 5 and α = 0.0027, which are commonly used in the field of quality control. Additionally, to analyze the sensitivity of the parameters, we set u equal to 0.1. The results remain consistent for other values of α, n, and u. We defined the base distribution as LG with parameters θ = 0.5 and p = 0.5. Using Algorithm 2, the percentile bootstrap control limits for this distribution were calculated as follows: MLCL = 0.02772982, MUCL = 1.241885, SDLCL = 0.002509196, and SDUCL = 0.03451446. For the sensitivity analysis of the θ parameter, we fixed p at 0.5 and varied θ from 0.05 to 0.95. Conversely, for the sensitivity analysis of the p parameter, we fixed θ at 0.5 and varied p from 0.05 to 0.95. We then assessed the effect of these changes on the ARL of the control chart. The results of sensitivity analysis for p and θ are presented in Figs 2 and 3 respectively.

ARL for different values of p based on the LGD with (θ = 0.5, p = 0.05).

ARL for different values of θ based on the LGD with (θ = 0.5, p = 0.05).

As seen in Fig 2, the value of ARL for p = 0.5 (ARL0) is approximately 370. As we move further away from pp, the ARL value (ARL1) decreases. However, the slope of this decrease is steeper on the left side, such that for equal distances from 0.5, the ARL value for lower p values is less. For example, for p = 0.25, the ARL value is approximately 5, while for p = 0.75, which is the same distance above 0.5, this value is nearly 108. This demonstrates that a decrease in the p parameter is detected by the control chart much earlier than an increase by the same amount in p. According to Fig 3, the situation is completely opposite for the θ parameter. That is, at equal intervals above and below 0.5, the ARL value for upper θ values are lower than for lower θ values. For example, when the θ is equal to 0.75, the ARL value is approximately 99, while this value is approximately 183 for θ = 0.25. This means that an increase in the θ parameter is detected by the control chart much earlier than a decrease of the same amount in this parameter.

To examine the effect of the selected percentile (u) on the implementation of the control chart, we again fixed the value of n at 5 and α at 0.0027. The base distribution for the qualitative characteristic is considered to be LG with θ = 0.5 and p = 0.5. When this distribution is changed to LG with parameters θ = 0.75 and p = 0.25, we calculated the ARL of the control chart for u values ranging from 0.05 to 0.95. The results are presented in Fig 4.

ARL for different values of u.

According to Fig 4, it can be seen that for each of the 0.05 to 0.65 percentiles, if a change in the LG distribution parameters occur as mentioned, this change is discovered on average after 20 samples. However, choosing higher percentiles will detect this change a little earlier. Specifically, the ARL value for u = 0.95 is almost equal to 15, which is smaller than that of the other percentiles. Therefore, if we want to detect the change in the LG distribution parameters from θ = 0.5,p = 0.5 to θ = 0.75,p = 0.25 earlier, it is better to use the 0.95 percentile.

In general, the issue of selecting the appropriate percentile for the bootstrap control chart for the LG distribution depends on which parameter changes are most significant to us. Once we determine the extent of increase or decrease in the parameters of interest, we can decide, through simulations similar to this section, which of the percentiles would be more beneficial for the specific issue at hand.

6. Illustrative examples

In this section, we analyze the application of the bootstrap percentile control chart for the of the LG distribution. This will be demonstrated through a simulated example and a real-world dataset (S1 and S2 Tables).

6.1 A simulated example

To create the simulated example in the R software, we initially generated 25 sets of 5 samples each from the LG distribution with parameters (θ = 0.25, p = 0.25) using algorithm 1. Subsequently, with α = 0.0027, and following steps 2 to 7 of Algorithm 2, we derived the values for LCL and UCL for u = 0.05 and u = 0.1. The average of these values was then used as the LCL and UCL for phase I. The computed values are as follows: for u = 0.05, the LCL = 0.06308313 and the UCL = 2.079311; and for u = 0.1, the LCL = 0.1308916 and the UCL = 3.448046.

For the implementation phase (phase II), we initially generated 10 subgroups of size 5 from the LGD with parameters (θ = 0.25,p = 0.25) to demonstrate samples under control (samples 26 to 35). Additionally, we created five subgroups of size 5 from the LGD with parameters (θ = 0.75,p = 0.75) to represent samples out of control (samples 36 to 40). Figs 5 and 6 illustrate the quantile-based bootstrap control charts for the LGD subgroups associated with u = 0.5 and u = 0.1, respectively, covering subgroups 1 to 40 in the simulation. Each point in Figs 5 and 6 represents $[eqn]$ calculated from the estimation of the LGD parameters ( $[eqn]$ and $[eqn]$ ) in the corresponding subgroup, which is calculated from Eq 4.

Bootstrap percentile control chart with u = 0.5 for randomly generated subgroups of size 5 from the PLD.

Bootstrap percentile control chart with u = 0.1 for randomly generated subgroups of size 5 from the PLD.

As shown in Figs 5 and 6 it is clear that samples 36, 37, 39, and 40 exhibit out-of-control behavior in both the bootstrap percentile control charts with u = 0.05 and u = 0.1. It is important to note that the percentiles of process subgroups have shifted starting from sample 36 due to changes in the LGD parameters. Both of our control charts effectively detected this shift in sample 36 with an ARL of 1, which is a significant achievement. This scenario demonstrates the strong performance of the bootstrap percentile control charts presented in this article and highlights that the choice of percentile u (whether 0.05 or 0.1) has no impact on the performance of the control charts introduced.

We have also presented the $[eqn]$ Shewhart control chart for the simulated data in this example, created using Minitab software, as shown in Fig 7. In this control chart, each point represents the average of five subgroups. The lower control limit (LCL) and upper control limit (UCL) were calculated using the first 25 subgroups as phase 1 data, and the averages of subgroups 26 to 40 were added as phase 2 points. As can be seen in Fig 7, all the points are under control, and this chart could not identify the change in the LG distribution parameters from subsamples 36 to 40. Therefore, this example demonstrates that the $[eqn]$ Shewhart control chart is not suitable for such data. In contrast, according to Figs 5 and 6, our proposed control chart detected this change immediately, indicating its high efficiency.

The X¯ Shewhart control chart for randomly generated subgroups of size 5 from the PLD.

6.2 A real data set example

To evaluate the application of the control charts presented in Section 2 to a real dataset, we examined the dataset "survival times in years of gastric cancer patients" obtained from [64]. Given that the data were originally ordered, we randomized them into 9 samples of size 5 to enable their use in the control charts. Table 10 shows the observed subgroups of the dataset. Using the "fitdistrplus" package in R software, we modeled the dataset with the LG distributions. The estimated LGD parameters were determined as $[eqn]$ and $[eqn]$ using the maximum likelihood (ML) method (stage 2 of Algorithm 2). Additionally, the Kolmogorov-Smirnov test statistic for fitting the LG distribution to the dataset is 0.09487884 with a p-value of 0.8124, indicating that the LGD fits the dataset well. Furthermore, we calculated the Cramer-von Mises statistic to be 0.06430339, the Anderson-Darling statistic to be 0.46977558, the Akaike Information Criterion (AIC) to be 120.345, and the Bayesian Information Criterion (BIC) to be 123.9583, once again demonstrating a good fit of the LGD to the dataset. Fig 8 shows the histogram, Q-Q plot, CDF plot, and P-P plot of the survival time dataset, along with the LG distribution fitted to it.

Histogram (top left), Q-Q plot (top right), CDF plot (bottom left), and P-P plot (bottom right) of survival time dataset along with the LG distribution fitted to it.

Table 10: Observed subgroups of the survival times dataset, along with the Q^(u = 0.05) for each subgroup.

Applying stages 3 to 7 of Algorithm 2 with α = 0.0027, u = 0.05, n = 5, B = 10000, and k = 100 yielded the following results: MLCL = 0.003753354, SDLCL = 0.0003751751, MUCL = 0.2262595, and SDUCL = 0.00220619. Table 10 shows the $[eqn]$ for each subgroup, where u is set to 0.05, and Fig 9 illustrate the bootstrap percentile control chart for the subgroups of size n = 5 for the survival times in years of gastric cancer patient’s dataset with u = 0.05.

Bootstrap percentile control chart with u = 0.05 for subgroups of size 5 of the survival times data set.

As seen in Fig 9, the 5th percentile of all observations of the subgroups is in control. Therefore, this control chart can be used for phase II monitoring of the 5th percentile subgroups with a sample size of 5 for future data.

6.3 A real data set example for phase II

Continuing from the previous example to demonstrate how to implement the control chart introduced for Phase II, we have considered the dataset number 6 presented in [34]. In the mentioned article, the LG distribution was fitted to this dataset with parameters θ = 0.5455 and p = 0.6348. We want to determine whether we can detect this change in the parameters of the LG distribution using the control chart presented in Phase I for the survival time dataset. For this purpose, we have taken the initial 30 data points of dataset number 6 from [34] and considered them as phase II data in subgroups of 5 from subgroups 10 to 15. Table 11 shows the subgroups along with $[eqn]$ corresponding to each subgroup for u = 0.05.

Table 11: Observed subgroups of the phase II dataset, along with the Q^(u = 0.05) for each subgroup.

Fig 10 shows the bootstrap percentile control chart for the of data from phases I and II simultaneously, where the LCL and UCL are calculated from the data of phase I, and the points on the chart represent the $[eqn]$ calculated for each subgroup. In this figure, points 1 to 9 relate to the phase I dataset, while points 10 to 15 pertain to the phase II dataset. As seen in Fig 10, subsample number 13 is out of control. This indicates that the bootstrap control chart has detected a change in the parameters of the LG distribution after four sampling iterations.

Bootstrap percentile control chart with u = 0.05 for phases I and II data sets.

Fig 11 shows the $[eqn]$ Shewhart control chart for the data sets of phases I and II, where the UCL and LCL are calculated from the phase I data set, and the points on the chart represent the subgroup means. As seen in this figure, all points are under control, and if we were using thiss control chart, we would not be able to detect changes in the parameters of the LG distribution.

The X¯ Shewhart control chart for phases I and II data sets.

7. Discussion and conclusions

In many industrial settings, quality characteristics often follow skewed distributions, making traditional methods like Shewhart control charts unsuitable. To address this, we can use more flexible distributions, such as the Lindley geometric (LG) distribution, which is valuable for modeling material strength and failures. Understanding the strength distribution is crucial for structural design, especially when lower percentiles indicate reduced tensile strength.

This research introduces a new control chart that employs parametric bootstrap techniques for LG distribution percentiles. We first outlined key features of the LG distribution, including quantile calculation, parameter estimation, and random number generation. Then, we developed an algorithm to determine the upper and lower limits of the bootstrap percentile control chart.

We evaluated the effectiveness of the proposed control chart through simulations, calculating average run lengths (ARLs) and their standard deviations in both in-control and out-of-control scenarios. The simulations, conducted using R software, revealed several important findings:

Larger subgroup sizes result in narrower control limits.Higher percentiles lead to wider control limits.A lower significance level (α) results in wider control limits.Increases in LG distribution parameters narrow the control limits.Changes in the parameter p are detected more quickly than changes in the parameter θ.Detection is faster when both θ and p are altered.Increasing the sample size (n) decreases the ARL.The choice of percentile (0.05, 0.1, 0.75 or 0.95) has minimal impact on the control chart’s performance.

In conclusion, the simulation results demonstrate the high efficiency of the proposed control chart. Its performance was further validated through both simulated and real-world examples, both of which highlighted its exceptional effectiveness.

8. Appendix

8.1 The R code of Algorithm 1

Generating random numbers from LGDlibrary(lamW)qlingeo<-function(p,Theta, pp){ stopifnot(p < 1, p > 0, Theta> 0) x = -lamW::lambertWm1(-(p-1)*(Theta+ 1)exp(-(Theta+ 1))/(ppp-1))/Theta- 1/Theta- 1 return(x)}rndlingeo<-function (n,Theta, pp){ y = stats::runif(n, min = 0, max = 1) randdata = qlingeo(y, pp, Theta) return(randdata)}

8.2 The R code of Algorithm 2

Constructing the bootstrap percentile control chart for the **LGD.fit_distribution <- function(sample){y<-samplelogL<-function(x){n<-length(y)th<-x[1]p<-x[2]lnL<-2nlog(th)-nlog(th+1)+nlog(1-p)+sum(log(1+y))-thsum(y)-2sum(log(1-p(1+thy/(th+1))exp(-thy)))return(-lnL)}opt<-optim(c(0.5,0.5),lower = c(0.01,0.01),upper = c(10,0.999),logL, method = "L-BFGS-B")phat <-opt$par[1]thetahat <-opt$par2return(c(phat, thetahat))}LGBootCC <- function(u, alpha, theta, p, n = 5, BootStrap_iterations = 1E+4,MontiCarlo_iter = 100){LCL <-0UCL <-0for(itr in 1:MontiCarlo_iter) { QS<-0 Theta<-0 P<-0  for(i in 1:BootStrap_iterations)   {   sample <- rndlingeo(n, theta, p)   result <- fit_distribution(sample)   if (!is.null(result))   {    Theta[i] <- result[1]    P[i] <- result[2]    QS[i] <- qlingeo(u, Theta[i], P[i])   }  }QS <- sort(QS)LCL[itr] <- QS[as.integer(BootStrap_iterations * alpha / 2)]UCL[itr] <- QS[as.integer(BootStrap_iterations * (1 ‐ alpha / 2))]cat("LCL", LCL[itr], "UCL", UCL[itr], "\n")}MLCL <- mean(LCL)MUCL <- mean(UCL)SDLCL <- sd(LCL)SDUCL <- sd(UCL)cat("Mean of LCL (MLCL): ", MLCL," Mean of UCL (MUCL): ", MUCL, "\n")cat("Standard deviation of LCL (SDLCL): ", SDLCL," Standard deviation of UCL (SDUCL): ", SDUCL, "\n")}

8.3 The R code of fitting the LG distribution to the survival times dataset

dlingeo<-function (x, pp, Theta){ pdf = (Theta^2/(Theta + 1))(1 ‐ pp)(1 + x)* (exp(-Thetax))(1 ‐ pp*(1 + (Thetax)/(Theta + 1))(exp(-Thetax)))^(-2) return(pdf)} plingeo<-function (q, pp, Theta){ cdf = (1 - (1 + (Thetaq)/(Theta + 1))exp(-Thetaq))/(1 ‐ pp*(1 + (Theta*q)/(Theta + 1))exp(-Thetaq))return(cdf)}# Run qlingeo function of the Algorithm 1 codelibrary(fitdistrplus)servtime<-c(1.326, 0.841, 0.282, 2.830, 0.121,0.644, 0.197, 1.581, 2.178, 1.553,1.447, 2.343, 0.863, 3.658, 0.132,4.033, 3.978, 2.416, 0.534, 0.501,0.458, 4.003, 0.260, 1.099, 0.696,0.164, 1.271, 0.641, 3.743, 0.395,0.203, 0.296, 0.529, 1.485, 2.825,0.115, 1.589, 2.444, 1.219, 3.578,0.540, 0.507, 0.466, 0.047, 0.334)fitlingeo <- fitdist(servtime,"lingeo", start = list(pp = 0.25,Theta = 0.75),lower = c(0.1,0.2),upper = c(0.99,10))as.vector(summary(fitlingeo)$estimate)gofstat(fitlingeo)plot(fitlingeo)

Supporting information

S1 TableObserved subgroups of the survival times dataset.(DOCX)

S2 TableObserved subgroups of the phase II dataset.(DOCX)