CC BY
1
C.A. Hesse, D. B. Boyetey, E.D. Kpeglo, A.A. AshiagborEVALUATION OF PARAMETRIC ESTIMATION METHODRT&A,No1 (82) FOR THE GAMMA DISTRIBUTION...Volume 20, March2025
EVALUATION OF PARAMETRIC ESTIMATION METHODS FOR THE GAMMA DISTRIBUTION USING MAXIMUM LIKELIHOOD AND BAYESIAN APPROACHES IN A CENSORED LIFE-TESTING STRATEGY WITH MARKOV CHAIN MONTE CARLO SIMULATIONS
Christian AkrongHesse1,DominicBuerBoyetey2, Emmanuel DodziKpeglo3*Albert AyiAshiagbor4
�
1,3*Department of Economics and Actuarial Science, University of Professional Studies, Accra, Ghana2Department of Built Environment, University of Environmentand Sustainable Development, Somanya,Ghana4 Department of Business Management, University of Pretoria, South Africa1,[email protected], [email protected],3*[email protected]@up.ac.za
Abstract
The goal of this study was to address the computational challenges associated with parametric estimation of the gamma distribution by evaluating the performance of the maximum likelihood and maximum a-posteriori estimation methods within the framework of Markov Chain Monte Carlo simulations. This was done by first assuming a censored life-testing strategy that terminates on the rth failure from a given sample of n electronic devices. Second, we obtained the joint distribution function of the first r-order statistic by arranging the r values in order of magnitude. Finally, we explored through the Markov Chain Monte Carlo framework using the maximum likelihood and maximum a-posteriori to estimate the gamma distribution parameters. The findings of this study suggest that both estimation methods were not significantly different from the actual hypothesized parameter values. Further, we observed that irrespective of the prior distribution used for the Bayesian maximum a-posteriori Markov Chain Monte Carlo estimation, the resulting parametric estimates of the gamma distribution remain the same, confirming the assertion that the Bayesian maximum a-posteriori Markov Chain Monte Carlo approach is a valuable tool for informative posterior analysis. The study�s uniqueness lies in adopting a censored life-testing strategy centered on the joint distribution function of the first r-order statistic.
Keywords: bayesianinference, gammadistribution, maximumlikelihoodestimation,maximum a-posteriori,reliabilityanalysis
C.A. Hesse, D. B. Boyetey, E.D. Kpeglo, A.A. AshiagborEVALUATION OF PARAMETRIC ESTIMATION METHODRT&A,No1 (82) FOR THE GAMMA DISTRIBUTION...Volume 20, March2025
I.Introduction
The Gammadistributionhasbeenextensivelystudied and developedinstatisticalinference.Itispreferred over other probability distributions for its superior applications in insurance and finance.Dickson[1]andVeazie etal. [2]projectthe twoorthree parameter Gammadistributionfor itsability tomodel highlyskewedpositiveandornegativedatapoints.Gammadistributionapplicationismostevidentinstudyingrandomvariablessuchaswaitingtimes, claimsize orfrequencyandinvestmentreturns. While everyprobabilitydistributionofferssomedistinctadvantageinspecificcontexts, Deyetal.[3]support thegammadistributionforitscomputationalefficiencyandmemorylessproperty. Memorylesspropertymeansthatthe occurrencesofpasteventsdonotinfluencetheprobabilityoffutureoccurrencesofthe event(Noguchi&Robles[4]; Shore [5];Tao[6]).Inother words,the memorylesspropertyofthe gammadistributionmakesiteasytostudythe probabilityofthe occurrence ofaneventindependentlyofthe probabilityoffuture occurrencesof theevent.Thispropertydemonstratestheusefulnessandversatilityofthegammadistributioninseveral survival andreliabilitystudies.
Theuseofthegammadistributioninreliability analysis extends tootherfields suchas engineering, manufacturingand biomedicalresearch. For instance, inmodelingtime-to-failure ofmanufacturedcomponents,Elsayed[7]employedthegammadistributionalongsidetheMaximumLikelihood Estimation(MLE)method toestimate the parameters. Similarly, Shipesetal. [8]usedthePoisson-GammaModel intheirsurvival analysisoftime-to-eventclinicaltrialdata. These studiesunderscore the flexibilityofthe gammadistributionincapturingdiverse eventpatterns. More importantly,theparametricestimationofthegammadistributioninthesestudiesoffersroomforfurtherexploration. Meeker and Escobar[9]explainedparametricestimationinbothreliabilityandsurvival analysisasfittingaspecificprobabilitydistribution(e.g.,gamma,exponential,Weibull)totheobserved failure dataand estimating its parameters (location, shape and scale parameters).
Literatureaboundswithdifferent combinationsofestimationandsimulationmethodstoestimate the gammadistributionparameters. Severalstudies(Ghosh&Hamedani[10];Junmei& Liqin[11]) seemtoopt fortheMLEmethodduetoitsoptimalandconsistent parameterestimators.Ghosh andHamedani[10]provided detailed propertiesofthe two-parameter gammadistributionusingthe MLE method toinvestigate itsmoments, hazard functionand reliabilityparameters. Whenapplied toalifetime dataset, the gammadistributionproduced asuperior fitcompared toothermodels.Whilethegammadistributionisapplaudedforitsflexibilityinmodel fitting,Ozsoy,Unsal andOrkcu[12]warnedofapotentialcomputationalcomplexitywhenMLEisusedforitsparameterestimation.Theyexplainedthat thedistributionfunction(orsurvivalfunction)ofthegammadistributionisnotavailableinaclosedformiftheshapeparameterisnotaninteger,therebymakingthe use ofMLE anear futile exercise. Thisnotwithstanding, studies(Hamadaet al.[13];Rubinstein&Kroese [14])have adapted numericalmethodstoevaluate the parametersofthe gammadistributionbyexploringacombinationofBayesianestimationand simulationprocedures. Hamada etal.[13]foundtheBayesianestimationuseful intheirprobabilisticframeworkforreliabilityestimationasit incorporatesadditional informationaboutthedistributionknownasaprior.TheBayesianframeworkentailscarefulelicitationsofprior expertinformationtoenhance the data,leadingtoimproved predictionofextreme cases(Coles&Tawn[15]). Hussainetal.[16]andKoholeetal.[17]proffer theBayesianmethod, asitatleastoffersawayaroundthecomplexityoftherootofthemaximum likelihoodequation known toexistin MLE.The Bayesianapproach, therefore,appears more flexible and informative through its posterior analysis.
Inrecentdecades, the surge instatisticalapplications has sparkedagrowing interestinBayesianparametricsimulation, givingrise tothe efficientconceptofMaximuma-Posteriori(MAP).Servingasthe BayesiancounterparttoMLE,MAPestimationentailsidentifyingparametervalues
C.A. Hesse, D. B. Boyetey, E.D. Kpeglo, A.A. AshiagborEVALUATION OF PARAMETRIC ESTIMATION METHODRT&A,No1 (82) FOR THE GAMMA DISTRIBUTION...Volume 20, March2025
that maximizetheposteriordistributionandact asestimatesfortheunknownparameters(Hesseetal., 2016). Whenthere isanoninformative prior inthe Bayesiananalysis, the MAPestimate isthesameasthat oftheMLE.DuetothecomputationalintensityofMAPresultingfromthe incorporationofprior information, Hesse etal.[18]turnedtoMarkovchainMonteCarlo(MCMC)toobtainsamplesfromthe posterior distribution, enablingthe estimationofregressionparameters. The conceptofMCMC ispopularinfieldssuchasmanufacturing,physicsandfinance,anditusesprobabilitydistributionstomake selections(Benson&Kellner [19]).Inreliability assessment,Naess,Leira andBatsevych[20])noted thatthe MCMCcancheckfailure criterion, regardlessofthedistributionorsystemcomplexity.Fauzietal.[21])reliedontheMCMCalgorithms forsampling fromaposteriordistribution,essentially,tosimulate systembehavior and estimate reliabilitymetrics.
Thisstudy aims totacklethecomputational challengesassociatedwithparametricestimationofthe gammadistributionbyevaluatingthe performance oftwoestimationprocedures: MaximumLikelihoodEstimation(MLE) andMaximuma-Posteriori (MAP)estimation.Theseassessments will be conducted withinthe frameworkofMarkovChainMonteCarlo(MCMC)simulations.Inthis paper, the twoMCMC-basedestimationtechniquesaredenotedasMLE_MCMCandMAP_MCMC.Thestudy�s uniqueness lies intheadoptionofacensoredlife-testingstrategy,terminatingupontheoccurrence ofthe rth(where r < n)failure. Thisapproachdivergesfromclassicallife testing, whichrequiresthe complete failure ofalln samples.Thestudyconcentratesonthejoint distributionfunctionofthefirst r-order statistic, preciselythe smallestr values,asanalternativetoutilizingthecompletedatasetforestimation.Additionally,welookintothesensitivityofMAP_MCMCbyapplyingvariouspriordistributions.Thisstudyisrelevantasitillustratesthatthejointdistributionfunctionofthefirst r-orderstatisticprovesmore suitable for estimatingthe parameter(s)oftheprobabilitydensityfunction(pdf) of thetime-to-failure randomvariable for any engineered device.
II.Methods
Consider n samplesofmanufactured componentsthatwere subjected toreliabilitylifetestsfromacertainpopulationofinterest. The randomvariable T ofinterest isthetimeit takesuntilthecomponentfails.Supposethe underlyingfailure timesare....(1), ., ....(....) where ....(....) .....(....+1), .... = 1, ., .....1.Andlet ........(....) bethedistributionfunctionofT andlet........(....) beitsprobabilitydensityfunction(pdf).Assumingfurtherthat thereliabilitylifetestsconcludeat therthfailure, where r is lessthanorequalton, the number offailuresistreated asafixed value, while the failure timesare regarded asrandomvariables. We employthe gammadistributiontomodelthe time-to-failure randomvariablesin thisscenariooflifetesting,assumingthatthefailurerateisnotconstant.Thegammadistributionispreferred inthisinstance because itexhibitsafailure rate thatfollowsabathtub-shaped curve (decreasingfailure ratesatthe initialphase and increasingfailure ratesatalaterphase).Blanksby andLyons [22]assertthatthegammadistributionallowsforflexibilityincapturingdiverse failure rate behavioursand iswell-suited forscenarioswhere the hazardfunction varies over time. The next subsection presentsa synopsis of the gamma distribution.
I.GammaTime-to-FailureRandom Variable
Thecontinuousrandom variableT,is said to have the gamma distribution with parameters . > 0and. > 0if its pdfis given by:
.................1.............
........(....)= , ....>0 (1)
....(....).
where ....(....)= . .........1................. isthegammafunction.Thecumulativedistributionfunctionisthe
0
C.A. Hesse, D. B. Boyetey, E.D. Kpeglo, A.A. AshiagborEVALUATION OF PARAMETRIC ESTIMATION METHODRT&A,No1 (82) FOR THE GAMMA DISTRIBUTION...Volume 20, March2025
regularizedgamma function:
....(...., ........)
........(....)= , ....> 0. (2)
....(....) ........
where ....(...., ........) isthelowerincomplete gammafunction, thatis....(...., ........)= . .........1..................This
0
gammadistributionisinthetwo-parameterfamily ofcontinuous probability distributions.Theseare:
�
Shape Parameter (.):Thisparameterdeterminesthe shape ofthedistribution.It isapositive real number.
�
Scale Parameter (.):Thisparameter isassociated withthe rate ofevents. Itisalsoapositive realnumber. The densityand cumulativedistributionfunctionsare sometimesexpressed in terms of the scale parameter, .... =1/.....
Table 1 shows some characteristicsof the gamma distribution(Mannetal.[23]).
Table 1: Properties of the gamma distribution
.... ....
Variance
....2
MedianNosimpleclosed
.....1
Modefor..1
.... .....1
Reliabilityfunction(........)....
............. . ....!
....=1
Inthe nextsubsection, we exploretwoapproachesfor estimatingthe parametersofthe gammadistributionthroughtheMarkovChainMonte Carlo(MCMC)simulationframework: (a)MaximumLikelihoodEstimation(MLE) and(b)Maximuma-Posteriori(MAP).
II.ParametricEstimationoftheGammaDistribution
According toOfosuandHesse[24],thelikelihoodfunction L ofthefirstr-orderstatistics,....(1) . ....(2) .......(....),of the random variable ofinterest in this study canbe specified as:
....= ........(1),.,....(....)(....1, ., ........) ....
....!
[1 .........(........)]......... .........(........)
=
(.........)!
....=1
....
......... .....1.................
....! ....(...., ............) ................
= .1 . ... .
(.........)! ....(....) ....(....)
....=1 .... ......... ....
....! 11
............)]......... . ............ ...............1..................
=[....(....) .....(...., .. .
(.........)! ....(....) ....(....)
....=1 ....
....!1 .... ....
............)]......... . ............ ...............1...............=1........
=[....(....) .....(...., . .
(.........)! ....(....)
....=1
C.A. Hesse, D. B. Boyetey, E.D. Kpeglo, A.A. AshiagborEVALUATION OF PARAMETRIC ESTIMATION METHODRT&A,No1 (82) FOR THE GAMMA DISTRIBUTION...Volume 20, March2025
(3)
Thenaturallogarithm ofthelikelihoodfunction gives:
....
ln ....= ln ....! .ln(.........)! + (.........) ln[....(....) .....(...., ............)] .....ln ....(....)+ ........ln ....+(.....1) . ln .........
....=1
....
.........=1......... (4)Thisfunctionyieldsthefollowinglogarithmiclikelihoodequations:
.....................
....ln.... ........ .... .........1(.........)........
= .. ln ........ + =0 (5)
....=1
........ .... [....(....).....(...., ............)]
....ln.... (.........)..... ' (....)..... ' (...., ............). ........ ' (....) ....
= ++ ....ln ....+ . ln ........ =0 (6)
....=1 ........ [....(....).....(...., ............)] ....(....)
SolvingEquations(5)and(6)isnotablychallenging.When astraightforwardsolution tothelikelihood equationsiselusive, variousproceduresare available forMLE. CommonmethodsencompassIterative Methods, the Expectation-Maximization(EM)Algorithm,GradientDescent,Quasi-NewtonMethods, Monte CarloMethods, Profile Likelihood, Bootstrapping, and NumericalOptimization(Dempster et al.[25];Gilks etal.[26];Nocedal&Wright[27];Pressetal.[28]).In caseswhereobtainingasolution tothelog-likelihood equations proves to be difficult, we turn toMarkovChainMonte Carlo(MCMC)samplingtechniquestogenerate samplesfromthe likelihoodfunction. Consequently,weappliedtheMCMC estimation techniqueandreferredtoit asMLE_MCMC. The primaryobjective istodetermine parameter estimatesthatmaximize thelikelihood functiongiventhe sample data. The MLE_MCMCapproachidentifiesthe mode ofthe simulatedMCMCsamplefromthebivariatelikelihoodfunction in Equation (3),representingthepoint estimate of the parameter vector.... =(...., ....). Thatis,
....
............)]......... . 1 ............(..... .....1).......... ........
................ = arg .................[....(....) .....(...., . ........ ....=1 .... . (7)
....=1 ....(....)
ThefollowingalgorithmisthedescriptionforthemultivariateMetropolisHastingsprocedure (Steyvers [36]):
�
Set ....=1
�
Generatean initialvaluefor....~....(....1, ....2).
�
Repeat
....= ....+1
DoaMH stepon.,Generate a proposal..... ~....(...., ....2);
....(...../....)
Evaluatetheacceptanceprobability.... = .............1, .;
....(..../....)
Generate au fromaUniform(0,1)DistributionIf........., accept the proposal and set.... = .....
� Until ....= .....
Maximum APosteriori(MAP)estimationisthe BayesiancounterparttoMaximumLikelihoodEstimation(MLE),incorporatingadditionalinformation through thepriordistribution.Now,thejointpdf of ....(1), ., ....(....) and.... =(...., ....) is givenby....(....1, ., ........, ....)= ........(1),.,....(....)(....1, ., ........|....)....(....),
where ....(....) isthe priordistributionoftheparameter..Weassume. and. areindependentand1 .....(..../....+..../....)
exponentiallydistributedwithmeansa andb, respectively. Thus,....(....)= , .... >0,.... >
........
0.
C.A. Hesse, D. B. Boyetey, E.D. Kpeglo, A.A. AshiagborEVALUATION OF PARAMETRIC ESTIMATION METHODRT&A,No1 (82) FOR THE GAMMA DISTRIBUTION...Volume 20, March2025
.... ....
....(....1, ., ........, ....)= ....![....(....) .....(...., ............)]......... . 1 . ............ ..... ..............1................................+..... . (8)
....=1
........(.........)! ....(....)
Thus, the marginal p.d.f. of....(1), ., ....(....) is
.. ........(1),.,....(....)(....1, ., ........)= . . ....(....1, ., ........, ....)................
00
.... .... ....! ..
....+
= .. [....(....) .....(...., ............)]......... . 1 . ............ ..... ..............1............................ ......................
....=1
........(.........)! 0 0 ....(....)
(9)The conditional p.d.f. of., given....(1), ., ....(....), is therefore defined by
....(....1, ., ........, ....)....(....|....1, ., ........)=
........(1),.,....(....)(....1, ., ........) ....
....!1 ....
[....(....) .....(...., ............)]......... . ............ ..... ..............1................................+.....
....=1
........(.........)! ....(....).
=
....
....! .. 1 ....
............)]......... . ............ ..... .....1.................}..............+.....
....=1{........ ................ ........(.........)! .0 .0 [....(....) .....(...., ....(....).
.... .... ....
= ....[....(....) .....(...., ............)]......... . 1 . ............(.....=1 .... .............1 )..............+....+.........=1......... . (10)
....(....)
where K isindependentof. and..ThetypicalapproachinBayesianestimationistoemploytheposterior mean, ....(....|....1, ., ........), asa pointestimatefor. (Hesse etal. [18]).TheMaximumaposteriori(MAP)estimatorof. isthevaluethat maximizestheposteriordistribution.ThisstudyusestheMarkovChainMonte Carlo(MCMC)samplingapproachtodrawsamplesfromthe posteriordistribution. Thisspecificmethod ofestimation, denoted asMAP_MCMCfor the purpose ofthisstudy, identifiesthe mode ofthe posterior distribution, representingthe pointestimate for the parameter ..Thus,
.... .... ....
....+.....
................ = arg .............= ....[....(....) .....(...., ............)]......... . 1 . ............(.....=1 .... .............1 )..............+ ....=1........... (11)
....(....)
III.Results andDiscussion
Consider ascenariowhere acompanyspecializesinmanufacturingaspecificdevice, itscomponents,or equipment. The companyisdedicated toassessingthe reliabilityofthese items. Inthiscontext,reliabilitydenotestheprobabilityofadevicesuccessfullyfulfillingitsintendedfunction.Asample of n suchdevicesfromapopulationofinterestisplaced inanenvironmentcloselyresemblingthe conditions the itemswill encounter in actual use. One or more stresses of constant severityare thenapplied tosimulate real-world scenarios. Giventhat,200ofthese devicesare programmed toundergoreliabilitylifetestingwiththetest truncatingwhenthe100thfailed device wasobserved. The first 100 sampled time-to-failureunitswereorderedandfittedtothegammadistribution.
Given thatT isa continuous randomvariablewithadistributionfunction........(....), then, accordingtothe probabilityintegraltransformationconcept, ....= ........(....),followsthecontinuousuniformdistributionovertheinterval(0,1) (Ofosu&Hesse[24]). The inverse transform sampling method isemployedsubsequently togeneratesamples fromthegammadistributionusing thefollowing steps:
�
Generate a random numberu fromthestandarduniformdistributionintheinterval(0,1).
�
Find the inverse of the desired cumulative density functions (CDF), denoted as.....1(....).
�
Compute....= .....1(....). The computed values of the random variable X correspondto thedesireddistributionwithprobabilitydensityfunction........(....).
C.A. Hesse, D. B. Boyetey, E.D. Kpeglo, A.A. AshiagborEVALUATION OF PARAMETRIC ESTIMATION METHODRT&A,No1 (82) FOR THE GAMMA DISTRIBUTION...Volume 20, March2025
There isnoclosed-formexpressionforthe gammadistribution�sinversecumulative distributionfunction(CDF). However, variousnumericalmethodsand statisticalsoftware packagesare available tocalculate quantiles(inverse CDFvalues)for the gammadistribution. Commonlyemployed numericalmethodsinclude the Newton-Raphson MethodandBrent'sMethod(Burden&Faires [29];Kincaid&Cheney[30];Dahlquist&Bjorck[31]). Inaddition, other statisticalsoftware packageslike R, Python(SciPy), MATLAB and MS. Excelprovide functionstocompute quantilesforthegammadistribution(Eaton[32];Hanselman &Littlefield[33]).
We, therefore,usethe'gaminv'function in MATLAB(MathWorks[34])tocalculatequantilesforagivenprobability, simulated fromthe uniformdistributionover the interval(0, 1). Table 2displaysthe first100outof200ordered datapointssimulatedfromthegammadistributionwithparametersalpha (.) = 10 andbeta (.) = 0.05, resulting in a mean (�) of 200. Itis assumed that theseobservationsrepresentthe outcomesofareliabilitylifetestinvolving200devicesuntilthe failure ofthe100thdevice.
Table 2: Ordered data simulated from the gamma distribution
73.23076.64381.30885.70988.57890.10391.54591.88495.09197.130
98.56499.319101.788105.722107.333107.880110.057111.884113.007115.471
117.521121.779125.079125.766128.059128.199130.834132.564134.648135.505
135.923136.019138.582141.671142.111143.420144.675146.128147.288147.939
152.135153.413153.475155.564155.716156.338156.353157.133157.848157.903
159.273159.895159.900160.563160.563161.096162.658162.837168.294168.322
169.501169.539170.437170.723170.831172.767173.089173.164173.700175.040
175.844176.285176.485177.148177.296177.426177.463178.349178.942179.484
181.027181.703181.777182.658182.677182.807185.919188.967189.609189.759
189.848190.031191.716192.004192.051192.169192.448194.546194.670195.500
On the MLE_MCMC parameter estimates, we deduce from Equation (3), given thatn = 200,r = 100,........ = 195.500 fromTable 2, together withthe Metropolis-Hastings algorithm,as describedabove,were employedtodrawsamplesfromthe likelihood function. The MATLAB code for implementingthecomponent-wiseMetropolis samplerforthelikelihoodfunctionis providedinListings A1andA2, inthe appendix. The mode ofthisbivariate sample providesthe maximumlikelihood estimate for the parametersalpha(.) andbeta (.) ofthegamma distribution,whicharegivenby................ = 10.4169 and................ = 0.0487, respectively. Note thatthemaximumlikelihood estimatesare notsignificantlydifferentfromthe actualvalue ofthe parameters, thatis. =10 and. =0.05. Figure 1showsthegraph oftheactualandtheestimatedgammafunction with thegiven parameters.
The closenessofthese estimated parametersofthe gammadistributionisconsistentwithobservations madeby Sauloetal.[35], whoobserved thatthe generalized gammadistributiongenerated veryclose valuesofthe log-likelihoodfunction when comparedtotheDagumdistributioninanMLE procedure forcensored remissiontimesofcancer patients. The closenessofthe valuesaffirms thattheMLE_MCMCapproachis as trustworthy as theclassicalMLE.
C.A. Hesse, D. B. Boyetey, E.D. Kpeglo, A.A. AshiagborEVALUATION OF PARAMETRIC ESTIMATION METHODRT&A,No1 (82) FOR THE GAMMA DISTRIBUTION...Volume 20, March2025
Figure 1: pdf of the gamma distribution
Similartothe MLE_MCMCestimates, we utilize the datafromTable 2and the Metropolis-Hastingsalgorithmtosimulate asample fromthe Bayesianposterior distribution. We further assumedthat . and. are independentand exponentiallydistributed withmeansa = 15 andb = 0.08.The MATLAB code for implementingthe Metropolis-Hastingssamplerforthe posteriordistributionin Equation (9)isprovidedin ListingsA3andA4,in theappendixtoobtain theMAP_MCMCparameter estimatesofthegammadistribution.
The resultsrevealthat theMAP_MCMCestimatesof. and. are ................ = 10.4169 and................ = 0.0487 whichpreciselymatchthe resultsobtained whensamplingisdone directlyfromthelikelihoodfunction.Theseresultsremain consistenteven when varyingthevalues ofa andb,such asa = 20 andb =1. Inother words, the parameter estimatesofthe gammadistributionfromthe MAP_MCMCapproachproduce preciselythe same estimatesastheMLE_MCMCapproach. Theseeminglynodifference betweenthe twoestimationmethodsconfirmstheirflexibilityinresolvingestimationcomplexitiesassociatedwiththegammadistributionfortime-to-failure variables.
Further, we explore the sensitivityofthe MAP_MCMCestimator todifferentpriordistributionsbyconductingarepeatedMCMCsimulationassumingthat theparameters. and. are independentand followLognormal, Pareto, Weibull, and Gumbelpriors. The resultsofthe MAP_MCMCestimatefor. and. acrossthesefourpriordistributionsareconsistent withthoseobtainedusingtheexponentialprior.The implicationisthatirrespective ofthe distributionused, the outcome forthe twoestimationproceduresormethodsproduced thesame parameter estimates. We, therefore,concurwithpreviousstudies(Hussainet al.[16];Kohole etal.[17]toconcludethattheBayesianMCMCapproachoffersamore flexible andinformativeposterior analysis.Inaddition, ourestimationprocedure affirmsFauzietal. [21]assertionthatBayesianMCMCalgorithms forsampling fromaposteriordistributionessentiallysimulatesystembehaviorandestimatereliabilitymetrics.
IV.Conclusions
Theapplicabilityofthegammadistributiontostudythedistributionofrandomvariables,suchastime-to-failure ofanevent, hasproventobe effective andversatileinmanyreliabilitystudies.Thepreoccupationofthese studiescentersonthe parametricestimationand derivationofother propertiesofthe gammadistribution. The challenge, however, isthe computationalcomplexitiesencounteredwhen theclassical maximumlikelihoodestimation(MLE)methodisusedforthe
C.A. Hesse, D. B. Boyetey, E.D. Kpeglo, A.A. AshiagborEVALUATION OF PARAMETRIC ESTIMATION METHODRT&A,No1 (82) FOR THE GAMMA DISTRIBUTION...Volume 20, March2025
parameterestimation.It iscontestedthat MLEdoesnot provideastraightforwardsolutiontothelog-likelihoodequations sincetheinversecumulativedistributionfunction(CDF)ofthegamma distributioncannotbe expressed inaclosedform. Consequently, manystudieshave resorted tonumericalmethodstoevaluate gammadistributionparametersbyexploringacombinationofBayesianestimationand other simulation procedures.
Therefore, inthisstudy, we relied onone ofthe BayesianparametricsimulationmethodscalledMaximuma-Posteriori(MAP)toestimate the gammadistributionparameters. The MAPisconsidered the BayesiancounterparttoMLE. Itentailsidentifyingparametervaluesthatmaximizethe posterior distributionand actasestimatesfor the unknownparameters. Subsequently, we used the MarkovChainMonte Carlo(MCMC)simulationprocedure toobtainsamplesfromthe posteriordistributiontoobtainthegammaparameters.TheMCMCwas necessary becauseofthecomputationalintensitythe MAPrequirestoincorporate the prior information. Our objective wastoevaluatetheperformanceof theMaximumLikelihoodEstimation(MLE) andMaximuma-Posteriori(MAP)estimationmethodsunder the frameworkofMarkovChainMonteCarlo (MCMC)simulations.Thisevaluationwasdonebyfirstassumingacensoredlife-testingstrategythat terminatesonthe rth(where r < n)failure fromagivensample ofn electronicdevices. Second, weobtainedthejointdistributionfunctionofthefirst r-order statisticbyarrangingthe valuesofr in order ofmagnitude. Finally, we explored, throughthe MCMCframework, the parametricestimationofthegammadistribution usingtheMLEandtheMAP.
Thefindingofthestudysuggeststhat bothestimationmethodsyieldedtheexact parameterestimatesof10.4169 and0.0487,respectively,foralpha(.) andbeta (.) ofthegammadistribution.These estimateswere notsignificantlydifferentfromthe actualhypothesized value of . =10 and. =0.05. The seeminglynodifference betweenthetwoestimationmethodsconfirmstheir flexibilityinresolvingestimationcomplexitiesassociatedwiththegammadistributionfortime-to-failure variables. Further, we observed that irrespectiveofthepriordistributionusedfortheMAP_MCMCestimation, the resultingparametricestimatesofthe gammadistributionremainthe same (unchanged), confirmingthe assertionthatthe BayesianMCMCapproachisavaluable toolforinformative posterioranalysis.
Declarations
Funding
None
Conflict of interest
The authors declare that they have no competing interests.
Availability of data and materials
The dataset used in the analysis are contained in the study.
AcknowledgementsNotapplicable
References
[1]Dickson,D.C.(2012).Risk andrelativity:Estimating andunderstandinginsuranceclaimdistributionsusingstatistical models(2nded.).John Wiley & Sons.
[2]Veazie, P., Intrator, O., Kinosian, B., &Phibbs, C. S. (2023). Better performance for rightskeweddata usinganalternativegamma model.BMC Medical Research Methodology,23(1),298.
[3]
Dey, S., Elshahhat, A., &Nassar, M. (2023).Analysisofprogressive type-IIcensored gammadistribution.Computational Statistics,38(1),481-508.
C.
A. Hesse, D. B. Boyetey, E.D. Kpeglo, A.A. AshiagborEVALUATION OF PARAMETRIC ESTIMATION METHODRT&A,No1 (82) FOR THE GAMMA DISTRIBUTION...Volume 20, March2025
[4]Noguchi, K., &Robles, K.F. (2022). OnGeneratingDistributionswiththe MemorylessProperty.The American Statistician,76(3),280-285.
[5]Shore, H. (2023). Anovelapproachtomodelingsteady-state process-timewithsmoothtransitionskeweddatausing analternative gammamodel.BMC Medical Research Methodology,23(1), 298.
[6]Tao, Y. (2023). Memorylesspropertyofthe income distributionasanindicationfor testing
theequalityofopportunity:EvidencefromChina(1978-2015).https://www.preprints.org/manuscript/202303.0279/v1
[7]Elsayed, E. A. (2012). Reliability engineering (2nd ed.). John Wiley & Sons.
[8]Shipes, V. B., Meinzer, C., Wolf, B. J., Li, H., Carpenter, M. J., Kamel, H., &Martin, R. H.(2023).Designing aphase-III time-to-eventclinicaltrial usingamodifiedsamplesizeformulaandPoisson-Gammamodelforsubject accrualthat accountsforthelaginsiteinitiationusingthePERTdistribution.Statistics in Medicine,42(30),5694-5707.
[9]Meeker, W. Q., &Escobar, L. A. (1998). Statistical methods forreliability data.JohnWiley &Sons.
[10]Ghosh, I. &Hamedani, G. (2017). Gamma-kumaraswamy distribution in reliability analysis: properties and applications. Intechopen, http://dx.doi.org/10.5772/66821
[11]Junmei, Z., & Liqin, L. (2023).Estimatingparametersofthegamma distributioneasilyandefficiently.Communications in Statistics-Theory and Methods, 1-9.
[12]Ozsoy, V.S., Unsal, M.G. &Orkcu, H.H. (2020). Use ofthe heuristicoptimizationintheparameter estimationofgeneralized gammadistribution: ComparisonofGA, DE, PSOand SAmethods. Comput. Stat. 35,1895�1925.
[13]Hamada, M. S., Wilson,A. G., Reese, C. S., &Martz, H. F. (2008).Bayesian reliability. SpringerScience & Business Media.
[14]Rubinstein, R. Y., & Kroese, D. P. (2016). Simulation and the Monte Carlo method.John Wiley& Sons.
[15]Coles, S. G. &Tawn, J. A. (1996). ABayesianAnalysisofExtreme RainfallData. Appl Statistics, 45,463-478.
[16]Hussain, I., Haider, A.,Ullah, Z., Russo, M., Casolino, G. M., &Azeem, B. (2023). Compa�rativeanalysisofeightnumerical methodsusingWeibull distributiontoestimatewindpowerdensity forcoastalareas inPakistan.Energies,16(3),1515.
[17]Kohole, Y. W., Fohagui, F. C. V., Djiela, R. H. T., &Tchuen, G. (2023). Wind energypotentialassessmentforco-generationofelectricityand hydrogeninthe farNorthregionofCameroon.Energy Conversion and Management,279,116765.
[18]Hesse, C. A., Oduro,F.T, Ofosu, J. B.&Kpeglo, E. D. (2016). Assessingthe RiskofRoadTrafficFatalities Across Sub-PopulationsofaGiven GeographicalZone,UsingaModifiedSmeed�s Model.International Journal of Statistics and Probability, 5(6),121 �132.
[19]Benson, R. &Kellner, D. (2020).Monte Carlo simulation for reliability,in2020 AnnualRelia�
bilityandMaintainabilitySymposium(RAMS),2020,1�6. doi:10.1109/RAMS48030.2020.9153600.
[20]Naess, A., Leira, J. B.&Batsevych, O. (2009). Systemreliabilityanalysisbyenhanced MonteCarlosimulation,Structural Safety, 31(5),349�355.doi:10.1016/j.strusafe.2009.02.004.
[21]Fauzi, N. F. M., Roslan, N. N. R. &Ridzuan, M. I. M. (2023). Reliabilityperformance ofdistributionnetworkbyvariousprobabilitydistributionfunctions.IJECE, 13(2),2316-2325.
[22]Blanksby,P. E., & Lyons, B. (2023).Reliability and life testing (3rd ed.). Springer Nature.
[23]Mann, N. R., Schafer, R. E &Singpurwalla, N. D. (1974). Methods for Statistical Analysis of Reliability and Life Data.John Wiley&Sons,Inc.
[24]Ofosu,J.B.andHesse, C. A. (2010). Introduction to probability and probability distributions.EPP Books Services,Accra.
C.A. Hesse, D. B. Boyetey, E.D. Kpeglo, A.A. Ashiagbor
EVALUATION OF PARAMETRIC ESTIMATION METHOD RT&A,No1 (82)
FOR THE GAMMA DISTRIBUTION... Volume 20, March2025
[25]Dempster, A. P., Laird, N. M., &Rubin, D.B. (1977). Maximum likelihood from incomplete data via the EM algorithm.JournaloftheRoyalStatistical Society. Series B (Methodological), 39(1),1-38.
[26]Gilks, W. R., Richardson, S., &Spiegelhalter, D.J. (2004). Markov chain Monte Carlo in practice.ChapmanandHall/CRC.
[27]Nocedal, J., & Wright, S. J. (2006). Numerical optimization (2nded.).Springer.
[28]Press, W. H., Teukolsky, S. A., Vetterling, W. T., &Flannery, B. P. (2007). Numerical recipes in C++: The art of scientific computing (2nd ed.). Cambridge University Press.
[29]Burden, R. L., & Faires, J. D. (2000). Numerical Analysis. Brooks Cole.
[30]Kincaid, D., &Cheney, W. (2002). Numerical Analysis:MathematicsofScientificComputing.AmericanMathematicalSoc.
[31]Dahlquist, G., & Bjorck, A. (2008). Numerical Methods in Scientific Computing.SIAM.
[32]Eaton,J.W.(2002).GNU Octave and MATLAB in practice. Prentice Hall.
[33]Hanselman, D. C., & Littlefield, B. L. (2005). Mastering MATLAB 7. Prentice Hall.
[34]MathWorks.(2020).MATLAB documentation: Statistics Toolbox-Quantile.https://www.mathworks.com/help/matlab/ref/quantile.html
[35]Saulo, H., Vila, R., Borges, G. V., Bourguignon, M., Leiva, V., &Marchant, C. (2023). Mode�lingincome datavianewparametricquantile regressions: Formulation, computationalstatistics, andapplication.Mathematics,11(2),448.
[36]Steyvers,M.(2011).Computational statistics with MATLAB. UniversityofCalifornia, Irvine,psiexp.ss.uci.edu/research/teachingP205C/205C.pd.
Appendix
Listing A1: Likelihood function for . and .
1.functiony = MLE_gamma(alpha,beta,x,n,xr,r)
2.y=(gamma(alpha)-gammainc(alpha,beta*xr))^(n-r)*(1/gamma(alpha))^n*beta^(r*alpha)*prod(x.^(alpha-1))*exp(-beta*sum(x));
Listing A2: Implementation of Metropolis Hastings algorithm in MATLAB using the likelihood function
1.
% % Metropolis procedure to sample from the posteriordistribution
2.
%Component-wise updating. Use a normal proposaldistribution
3.
%%Set uptheImport Optionsandimport thedata
4.
opts = spreadsheetImportOptions("NumVariables", 1);
5.
% Specify sheet and range
6.
opts.Sheet ="Sheet1";
7.
opts.DataRange ="A1:A100";
8.
%Specifycolumn namesandtypes
9.
opts.VariableNames ="x";
10.
opts.VariableTypes ="double";
11.
%Import thedata
12.
x = readtable("C:\Users\USER\OneDrive\New Papers\Life Testing\Gamma.xlsx", opts,
"UseExcel", false);13.x=table2array(x);
14.r=length(x);15.xr=x(100);16.n=200;
17.
% % Initialize the Metropolis sampler
18.
T=5000;% Set the maximum number of iteration
19.
propsigma=[0.01,0.0001];%standarddeviationofproposal distribution
20.
parametermin=[9,0.04];%defineminimumforalphaandbeta
21.
parametermax=[11,0.06]; % define maximumforalphaandbeta
22.
seed=1; rand('state', seed ); randn('state',seed ); %#ok<RAND> % set the randomseed
23.
state=zeros(2,T); % storage space for the state of the sampler
24.
alpha=unifrnd(parametermin(1),parametermax(1)); %Startvalueforalpha
25.
beta=unifrnd(parametermin(2),parametermax(2)); %Start valueforbeta
26.
t=1; %initializeiterationat1
27.
state(1,t)=alpha; %savethecurrent state28.state(2,t)=beta;
C.A. Hesse, D. B. Boyetey, E.D. Kpeglo, A.A. Ashiagbor
EVALUATION OF PARAMETRIC ESTIMATION METHOD RT&A,No1 (82)
FOR THE GAMMA DISTRIBUTION... Volume 20, March2025
29.
%%Start sampling
30.
whilet<T% Iterate until we have T samples31.t=t+1;
32.% % Propose a new value for alpha33.new_alpha=normrnd(alpha,propsigma(1));
34.
pratio=MLE_gamma(new_alpha,beta,x,n,xr,r)/MLE_gamma(alpha,beta,x,n,xr,r);
35.
a=min([1pratio]);%Calculatetheacceptanceratio
36.
u=rand; % Draw a uniform deviate from [0 1]
37.
ifu<a %Doweacceptthisproposal?
38.
alpha=new_alpha;% proposal becomes newvalue for alpha
39.
end
40.
% % Propose a new value for beta41.new_beta=normrnd(beta,propsigma(2));
42.
pratio=MLE_gamma(alpha,new_beta,x,n,xr,r)/MLE_gamma(alpha,beta,x,n,xr,r);
43.
a=min([1pratio]);%Calculatetheacceptanceratio
44.
u=rand; % Draw a uniform deviate from [0 1]
45.
ifu<a %Doweacceptthisproposal?
46.
beta=new_beta; % proposal becomes newvalue for beta
47.
end
48.
%%Savestate
49.
state(1,t) = alpha;
50.
state(2,t) = beta;
51.
end52.Mean=mean(state,2)53.Mode=mode(state,2)
Listing A3: Likelihood function for . and .
1.functiony = MLE_gamma(alpha,beta,x,n,xr,r,a,b)
2.y=(gamma(alpha)-gammainc(alpha,beta*xr))^(n-r)*(1/gamma(alpha))^n*beta^(r*alpha)*prod(x.^(alpha-1))*exp(-beta*sum(x)-alpha/a-beta/b);
Listing A4: Implementation of Metropolis Hastings algorithm in MATLAB using the likelihood function
1.% % Metropolis procedure to sample from the posterior distribution
2.
%Component-wise updating. Use a normal proposaldistribution
3.
%%Set uptheImport Optionsandimport thedata
4.
opts = spreadsheetImportOptions("NumVariables", 1);
5.
% Specify sheet and range
6.
opts.Sheet ="Sheet1";
7.
opts.DataRange ="A1:A100";
8.
%Specifycolumn namesandtypes
9.
opts.VariableNames ="x";
10.
opts.VariableTypes ="double";
C.A. Hesse, D. B. Boyetey, E.D. Kpeglo, A.A. AshiagborEVALUATION OF PARAMETRIC ESTIMATION METHODRT&A,No1 (82) FOR THE GAMMA DISTRIBUTION...Volume 20, March2025
11.
%Import thedata
12.
x = readtable("C:\Users\USER\OneDrive\New Papers\Life Testing\Gamma.xlsx", opts,
"UseExcel", false);13.x=table2array(x);
14.r=length(x);15.xr=x(100);16.n=200;17.a=1518.b=0.08
19.
% % Initialize the Metropolis sampler
20.
T=5000;% Set the maximum number of iteration
21.
propsigma=[0.01,0.0001];%standarddeviationofproposaldistribution
22.
parametermin=[9,0.04];%defineminimumforalphaandbeta
23.
parametermax=[11,0.06]; % define maximum for alpha and beta
24.
seed=1; rand('state', seed ); randn('state',seed ); %#ok<RAND> % set the randomseed
25.
state=zeros(2,T); % storage space for the state of the sampler
26.
alpha=unifrnd(parametermin(1),parametermax(1)); %Startvalueforalpha
27.
beta=unifrnd(parametermin(2),parametermax(2)); %Start valueforbeta
28.
t=1; %initializeiterationat1
29.
state(1,t)=alpha; %savethecurrent state30.state(2,t)=beta;
31.
%%Start sampling
32.
whilet<T% Iterate until we have T samples33.t=t+1;
34.% % Propose a new value for alpha35.new_alpha=normrnd(alpha,propsigma(1));
36.
pratio=MLE_gamma(new_alpha,beta,x,n,xr,r,a,b)/MLE_gamma(alpha,beta,x,n,xr,r,a,b);
37.
a=min([1pratio]);%Calculatetheacceptanceratio
38.
u=rand; % Draw a uniform deviate from [0 1]
39.
ifu<a %Doweacceptthisproposal?
40.
alpha=new_alpha;% proposal becomes newvalue for alpha
41.
end
42.
% % Propose a new value for beta43.new_beta=normrnd(beta,propsigma(2));
44.
pratio=MLE_gamma(alpha,new_beta,x,n,xr,r,a,b)/MLE_gamma(alpha,beta,x,n,xr,r,a,b);
45.
a=min([1pratio]);%Calculatetheacceptanceratio
46.
u=rand; % Draw a uniform deviate from [0 1]
47.
ifu<a %Doweacceptthisproposal?
48.
beta=new_beta; % proposal becomes the newvalue for beta
49.
end
50.
%%Savestate
51.
state(1,t) = alpha;
52.
state(2,t) = beta;
53.
end54.Mean=mean(state,2)55.Mode=mode(state,2)
