Metalog distribution

The history of probability distributions can be viewed, in part, as a progression of developments towards greater flexibility in shape and bounds when fitting to data. The normal distribution was first published in 1756,[3] and Bayes’ theorem in 1763.[4] The normal distribution laid the foundation for much of the development of classical statistics. In contrast, Bayes' theorem laid the foundation for the state-of-information, belief-based probability representations. Because belief-based probabilities can take on any shape and may have natural bounds, probability distributions flexible enough to accommodate both were needed. Moreover, many empirical and experimental data sets exhibited shapes that could not be well matched by the normal or other continuous distributions. So began the search for continuous probability distributions with flexible shapes and bounds.

Early in the 20th century, the Pearson[5] family of distributions, which includes the normal, beta, uniform, gamma, student-t, chi-square, F, and five others,[6] emerged as a major advance in shape flexibility. These were followed by the Johnson[7][8] distributions. Both families can represent the first four moments of data (mean, variance, skewness, and kurtosis) with smooth continuous curves. However, they have no ability to match fifth or higher-order moments. Moreover, for a given skewness and kurtosis, there is no choice of bounds. For example, matching the first four moments of a data set may yield a distribution with a negative lower bound, even though it might be known that the quantity in question cannot be negative. Finally, their equations include intractable integrals and complex statistical functions, so that fitting to data typically requires iterative methods.

Early in the 21st century, decision analysts began working to develop continuous probability distributions that would exactly fit any specified three points on the cumulative distribution function for an uncertain quantity (e.g., expert-elicited ${\displaystyle 0.10,0.50}$, and ${\displaystyle 0.90}$ quantiles). The Pearson and the Johnson family distributions were generally inadequate for this purpose. In addition, decision analysts also sought probability distributions that would be easy to parameterize with data (e.g., by using linear least squares, or equivalently, multiple linear regression). Introduced in 2011, the class of quantile-parameterized distributions (QPDs) accomplished both goals. While being a significant advance for this reason, the QPD originally used to illustrate this class of distributions, the Simple Q-Normal distribution,[9] had less shape flexibility than the Pearson and Johnson families, and lacked the ability to represent semi-bounded and bounded distributions. Shortly thereafter, Keelin[1] developed the family of metalog distributions, another instance of the QPD class, which is more shape-flexible than the Pearson and Johnson families, offers a choice of boundedness, has closed-form equations that can be fit to data with linear least squares, and has closed-form quantile functions, which facilitate Monte Carlo simulation.

The metalog distribution is a generalization of the logistic distribution, where the term "metalog" is short for "metalogistic". Starting with the logistic quantile function, ${\displaystyle x=Q(y)=\mu +s{\mbox{ }}\ln {\Bigl (}{y \over {1-y}}{\Bigr )}}$, Keelin substituted power series expansions in cumulative probability ${\displaystyle y=F(x)}$ for the ${\displaystyle \mu }$ and the ${\displaystyle s}$ parameters, which control location and scale, respectively.[10]

${\displaystyle \mu =a_{1}+a_{4}(y-0.5)+a_{5}(y-0.5)^{2}+a_{7}(y-0.5)^{3}+a_{9}(y-0.5)^{4}+\dots }$

${\displaystyle s=a_{2}+a_{3}(y-0.5)+a_{6}(y-0.5)^{2}+a_{8}(y-0.5)^{3}+a_{10}(y-0.5)^{4}+\dots }$

Keelin's rationale for this substitution was fivefold.[10] First, the resulting quantile function would have significant shape flexibility, governed by the coefficients ${\displaystyle a_{i}}$. Second, it would have a simple closed form that is linear in these coefficients, implying that they could easily be determined from CDF data by linear least squares. Third, the resulting quantile function would be smooth, differentiable, and analytic, ensuring that a smooth, closed-form PDF would be available. Fourth, simulation would be facilitated by the resulting closed-form inverse CDF. Fifth, like a Taylor series, any number of terms ${\displaystyle k}$ could be used, depending on the degree of shape flexibility desired and other application needs.

