Applying User Surveys and Accelerated Tests Data to Estimate Reliability of New Consumer Products Using a Discrete Life Distribution Model

Assessing the reliability of consumer products that are subject to random discrete damage is critical during their design. Previous studies have approximated a continuous life distribution for consumer products by treating the occurrence (cycles) of accumulating damage as a continuous random variable. However, when the lifetime of a product is only a few damage cycles (e.g., ten drop cycles of a laptop), using a discrete lifetime distribution is more accurate. Using a discrete lifetime distribution is challenging because it contains a summation term with an unknown upper bound, which makes calculating its likelihood function cumbersome. This paper proposes a method to address this issue. First, the upper bound of the summation term is approximated through a gradient descent algorithm and the Maximum Likelihood Estimation method. Then, the upper bound is fixed, and the other parameters of the reliability model are estimated using Bayesian analysis. The paper presents a hypothetical case study that shows that using the discrete model leads to a more accurate estimation of product life when dealing with a small number of cycles.


I. INTRODUCTION
The reliability of a product is the probability that a product performs its intended function adequately for a specified time under specified use conditions [1]. The reliability of a product, thus, depends on its usage time, applied stresses (i.e., agents that cause damage to the product), and stress adjustors (i.e., conditions that increase or decrease the stress magnitude or stress absorption).
To estimate product reliability, some previous studies used damage-based reliability estimation approaches to consider the effect of various stress and stress adjustors in reliability analysis [2,3,4]. These methods are applicable to reliability test data or sensor data that are costly to collect.
Our previous studies [5,6,7] developed a stress-based reliability estimation approach that utilized user survey data as a cost-effective source for estimating reliability.
That approach assumed a continuous lifetime distribution and, through Bayesian analysis, estimated the reliability model of a consumer product. In addition, the concept of the equivalent cycle was defined to consider the effect of various stresses and stress adjustors.
A continuous distribution is applicable when the random variable representing life is continuous (e.g., expressed in terms of calendar time), or the product undergoes a large number of damage-accumulating usage cycles (e.g., hundreds of cycles of a device dropping on a hard floor). On the other hand, when the performance of a product is measured occasionally (e.g., every week or every month), or the product undergoes a few cycles of damage-producing usage before failing, a discrete random variable and discrete forms of parametric distributions can more accurately estimate the life distribution or reliability model [8,9]. This paper extends our previous approach for discrete lifetime distributions. Some discrete life distributions have a summation term whose upper bound can be expressed as an unknown equivalent cycle. Such upper bounds make reliability estimation through a Bayesian or Maximum Likelihood Estimation (MLE) method (i.e., a method for estimating the parameters of a parametric distribution) very difficult. The discrete Weibull distribution type III has this issue. Ideally, in reliability analysis, a method for selecting the best distribution should be used. However, to demonstrate the proposed generic approach in this paper, the case study considers the discrete Weibull type III distribution as a complicated discrete distribution with the summation term.
The Weibull distribution is a well-known lifetime distribution that is widely used in reliability engineering. The Weibull distribution is available in continuous and discrete forms. Three types of discrete Weibull distributions have been proposed: (1) the type I has a reliability function that mimics the reliability function of a continuous Weibull distribution, (2) the type II has a hazard function that mimics the hazard function of a continuous Weibull distribution, and (3) the type III is more generic and does not follow any function of a continuous Weibull distribution [10].
Various extensions of discrete Weibull distributions have also been introduced, and different methods for estimating their parameters have been suggested. For instance, Jia et al. [11] proposed a discrete extended Weibull distribution and used MLE for estimating its parameters. Barbiero [8] proposed three methods, including the method of proportion, MLE, and the method of moments for estimating the parameters of a discrete Weibull type III distribution. The above studies assumed that the product is used under a fixed stress level, and the number of uses (damage producing) cycles is equivalent to the number of times the stress is applied to the product. On the other hand, the approach presented in this paper assumes that consumer products are used under varying stress levels, and the concept of equivalent cycles is used to calculate the number of cumulatively damaging usage cycles under a reference stress condition. As discussed earlier, when using the idea of equivalent cycles, the upper bound of the summation term becomes an unknown value that is hard to estimate. This paper considers a consumer product that is used under various stress conditions, and its failure time is expected to occur within a few applied stress cycles (loads).
Thus, the discrete forms of lifetime distributions are proposed for estimating the reliability model. The reliability estimation approach developed in our previous studies [6,7] is generic and can be used with a discrete lifetime distribution, but the Bayesian method alone may not be enough to estimate the reliability model's parameters, and further elaboration is needed.
The reliability model of a discrete distribution (such as the Weibull Type III) might have a summation term with a time variable as its upper bound [8,10]. In the reliability estimation approach proposed in our previous studies [5,6], this upper bound becomes an equivalent cycle, which is a function of one or more unknown parameters. This makes calculating the summation term and consequently estimating the parameters difficult. Therefore, this study suggests initializing the upper bound and updating it using MLE and a gradient descent algorithm until the upper bound has converged. The upper bound is then fixed, and Bayesian analysis is used to estimate the parameters of the reliability model as discussed in previous studies.
This paper also presents a case study that shows the application of the proposed approach using simulated survey and test datasets for a consumer product with the failure mode of cracking caused by accidental drops. It is assumed that the equivalent cycles of the devices follow a discrete Weibull distribution type III. Then, the parameters of the reliability model are estimated using the proposed approach. The final reliability model of the product is compared with the reliability model estimated in the previous study using the continuous Weibull distribution [5].
The rest of this paper is organized as follows: Section II explains the reliability estimation approach assuming a discrete lifetime distribution. Section III presents the case study. Finally, Section IV concludes the paper.

