The MLP Distribution: A Modified Lognormal Power-Law Model for the Stellar Initial Mass Function
This work explores the mathematical properties of a distribution introduced by Basu & Jones (2004), and applies it to model the stellar initial mass function (IMF). The distribution arises simply from an initial lognormal distribution, requiring that each object in it subsequently undergoes exponential growth but with an exponential distribution of growth lifetimes. This leads to a modified lognormal with a power-law tail (MLP) distribution, which can in fact be applied to a wide range of fields where distributions are observed to have a lognormal-like body and a power-law tail. We derive important properties of the MLP distribution, like the cumulative distribution, the mean, variance, arbitrary raw moments, and a random number generator. These analytic properties of the distribution can be used to facilitate application to modeling the IMF. We demonstrate how the MLP function provides an excellent fit to the IMF compiled by Chabrier (2005) and how this fit can be used to quickly identify quantities like the mean, median, and mode, as well as number and mass fractions in different mass intervals.
keywords:accretion—stars: clusters—stars: formation—stars: mass function
The distribution of stellar and substellar masses at birth, the Initial Mass Function (IMF), is a key feature of star formation. It has been studied intensely since first estimated by Salpeter (1955), who measured a power-law tail for high masses of the approximate form . Subsequent work has established a shallower slope at masses less than and a turnover at approximately when the masses are placed in logarithmically spaced bins. These compilations of the IMF have tended to identify a power-law profile at high masses (e.g., Scalo, 1998; Kroupa, 2001, 2002), although earlier work (Miller & Scalo, 1979) did fit the IMF with a lognormal distribution. Chabrier (2003) (see also Chabrier (2005)) has compiled an IMF in the substellar and low mass stellar regime and advocates a lognormal fit for masses up to and a power-law fit for . By appealing to the Central Limit Theorem (CLT), Chabrier (2003) claims a better rationale for the lognormal fit at low masses than the approach of using broken three-component power-laws (Scalo, 1998; Kroupa, 2001, 2002). However, Chabrier’s overall fit also requires joining a lognormal with a power law at high masses, so has one joining condition instead of two. An ideal next step is to find a single function with no joining conditions that has a rationale that is at least on par with appeals to the CLT.
There are many disciplines in which a desired distribution is one that is like a lognormal at low and intermediate values, with a characteristic peak and turnover, but transitions to a power law distribution at high values. Besides astronomy, this need has arisen in fields as diverse as biology (Limpert et al., 2001), computer science (Mitzenmacher, 2006), ecology (Allen et al., 2001), and finance (Mandelbrot, 1997), to name several. A review of common statistical resources reveals that very few analytic functions of this type exist, hence the attempts to fit empirical distributions by patching together different functions over different domains. Power-law distributions were first introduced by Pareto (1896) to explain the distribution of incomes seen in data from many different countries. The distribution of city sizes also shows a power-law character and regularity across many countries. It is referred to as Zipf’s Law (Zipf, 1949). Henceforth, we refer to pure power-law distributions as Pareto distributions, in conformity with much of the statistical literature.
In this paper, we analyze and characterize the properties of a hybrid three-parameter probability density function (pdf) introduced by Basu & Jones (2004). We feel that it can be used to fruitfully model data sets that exhibit both lognormal-like and power-law behavior. Indeed, the MLP distribution illustrates the fact that many generative processes that lead to a lognormal distribution, can with some modification, yield a power-law tail instead. As a result of its origin as a modified lognormal, two parameters of the MLP distribution are identified as and , preserving the notation of the lognormal distribution, while the third parameter is , the power-law index that also characterizes the Pareto distribution. However, the parameters and no longer represent the mean and variance, respectively, of the logarithm of the random variable, as they do for the lognormal distribution. This three-parameter function is a simpler and more readily usable form of a more general four-parameter pdf derived by Reed in several papers (Reed, 2002, 2003). The latter function arises from a stochastic growth law, Geometric Brownian Motion, rather than the pure exponential growth used by Basu & Jones (2004). An advantage of the pdf we use here is that it introduces only one additional parameter beyond that in the lognormal. As a result, it is a natural first step when fitting data that may look like a modified lognormal, and its relatively compact analytic closed-form expression makes it easy to use with common fitting techniques.
The model of Basu & Jones (2004) falls into one of two major categories of IMF models in modern day star formation theory, that of an accretion based scenario in which temporal effects are important. The basic idea is that star formation is a killed (or stopped) process, with mass accretion terminated by events such as stellar outflows, ejection from the mass reservoir, or any other reason for emptying the mass reservoir. The distribution of accretion lifetimes then plays a key role in setting the shape of the IMF, while the Jeans mass does not. Other models in this category include those of Adams & Fatuzzo (1996), Bate & Bonnell (2005), and Myers (2009, 2014), with the latter developing a model that can fit the observed IMF by tuning several parameters. The alternate scenario is one in which the mass distribution is determined by the spatial properties of the gas distribution and by gravitational instability, so that some combination of the turbulent spectrum and the Jeans mass set the IMF, e.g., Padoan & Nordlund (2002); Hennebelle & Chabrier (2009), with the latter developing an analytic approach. Although the fitting of analytic functions to the IMF cannot alone settle which models may be most suitable, they can however greatly facilitate analysis in fields like galaxy studies, where the conclusions depend strongly on an adopted IMF model. Analytic functions also allow a simpler analysis of the effect of varying IMF parameters.
The paper is organized as follows. In Section 2 and Section 3 we introduce the lognormal and Pareto distributions, respectively, and present some of their relevant properties. This is done for completeness of the presentation and to add context when reading the new results about the MLP distribution. In Section 4 we discuss the formulation of Basu & Jones (2004) that leads to the MLP distribution. We then examine some relevant mathematical properties of this distribution in Section 5, which includes expressions for its cumulative distribution function, mean and variance, arbitrary raw moments, and an approximation to its mode. These expressions, excluding the approximate mode, are shown to reduce to the corresponding lognormal expressions in the appropriate limit. In Section 6, we fit the MLP distribution to the IMF of Chabrier (2005) and then use the analytic function to quickly estimate some of the above described IMF properties, as well as a cumulative mass fraction. Some closing remarks are given in Section 7. The derivations of the expressions given in Section 5 are included in Appendix B, which also contains relevant properties of the error and complementary error functions and related integrals used in this paper.
2 The Lognormal Distribution
According to the CLT of probability theory (Gut, 2005), if are identically distributed, independent random variables with mean and standard deviation , then converges in distribution to the standard normal variable with zero mean and variance of unity. In addition, the identical distribution assumption for may be dropped, and the result will also follow provided certain conditions (Lindeberg’s) are satisfied (Gut, 2005). As pointed out in Golberg (1984) the CLT is used to partially justify why so many observable phenomena appear to be normally distributed: if the variable of interest is thought to be influenced by the sum of a large number of independent factors, then the CLT can be invoked to explain the apparent normality. However, it is frequently the case that the variables of interest are non-negative, whereas normal variables are not. One way to circumvent this complication is to have a distribution with similar properties but which is always positive. Suppose that above was instead assumed to be of the form . Then , where . Note that satisfy the conditions of the CLT, so that converges to a normal random variable, . Moreover, if is normal, then has a lognormal density (Golberg, 1984; Aitchison & Brown, 1957; Crow & Shimizu, 1988),
Thus, under the assumptions above, converges to a lognormal. We list some of the relevant properties of the lognormal distribution below:
Cumulative Distribution Function (cdf)
where is the error function (see Equation 29 in the appendix).
3 The Pareto Distribution
Pareto distributions have been used extensively to model a wide variety of phenomena in the sciences and social sciences, such as the size of forest fires, the intensity of earthquakes, the citations to papers, and the population of cities. For a recent review, see (Clauset et al., 2009).
A random variable has a Pareto distribution if its probability density function is
Cumulative Distribution Function (cdf)
4 The MLP distribution
The key element in the derivation of the MLP distribution is that even though the initial values of a random variable may have a lognormal density, later time evolution could in fact skew their distribution. In other words, the subsequent time of growth of the random variable (representing a physical quantity) is itself another random variable. Our pdf can be derived analytically on this basis using a few simplifying assumptions.
Imagine that initial values of a quantity are drawn randomly from a lognormal distribution with parameters and . If the subsequent growth of an object with is characterized by exponential growth
with fixed growth rate , then the multiplicative relation between values of individual objects is preserved. A lognormal distribution is then maintained with the same but the mean of the logarithmic values shifts to after a fixed time . However, we can treat the time as a random variable and draw it from an exponential pdf
where is a stopping rate. In this case, we can derive the probability density function of final masses as
Using the identity in subsection A.2 we find the closed form
Figure 1 shows the MLP density (Equation 14) and a lognormal distribution for similar values of parameters. The specific values of the parameters used here are somewhat arbitrary, but correspond approximately to a best fit lognormal for the low mass end of a stellar mass distribution, as modeled in Basu & Jones (2004). While is largely a scale-dependent parameter, the other two parameters are expected to fall in the approximate range and based on fits of lognormal and power-law distributions to a wide range of phenomena in the sciences and social sciences (Limpert et al., 2001; Clauset et al., 2009). The function derived here is related to a four-parameter distribution (Reed, 2002, 2003) that can be derived under the assumption of geometric Brownian motion, , which is a stochastic growth law where represents , the uniform random variate in the interval , and is an amplitude of the fluctuations. The resulting pdf is a double-tailed Pareto distribution, with coefficients that are roots of a quadratic equation. The pdf we derive here results from a simpler model and has the advantage of a single-expression closed form. It may also more easily correspond to many data sets that seemingly warrant only one additional parameter (beyond the original two of the lognormal distribution) in order to get a reasonable fit and quantify the data.
5 Relevant Mathematical Properties of the MLP Distribution
A few relevant mathematical properties of the MLP distribution are examined below, leaving most of the necessary calculations in Appendix B. Let denote a random variable with an MLP distribution of density .
5.1 Cumulative Distribution
5.2 Mean, Variance, and Raw Moments
The th raw moment of , defined as the expectation value of ,
exists if and only if , in which case it is given by
Note that this expression is exactly the formula for the raw moments of a lognormal distribution (Equation 2) with parameters and , scaled by the factor , and in the limit as , the expressions are identical. This is consistent with the derivation of the MLP distribution in Section 4. The limit corresponds to for finite death rate , so that the drift term vanishes, and the distribution remains a lognormal with mean . Assuming , we can obtain the following expressions, including the mean and variance of the distribution:
Higher moments around the mean can be computed using Equation 18 with the identity
To find the mode, that is, the value that maximizes the MLP pdf in Equation 14, we must solve the transcendental equation
Although the solution to Equation 23 will generally require the implementation of numerical methods, we note that if then provides an approximate solution to Equation 23, which in terms of the original parameters results in
The approximation in Equation 25 is useful only when the assumption is closely met, and behaves poorly otherwise. However, when a precise numerical solution is required, one can use this approximation as a starting point in the iteration procedure being implemented. It is worth noting that even if one tried to find the peak in the space vs , the resulting equation to be solved is in fact the same, owing to the fact that is just with a different value of , and .
5.4 Random Number Generation
For practical purposes of comparing data sets with a model distribution, it is valuable to be able to draw a random sample from the model. Using the definition of the lognormal random variable (see Section 2), a random number drawn from a lognormal pdf (Equation 1) is an element of the lognormal random variate
where is the normal random variate with zero mean and variance of unity. For drawing from an exponential pdf (Equation 12), the exponential random variate is
where is the uniform random variate in the interval . The above formula can be obtained from Equation 12 through the general method of calculating the cdf, inverting it, and then drawing the argument as an element of the uniform random variate . We can use Equation 26 and Equation 27 to derive a random variate for the MLP distribution. We note that the MLP distribution (Equation 14) is formally obtained from an initial lognormal pdf with mean under the transformation , in which is chosen randomly from an exponential distribution. Therefore, we can use Equation 26 and Equation 27 to write that the MLP random variate is
6 Fitting the MLP function to the IMF
Here, we find a set of parameters for the MLP distribution that fit the updated IMF of Chabrier (2005). We use a normalized form of Equation (1) in Chabrier (2005). It is divided into logarithmic bins of mass normalized by with intervals of 0.05 over the range of and we use a Levenberg-Marquardt least-squares minimization method to fit the MLP function. The best fit parameters are , , and . Figure 2 shows the best fit MLP function along with the Chabrier function, both in normalized form. Figure 3 shows three sets of random samples drawn from this best fit MLP function, using Equation 28 and sample sizes of 100, 1000, and 10000, respectively. The samples are made into histograms with logarithmic bin width of , and each histogram is represented by a different symbol. We find that samples of well over 100 are needed to adequately populate the power-law tail. For MLP samples of about 100 or less, fitting a lognormal pdf typically yields a better goodness-of-fit statistic than the MLP distribution, given that it also has one less fitting parameter.
Some sample routines for fitting the MLP function or drawing random samples from it can be found at http://www.astro.uwo.ca/basu/mlp.html.
With a best fit MLP function in hand, it is relatively easy to quantify several important properties of the IMF. Equation 19 can be used to calculate the mean mass in the distribution; it is . If we were to truncate the distribution at , the mean mass is about 10% lower at . The mode of the distribution is calculated numerically to be . However, the usual practice is to plot data as rather than , so the usually quoted IMF peak is that of . We find this to be for our best fit parameters. The mode of is found analogously to that for described in subsection 5.3, but with lowered by one integer value.
Figure 4 shows the cumulative MLP distribution for the best fit parameters. The legend panel in Figure 4 illustrates the specific values for for several relevant mass values, e.g., one solar mass, as well as the minimum stellar mass, (Chabrier & Baraffe, 2000). Although the data used to generate the IMF quoted by Chabrier (2005) does not span the entire hypothetical mass range , one can use the numbers in Figure 4 to quickly estimate some highlights of the mass function, with interesting pairs of marked by symbols in Figure 4: about one quarter of objects are substellar, about one third of objects have masses less than , the median mass is , just over 90% of objects have masses less than , and only 9% of objects have masses in the range . Figure 5 shows the mass fraction (with closed form given in Equation 55) for the best fit parameters, normalized by the value at . Symbols denote interesting pairs of , and many provide an interesting contrast to their counterparts in Figure 4. For example, only about 2% of the total mass is in substellar objects, a majority of mass is tied up in objects more massive than , and half of all mass is in objects above mass .
The ability to model stellar populations with an IMF that has clearly identifiable number or mass fractions in different mass ranges should be a key aid to galaxy studies (e.g., Kennicutt, 1998), for example in the determination of star formation rates, chemical evolution, and mass-to-light ratios.
We have derived several important properties of the modified lognormal power-law (MLP) probability distribution function that has been recently introduced in the literature. The three-parameter MLP function has the salutary properties of a main body that resembles a lognormal, including a peak value and a decline toward low values, as well as a power-law tail at high values. This function can potentially be applied in a variety of fields where empirical distributions may have a power-law tail that coexists with a peak at lower values. The MLP distribution can also help to settle a frequent contentious question: is a data set consistent with a lognormal distribution or does it also show evidence for a power-law tail? Simultaneous fitting of the lognormal and MLP distributions to the same data sets can help to answer this question in many fields.
Comparison of real data sets with the MLP distribution is facilitated by the results presented in this paper. We have derived analytic expressions for the cumulative distribution function, the mean, variance, and higher moments of the distribution. We also derived an approximation for the mode that can serve as an initial guess for nonlinear solvers that can iterate to a more exact solution. The random variate of the MLP distribution has also been introduced. Together, these results can put the MLP distribution on a more equal footing with many classical distributions (e.g., lognormal, exponential, Pareto, Rayleigh) that are frequently used to fit empirical data. The use of the MLP distribution to fit an empirical data set has been demonstrated and we show how useful information about the mean, mode, median, and distributions of number or mass fractions can be easily calculated.
SB acknowledges the hospitality of the Isaac Newton Institute of Mathematical Sciences at Cambridge University during the initial stages of writing of this paper. This work was supported by an NSERC Discovery Grant to SB.
- Abramowitz & Stegun (1964) Abramowitz M., Stegun I. A., 1964, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, Dover, New York, NY
- Adams & Fatuzzo (1996) Adams F. C., Fatuzzo M., 1996, ApJ, 464, 256
- Aitchison & Brown (1957) Aitchison J., Brown J. A. C., 1957, The Lognormal Distribution, Cambridge University Press, London
- Allen et al. (2001) Allen A. P., Li B. L., Charnov E. L., 2001, Ecology Letters, 4, 1
- Basu & Jones (2004) Basu S., Jones C. E., 2004, MNRAS, 347, L47
- Bate & Bonnell (2005) Bate M. R., Bonnell I. A., 2005, MNRAS, 356, 1201
- Chabrier (2003) Chabrier G., 2003, PASP, 115, 763
- Chabrier (2005) Chabrier G., 2005, in Corbelli E., Palla F., Zinnecker H., eds, Astrophysics and Space Science Library, Vol. 327, The Initial Mass Function 50 Years Later, Springer-Verlag, Dordrecht, p. 41
- Chabrier & Baraffe (2000) Chabrier G., Baraffe I., 2000, ARA&A, 38, 337
- Clauset et al. (2009) Clauset A., Shalizi C. R., Newman M. E. J., 2009, SIAM Review, 51, 661
- Crow & Shimizu (1988) Crow E. L., Shimizu K., 1988, Lognormal Distributions: Theory and Applications, CRC Press, New York, NY
- Golberg (1984) Golberg M. A., 1984, An Introduction to Probability Theory with Statistical Applications, Plenum Press, New York, NY
- Gut (2005) Gut A., 2005, Probability: A Graduate Course, Springer, New York, NY
- Hennebelle & Chabrier (2009) Hennebelle P., Chabrier G., 2009, ApJ, 702, 1428
- Kennicutt (1998) Kennicutt Jr R. C., 1998, in Gilmore, G., Howell, D., eds., ASP Conf. Ser. Vol. 142, The Stellar Initial Mass Function, Astron. Soc. Pac., San Francisco, p. 1
- Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
- Kroupa (2002) Kroupa P., 2002, Science, 295, 82
- Limpert et al. (2001) Limpert E., Stahel W. A., Abbt M., 2001, BioScience, 51, 341
- Mandelbrot (1997) Mandelbrot B., 1997, Fractals and Scaling in Finance, Springer-Verlag, New York, NY
- Miller & Scalo (1979) Miller G. E., Scalo J. M., 1979, ApJS, 41, 513
- Mitzenmacher (2006) Mitzenmacher M., 2006, Internet Mathematics, 1, 226
- Myers (2009) Myers P. C., 2009, ApJ, 706, 1341
- Myers (2014) Myers P. C., 2014, ApJ, 781, 33
- Padoan & Nordlund (2002) Padoan P., Nordlund Å., 2002, ApJ, 576, 870
- Pareto (1896) Pareto V., 1896, Cours d’ï¿½conomie Politique, Droz, Geneva
- Reed (2002) Reed W. J., 2002, Journal of Regional Science, 42, 1
- Reed (2003) Reed W. J., 2003, Phys. A, 319, 469
- Salpeter (1955) Salpeter E. E., 1955, ApJ, 121, 161
- Scalo (1998) Scalo J. M., 1998, in Gilmore, G., Howell, D., eds., ASP Conf. Ser. Vol. 142, The Stellar Initial Mass Function, Astron. Soc. Pac., San Francisco, p. 201
- Seggern (1990) Seggern D. H. V., 1990, CRC Handbook of Mathematical Curves and Surfaces, CRC Press, Boca Raton, FL
- Zipf (1949) Zipf G. K., 1949, Human Behavior and the Principle of Least Effort, Addison-Wesley, Reading, MA
Appendix A Error Function
a.1 Definition and Basic Properties
The Error function, denoted by , and the Complementary Error Function, denoted by , are defined as
Note from the definition that is an odd function. Also,
(Abramowitz & Stegun, 1964). The derivative and integral of are:
where is an integration constant. Finally, the Taylor expansion for is given by
A plot of is presented in Figure 6.
a.2 Related Integrals used in our analysis
Consider the integral
By completing the square and letting it can be brought to
where is an integration constant. In particular, using the properties of and its relation to , we have
The result in Equation 36 can also be used to find integrals of the form
where and . By using the substitution , we can bring it to the form
which we can evaluate to
Appendix B Derivation of the Results in section 5
b.1 The MLP as a Probability Density Function
We can explicitly verify that the MLP density (Equation 14) is indeed a valid probability density function for all , i.e., it is positive and integrates to 1. Positivity is clear from the definition. To verify the normalization condition, let , , and . We then have
Letting and , we can write this integral as
Consider the first term on the right hand side in Equation 42. As , (since and is odd), hence for the term vanishes in this limit. As , the limit takes on an indeterminate form. Applying L’Hospital’s rule we can see that the limit is also zero in this case:
where . We can evaluate the second (integral) term with Equation 40:
b.2 Cumulative Distribution
The cumulative distribution function,
and for our case is defined on . To find the closed form we first we apply integration by parts, using also the same definitions for , , and as in subsection B.1:
We again use the identity in Equation 40 to evaluate the integral term and obtain
which, upon returning to the original parameters, becomes
b.3 Raw Moments
Next we derive a closed form for arbitrary raw moments of the distribution, as well as an expression for its variance. Let be an MLP random variable with probability density function . The th raw moment of , defined as the expectation value of , is given by
with , and as defined in subsection B.1. Before we arrive at a closed form expression for the moments we consider the convergence of the integral in Equation 49. This integral diverges for . To see this, note that if then Now write Equation 49 as
where such that . The existence of such is ensured by the fact that for , is a continuous strictly increasing positive function having an upper limit of 2. Then
where the integral on the right hand side diverges for . Thus, the integral in Equation 49 diverges for . Conversely, the moments are finite for . For any ,
b.4 Mass Fraction
We can obtain the expression for the mass fraction (up to a mass ),
using the same approaches employed in the computation of the cumulative distribution and the raw moments in the previous sections. Since the mean exists only if , it is also natural to consider the closed-form expression for the mass fraction for . This can be achieved by replacing all explicit appearances of with in Equation 47, while maintaining the dependencies in the parameters , and in the original form, as we did for the calculation of the raw moments. Thus,
which, upon returning to the original parameters, becomes
Note that, as expected, becomes the expression for the expectation value as (for ).