Note that the subscripts of the ${\displaystyle a}$-coefficients are such that ${\displaystyle a_{1}}$ and ${\displaystyle a_{4}}$ are in the ${\displaystyle \mu }$ expansion, ${\displaystyle a_{2}}$ and ${\displaystyle a_{3}}$ are in the ${\displaystyle s}$ expansion, and subscripts alternate thereafter. This ordering was chosen so that the first two terms in the resulting metalog quantile function correspond to the logistic distribution exactly; adding a third term with ${\displaystyle a_{3}\neq 0}$ adjusts skewness; adding a fourth term with ${\displaystyle a_{4}\neq 0}$ adjusts kurtosis primarily; and adding subsequent non-zero terms yields more nuanced shape refinements.[10]^: p.252

Rewriting the logistic quantile function to incorporate the above substitutions for ${\displaystyle \mu }$ and ${\displaystyle s}$ yields the metalog quantile function, for cumulative probability ${\displaystyle 0<y<1}$.

${\displaystyle M_{k}(y)=\left\{{\begin{array}{ll}a_{1}+a_{2}\ln {\Bigl (}{y \over {1-y}}{\Bigr )}&{\mbox{for }}k=2\\a_{1}+a_{2}\ln {\Bigl (}{y \over {1-y}}{\Bigr )}+a_{3}(y-0.5)\ln {\Bigl (}{y \over {1-y}}{\Bigr )}&{\mbox{for }}k=3\\a_{1}+a_{2}\ln {\Bigl (}{y \over {1-y}}{\Bigr )}+a_{3}(y-0.5)\ln {\Bigl (}{y \over {1-y}}{\Bigr )}+a_{4}(y-0.5)&{\mbox{for }}k=4\\M_{k-1}(y)+a_{k}(y-0.5)^{k-1 \over 2}&{\mbox{for odd }}k\geq 5\\M_{k-1}(y)+a_{k}(y-0.5)^{{k \over 2}-1}\ln {\Bigl (}{y \over {1-y}}{\Bigr )}&{\mbox{for even }}k\geq 6\\\end{array}}\right.}$

Equivalently, the metalog quantile function can be expressed in terms of basis functions: ${\displaystyle M_{k}(y)=\sum _{i=1}^{k}a_{i}g_{i}(y)}$, where the metalog basis functions are ${\displaystyle g_{1}(y)=1,{\mbox{ }}g_{2}(y)=\ln {\Bigl (}{y \over {1-y}}{\Bigr )},}$ and each subsequent ${\displaystyle g_{i}(y)}$ is defined as the expression that is multiplied by ${\displaystyle a_{i}}$ in the equation for ${\displaystyle M_{k}(y)}$ above. Note that coefficient ${\displaystyle a_{1}}$ is the median, since all other terms equal zero when ${\displaystyle y=0.5}$. Special cases of the metalog quantile function are the logistic distribution (${\displaystyle k=2}$) and the uniform distribution (${\displaystyle k\geq 4,a_{1}=0.5,a_{4}=1,a_{i}=0}$ otherwise).

Differentiating ${\displaystyle x=M_{k}(y)}$ with respect to ${\displaystyle y}$ yields the quantile density function[11] ${\displaystyle q(y)=dx/dy}$. The reciprocal of this quantity, ${\displaystyle (q(y))^{-1}=dy/dx=f(Q(y))}$, is the probability density function expressed as a p-PDF,[12]

${\displaystyle m_{k}(y)=\left\{{\begin{array}{ll}{y(1-y) \over {a_{2}}}&{\mbox{for }}k=2\\{\Bigl (}{a_{2} \over {y(1-y)}}+a_{3}{\Bigl (}{y-0.5 \over {y(1-y)}}+\ln {y \over {1-y}}{\Bigr )}{\Bigr )}^{-1}&{\mbox{for }}k=3\\{\Bigl (}{a_{2} \over {y(1-y)}}+a_{3}{\Bigl (}{y-0.5 \over {y(1-y)}}+\ln {y \over {1-y}}{\Bigr )}+a_{4}{\Bigr )}^{-1}&{\mbox{for }}k=4\\{\Bigl (}{1 \over {m_{k-1}(y)}}+a_{k}{k-1 \over 2}(y-0.5)^{(k-3)/2}{\Bigr )}^{-1}&{\mbox{for odd }}k\geq 5\\{\Bigl (}{1 \over {m_{k-1}(y)}}+a_{k}{\Bigl (}{(y-0.5)^{k/2-1} \over {y(1-y)}}+({k \over 2}-1)(y-0.5)^{(k/2-2)}\ln {y \over {1-y}}{\Bigr )}{\Bigr )}^{-1}&{\mbox{for even }}k\geq 6,\\\end{array}}\right.}$

which may be equivalently expressed in terms of basis functions as

${\displaystyle m_{k}(y)=\left(\sum _{i=1}^{k}a_{i}{{dg_{i}(y)} \over {dy}}\right)^{-1}}$ where ${\displaystyle 0<y<1}$.

Note that this PDF is expressed as a function of cumulative probability, ${\displaystyle y}$, rather than variable of interest, ${\displaystyle x}$. To plot the PDF (e.g., as shown in the figures on this page), one can vary ${\displaystyle y\in (0,1)}$ parametrically, and then plot ${\displaystyle x=M_{k}(y)}$ on the horizontal axis and ${\displaystyle m_{k}(y)}$ on the vertical axis.

Based on the above equations and the following transformations that enable a choice of bounds, the family of metalog distributions is composed of unbounded, semibounded, and bounded metalogs, along with their symmetric-percentile triplet (SPT) special cases.

As defined above, the metalog distribution is unbounded, except in the unusual special case where ${\displaystyle a_{i}=0}$ for all terms that contain ${\displaystyle \ln {\Bigl (}{y \over {1-y}}{\Bigr )}}$. However, many applications require flexible probability distributions that have a lower bound ${\displaystyle b_{l}}$, an upper bound ${\displaystyle b_{u}}$, or both. To meet this need, Keelin used transformations to derive semi-bounded and bounded metalog distributions.[1] Such transformations are governed by a general property of quantile functions: for any quantile function ${\displaystyle x=Q(y)}$ and increasing function ${\displaystyle z(x),x=z^{-1}(Q(y))}$ is also a quantile function.[13] For example, the quantile function of the normal distribution is ${\displaystyle x=\mu +\sigma \Phi ^{-1}(y)}$; since the natural logarithm, ${\displaystyle z(x)=\ln(x-b_{l})}$, is an increasing function, ${\displaystyle x=b_{l}+e^{\mu +\sigma \Phi ^{-1}(y)}}$ is the quantile function of the lognormal distribution. Analogously, applying this property to the metalog quantile function ${\displaystyle M_{k}(y)}$ using the transformations below yields the semi-bounded and bounded members of the metalog family. By considering ${\displaystyle z(x)}$ to be metalog-distributed, all members of the metalog family meet Keelin and Powley's[9] definition of a quantile-parameterized distribution and thus possess the properties thereof.

${\displaystyle {\begin{array}{l|c|c|c|c|c}&&&&&{\mbox{shape}}\\&{\mbox{transformation}}&{\mbox{quantile function}}&{\mbox{PDF}}&{\mbox{where}}&{\mbox{parameters}}\\\hline {\mbox{Metalog (unbounded)}}&z(x)=x&M_{k}(y)&m_{k}(y)&0<y<1&k-2\\\hline {\mbox{Log metalog}}&z(x)=\ln(x-b_{l})&M_{k}^{\log(}y)=b_{l}+e^{M_{k}(y)}&m_{k}^{\log(}y)=m_{k}(y)e^{-M_{k}(y)}&0<y<1&k-1\\{\mbox{(bounded below)}}&&=b_{l}&=0&y=0\\\hline {\mbox{Negative-log metalog}}&z(x)=-\ln(b_{u}-x)&M_{k}^{\operatorname {nlog} }(y)=b_{u}-e^{-M_{k}(y)}&m_{k}^{\operatorname {nlog} }(y)=m_{k}(y)e^{M_{k}(y)}&0<y<1&k-1\\{\mbox{(bounded above)}}&&=b_{u}&=0&y=1\\\hline {\mbox{Logit metalog}}&z(x)=\ln {\Bigl (}{x-b_{l} \over {b_{u}-x}}{\Bigr )}&M_{k}^{\operatorname {logit} }(y)={b_{l}+b_{u}e^{M_{k}(y)} \over {1+e^{M_{k}(y)}}}&m_{k}^{\operatorname {logit} }(y)=m_{k}(y){{\bigl (}1+e^{M_{k}(y)}{\bigr )}^{2} \over {(b_{u}-b_{l})e^{M_{k}(y)}}}&0<y<1&k\\{\mbox{(bounded)}}&&=b_{l}&=0&y=0\\&&=b_{u}&=0&y=1\end{array}}}$

Note that the number of shape parameters in the metalog family increases linearly with the number of terms ${\displaystyle k}$. Therefore, any of the above metalogs may have any number of shape parameters. By contrast, the Pearson and Johnson families of distributions are limited to two shape parameters.

Bounded SPT metalog parameterized with CDF data ${\displaystyle (20,0.1),(30,0.5),}$ and ${\displaystyle (50,0.9)}$ and with lower and upper bounds ${\displaystyle 0}$ and ${\displaystyle 100}$ respectively.

The symmetric-percentile triplet (SPT) metalog distributions are a three-term ${\displaystyle (k=3)}$ special case of the unbounded, semi-bounded, and bounded metalog distributions.[14] These are parameterized by the three ${\displaystyle (x,y)}$ points off the CDF curve, of the form ${\displaystyle (x_{1},\alpha )}$, ${\displaystyle (x_{2},0.5)}$, and ${\displaystyle (x_{3},1-\alpha )}$, where ${\displaystyle 0<\alpha <0.5}$. SPT metalogs are useful when, for example, quantiles ${\displaystyle (x_{1},x_{2},x_{3})}$ corresponding to the CDF probabilities (e.g. ${\displaystyle 0.1,0.5,0.9}$) are elicited from an expert and used to parameterize the three-term metalog distributions. As noted below, certain mathematical properties are simplified by the SPT parameterization.

The metalog family of probability distributions has the following properties.

Fitting to data

10-term log metalog distribution over maximum annual river gauge height (ft) from 1920 to 2014 for the Williamson River below Sprague River confluence, Chiloquin, Oregon. Data source: USGS.

The coefficients ${\displaystyle {\boldsymbol {a}}}$ can be determined from data by linear least squares. Given ${\displaystyle n}$ data points ${\displaystyle (x_{i},y_{i})}$ that are intended to characterize a metalog CDF, and an ${\displaystyle n\times k}$ matrix ${\displaystyle {\boldsymbol {Y}}}$ whose elements consist of the basis functions ${\displaystyle g_{j}(y_{i})}$, then as long as ${\displaystyle {\boldsymbol {Y}}^{T}{\boldsymbol {Y}}}$ is invertible, the column vector ${\displaystyle {\boldsymbol {a}}}$ of the coefficients ${\displaystyle a_{i}}$ is given by ${\displaystyle {\boldsymbol {a}}=({\boldsymbol {Y}}^{T}{\boldsymbol {Y}})^{-1}{\boldsymbol {Y}}^{T}{\boldsymbol {z}}}$, where ${\displaystyle n\geq k}$ and column vector ${\displaystyle {\boldsymbol {z}}=(z(x_{1}),\ldots ,z(x_{n}))}$. If ${\displaystyle n=k}$, this equation reduces to ${\displaystyle {\boldsymbol {a}}={\boldsymbol {Y}}^{-1}{\boldsymbol {z}}}$, where the resulting metalog CDF runs through all data points exactly. For SPT metalogs, it further reduces to expressions in terms of the three ${\displaystyle (x_{i},y_{i})}$ points directly.[14]

An alternate fitting method, implemented as a linear program, determines the coefficients by minimizing the sum of absolute distances between the CDF and the data, subject to feasibility constraints.[19]

How metalogs converge to standard normal distribution as ${\displaystyle k}$ increases from 2 to 10

Weibull distributions (blue) closely approximated by nine-term semi-bounded metalog distributions (dashed, yellow)

For 3,474 steelhead trout caught and released on the Babine River in British Columbia during 2006-2010, empirical weight data (histogram) and 10-term log metalog PDF (blue curve) fit to this data by least squares.

Due to their shape and bounds flexibility, metalogs can be used to represent empirical or other data in virtually any field of human endeavor.

• Astronomy. Metalogs were applied to assess the risks of asteroid impact.[26]
• Cybersecurity. Metalogs were used in cyber security risk assessment.[19][27]
• Eliciting and combining expert opinion. Statistics Canada elicited expert opinions on future Canadian fertility rates from 18 experts, which included the use of spreadsheet-based real-time PDF feedback based on five-term metalogs. The individual expert opinions were then weighted and combined into an overall metalog-based forecast.[24]
• Empirical data exploration and visualization. In fish biology, a 10-term log metalog distribution (bounded below at 0) was fit to the weights of 3,474 steelhead trout caught and released on the Babine River in British Columbia during 2006–2010. The bimodality of the resulting distribution has been attributed to the presence of both first-time and second-time spawners in the river, the latter of which tend to weigh more.[28]
• Hydrology. A 10-term semi-bounded metalog was used to model the probability distribution of annual river gauge heights.[29]
• Oil field production. Semi-bounded SPT metalogs were used to analyze biases in projections of oil-field production when compared to observed production after the fact.[30]
• Portfolio management. SPT metalogs have been used to model commercial value of new products and product portfolios.[31]
• Simulation input distributions. To support a bidding decision, uncertainty about the future value of each of 259 financial assets was represented as an SPT metalog. A simulation of total portfolio value was shown to yield more realistic results than a corresponding simulation based on discrete low, median, and high values for each asset.[32]
• Simulation output distributions. Metalogs have also been used to fit output data from simulations in order to represent those outputs as closed-form continuous distributions (both CDFs and PDFs). Used in this way, they are typically more stable and smoother than histograms.[32]
• Sums of lognormals. Metalogs enable a closed-form representation of known distributions whose CDFs have no closed-form expression. Keelin et al. (2019)[18] apply this to the sum of independent identically distributed lognormal distributions, where quantiles of the sum can be determined by a large number of simulations. Nine such quantiles are used to parameterize a semi-bounded metalog distribution that runs through each of these nine quantiles exactly. Quantile parameters are stored in a table, which can then be interpolated to yield in-between values; these values are guaranteed to be feasible by the convexity property above.

Metalog panel for steelhead weight data

For a given application and data set, choosing the number of metalog terms ${\displaystyle k}$ depends on context and may require judgment. For expert elicitation, three to five terms is usually sufficient. For data exploration and matching other probability distributions such as the sum of lognormals, eight to 12 terms is usually sufficient. A metalog panel, which displays the metalog PDFs corresponding to differing numbers of terms ${\displaystyle k}$ for a given data set, may aid this judgment. For example, in the steelhead weight metalog panel,[1] using less than seven terms arguably underfits the data by obscuring the data's inherent bimodality. Using more than 11 terms is unnecessary and could, in principle, overfit the data. The case with 16 terms is infeasible for this data set, as indicated by the blank cell in the metalog panel. Other tools (such as regularization, Akaike information criterion, and Bayesian information criterion) may also be useful. For example, when applied to the steelhead weight data, the AIC ranking of metalog distributions from 2-16 terms along with a wide range of classical distributions identifies the 11-term log metalog as the best fit to this data. A similar BIC ranking identifies the 10-term log metalog as the best fit. Keelin (2016)[1] offers further perspectives on the distribution selection within the metalog family.[33]

The metalog distributions belong to the group of distributions defined in terms of the quantile function, which include the quantile-parameterized distributions, the Tukey lambda distribution, its generalization, GLD,[34] the Govindarajulu distribution[35] and others.[13] The following distributions are subsumed within the metalog family:

• The logistic distribution is a special case of the unbounded metalog where ${\displaystyle a_{i}=0}$ for all ${\displaystyle i>2}$.
• The uniform distribution is a special case of: 1) the unbounded metalog where ${\displaystyle k\geq 4}$, ${\displaystyle a_{1}=0.5}$, ${\displaystyle a_{4}=1}$ and ${\displaystyle a_{i}=0}$ otherwise; and 2) the bounded metalog where ${\displaystyle k\geq 2}$, ${\displaystyle b_{l}=0}$, ${\displaystyle b_{u}=1}$, ${\displaystyle a_{2}=1}$, and ${\displaystyle a_{i}=0}$ otherwise.
• The log-logistic distribution, also known as the Fisk distribution in economics, is a special case of the log metalog where ${\displaystyle b_{l}=0}$, and ${\displaystyle a_{i}=0}$ for all ${\displaystyle i>2}$.
• The log-uniform distribution is a special case of the log metalog where ${\displaystyle k\geq 4}$, ${\displaystyle a_{1}=0.5}$, ${\displaystyle a_{4}=1}$, and ${\displaystyle a_{i}=0}$ otherwise.
• The logit-logistic distribution[36] is a special case of the logit metalog where ${\displaystyle a_{i}=0}$ for all ${\displaystyle i>2}$.