II. THE RELIABILITY ESTIMATION APPROACH
This section describes the proposed approach for estimating the reliability model of a new consumer product with a few life cycles using reliability test data of the new product coupled with user survey data and reliability test data of a similar product (e.g., the older versions of the new product). The current approach extends the approach proposed in our previous study by assuming a discrete lifetime distribution [6]. The outline of the approach is shown in Figure 1. The approach has three steps: (1) calculating the equivalent cycles, (2) removing bias from the survey responses, and (3) estimating the parameters of the reliability model of the new product using a sequential Bayesian approach. The steps are described in the following sections.

A. CALCULATING THE EQUIVALENT CYCLES
Users of consumer products cycle their devices under various stress profiles (i.e., a combination of stress and stress adjustors). The amount of damage after each cycle depends on the stress profile of that cycle. If only the number of cycles is used for reliability assessment (i.e., a simple cycle-to-failure analysis), the effect of different stress profiles will be disregarded. A reference stress profile is selected to consider the actual stress profiles in the reliability analysis. Then the number of cycles under various stress profiles is converted into the number of cycles under the reference profile through the stress-life model of the product. The number of cycles under the reference stress profile is called the equivalent cycles. Figure 2 schematically shows the concept of equivalent cycles. Assume a user who drops his handheld electronic device various times under various stress profiles. The equivalent cycles are the number of cycles under the reference stress profile that causes the same damage as the actual cycles under actual stress profiles.
Converting the actual cycles into equivalent cycles through the stress-life model of the product requires quantitative stresses and stress adjustors.
Nevertheless, surveys of device users usually collect subjective data. Therefore, our approach scores each actual usage cycle's stresses and stress adjustors and combines the values through an additive or a multiplicative stress-index (SI) model to find a unique SI value for each usage cycle.
Actual Use Conditions Reference Condition  This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication.  (1) to (4), where is the j-th stress value of the i-th device and is its corresponding weight, and is the k-th stress adjustor value of the i-th device and is its corresponding weight. The weights of the additive models are unknown and must be estimated with the other parameters of the reliability model in a Bayesian analysis. A decision on the choice of the function depends on the nature of the data and can be made judgmentally or formally optimized. More information about the SI models is provided in [5]. Then, the approach selects an appropriate stress-life model and calculates the equivalent cycles. The process of calculating the equivalent cycles is shown in Figure 3. The calculated equivalent cycles at this step are parametric because they contain the unknown parameters of the stress-life model. These parameters are later estimated through our proposed approach.

