Exploring the Potential Use of the Birnbaum-Saunders Distribution in Inventory Management
Choosing the suitable demand distribution during lead-time is an important issue in inventory models. Much research has explored the advantage of following a distributional assumption different from the normality. The Birnbaum-Saunders (BS) distribution is a probabilistic model that has its genesis in engineering but is also being widely applied to other fields including business, industry, and management. We conduct numeric experiments using the R statistical software to assess the adequacy of the BS distribution against the normal and gamma distributions in light of the traditional lot size-reorder point inventory model, known as (, ). The BS distribution is well-known to be robust to extreme values; indeed, results indicate that it is a more adequate assumption under higher values of the lead-time demand coefficient of variation, thus outperforming the gamma and the normal assumptions.
1. Introduction and Bibliographical Review
Inventory management permeates decision-making in countless firms. The topic has been extensively studied in academic and corporate spheres, for example, Braglia et al.  and Cai et al. . The key questions which the inventory management seeks to answer, usually influenced by a variety of circumstances, are as follows: when to order, determining an economic order quantity (EOQ) or lot size, and how much safety stock (SS) to keep, establishing a reorder point (ROP); see Namit and Chen  and Porras and Dekker .
According to Wanke , inventory management involves a set of decisions whose objective is to match existing demand with the supply of products and materials over space and time. This objective allows us to achieve specified costs and service levels, considering product, operation, and demand characteristics. It is known that the inventory total cost (TC) is a function of ordering, holding, and shortage costs; see Hillier and Lieberman .
The importance attached by firms to inventory management can be attributed to the following: first and foremost, it is the need to ensure that products, given the competitive pressure exercised by markets, are always supplied to customers at the least possible cost; see Eaves . Second, some other factors contribute to a high concern with inventory management, such as product diversity or behavior; see Huiskonen . High opportunity costs also contribute to this concern, thus affecting the financial indicators on which assessments of firm performance are based; see Wanke .
The inventory management models are frequently classified in two types: pull and push. On the one hand, according to Ballou and Burnetas , pull-type planning models range from those that set inventory levels based on the EOQ to those fixed in proportion to forecasted demand. The EOQ model is the simplest and most fundamental of all inventory models because it describes important trade-offs between fixed ordering and holding costs; see Nahmias . Despite its shortcomings, the basic EOQ model is the cornerstone of several software packages for inventory control; see Lee and Nahmias . Interested readers can refer to Yan and Wang  and Min et al.  for more details about the EOQ model. Today EOQ is used in conjunction with ROP in inventory control models to determine cycle and SS under demand per unit of time (DPUT) and lead-time (LT) uncertainty. These models are well described in most logistics and operations textbooks, as are their underlying assumptions; see Nahmias . Models that set inventory levels in proportion to forecasted demand constitute a particular form of the periodic review model, except that replenishment quantities are not based on the EOQ model; see Ballou . On the other hand, push-type planning models take place when inventory decisions are based on the demand or its forecast at multiple downstream stocking locations, similar to the resource planning system of logistics.
According to Silver et al. , Eaves and Kingsman , Syntetos et al. , and Boylan et al. , demand is one of the main factors in inventory management models. In general, demand may be classified from two settings. First, it can be deterministic or random; if random, the demand follows a statistical distribution (also known as probabilistic model); otherwise, the demand is constant, which implies a degenerate statistical distribution; that is, its variance is equal to zero. Second, demand might be independent or dependent. For more details, see Disney et al. , Porras and Dekker , Wanke , and Rojas et al. .
Demand uncertainties directly affect the operation of the physical system of logistics. Moreover, to be closer to reality, single or multiple period inventory models must take into account that demand is occurring in a random fashion, which is explained by several factors. Thus, DPUT is taken to be a random variable (RV). Furthermore, during LT, due to the mentioned randomness, the corresponding demand (LTD) is also a RV; therefore, the behavior of DPUT and LTD must be described by statistical distributions; see Johnson et al. [21, 22]. The Gaussian (or normal) distribution is often used for describing the data of these two RVs (DPUT, LTD) involved in inventory models. However, it is well-known that the normal distribution is validly used for RVs that take negative and positive values with a symmetrical behavior. Hence, first, quantities less than zero could be admitted when the modeling is carried out under the normal distribution, which is not possible in real-world situations for DPUT and LTD, because they only admit values greater than zero; see Nahmias . Second, another drawback using the normal model is that DPUT and LTD data often follow asymmetric distributions; see Moors and Strijbosch . Mentzer and Krishnan  studied the nonnormality effect on inventory models and found that the normal distribution is appropriate in few practical cases; see also Eppen and Martin . A recent case study with DPUT data of 89 food products supports such nonnormality; see Leiva et al.  and Rojas et al. . In any case, the normality assumption must be checked by goodness-of-fit methods; see Barros et al. . Thus, the use of the normal distribution to model DPUT and LTD and then to determine the ROP and SS can lead to wrong results, resulting in shortages or excess inventories. Nonnormal distributions with positive support that have been used for describing DPUT in inventory management include models such as gamma or Erlang, inverse Gaussian, log-normal, Pearson, Poisson, uniform, and Weibull; see Burgin , Tadikamalla , Lau , Wanke , Cobb et al. , and Pan et al. .
A unimodal, two-parameter probability model with positive support and asymmetry to the right that is receiving considerable attention is the Birnbaum-Saunders (BS) distribution; see Birnbaum and Saunders  and Johnson et al. [22, pages 651–663]. The BS distribution has good properties and is related to the normal distribution and implemented in the R statistical software (http://www.r-project.org) via a package called gbs; see R Team . Although the BS distribution has its genesis from engineering, its applications range across diverse fields as business, industry, and management, which have been conducted by an international, transdisciplinary group of researchers; see, for example, Jin and Kawczak , Podlaski , Bhatti , Lio et al. , Paula et al. , Marchant et al. , and Leiva et al. [42, 43]. In addition, although originally conceived as a count model, the BS distribution includes the duration of the counting period (daily or weekly), which obviates having to collect additional data, among other properties; see Fox et al. . In sum, the BS distribution is a good candidate for describing demand data in inventory models; see Leiva et al.  and Rojas et al. .
Our main objective is to explore the use of the BS distribution in inventory management. Differently from previous studies that exclusively considered the effects of one given distribution on inventory decision-making, we also analyze its adequacy in light of different operating characteristics and costs. Specifically, we assess how the BS, gamma, and normal LTD distributions interact with relevant product characteristics and affect the optimal EOQ and SS inventory indicators in terms of the optimization of the TC function. We minimize this function using stochastic programming, a technique where constraints and/or objective function of the problem to be optimized contain RVs that can follow any distribution; see Shapiro et al.  and Thangaraj et al. . We solve the problem of stochastic programming with a search heuristic called differential evolution (DE), which is a global numerical optimization approach based on genetic algorithm concepts; see Storn and Price  and Price et al. . We implement our results in R code, which is available upon request from the authors.
Section 2 reviews general aspects of inventory management models and the statistical distributions used in this study. Section 3 explores the effect of different LTD distributions in inventory management, introducing the simulation scenario, formulating the stochastic programming model, discussing the DE algorithm, and providing a numerical study. Section 4 concludes the study making some considerations on the management of our findings and on future research.
In this section, we discuss general aspects of inventory management models and demand statistical distributions used for the methodology presented in Section 3.
2.1. Inventory Management Models
The model is based on the ROP () and the EOQ model given bywhere is the DPUT rate in units of the product and , are the ordering and holding costs, respectively; see Yan and Wang  and Min et al. . However, as mentioned, DPUT is a RV. Then, given in (1) must be calculated as the mean (expected value) of the DPUT distribution that adequately fits the data. Specifically, let be a RV corresponding to the DPUT at time , forming a sequence of independent and identically distributed RVs with mean and variance . In addition, let be a RV corresponding to the LT between the ordering of a product and its delivery (expressed in time units) with mean and variance . Then, the LTD is given bywith probability density function (PDF) and whose expectation and variance are, respectively, defined asThe ROP can be computed from expressed in (3). However, to be protected from randomness of the LTD, it is necessary to include a SS, which allows the ROP to becomewhere , are defined in (3) and is the safety factor (SF) or number of standard deviations (SDs) of the LTD. Note that although the model is based on (1) and (4), it is possible to see that is obtained by . Thus, we refer to the model as thereafter.
The expected TC of the inventory is given bywhere (in units of the product) is given in (1) and , , and are given in (4); is the holding cost (in $ per $ per unit of time); is the ordering cost (in $ per each replenishment order placed); is the shortage cost (in $ incurred whenever a stock-out occurs); and is the PDF of the LTD given in (2). In order to minimize the expected TC defined in (5), we optimize the indicator given in (1) altogether with given in (4).
2.2. Demand Distributions
Notice that it is necessary to specify the LTD distribution to determine the SS given in (4), which allows the SF to be established; see Porras and Dekker . In order to facilitate the calculation of the SF, the LTD has been traditionally modeled with the normal distribution; see Silver and Peterson . Thus, the SF for a specific service level can be obtained from a percentile of the standard normal distribution, denoted by . However, as mentioned, various studies criticize the normality assumption. Therefore, the use of the normal distribution to determine ROP and SS given in (4) is questionable, leading to possible stock shortage or excess.
Silver  pointed out that in most models leading to inventory management decisions some assumptions are made, often in an implicit way. The effects of these assumptions on costs and service levels should be taken into account. The most common ones are (i) to assume a demand distribution (e.g., normal) and (ii) to suppose that the distribution parameters are known (e.g., the mean and SD) or estimated from the demand data. Lau  presented a model for computing EOQs and SSs given in (1) and (4), respectively, using the first four moments, that is, mean, variance, third moment reflecting skewness, and fourth moment reflecting kurtosis of any given LTD distribution. Lau  also pointed out the risk of misleading decisions regarding ROP and customer service level when one considers a normally distributed LTD. The 95th percentile of the distribution is often used to set service levels. Next, we present some mathematical features for the three LDT distributions to be considered in this study, that is, the BS, gamma, and normal models.
The Normal Distribution. A RV following a normal distribution with mean and variance is denoted by , where “” means “distributed as”. In this case, PDF, cumulative distribution function (CDF), and (QF) quantile function of are, respectively, where and , for , with and being the inverse CDF or QF. In addition, the coefficients of variation (CV), skewness or asymmetry (CS), and kurtosis (CK) of are, respectively,
The BS Distribution. A RV following a BS distribution with shape and scale parameters is denoted by . In this case, the PDF, CDF, and QF of are, respectively,where , for , is the CDF, is the QF, and is the inverse CDF of . Note that ; that is, is also the median or 50th percentile of the distribution. The mean, variance, CV, CS, and CK of are, respectively, In addition, the RVs and are related by Also, note that follows a chi-squared distribution with one degree of freedom. The BS distribution holds the scale and reciprocation properties; that is, (i) , with , and (ii) , respectively.
The Gamma Distribution. A RV following a gamma distribution with shape and scale parameters is denoted by . In this case, the PDF and CDF of are, respectively, where and stand for the usual and incomplete gamma functions, respectively. The corresponding QF given by , for , must be obtained by solving this equation with an iterative numerical method. The mean, variance, CV, CS, and CK of are, respectively, The gamma distribution also shares the scale property; that is, , with .
3. Assessing the Impact of Different Distributions
In this section, we introduce our simulation scenario and formulate the stochastic programming used to optimize the expected TC associated with the model. Then, we discuss the DE algorithm, which allows us to solve the problem of stochastic programming, and provide a numerical study performed with the R software. We evaluate how different LTD distributions (BS, gamma, and normal) interact with different inventory indicators (demand, cost, and LT) and their underlying EOQs and SSs. We assess under what circumstances a distributional assumption is preferable to the other one in terms of the expected TC.
3.1. Scenario of the Simulation Study
Assume BS, gamma, and normal distributions for the LTD. Then, fix values for the parameters of these distributions by considering values for means and SDs of DPUT and LT generated from uniform distributions. Now, generate holding, ordering, and shortage costs also from uniform distributions. This allows us to establish the expected TC to be minimized. Ten thousand (10000) different simulated scenarios of means and SDs for DPUT and LT as well as holding, ordering, and shortage costs are generated using an R package called stats. The values of these uniformly distributed inventory indicators used to build the scenarios are presented in Table 1. They were chosen based on values proposed in selected papers focused on managerial and industrial applications compiled by Wanke ; see Table 2 in this reference for more details about values that are frequently used to generate simulation scenarios in inventory management problems. A discussion about this can be found in Wanke .
3.2. Stochastic Programming
Once the values for the inventory indicators are defined and the distributional assumptions (BS, gamma, and normal) for the LTD established, stochastic programming is performed on the expected TC function given in (5); see Namit and Chen . According to the values provided in Table 1, we assume values for means and SDs of the DPUT and LT, as well as for the holding, ordering, and shortage costs, to be uniformly distributed in the objective function corresponding to the expected TC given in (5). The decision variables of the programming are and given in (1) and (4), respectively, whereas , , , , , , and are given from (5). Therefore, our optimization problem of the expected TC can be visualized like a model of stochastic programming formulated as
The problem of stochastic programming formulated in (14) is aimed (i) at identifying the most adequate inventory policy for each of the distributional assumptions and (ii) at minimizing the expected TC. The optimized problem can provide useful information for academics and practitioners on how these assumptions interact with product, operation, and demand characteristics. To solve the problem of stochastic programming formulated in (14), we use the DE algorithm detailed in the next section.
3.3. Differential Evolution
DE is a member of the family of genetic algorithms, which mimic the process of natural selection in an evolutionary manner; see Holland . A genetic algorithm solves optimization problems with biology-inspired operators of crossover, mutation, and selection, generating successive populations of individuals (solutions or generations). Then, the DE algorithm optimizes problems by evolving a population of candidate solutions employing the mentioned operators. The DE algorithm uses floating-point techniques for obtaining the solutions and arithmetic operations in their mutation, in contrast to classic genetic algorithms. In addition, the DE algorithm finds the global optimum of the objective function, which is not required to be either continuous or differentiable; see Thangaraj et al.  and Mullen et al. .
The DE algorithm has also been used to optimize problems that arise in inventory management, such as joint replenishment, replenishment coordination, and inventory location-allocation; see Qu et al. [54, 55] and Wang et al. . In the present study, the problem of stochastic programming formulated in (14) is solved with the DE algorithm, which allows us to find the optimal values of and that minimize the expected TC function given in (5), for 10000 simulated scenarios.
In what follows, we first discuss the DE algorithm used in our research work; see Storn and Price  and Price et al.  for more details about the algorithm. Then, we sketch the content of an R package called DEoptim, which implements the DE algorithm and was first published on CRAN in 2005. Since becoming publicly available, it has been used by several authors to solve optimization problems arising in diverse domains. We refer interested readers to Ardia et al.  and Mullen et al.  for a detailed description of the package.
Let be the number of members (tuning parameter vector) in the population, where denotes the dimension of the vector , in our case given by . The DE algorithm needs a starting population, which is obtained by sampling the objective function at multiple randomly chosen initial points (generation or solution zero (0)) for and . Before the population is initialized, both lower and upper bounds for each tuning parameter must be specified. Parameter bounds establish the domain from which the vectors at the generation 0 are chosen. To establish the generation 0, guesses for the optimal value of must be provided, using either random values within a range defined by the practitioner or values fixed by him/her. Each generation creates a new population from the current population members , where indexes the tuning parameter vector that make up the population and indexes the generation. The new generation is obtained using differential mutation of the population members. In our research, denotes the number of simulated scenarios and denotes the tuning parameters for the lot size and the number of demand SDs. An initial mutant parameter vector is created by choosing three members of the population , , and , at random. Then, the elements of the initial mutant parameter vector are generated bywhere in our case and is a positive scale factor whose effective value is usually less than one (usual default: ). After the first mutation operation, it continues until mutations have been made or until certain crossover probability (CP) in is less than , where is a random number from the uniform distribution in . CP controls the fraction of the tuning parameter values that are copied from the mutant and approximates the probability that a parameter value is inherited from the mutant, since at least one mutation always occurs. Mutation is applied in this way to each member of the population.
Calculations in our simulation study were performed with the aid of the DEoptim package, which consists of the core function DEoptim() whose arguments are as follows:(i)fn is the function to be minimized, which must have as its first argument the vector of real-valued parameters to optimize and return a scalar real result.(ii)lower, upper correspond to two vectors establishing scalar real lower and upper bounds on each tuning parameter to be optimized; the th element of the lower and upper vectors corresponds to the th parameter; the implementation searches the global optimum of fn between lower and upper.(iii)… allows the user to pass additional arguments to the function fn.(iv)control is a list whose default value is the return value of DEoptim.control() but whose main elements are interpreted as follows:(a)VTR specifies the global minimum of fn if it is known or if you wish to cease optimization after having reached a certain value (default = -Inf).(b)strategy defines the differential evolution strategy used in the optimization procedure, described in detail by Mullen et al. .(c)NP is the number of population members ( or 50), in our case denoted by .(d)bs: if bs is FALSE, then every mutant is tested against a member in the previous generation, and the best value survives into the next generation; if bs is TRUE, then the old generation and NP mutants are sorted by their associated objective function values, and the best NP vector proceeds into the next generation (default = FALSE).(e)itermax is the maximum number of iterations (i.e., population generations) to be allowed (default = 200).(f)CR is the CP from the interval (default = 0.9).(g)F is the stepsize from interval , in our case denoted by (default = 0.8).(h)trace, initialpop, storepopfrom, storepopfreq, checkWinner, and avWinner are other elements of the list control, which are described in detail by Mullen et al. ; see also Price et al. .
The return value of the function DEoptim() is a member of the S3 class DEoptim. Members of this class have a plot and a summary that allow us to analyze the optimizer’s output. In our application of the DEoptim package, decision variables and were lower-bounded at zero. As regards the list of tuning parameters used in our research, the default values for the DEoptim package were considered.
3.4. Numerical Results and Discussion
The underlying idea of performing a sensitivity analysis on the testing variables related to product, demand, and operational characteristics is to discriminate between groups where the three distributional assumptions led to minimal TCs. Table 2 summarizes the numerical experiments conducted with the R software for 10000 simulated scenarios using the DE algorithm. Note that the BS distribution yields, on average, smaller TCs, EOQs, SSs, and, therefore, smaller inventory levels, in comparison to the normal and gamma assumptions.
Table 3 depicts the adequacy of the distributional assumption built upon the empirical evidence presented in Silver et al. . We conclude that the BS assumption for the LTD was more adequate in 81.60% of cases in terms of TCs and inventory levels; moreover, it prevails in circumstances of a higher LTD CV in comparison to gamma and normal assumptions, that is, a CV of 0.56 for the BS distribution against CVs of 0.44 and 0.34 for the gamma and normal distributions, respectively. It is worth mentioning that, due to the numeric integral for computing the BS CDF, which diverged in less than 1% of cases, the gamma distribution was a preferable assumption in terms of TCs and inventory levels in such cases.
This paper has assessed how different demand distributions during lead-time interact with relevant characteristics of the product. We have considered the choice of the optimal inventory policy in terms of total costs of the inventory management. Differently from most previous studies on the topic that exclusively explored the effects of one single distributional assumption, this paper also analyzed its adequacy in light of some inventory management key elements for decision-making, such as cost, demand, and lead-time. It was shown that the Birnbaum-Saunders distribution outperformed the normal and gamma assumptions with respect to demand uncertainty during the lead-time. The contributions of this paper are on both the academic and the practical sides. Departing from what was found in previous studies, the obtained results provide a guidance on the selection of the most appropriate distributional assumption for the demand during the lead-time. Specifically, while the Birnbaum-Saunders distribution is more adequate to handle demand during lead-time with high coefficients of variation, the gamma and the normal assumptions should be restricted to well-behaved patterns. This paper also can present a practical contribution by means of numerical analyses conducted with the aid of a computational code developed for such purposes in the R statistical software.
With respect to future research endeavors, one possible extension of this work could be to analyze the results presented here in terms of inventory management in decentralized systems. The basic idea would be to assess, via numerical studies, the extent to which the decision to pool inventories is affected by different distributional assumptions and by their levels of skewness and kurtosis. In addition, in practical case studies, the usage of diagnostic statistical methods to detect influential data of demand or LT can be studied for this type of inventory models. In that case, demand distributions that provide parameter estimators to be robust to atypical data of demand must be considered; see, for example, Paula et al. .
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
The authors thank the editors and referees for their constructive comments on an earlier version of this paper which resulted in this improved version. This research was supported by COPPEAD, Capes, and CNPq from Brazil and FONDECYT 1120879 from Chile.
K. Namit and J. Chen, “Solutions to the inventory model for gamma response time demand,” International Journal of Physical Distribution & Logistics Management, vol. 29, pp. 138–154, 1999.View at: Google Scholar
F. Hillier and G. J. Lieberman, Introduction to Operational Research, McGraw Hill, New York, NY, USA, 2005.
A. Eaves, Forecasting for the ordering and stock-holding of consumable spare parts [Ph.D. dissertation], Lancaster University, Lancaster, UK, 2002.
S. Nahmias, Production and Operations Analysis, McGraw-Hill, New York, NY, USA, 2001.
H. L. Lee and S. Nahmias, “Single-product, single-location models,” in Logistics of Production and Inventory, S. C. Graves, A. H. G. Kahn, and P. H. Zipkin, Eds., pp. 3–55, North-Holland, Amsterdam, The Netherlands, 1993.View at: Google Scholar
E. A. Silver, R. Peterson, and D. F. Pyke, Inventory Management and Production Planning and Scheduling, Wiley, New York, NY, USA, 1998.
N. L. Johnson, S. Kotz, and N. Balakrishnan, Continuous Univariate Distributions, vol. 1, Wiley, New York, NY, USA, 1994.View at: MathSciNet
N. L. Johnson, S. Kotz, and N. Balakrishnan, Continuous Univariate Distributions, vol. 2, Wiley, New York, NY, USA, 1995.
J. T. Mentzer and R. Krishnan, “The effect of the assumption of normality on inventory control/customer service,” Journal of Business Logistics, vol. 6, no. 1, pp. 101–120, 1985.View at: Google Scholar
P. R. Tadikamalla, “The inverse Gaussian approximation to the lead-time demand in inventory control,” International Journal of Production Research, vol. 19, pp. 213–219, 1981.View at: Google Scholar
H. S. Lau, “Toward an inventory control system under non-normal demand and lead-time uncertainty,” Journal of Business Logistics, vol. 10, pp. 88–103, 1989.View at: Google Scholar
R Team, R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria, 2014, https://www.r-project.org.
X. Jin and J. Kawczak, “Birnbaum-Saunders and lognormal kernel estimators for modelling durations in high frequency financial data,” Annals of Economics and Finance, vol. 4, no. 1, Article ID 103124, pp. 103–124, 2003.View at: Google Scholar
A. Shapiro, D. Dentcheva, and A. Ruszczynski, Lectures on Stochastic Programming: Modeling and Theory, SIAM, Philadelphia, Pa, USA, 2009.
R. Thangaraj, M. Pant, P. Bouvry, and A. Abraham, “Solving multi objective stochastic programming problems using differential evolution,” in Swarm, Evolutionary, and Memetic Computing, B. K. Panigrahi, S. Das, P. N. Suganthan, and S. S. Dash, Eds., vol. 6466 of Lecture Notes in Computer Science, pp. 54–61, Springer, New York, NY, USA, 2010.View at: Publisher Site | Google Scholar
K. V. Price, R. M. Storn, and J. A. Lampinen, Differential Evolution: A Practical Approach to Global Optimization, Springer, Berlin, Germany, 2006.
E. A. Silver and R. Peterson, Decision Systems for Inventory Management and Production Planning, Wiley, New York, NY, USA, 1985.
J. H. Holland, Adaptation in Natural Artificial Systems, University of Michigan Press, Ann Arbor, Mich, USA, 1976.
K. M. Mullen, D. Ardia, D. L. Gil, D. Windover, and J. Cline, “DEoptim: an R package for global optimization by differential evolution,” Journal of Statistical Software, vol. 40, no. 6, pp. 1–26, 2011.View at: Google Scholar
L. Wang, H. Qu, S. Liu, and C. Chen, “Optimizing the joint replenishment and channel coordination problem under supply chain environment using a simple and effective differential evolution algorithm,” Discrete Dynamics in Nature and Society, vol. 2014, Article ID 709856, 12 pages, 2014.View at: Publisher Site | Google Scholar
D. Ardia, K. Boudt, P. Carl, K. Mullen, and B. Peterson, “Differential evolution with DEoptim: an application to non-convex portfolio optimization,” R Journal, vol. 3, pp. 27–34, 2011.View at: Google Scholar