Freely available software tools can be used to work with metalog distributions:

• Excel workbooks. By pasting or typing in CDF data, metalogs (with choice of bounds) are instantly displayed.
  • SPT metalogs workbook[37] calculates 2–3 term metalogs determined by three ${\displaystyle (x_{i},y_{i})}$ CDF data.
  • Metalogs workbook[38] calculates 2–16 term metalogs (including metalog panel) determined by 2-10,000 ${\displaystyle (x_{i},y_{i})}$ CDF data.
  • ELD (equally likely data) Metalog workbooks[39] calculate 2–16 term metalogs determined by 2–10,000 ${\displaystyle (x_{i})}$ CDF data, where ${\displaystyle y_{i}}$'s and metalog panel are automatically calculated.
• R. rmetalog[40] (on Complehensive R Archive Network, CRAN).
• Python. Pymetalog[41] closely mirrors the R package. Metalogistic[42] takes advantage of the SciPy platform.
• MakeDistribution.com[43] facilitates experimentation with metalogs parameterized by several CDF data points. The SPT metalog calculator,[44] metalog calculator[45] and ELD metalog calculator[46] are online versions of the Excel Workbooks.
• SIPmath Modeler Tools[47] support metalog distributions in an Excel add-In for simulation.
• Lumina's Analytica Free 101 software[48] for modeling and aiding difficult decisions.
• BayesFusion's Metalog Builder [49] allows for interactive building of metalog distributions. BayesFusion's GeNIe [50] (academic version of the software is free for academic research and teaching) implements the metalog distributions.

Commercially available packages also support the use of metalog distributions:

• FrontLine Solvers: Analytic Solver, RASON, and Solver SDK,[51] software for optimization. Automatically fits user data to the full range of (bounded and unbounded, multi-term) metalog distributions and provides option to compare metalog distributions with classical distributions based on user-selected goodness of fit criteria.
• Lone Star Analysis: TruNavigator and AnalyticsOS software [52] for predictive and prescriptive analytics.

• The Metalog Distributions website, www.metalogs.org
• The Metalog Distributions YouTube channel, educational videos