B. ESTIMATING AND REMOVING BIAS FROM USERS' RESPONSES
When the equivalent cycles of the surveyed devices are calculated, the results are biased (that is, inaccurate) because the users' responses to the applied stresses, stress adjustors, and the number of usage cycles are biased. To measure the amount of bias in the user responses, the lifetime distribution of a similar product (e.g., an older version of the new product) is estimated twice: the first time using its user survey data, and the second time using its reliability test data. The reliability test data is assumed to be the most accurate representative of life because the test conditions are recorded in the laboratory. Thus, the differences between the two life distributions occur because of the bias in the survey responses, and, if the bias is removed, the distributions will become similar. We assume that, if all equivalent cycles in the life distribution of the surveyed devices are multiplied by ψ (i.e., the bias value), the shifted distribution will be close to the life distribution of the tested devices. This is schematically shown in Figure 4. The distance between the shifted life distribution of the surveyed devices and the life distribution of the tested devices is then calculated using the KL divergence method and minimized by using a gradient descent algorithm to adjust the value of ψ. The KL divergence is a distance or similarity measure between two probability distributions [12]. The KL divergence equation is shown in (5) [13], where ∆ is the measured distance, illustrates the life distribution of the tested devices, ̅ is the vector of the parameters of , is the reference SI value, shows the life distribution of the surveyed devices, and represents the equivalent cycles of the devices. We suggest using a portion (x%) of the survey and test data about the similar product to estimate the ψ value and using the remaining portion (1 -x%) in the sequential Bayesian analysis to build the prior distribution (i.e., the intermediate  This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication.  posterior distribution in Figure 1) for the parameters of the reliability model of the new product.

C. ESTIMATING THE PARAMETERS OF THE RELIABILITY MODEL OF THE NEW PRODUCT
The parameters of the reliability model of the new product are estimated using a sequential Bayesian analysis [10]. The analysis involves three steps, as shown in the box at the bottom of Figure 1. All steps assume a parametric discrete life distribution for the equivalent cycles.
In the first step, non-informative or weakly-informative prior distributions are selected for the parameters of the reliability model. Next, the equivalent cycles of the second portion of the surveyed devices are multiplied by ψ, which now is a known value. This removes the bias from the equivalent cycles. The adjusted (biased-removed) equivalent cycles are then used as observations (i.e., likelihood data) in the first Bayesian analysis. The primary joint posterior distribution of the parameters is then estimated.
The primary joint posterior distribution is used as the prior distribution for the second Bayesian analysis in the second step. Then, the equivalent cycles of the second portion of the tested similar devices are used as observations, and the intermediate joint posterior distribution is estimated.
In the third step, the intermediate posterior distribution is used as the prior distribution, and the equivalent cycles of the new device are used as the observation. Finally, the final joint posterior distribution of the parameters is estimated in the Bayesian analysis and used to approximate the reliability model of the new product.
When assuming continuous life distributions, Bayesian analysis is straightforward. But for some discrete life distributions, the analysis is complicated. For instance, in the probability mass function (pmf) of a discrete Weibull distribution type III, shown in (6) [4], the upper limit of the summation is a time value, which, in our approach, becomes an equivalent cycle that depends upon the unknown parameters of the stress-life model, which are unknown. This makes calculating the summation term and completing the Bayesian analysis difficult. In the following, we propose a method to deal with this difficulty.
The proposed method utilizes a gradient descent algorithm to fix the upper bound of the summation term. Gradient descent is an optimization technique used to estimate the parameters of a linear model by minimizing its cost function [14]. The cost function measures the difference between a variable's predicted and true value, and these differences are used to update the estimates in each iteration. For example, consider a linear model with the slope of W and the intercept of b (i.e., ( ) = . + ). The gradient descent estimates for W and b are shown by (7) and (8), respectively, where l indicates the number of iterations in the gradient descent technique, is the learning rate, and ( ) are the variables of the model, is the true value of ( ), and W and b are the parameters [14].  (7) and (8). Each iteration calculates the upper-bound of the summation term of (6) using the latest gradient descent estimates and, through the MLE method, estimates the true value of ( ) (i.e., ). This process continues until the upper bound has converged. Then the upper bound is fixed, and the other parameters of (6) are estimated in the sequential Bayesian analysis. The outline of the process of estimating the upper bound is shown in Figure 5.   This section shows the application of the proposed reliability estimation approach using simulated reliability data for an electronic product (e.g., a laptop) with the failure mode of cracking caused by a few (e.g., 10 to 20) accidental drops. The approach for simulating the test data of the new and old products and the survey data of the old product is explained. Then, the proposed reliability estimation approach is used to estimate the reliability model of the new product.

