Frontiers in Computing and Intelligent Systems ISSN: 2832-6024 | Vol. 12, No. 1, 2025 199 Study on Parameter Estimation of Load‐Sharing Parallel Systems under Log‐Logistic Component Lifetimes Shu Xiang Tianjin University of Commerce, Tianjin, China Abstract: This article addresses the problem of parameter estimation in load‑sharing parallel systems with component lifetimes following log‑logistic distribution. By formulating the joint likelihood function for system lifetimes, we obtain robust estimates via maximum likelihood estimation based on a Gauss-Seidel iterative optimization algorithm and then employ bootstrap techniques to construct confidence intervals for the parameters, thereby quantifying the uncertainty of the estimates. Simulation studies demonstrate the high accuracy and stability of the proposed method in handling complex dependency structures, which provides a solid theoretical foundation for the reliability analysis of load-sharing mechanisms in complex engineering systems. Keywords: Load-sharing Parallel System; Log-Logistic Distribution; Parameter Estimation; Gauss-Seidel Optimization Algorithm; Bootstrap Method. 1. Introduction Load-sharing mechanisms are crucial in reliability analysis for modern engineering systems. When a component fails, the remaining components dynamically redistribute the load, which significantly affects the system's failure rate. These mechanisms are common in aerospace, power transmission, textile engineering, and materials science. For example, Rosen (1964) noted non-monotonic failure rates in fiber composite parallel systems after a single fiber failure. Carlson and Kardomateas (2005) linked the fiber bundle model in textile engineering and fatigue crack propagation in materials science to load-sharing mechanisms. Amadi (2016) found that transformer bank failures in power systems cause dynamic load redistribution. Traditional reliability models often assume that component lifetimes have constant or monotonic failure rates. However, these assumptions are inadequate for capturing the complex, non-monotonic failure rate characteristics observed in load- sharing systems. For instance, Lawless (2003) highlighted that the exponential distribution's constant failure rate assumption cannot model the dynamic, non-constant failure rates in load-sharing systems. Khan (2018) highlighted that the Weibull distribution's limitation with its monotonic hazard rate function assumption. Li et al. (2020) mentioned that the Kumaraswamy distribution's restricted parameter ranges limit its flexibility, despite describing non-monotonic failures. Zhong (2024) observed that the Gompertz distribution's high sensitivity to parameters will cause unstable estimation results. Thus, there's an urgent need for new reliability models which are flexible in mathematics and have practical applicability in reliability engineering. The log-logistic (LL) distribution has attracted attention in reliability analysis for its explicit probability density function and non-monotonic hazard function. Al-Shomrani et al. (2016) employed Markov chain Monte Carlo (MCMC) methods to demonstrate the LL distribution's superiority in survival data analysis, particularly in characterizing a two-stage process of accelerated failure followed by mitigated damage accumulation. Kariuki et al. (2024) introduced the extended log-logistic (ELL) distribution, which overcomes traditional LL distribution shape constraints through a synergistic mechanism of shape and scale parameters, thereby better accommodating complex survival data characteristics. Although the above studies do not directly target load-sharing systems, their theoretical contributions offer valuable insights for modeling such systems with non-monotonic failure rate characteristics. Maximum likelihood estimation (MLE) is key for analyzing load-sharing systems, but Zhou et al. (2021) found that MLE can produce biases and convergence issues with small sample sizes. To address these limitations, researchers have proposed various strategies. Kong and Ye (2016) improved k-out-of-n system reliability assessment by using interval estimation but assumed known load allocation rules, limiting engineering applicability. Park et al. (2020) proposed an improved EM algorithm framework, decomposing complex likelihood functions to enhance parameter estimation stability. Chang et al. (2019) integrated MCMC methods with Bayesian inference, significantly improving estimation accuracy while considering system uncertainties. However, balancing computational efficiency with estimation accuracy in complex load-sharing systems remains a pressing challenge. In response to the limitations of existing models in characterizing load-sharing mechanisms, this study introduces the log-logistic distribution with a non-monotonic hazard function to model the failure behavior of components in parallel systems. The Gauss-Seidel iterative optimization method is incorporated, which ensures reliable modeling and robust parameter estimation under complex dependency structures. The effectiveness of the proposed method is demonstrated through simulation experiments. The paper is organized as follows: Section 2 establishes a load-sharing parallel system model based on the LL distribution and derives the joint likelihood function for system lifetimes. Section 3 presents a parameter estimation method combining maximum likelihood estimation with the Gauss-Seidel iterative method and analyzes algorithm convergence. Section 4 evaluates the estimation performance of MLE under different sample sizes through simulation studies and examines its stability by constructing confidence intervals using the bootstrap method. Section 5 presents the conclusion. 200 2. Model Construction 2.1. Model Description To develop a load-sharing parallel system model and estimate its parameters, the following assumptions are made. 1. The system comprises n independent and identically distributed (IID) load-sharing parallel subsystems, each consisting of k components. All components collectively bear a constant total load, which is uniformly distributed among the operational components. These n subsystems are constructed under identical experimental conditions and are independent of each other. Within each subsystem, the k components are interconnected through load sharing, and all components follow the same log-logistic distribution. 2. The independent and identically distributed systems are tested, and their failure times are recorded. Let tij (for i=1, 2, ..., n; j=1, 2, ..., k) denotes the time interval between the failure of the (j-1)th component and the jth component in the ith system. 3. When a component in the system fails, the load it previously bore is redistributed among the remaining operational components, leading to an increase in the hazard rate of those surviving components. 4. The lifetime of each component initially follows a log- logistic distribution with parameters α1 and β1. As components fail successively, the parameters of the lifetime distribution for the remaining components are dynamically updated. For example, after the first component fails, the parameters change to α2 and β2; after the second failure, they become α3 and β3, and so on, until the last component is characterized by parameters αk and βk. To facilitate the construction of the likelihood function and parameter estimation in subsequent analysis, this paper adopts the survival function of the log-logistic distribution 1 ( ) = , > 0, > 0, > 0, 1+ ( )β F x x α β αx (1) and the probability density function     1 2 ( ; , ) = , > 0, > 0, > 0. 1+ β β βα α f x α β x α β αx x      (2) Here, the parameters and represent the scale and shape parameters, respectively. 2.2. Likelihood Construction Based on the above assumptions, we construct the likelihood function for parameter estimation. First, consider the case of the first failure time in a single system. Since all components are independent, the earliest failure time among the k components is the minimum of their lifetimes, and its survival function is  1 1 1 =1 (min ..., > ) = ( > ) = [ ( ; , )] ., k k k j j P F α βX X x P X x x Taking the derivative of the above survival function yields the density function of the first failure time 1 min 1 1 1 1( ) = [ ( ; , )] ( ; , ).kα βf x k F x α βf x Thus, the likelihood contribution of the first failure time ti1 in the i-th system is: 1 1 1 1 1 1 1 1 1 1 1( ) ; )] ; ).[ ( (k i i iL , t = k F t fα β β α, ,βtα  Similarly, for the j-th failure, the number of remaining components is k−j+1, and the density function of the failure time can be expressed using parameters (αj, βj). Its likelihood contribution is [( ) )] ).( 1) ( ; ( ;k j j ij ij jj j j jij jL , t = k j + F αt , f ,α β α β βt where j=1, 2, …, k. Then, let α = (α1, α2, …, αk), β = (β1, β2, …, βk), both belonging to the space (0,∞)k, t = {tij, i=1, 2, …, n; j=1, 2, …, k}, and αj>0, βj>0. For a single system, the overall likelihood function is the product of the likelihood contributions at each stage, i.e.,             =1 2 1 2 1 1 1 ( ) = ( +1)[ ( ; , )] ( ; , ) 1 = , ! 1 1 ! 1 j j j j j j j j j β j j j ij β j ij β j j j ij β k k k j ij ij j j k β j= j ij k k j j + j= j i L . α β α β β α α t k j F t f t k + α β t α t+ k α α α t+ t =                           βα t For n independent systems, the joint likelihood function is           2 1 1 1 1 1 , ! [ ( ( ; ! . ] 1 ; ) ) j j j j j jij β j j j ij β j n k n k j i j i= j= n k n k i j = j + i j= L α β α β= k F t , f t , = k β + α α t α t           t  3. Parameter Estimator Taking the logarithm of the above formula yields the log- likelihood function as follows                1 1 1 1 ( , ) = ln ! + ln ln 1 ln 2 ln 1 j n k j j j j ij i= j= n k β j ij i= j= l n k α + β + β α t k j + + α t .                α t To obtain the maximum likelihood estimators of (α, β), we differentiate the log-likelihood function ( , )l α t , resulting in the following likelihood equations with respect to the parameters αj and βj    1 2 0 ( , ) 1 j n j j β i=j j j j ij nβ β = k j + α α α + α t l             α t (3)        1 1 ln ln 2 0 ( , ) 1 j n n j ij j ij β i= i=j j j ij α tn = l + α t k j + β β + α t         α t (4) 201 It is evident that an explicit solution for the parameters (α, β) is not attainable. Therefore, the above likelihood equations must be solved numerically by using appropriate iterative techniques for the 2 parameters. As described in Ortega & Rheinboldt (2000), the Gauss-Seidel method is particularly suitable for the likelihood equations in this context. The specific steps are as follows 1. For j=1, 2, …, k, choose appropriate initial values αj (0), and substitute them into the equation (4) to solve for βj, denoted as βj (1); 2. Substitute βj (1) into the equation (3) to solve for αj , denoted as αj (1) ; 3. Repeat steps (1) and (2) until αj and βj (for j=1, 2, …, k) satisfy the convergence condition:         1 1 ,m+ m m+ m j j j jmax α α , β β <   where ϵ is a pre-set convergence threshold, typically chosen as ϵ=10-6 or smaller. Next, we further derive the asymptotic distribution of the proposed MLE estimators. For this purpose, we need to compute the second-order partial derivatives of the log- likelihood function with respect to αj and βj, denoted respectively as:      , , 1, ,2 2 2 2 2 j jj j , , , j = ...,k. α l β l βα l      α t α t α t   However, due to the analytical complexity of these second- order partial derivatives, it is difficult to derive their explicit forms. Therefore, this paper only describes the procedure. According to the asymptotic theory of maximum likelihood estimation, under standard regularity conditions (such as consistency and non-singularity of the information matrix), the MLE of the parameter vector Θ=(α1, β1, α2 , β2,…, αk, βk) ′ denoted as  1 1 ˆ ˆˆ ˆ ˆ ,k kΘ α ,β ,...,α ,β  follows an asymptotic normal distribution as n→∞, i.e.,     1ˆ ,,0d 2kn Θ Θ N ,   α F where  ,α F is the Fisher information matrix. This matrix can be expressed in a block form               1 2 , 0 0 0 0 , , 0 , ,0 k                       α α α α     F F F F where       ,m m ijFα F is a 2×2 matrix for m=1, 2, …, k, with its elements defined as follows       11 12 21 22 , , . , , , 2 m 2 m 2 m m m m 2 m 2 m F = E α F = F = E α β F l = E β l l                          α t α t α t    Since the second-order partial derivatives mentioned above are difficult to obtain explicitly, the observed information matrix is commonly used in practical applications by evaluating the second derivatives of the log-likelihood function at the estimated parameter values ˆ ˆ m m m mα βα= , = β . Next, for simplicity, we study the MLE estimation of the log-logistic distribution when β1= β2=…=βk = 1. In this case, for n identically and independently distributed load-sharing parallel systems, the log-likelihood function is         1 1 ( ) ln ! ln 2 ln 1 . n k j j ij i= j= l = n k + α k j + +α t   α t Differentiate the above log-likelihood function with respect to αj, and set the derivative equal to zero. The resulting equation is    1 2 0 1, ) 1 ( n ij i=j j j ij tn = k j + = , j = ...,k. α α + α t l     α t (5) Obviously, it is not possible to obtain a closed-form solution for the parameter α directly from the likelihood equation mentioned above. Therefore, the Gauss-Seidel method mentioned earlier can be used to obtain its solution, and thus obtain the MLE of α . Based on the result of the MLE, we further construct a joint confidence interval. Under the large sample conditions, according to asymptotic theory, the MLE α̂ of the parameter vector α approximately follows a multivariate normal distribution. Its covariance matrix is provided by the inverse of the observed Fisher information matrix evaluated at α̂ , i.e.,      ,ˆ 0d 2kn N ,  α αα where    1 α αF . Therefore, under the large sample assumption, for a given confidence level 1-α (where0<α<1), the joint 100(1-α)% confidence region for the parameter vector α can be approximated as      1 ˆ ˆ 2 k,α- - , α α α α α where 2 k,α denotes the upper α quantile of the chi- squared distribution with degrees of freedom. This region is a confidence ellipsoid centered at the MLE α̂ , which visually reflects the joint uncertainty of the estimates αj. 4. Simulation Study In this section, we validate the effectiveness of the proposed parameter estimation method through simulation experiments. It is assumed that the lifetimes of the 202 components in the system follow a log-logistic distribution with a single unknown scale parameter. We set α1=0.1, α2=0.2 and α3=0.4, and generate 20 sample observations ( =20) from a system consisting of three components ( =3). The simulated data is shown in Table 1. To further assess the performance of the method under varying sample sizes, we repeat the experiment for =20, 50,75, and 100, and compute the parameter estimation bias and mean squared error (MSE) under each condition. The results are presented in Table 2. The experimental results indicate that both the bias and MSE of the parameter estimates decrease with increasing sample size, demonstrating that the proposed method achieves higher estimation accuracy under larger sample sizes. Table 1. Failure time samples i ti1 ti2 ti3 1 2.9880591 1.4940295 0.74701476 2 17.0741586 8.5370793 4.26853966 3 4.8628502 2.4314251 1.21571254 4 24.0615608 12.0307804 6.01539020 5 30.3834534 15.1917267 7.59586335 6 0.3782369 0.1891185 0.09455923 7 7.3155493 3.6577746 1.82888732 8 24.9571190 12.4785595 6.23927974 9 7.8938256 3.9469128 1.97345639 10 5.7552712 2.8776356 1.43881780 11 32.6378642 16.3189321 8.15946604 12 5.6903928 2.8451964 1.42259819 13 11.8367554 5.9183777 2.95918884 14 8.4538181 4.2269091 2.11345453 15 0.8972792 0.4486396 0.22431980 16 25.6964345 12.8482172 6.42410862 17 2.4512873 1.2256436 0.61282181 18 0.3481921 0.1740961 0.08704803 19 3.5563181 1.7781590 0.88907951 20 32.3017068 16.1508534 8.07542670 Table 2. Parameter estimation bias and mean squared error (MSE) 20 50 75 100 0.108405 0.117066 0.127464 0.106234 Bias 0.008405 0.017066 0.027464 0.006234 MSE 0.000071 0.000291 0.000754 0.000039 0.230243 0.196246 0.203168 0.210520 Bias 0.030243 -0.003754 0.003168 0.010520 MSE 0.000915 0.000014 0.000010 0.000111 0.262784 0.322616 0.365505 0.392723 Bias -0.137216 -0.077384 -0.034495 -0.007277 MSE 0.018828 0.005988 0.001190 0.000053 Furthermore, to analyze the correlation between parameters, we construct confidence ellipses for the parameter pairs (α1,α2), (α1,α3), and (α2,α3) under different confidence levels 1− ={80%,90%,95%,99%}, as shown in Figures 1–3. The results suggest a positive correlation among α1, α2, and α3. Fig 1. Confidence region for (α1, α2) After obtaining the parameter point estimates, the Bootstrap method is further employed to construct confidence intervals. By adopting both the boot-p and boot-t methods, the accuracy and uncertainty of the parameter estimates can be evaluated from different perspectives, thereby enhancing the reliability of the results (Efron & Tibshirani, 1986; DiCiccio & Efron, 1996). Fig 2. Confidence region for (α2,α3) Fig 3. Confidence region for (α1,α3) 203 4.1. Boot-p Algorithm 1. Based on the known parameter estimates, construct a sample dataset tij (where i=1, 2, …, n and j=1, 2, …, k) of size n. 2. Randomly draw a bootstrap sample * ijt of the same size n from the original dataset tij , and compute the corresponding bootstrap point estimate * 1 2ˆ ˆ ˆ ˆ( )* * * kα ,α ,...,α α . 3. Repeat the sampling and estimation process B times to obtain B sets of bootstrap point estimates * 1 2ˆ ˆ ˆ* * B, ,...,α α α . 4. Sort the B bootstrap point estimates to form an ordered sequence * * * (1) (2) ( )ˆ ˆ ˆ B  α α α . Based on the given confidence level 1-γ (0<γ<1), the 100(1−γ)% confidence interval for parameter αj is expressed as * * ( /2) ((1 /2) )ˆ ˆB B,   α α . 4.2. Boot-t Algorithm 1. Based on the known parameter estimates, construct a sample dataset tij (where i=1, 2, …, n and j=1, 2, …, k) of size n. 2. Randomly draw a bootstrap sample * ijt of the same size n from the original dataset tij , and compute the corresponding bootstrap point estimate 1 2 *ˆ ˆ ˆ ˆ( )* * * kα ,α ,...,α α . 3. For each bootstrap sample, calculate the pivot quantity * * *ˆ ˆ ˆ( ) / jj j jR     , where *ˆ j denotes the standard variance of the j-th bootstrap estimate. 4. Repeat steps (2) and (3) a total of B times to obtain B sets of pivot quantities * * * 1 2, ,..., BR R R . 5. Sort the B pivot quantities to form an ordered sequence * * * (1) (2) ( )BR R R   . Based on the given confidence level 1-γ (0<γ<1), the 100(1−γ) % confidence interval for parameter αj is expressed as     ˆ ˆˆ ˆ* * * * j j B / 21-γ / 2j jB α α- σ R , - σ R      . To validate the effectiveness of the proposed method and demonstrate its performance under different sample sizes, experiments were conducted using the R programming language. The system parameters were set as k=3, α1=0.1, α2=0.2 and α3=0.4, with sample sizes n selected as 20, 50, 75, and 100. By running the Bootstrap procedure 5000 times (B=5000), parameter estimates, standard errors (SE), and the widths of 95% Bootstrap confidence intervals were calculated for each sample size. The results are presented in Table 3. Table 3. Ninety-five percent bootstrap CIs Estimates n 20 50 75 100 Bootstrap estimates [SE] 0.1008206 0.1109668 0.1157372 0.1200843 0.030421 0.026519 0.025917 0.026912 0.1830094 0.1929322 0.1959856 0.1990104 0.036741 0.028989 0.027044 0.026512 0.2990239 0.3270942 0.3436778 0.3602818 0.094801 0.064075 0.054943 0.058538 Bootstrap-p CI {widths} ∈ 0.0563,0.1765 ∈ 0.0682,0.1680 ∈ 0.0730,0.1704 ∈ 0.0757,0.1745 {0.1202} 0.0998 0.0974 0.0988 ∈ 0.1114,0.2427 ∈ 0.1373,0.2428 ∈ 0.1461,0.2485 ∈ 0.1504,0.2483 {0.1312} 0.1055 0.1023 0.0979 ∈ 0.1946,0.5353 ∈ 0.2343,0.4694 ∈ 0.2618,0.4803 ∈ 0.2799,0.5026 0.1353 0.1066 0.0921 0.1284 Bootstrap-t CI {widths} ∈ 0.0317,0.1670 ∈ 0.0627,0.1693 ∈ 0.0868,0.1790 ∈ 0.0175,0.1459 0.1353 0.1066 0.0921 0.1284 ∈ 0.2165,0.3616 ∈ 0.1561,0.2470 ∈ 0.1445,0.2770 ∈ 0.1716,0.2725 0.1451 0.0909 0.1325 0.1009 ∈ 0.0080,0.3305 ∈ 0.1361,0.4348 ∈ 0.2858,0.4375 ∈ 0.2862,0.5021 0.3385 0.2987 0.1517 0.2160 From the results, it can be observed that as the sample size increases, the standard errors of the bootstrap estimates gradually decrease, and the confidence intervals become narrower. This indicates a significant improvement in the precision and stability of the parameter estimates. Notably, the boot-p method demonstrates more consistent and reliable estimation ranges across different parameters, further validating its robustness and applicability in this context. When combined with Table 2, the Gauss-Seidel method demonstrates smaller estimation bias and mean squared error (MSE) with larger sample sizes, indicating its suitability for precise point estimation. In contrast, the Bootstrap method quantifies estimation uncertainty through confidence intervals (CIs). While it exhibits greater fluctuation with smaller samples, its estimates stabilize as the sample size increases. In summary, the Gauss-Seidel method is better suited for high-precision estimation, whereas the Bootstrap method has a comparative advantage in evaluating estimation uncertainty. 5. Summary This study proposes a novel modeling approach for load- sharing parallel systems based on the log-logistic distribution and combines it with the Gauss-Seidel iterative optimization algorithm to achieve robust parameter estimation in complex dependency structures. By employing log-logistic 204 distribution to characterize the distribution of the component lifetimes, this paper overcomes the limitations of traditional methods in handling non-independent failures and non- monotonic hazard rate characteristics. Simulation results demonstrate that the proposed method significantly improves the accuracy and stability of parameter estimation and explores the impact of different sample sizes on estimation precision. The model not only enhances computational efficiency but also provides effective theoretical support for reliability analysis of load-sharing mechanisms in complex engineering systems. References [1] Rosen B. W. Tensile failure of fibrous composites. American Institute of Aeronautics and Astronautics, 1964, 2(11): 1985- 1991. [2] Amadi H. N., Izuegbunam F. I. Analysis of transformer loadings and failure rate in Onitsha electricity distribution network. American Journal of Electrical and Electronic Engineering, 2016, 4(6): 157–163. [3] Carlson R. L., Kardomateas G. A. An introduction to fatigue in metals and composites. Chapman & Hall, 1996. [4] Lawless J. F. Statistical Models and Methods for Lifetime Data (2nd ed.). Wiley, 2003. [5] Khan S. A. Exponentiated Weibull regression for time-to-event data. Lifetime Data Analysis, 2018, 24(2): 328-354. [6] Li Danqing, Xiong Wenjie, Zhang Zhengcheng. Parameter Estimation of Load Sharing Parallel System under Kumaraswamy Distribution [J] Journal of Anhui Normal University (Natural Science Edition), 2020, 43(04): 307-314. [7] Zhong Xin 'ao. Parameter Estimation of Load-sharing Systems under Gompertz Distribution and Exponential Pareto Distribution [D] Hainan Normal University, 2024. [8] Al-Shomrani A., et al. Log-logistic distribution for survival data analysis using MCMC. SpringerPlus, 2016, 5: 1774. [9] Kariuki P., et al. Properties, estimation, and applications of the extended log-logistic distribution. Scientific Reports, 2024, 14: 68843. [10] Zheng X., Chiang J.Y., Tsai T.R., Wang S. Estimating the failure rate of the log-logistic distribution by smooth adaptive and bias-correction methods. Computers & Industrial Engineering, 2021, 156: 107188. [11] Kong Y., Ye Z. Interval estimation for k-out-of-n load-sharing systems. IIE Transactions, 2016, 49(3): 344-353. [12] Park C., Wang M., Alotaibi R., Rezk H. R. Load-Sharing Model under Lindley Distribution and Its Parameter Estimation Using the Expectation-Maximization Algorithm. Entropy, 2020, 22(11): 1329. [13] Chang F., Zhe Y., Di S., Haifeng L., Caisheng W., Zhiwei W., Jie L. A classical and Bayesian estimation of a k-components load-sharing parallel system. IEEE Transactions on Reliability, 2019, 68(3): 1050-1062. [14] Ortega J. M., Rheinboldt W. C. Iterative Solution of Nonlinear Equations in Several Variables. SIAM, 2000. [15] Efron B., Tibshirani R. J. Bootstrap Methods for Standard Errors, Confidence Intervals, and Other Measures of Statistical Accuracy. Statistical Science, 1986, 1(1): 54-75. [16] DiCiccio T. J., Efron B. Bootstrap confidence intervals. Statistical Science, 1996, 11(3): 189-212.