A. THE SIMULATED DATASETS
We simulated an accelerated life test (ALT) dataset for the older version of the new product. This dataset has 100 samples (80 failed and 20 were right-censored). It includes information about the applicable stress (i.e., drop), stress adjustors (i.e., drop height, surface type, and the severity of drop), and the number of cycles to failure or right-censoring. This simulated ALT dataset is based on the following assumptions: 1. During the reliability test, samples were dropped from 0.5 m, 1 m, 1.5 m, and 2 m heights on a hard surface during strenuous or harsh activity.
2. Stress and stress adjustors were scored and combined through a multiplicative SI model. The scores for heights were 25, 50, 75, and 100, for the hard surface was 100 and for harsh and strenuous activity was 100.
3. The stress-life model was a known inverse power law (IPL), as shown in (9) [15], where N illustrates the number of cycles, SI is the stress-index value, A is 285, and n is 0.5. = − (9) 4. The reference stress profile was a drop from 2 m on a hard surface during a harsh activity.
5. The life distribution was a known discrete Weibull distribution type III, as shown in (6), where c=0.05 for the old product, c=0.01 for the new product, and β=0.5. Figure 6 shows the process of simulating the ALT dataset for the old product. We used the same approach for the new product but assumed that it was a better version of the old product and that its failure had been terminated (i.e., has greater cycles-to-failure).
We also simulated a user survey dataset for the old product. The survey dataset considered 500 users from different age groups with different heights, activities, and bias in the frequency of drops and the applicable stress profiles. This dataset is based on the following assumptions: 1. User dropped their devices under one or several stress profiles. For example, 1) knee height-semisoft surface-benign activity, 2) knee height-hard surface-harsh activity, 3) waist height-semihard surface-benign activity, 4) waist heighthard surface-harsh activity, 5) chest height-semisoft surfacebenign activity, 6) chest height-semihard surface-harsh activity, 7) head height-soft surface-hard activity, 8) head height-semihard surface-harsh activity.
Estimate the number of drops: Bias data: 2. Stress and stress adjustors were scored and combined through a multiplicative SI model. The scores for knee, waist, chest and head or higher heights were 25, 50, 75, and 100 for soft, semi-soft, semi-hard, and hard surfaces were 25, 50, 75, 100, and for benign and harsh activity were 50 and 100, respectively.
3. The stress-life model was a known inverse power law (IPL), as shown in (9), where N illustrates the number of cycles and SI is the stress-index value, A is 285, and n is 0.5.
4. The reference stress profile was a drop from the head or higher height on a hard surface during a harsh activity.
5. The survey data was initially generated from the discrete Weibull life distribution shown in (6). Then, the SI value and number of drops in the dataset were multiplied by random values taken from a N (0.80, 0.1) distribution to make the dataset biased. Taking random values from the normal distribution made it possible to assign different bias values to the users' responses (for the stress profile (SI value) and the number of drops). This bias value of each user remained constant during the ownership time. The process of simulating the user survey data is shown in Figure 7.

B. Estimating the reliability model of the new product
This section assumes that the stress-life model, life distribution, and reliability model of the old and new products are unknown. Therefore, we used the simulated datasets and the developed reliability estimation approach to estimate the reliability model of the new product. As discussed earlier, estimating the reliability model of the new product has three steps.
In the first step, the equivalent cycles of all surveyed and tested devices were calculated parametrically. For this purpose, we linearly scored the stress adjustors of the users' stress profiles and combined them through a multiplicative SI model to calculate the values of the devices. Also, the reference stress profile was selected as a drop from the head or higher height on a hard surface during a harsh activity. Its SI value was calculated through the multiplicative SI model ( =1000). Because drop was a mechanical stress (leading to damage), we assumed an IPL stress-life model and calculated the ratio between the number of equivalent drops and the actual number of drops, as shown in Figure 3. This resulted in (10), where ν is the equivalent cycle of the i-th device, d is the actual number of drops of the i-th device, is the SI value of the i-th device, and is the SI value of the reference stress profile equal to 1000.
In the second step, we removed bias from users' responses. We split the survey and test dataset of the similar product into two parts (30% in the first part and 70% in the second part [16]) and used the first portion of the two datasets to estimate the survey data bias. The parameters of the life distribution (i.e., the distribution of the equivalent cycles) of the tested and surveyed devices were estimated in an MLE analysis. The estimated parameters for the tested similar devices were = 0.43, = 0.046, = 0.50 and for the surveyed, similar devices were = 0.5, = 0.071, = 0.49. All equivalent cycles of the surveyed devices were then multiplied by an unknown number (i.e., ) to shift their life distribution. The distance between the shifted life distribution of the surveyed devices and the life distribution of the tested devices was then calculated using the KL divergence method, (5), and minimized through the gradient descent algorithm. This resulted in the bias value of =1.32, which was close to the inverse of the average bias originally introduced to the survey dataset ( =1.397). According to (12) and because the bias values for d and were randomly selected from (0.80, 0.1), the original mean bias value was 0.80 × 0.80 0.5 , which is approximately equal to 0.716, and its inverse is 1.397. To validate the estimated value, we multiplied all equivalent cycles of the surveyed devices by 1.32. We used the original value of n (n = 0.5) to estimate the numerical value of the equivalent cycles. Then, the reliability of the devices in the bias-reduced surveyed dataset was estimated using the Kaplan-Meier method, and we compared that with the true reliability model. As shown in Figure 8, the bias-reduced survey dataset follows the true reliability model (i.e., the solid black line). This provides evidence that the KL divergence method can remove the bias from users' responses. The estimated bias value was then multiplied by the equivalent cycles of the second portion of the surveyed devices to remove its bias.
In the third step, the second portion of the survey data (biasreduced), the second portion of the test data of the similar device, and the entire test data of the new device were used in a sequential Bayesian analysis to estimate the reliability model of the new product. It was assumed that the equivalent cycles of all datasets followed discrete Weibull distribution type III. The likelihood of the discrete Weibull distribution type III for an incomplete dataset (i.e., a dataset containing failed and right-censored devices) is shown in (11), where c and β are the distribution parameters, is the equivalent cycles of the i-th device, r is the number of failed devices, and m is the total number of devices.
(11) As discussed earlier, the upper bound of the summation term in the likelihood function was unknown and needed to be estimated before performing the Bayesian analysis. We distinguished between the in the upper-bound of the summation and in the exponential term of the likelihood function; the in the upper-bound was replaced by , where subscript g stands for guess, l shows the number of iterations, and i represents the device's number. The process of estimating is shown in Figure 9. This process had five steps. First, the equation of the equivalent cycles was linearized. Second, an initial guess for the parameter of the equivalent cycle was made. This parameter is shown by in Figure 9. Third, the equivalent cycles of devices were calculated using (12), where shows the actual number of drops for the i-th devices, is the SI value for the i-th device, is the reference SI value, and is the estimated parameter of the model after l iterations. Fourth, we put in the upper bound of the summation term and estimated the other parameters of the likelihood function through the MLE method. It was assumed that the MLE estimates were the true values of the parameters.
In the fifth step, we calculated ∆ through the gradient descent algorithm, selected a learning rate of 0.1, and updated the value of . The new value of was plugged into (12) again and the process of updating was continued until convergence. Then, the value of n and consequently the upperbound was fixed, and sequential Bayesian analysis was performed to estimate the joint distribution of the parameters of the reliability model of the new product.
The sequential Bayesian analysis consisted of three regular Bayesian analysis. The first Bayesian analysis assumed weakly informative prior distributions for the parameters of the reliability model. These distributions were Normal distributions with the mean values estimated by the MLE method and standard deviations that were half-normal distributions with a mean of zero and standard deviation of 10% of the MLE estimates. The second part of the biasreduced survey data was used as the observation in the first Bayesian analysis, and the primary posterior distributions of the parameters were estimated. The kernel density estimation (KDE) method was used for non-parametrically estimating the joint distribution of the parameters [16,17]. In the second Bayesian analysis, the primary joint posterior distribution of the parameters was used as a prior distribution. The second portion of the test data is about a similar product used as the likelihood data. The intermediate posterior distributions of the parameters were estimated, and the corresponding joint distribution was approximated using the KDE method. Finally, the third Bayesian analysis used the intermediate joint posterior distribution as the prior distribution and the test data of the new product as observation and estimated the final posterior distributions of the parameters. Figure 10 shows the estimated reliability models, the mean, and the true reliability model. The true reliability model is within the uncertainty region of the estimated models, and it is very close to the mean reliability model. Therefore, the proposed reliability estimation approach could adequately estimate the true reliability model.
We evaluated the sensitivity of our estimated reliability model to the selected reference stress profile. First, the reference stress profile was changed to a drop from a person's head or higher height on a hard surface during a benign (regular) activity. Then the new product's reliability model  was estimated. No noticeable change was observed in the estimated reliability model. We also evaluated the sensitivity of the reliability model to the selected SI model. An additive SI model replaced the multiplicative SI model, and the proposed approach estimated the reliability model. The true model still was reasonably inside the region of the estimated reliability models. However, the maximum difference between the true reliability model and the mean reliability model over time was about 0.1 higher than the maximum difference obtained using the multiplicative SI model. This was because the additive SI model had three more parameters (the weights of the three stress adjustors) than the multiplicative SI model, which resulted in more uncertainties.
We also estimated the reliability model of the new product using our previous approach, which assumes a continuous Weibull distribution [5]. Figure 11 shows the result of the continuous analysis, which used the same survey and test datasets. The difference between the continuous mean reliability model and the true reliability model and the difference between the discrete mean reliability model and the true reliability model over equivalent cycles are shown in Figure 12. As shown in this figure the maximum difference between the discrete mean reliability model and the true model is 2.5%. This difference is more than three times lower than the maximum difference between the continuous mean and true reliability models. Therefore, the discrete life distribution model resulted in a better reliability estimate and is a better model for products with few cycles to failure.
This study assumed a single failure mode for the hypothetical product. In the presence of multiple failure modes, the proposed approach should be used to estimate the FIGURE 11. The estimated reliability models (continuous distribution analysis) and the True reliability model. reliability model of each failure mode. Then, the total reliability model is measured as the events that the system is reliable under all failure modes. Thus, one should multiply the reliability models of all failure modes and estimate the total reliability model.

IV. CONCLUSION
Previous studies used the concept of equivalent cycles and estimated the reliability model of a new consumer product with large cycles-to-failure through a continuous lifetime distribution model using user survey data and reliability test data. This paper introduced and applied the concept of equivalent cycles and developed a mathematical method to estimate the reliability model of a new product with small cycles-to-failure through a discrete lifetime distribution model using user survey data and reliability test data.
This paper showed that a user survey can be a cost-effective and quick way to collect field data for estimating the reliability model of a product having a few cycles-to-failure. Besides, it was shown that using a discrete life distribution can more accurately estimate the reliability model of a product having a few discrete cycles-to-failure than using a continuous distribution. However, the discrete distribution analysis is computationally more involved because it requires additional steps related to the gradient descent algorithm and MLE.
The case study assumed a single failure mode. In the presence of multiple failure modes, one should multiply the reliability models of all failure modes and estimate the total reliability model. However, this case was not investigated in this study. This paper suggested using x% of the survey and test data of a similar product for estimating the bias value and This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and the remaining 1-x% for building prior distributions for the parameters of the reliability model of the new product. However, determining the value of x is out of the scope of this study and should be done through an optimization analysis in the future. The proposed approach is generic and can estimate the reliability model of a wide range of consumer products such as laptops and monitors.