Kaniadakis Functions beyond Statistical Mechanics: Weakest-Link Scaling, Power-Law Tails, and Modified Lognormal Distribution

Probabilistic models with flexible tail behavior have important applications in engineering and earth science. We introduce a nonlinear normalizing transformation and its inverse based on the deformed lognormal and exponential functions proposed by Kaniadakis. The deformed exponential transform can be used to generate skewed data from normal variates. We apply this transform to a censored autoregressive model for the generation of precipitation time series. We also highlight the connection between the heavy-tailed κ-Weibull distribution and weakest-link scaling theory, which makes the κ-Weibull suitable for modeling the mechanical strength distribution of materials. Finally, we introduce the κ-lognormal probability distribution and calculate the generalized (power) mean of κ-lognormal variables. The κ-lognormal distribution is a suitable candidate for the permeability of random porous media. In summary, the κ-deformations allow for the modification of tails of classical distribution models (e.g., Weibull, lognormal), thus enabling new directions of research in the analysis of spatiotemporal data with skewed distributions.

A feature of particular interest is the behavior of the tail(s) of a probability distribution, because the tails define the probabilities of extreme events. Distributions are characterized as sub-exponential (if their tail decays slower than the exponential) and super-exponential in the opposite case. The same function (e.g., the Weibull model) may exhibit transitions from sub-exponential to super-exponential by changing the value of a key parameter; in the Weibull case this is the modulus m: m < 1 leads to a sub-exponential and m > 1 to a super-exponential tail. Sub-exponential models are called heavy-tailed if the asymptotic • We show (Section 3) that the κ-deformed exponential and logarithmic functions (henceforth, κ-exponential and κ-logarithm) can be used to define normalizing transforms for non-Gaussian data, which extend the well-known (in statistics) Box-Cox family of transformations [46]. • We formulate an autoregressive, intermittent precipitation model based on the κmodified Box-Cox transform in Section 3.3. We show that the resulting precipitation time series has higher "peaks" than those obtained with the Box-Cox transform with the same parameter value. • We review the κ-Weibull distribution focusing on its connection with weakest-link theory (Section 4). This demonstrates that the κ-Weibull is a physically motivated generalization of the classical Weibull distribution for the mechanical strength of brittle materials, unlike modified Weibull distributions which fail to satisfy the weakestlink principle. • We show that for several physical quantities, including the thickness of magmatic sheet intrusions, the tensile strength of steel, earthquake waiting times, and precipitation amounts the κ-Weibull distribution provides a better fit than the Weibull according to model selection criteria. • We introduce the κ-lognormal distribution, which provides a deformation of the lognormal with lighter tails than the latter in Section 5. The κ-lognormal can be used to model asymmetric data distributions which concentrate more probability mass around the median than the lognormal. We discuss the importance of the generalized mean (power mean) of the lognormal distribution for estimating the effective permeability of heterogeneous porous media, and we calculate the generalized mean of the κlognormal distribution.

Mathematical Preliminaries
Let (Ω, F , P), where Ω is the sample space, F is a σ-field of subspaces of Ω, and P is a probability measure, define a probability space. A real-valued, scalar random variable X(ω) is defined by the mapping X : Ω → R, where R is the set of real numbers [47]. Furthermore, a stochastic process X(t; ω) indexed by the time t ∈ R is defined by the mapping X : R × Ω → R. In the following, the dependence on the state index ω is suppressed for convenience. In addition, the random variable X represents a load or the waiting time between consecutive "failure" events (e.g., earthquakes), or some other asymmetrically distributed variable.
The function F(x) : P(X ≤ x) : R → [0, 1] defines the cumulative distribution function (CDF) of X, or the marginal CDF of a stationary stochastic process {X(t)}. The expectation of X, assuming it is mathematically well defined, is given by means of E [X] = dF(x) x.
Assuming that F(x) is at least once differentiable, the probability density function (PDF) is given by the first derivative of the CDF, i.e., f (x) = dF(x)/dx.
The CDF is related to the so-called survival function, also known as the reliability function, which is given by S(x) = 1 − F(x). Whereas F(x) is a monotonically increasing function, S(x) is monotonically decreasing. The term "survival function" comes from reliability engineering: if X represents the strength or critical loading of a given system, F(x) is the probability that the system fails at loading level X ≤ x; then, S(x) is respectively the probability that the system remains intact (survives) at this loading level.
The quantile function, Q(p), where p ∈ [0, 1], returns the value x p ∈ R such that F(x p ) = p. Hence, Q(·) is the inverse of the CDF.
The hazard rate, also known as the hazard function, h(x) represents the conditional probability that the system fails for X ∈ (x, x + δx] where δx 1, conditioned on the survival of the system for X ≤ x. Let A denote the event that the system survives at level x and B denote the event of system failure in the interval [x, x + δx]. Then, by using the definition of conditional probabilities, h(x) = P(B|A) = P(B ∩ A)/P(A), the hazard rate is given by the following ratio: The asymptotic behavior of the hazard rate determines the probability of system failure with increasing load. The hazard rate for x → ∞ is determined by the tail of the probability functions f (x) and S(x). For certain probability models, e.g., the exponential and the gamma, h(x) tends to a constant as x → 0; for the Weibull model with modulus m > 1, h(x) tends to zero, whereas for models with power-law tails, the lognormal, and the Weibull model with m < 1, h(x) diverges as x → ∞. The hazard rate is an important factor in seismic risk assessment [48].

The κ-Exponential Function
The κ-generalized exponential is a one-parameter generalization of the exponential function, proposed by Kaniadakis [31,34]: with 0 ≤ κ < 1 and x ∈ R. The first few terms of the Taylor expansion of exp κ (x), reported in [49], are given by The emerging pattern persists for higher orders, i.e., terms of O(x 3 ) consist of the ordinary exponential expansion and a κ-dependent correction. The κ-exponential is expressed as the following power series [49]: The functions {ξ n (κ)} ∞ n=0 are polynomials of κ defined by the recurrence relations The polynomials ξ n (κ) for the first seven orders are given by It follows from Equation (3) that when x → 0 or κ → 0, exp κ (x) converges to the ordinary exponential, i.e., Equation (2) shows that the asymptotic behavior of exp κ (x) as x → ±∞ follows a power law [34,49], i.e., Based on the above, for x → +∞ the modified exponential exhibits a heavy tail, i.e. exp κ (−x) ∼ (2κx) −1/κ . Hence, exp κ (x) can be used to model subexponential probability distributions which are suitable for heavy-tailed data. The κ-exponential can also be introduced as the solution of a linear, first-order ordinary differential equation (ODE) with time-dependent rate [42]. Consider the ODE where r(x) is the following x-dependent rate function The solution of the ODE (15) is given by the function f (x) = exp κ (−βx). In case κ = 0, then r(x) = β, and the rate equation is solved by the standard exponential function f (x) = exp(−βx).

The κ-Logarithm Function
The inverse of the κ-exponential is the κ-logarithm, defined by the following function for x > 0: The κ-logarithm satisfies the equation ln κ exp κ (x) = x. In addition, it respects the κ-symmetry property ln κ (x) = ln −κ (x). A Taylor expansion of x κ and x −κ around κ = 0 leads to lim The first and second derivatives of the κ-logarithm are respectively given by Based on Equation (19a), the first derivative of the κ-logarithm is positive; therefore ln κ (x) is a monotonically increasing function. Based on Equation (19b), the second derivative of ln κ (x) is negative for 0 < κ < 1; therefore, ln κ (x) is a concave function for κ ∈ (0, 1).

Nonlinear Transformation of Data Based on the κ-Logarithm
Nonlinear, monotonic transformations are often applied to non-Gaussian data in order to restore normality [50][51][52]. This procedure, known as Gaussian anamorphosis, enables the use of data processing methods that are based on Gaussian assumptions. Various transforms are used in practice, including the Box-Cox [46] and Yeo-Johnson [53]. Such transforms can be generalized by means of the Kaniadakis functions. Below we focus on Box-Cox but the same arguments can be used for other transforms.

Box-Cox Transform and the Replica Trick
A widely used normalizing transformation in statistics is the so-called Box-Cox transform (BCT) [46]; the one-parameter version of BCT is given by the monotonic function The BCT is applied to skewed (non-Gaussian) data so that the transformed variable y = g BC (x) is better approximated by the Gaussian distribution. The BCT is applied to both time series and spatial data [52]. It is interesting to note that if λ < 0, g BC (x) takes positive (negative) values for x < 1 (x > 1), whereas if λ > 0, the g BC (x) takes positive (negative) values for x > 1 (x < 1). The BCT value for λ = 0 can be obtained by using either l'Hopital's rule or the Taylor expansion. The Taylor expansion of x λ around λ = 0 shows that Equation (21) shows that the logarithmic transform is a special case of the BCT for λ = 0. The inverse BCT is given by h BC (y) g −1 BC (y) = (λy + 1) 1/λ . In a different context, Edwards and Anderson [54] introduced the famous replica trick, which is also based on Equation (21), to study spin glasses. The replica trick is used to calculate the ensemble average (over the magnetic disorder) of the logarithm of the spin glass partition function, i.e., ln Z. By using the replica trick, ln Z is calculated by first evaluating g BC (Z), with λ = n ∈ N denoting the number of replicas (identical copies of the system), and then taking the limit n → 0.

The κ-Logarithmic Transform
The κ-logarithm transform (KLT) is a nonlinear, monotonic transformation from R + → R. Therefore, it can be used like the BCT for the Gaussian anamorphosis of positivevalued data. The KLT takes the form Equation (18) shows that the logarithmic transform is a special case of the KLT for κ → 0, as it is a special case of the BCT. To understand how the transformation works, consider the special cases κ = λ = 1 shown in Figure 1. The choice λ = 1 for the BCT is a simple linear shift of x to x − 1, whereas setting κ = 1 in the KLT leads to the nonlinear transformation x 2 −1 2x . The latter tends to the linear transformation x/2 for x 1. Figure 2 shows the two transformations for different values of λ = κ. Notice that the KL transformation is symmetric with respect to κ → −κ. The inverse of the KLT is given by h KL (y) = exp κ (y).

Application to Precipitation Modeling
Autoregressive (AR) models are used for different meteorological processes which exhibit memory. AR models exhibit short-term memory: the present depends on the past via the p most recent values of the process where m = E [y t ] is the expectation of y t , the set {φ i } p i=1 comprises the real-valued auto-regressive coefficients, { t } is the innovation process, and σ is its standard deviation. The innovation is typically considered to be standard Gaussian white noise, i.e., t ∼ N (0, 1), leading to a normally distributed time series {y t } T t=1 . Daily precipitation usually displays intermittency (i.e., intervals of zero precipitation) in addition to temporal variability of its intensity. These features can be modeled by using a censored autoregressive time series model [55]. Then, the amount of daily precipitation is given by x t = y t θ(y t − y c ), where θ(·) is the unit step function and y c is a censoring threshold. The threshold is selected so that F y (y c ) = p 0 , where p 0 is the probability of observing a dry day. The censored AR model leads to a truncated normal distribution for x t . Because this distribution does not adequately reflect the extreme values of precipitation, in practice the daily precipitation is given by means of the following censored and transformed autoregressive process, where h(·) : R → R is a monotonically increasing transform. The application of h(·) introduces skewness and increases the weight in the right tail of the distribution. The subtracted term h(0) restores the zero precipitation values after application of the transform. The function h(·) could, for example, represent the inverse BC or KL transforms, i.e., h(y) = (λy + 1) 1/λ or h(y) = exp κ (y), respectively. Figure 3 presents six realizations of the censored AR(1) model of Equation (24) using both the inverse BCT and KLT as the nonlinear transformation h(·), with equal values of κ and λ in each frame. The AR(1) model of Equation (24) is applied with φ 1 = 0.5 and σ = 0.6. The time series exhibit intermittent behavior due to censoring. The difference between BCT and KLT is negligible for small values of κ = λ. For κ ≈ 0, both transforms yield a censored lognormal process, because h KL (y) and h BC (y) converge to the normal exponential for κ = λ → 0. The KLT peaks become progressively higher compared to the BCT peaks as κ (and λ) increase. This behavior is due to the following inequality between the inverse transforms The increase in the relative height of KLT-based versus BCT-based peaks with increasing κ should not be confused with the fact that the peaks of x t are highest for κ = λ = 0, i.e., when h(·) is the exponential function (herein κ = 0 implies the limit κ → 0). This behavior is due to the inequality exp κ (y) < exp(y) for all y ≥ 0 and for 0 < κ ≤ 1, e.g., [49].

The κ-Weibull Distribution and Its Applications
Complex systems involve collections of interacting units and often exhibit probability distributions with power-law tails [23]. A long (power-law) right tail of the PDF implies that the occurrence probability for extreme events decays slowly.
One mechanism that can lead to the emergence of long tails is due to limited size of the observation window as illustrated in Figure 4. Consider an observation domain (indicated by the square domain in the center) which is part of a larger interconnected system (denoted by the oval-shaped area). Let us assume that in the entire system, the failure events occur at times t 1 < t 2 < t 3 < t 4 . Furthermore, assume that the interevent times t 2 − t 1 , t 3 − t 2 , and t 4 − t 3 follow a distribution which does not necessarily have a heavy tail. However, since the observed system involves only the square domain, the observed interevent times are t 2 − t 1 and t 4 − t 2 ; the latter results by adding the true interevent times t 3 − t 2 and t 4 − t 3 , leading to a larger period of quiescence than if the entire system were taken into account. The repeated occurrence of events outside the observed domain can thus inflate waiting times and transfer probability weight from the low and middle range of the probability distribution to the right tail. In such cases, the classical Weibull distribution may not be adequate because it has a right tail which decays at best (i.e., for m < 1) as a stretched exponential. In contrast, the κ-Weibull distribution has a flexible power-law tail with exponent equal to −m/κ.
The κ-Weibull distribution is a deformation of the Weibull model introduced in [39,40] to model the distribution of income in economy. The κ-Weibull has a power-law right tail which captures the observed Pareto law followed by income distributions. The κ-Weibull distribution was later applied to model the mechanical strength of materials and earthquake interevent times [43,44].
The κ-Weibull model admits explicit expressions for the main probability functions which are given by the following expressions [44]: Note that due to the asymptotic behavior of the κ-exponential given by (14), the survival function S κ (x) of the κ-Weibull follows a power law with tail exponent −m/κ: The power-law tail gives the κ-Weibull an advantage over the classical Weibull distribution for systems with algebraic decay of the right tail of the distribution. Figure 5 compares the tails of the survival functions for the Weibull and κ-Weibull models. Note that S κ (x) for κ = 0.1 is practically indistinguishable from the Weibull survival function.

Connection with Weakest-Link Scaling Theory
Weakest-link scaling theory underlies the classical Weibull distribution. Weakest-link scaling was proposed by Gumbel [65] and Weibull [66] in their works on the statistics of extreme values. This section provides (i) a brief review of weakest-link scaling in connection with the Weibull model and (ii) a demonstration that the κ-Weibull is also based on the same principle.
Weakest-link scaling treats a disordered system as a chain of critical clusters, also known as links or representative volume elements (RVEs). The term "weakest-link scaling" emphasizes the idea that the strength of an entire system is determined by the strength of its weakest link [67]. The concept of links implies the presence of critical subsystems. The failure of one such subsystem is presumably sufficient to cause failure of the entire system [21]. Thus, weakest-link scaling implies that the survival probability of the system is equal to the product of the link survival probabilities, i.e., where N eff is the number of links, S N eff (x) is the system survival function, and S The above dependence of the survival function is characteristic of brittle fibers and ceramic materials [19,20,22] and justifies the use of the Weibull strength distribution. It has also been shown that if the strength of the earth's crust follows the Weibull distribution, then the latter is also justified for the distribution of recurrence times between earthquakes under the conditions specified in [18].
Assuming a uniform link survival function S 1 (x), the system's survival function, Equation (27), becomes The number of links can also be expressed as the ratio of the system's volume, V over the link volume, V 0 , i.e., N eff = V/V 0 . To obtain the classical Weibull distribution, the link's survival function is assumed to have the exponential form , where x l is the link's scale parameter and m > 0 is the Weibull modulus or shape parameter [66]. This leads to the system survival function S N eff (x) = exp[−(x/x s ) m ], where x s = x l /N 1/m eff is the scale parameter for the entire system.
The κ-Weibull distribution can be obtained by simply replacing the exponential in the Weibull survival function with the deformed κ-exponential. However, this mathematically valid operation does not provide physical motivation for the κ-Weibull. The latter emerges in the framework of weakest-link scaling by using a modified link survival function. More precisely, let us assume that the link survival function depends on the parameter κ as shown in [43,44], i.e., that it satisfies wherex m l = x m l /κ. The parameter κ can be viewed as the inverse number of effective links, that is, N eff = 1/κ. The link survival S 1 (x) thus depends on the number of links, which implies a degree of interactivity in the system. In addition, the asymptotic behavior of S 1 (x) for x → ∞ is given by a power law, i.e., The link survival function for different system sizes N eff is plotted in Figure 6 against the variable z(x) = x m /x m l . The graphs exhibit the power-law asymptotic decline S 1 (x) ∼ 1/z(x) as well as slower decrease of S 1 (x) for increasing N eff .
Finally, the weakest-link scaling relation (28) in view of the link survival function (29) and the κ-exponential definition (2) leads to an S N eff (x) which is given by the κ survival function (25b).

κ-Weibull Plot for Graphical Testing
If we define the function Φ κ (x) = ln ln κ (1/S κ (x)), it follows from Equation (25b) that Φ κ (x) = m ln(x/x s ). This relation suggests a graphical approach to test if a given dataset {x n } N n=1 follows the κ-Weibull distribution: it suffices to test if the scatter plot of Φ κ (x n ) versus ln x n (for n = 1, . . . , N) is concentrated around a straight line with slope equal to m. This property allows a quick visual test of the fit between the data and the κ-Weibull distribution which is analogous to the widely used Weibull plot [68].
The linear dependence of Φ κ (x) on ln x is illustrated in the κ-Weibull plots of Figure 7. The graphs represent estimates of Φ κ (x) derived from six samples of 500 random numbers; the latter are generated from the κ-Weibull distribution by using the inverse transform sampling method [43]. The samples are drawn from κ-Weibull distributions with x s = 10, m = 0.9 and with x s = 10, m = 2; the samples with the same m value differ with respect to κ which takes values in {0.1, 0.5, 0.9}. The function Φ κ (x) is estimated by means of Φκ(x) = ln lnκ[1/Ŝ(x)], whereŜ(x) is the estimated survival function using the empirical staircase estimate of the CDF from the data, andκ is the maximum likelihood estimate of κ.
Although the use of graphical tools for estimating the tail exponent has an intuitive appeal, these tools can also be misleading. Thus, in the next section maximum likelihood estimation is used to determine the κ-Weibull parameters including the tail exponent. Statistical testing methods (e.g., Kolmororov-Smirnov test) can also be used to validate the hypothesis of a particular probability distribution model. Note that if the distribution parameters are not known a priori, but instead are estimated from the data (as is typically done in practice), Kolmororov-Smirnov testing must be implemented by using a Monte Carlo resampling approach as described in [27]. This testing approach was applied to probability models for earthquake recurrence times in [43].

Application to Real Data
We investigate the κ-Weibull as an alternative to the Weibull distribution for different data. These include a dataset comprising measurements of tensile strength of carbon fibres [69], daily averaged wind speed in Cairo (Egypt) [70], thickness of magmatic sheet intrusions (dykes) for different tectonic settings [71], tensile strength of steel [72], and earthquake recurrence times [73,74].
The aforementioned datasets are fitted to the Weibull and κ-Weibull distributions by using maximum likelihood estimation. The Matlab code used to estimate the model parameters is publicly available [75]. The results of the fits are presented in Table 1. The entries include the maximum likelihood estimates of the model parameters as well as the optimal negative log-likelihood (NLL) values for each fit. The lower the NLL of a given distribution model is, the better its fit to a particular dataset. For all the cases listed in Table 1, the NLL is lower for the κ-Weibull than for the Weibull distribution. However, the κ-Weibull involves three parameters whereas the Weibull model involves only two parameters. To account for the different number of parameters, the selection of the optimal model can be performed by means of the Akaike information criterion (AIC) [76], i.e., AIC = 2k + 2NLL, where k is the number of free parameters for each model (k = 2 for the Weibull and k = 3 for the κ-Weibull). The best model has the lowest AIC value. Because AIC κW − AIC W = 2(1 + NLL κW − NLL W ), AIC favors the κ-Weibull only if NLL κW − NLL W < −1. This condition is satisfied for all but the first two datasets.
A more stringent condition is provided by the Bayesian information criterion, BIC = k log N + 2NLL, where N is the data size [77]. For N > 8, the BIC imposes a bigger penalty on model complexity than AIC. The difference in BIC values for the κ-Weibull and the Weibull models is given by BIC κW − BIC W = log N + 2(NLL κW − NLL W ). Thus, under the BIC the κ-Weibull is optimal if NLL κW − NLL W < −(log N)/2. The BIC condition also favors the κ-Weibull distribution for the datasets 3-9. In order to allow easier comparison, the values of AIC and BIC (divided by the number of sample points) for both the Weibull and κ-Weibull distributions are listed for each dataset in Table 2.

The κ-Lognormal Distribution
The lognormal distribution is often used to model long-tailed processes [23]. In this section we derive a generalization of the lognormal which is based on the κ-deformation of the exponential function. Let us assume that the random variable Y follows the normal distribution with marginal PDF given by We then define the random variable X = g κ (Y) by means of the κ-exponential transformation To determine the PDF of the random variable X we use the standard integral relation for the PDF under the nonlinear transformation g κ (y) Because the nonlinear transform is monotonic and therefore invertible, the PDF of X is given by means of the following mapping [47,52]: Taking into account the normal PDF of Equation (30), the inverse transform which is given by g −1 κ (x) = ln κ (x), and the first derivative of ln κ (·) which is given by Equation (19), it follows that the PDF of the κ-lognormal distribution is given by It is clear that f x (x; κ = 0) recovers the lognormal distribution because lim κ→0 ln κ (x) = ln x and lim κ→0 x κ +x −κ 2 = 1. Moreover, the κ-lognormal PDF given by Equation (32) is symmetric with respect to κ.
The Box-Cox transform of the Gaussian PDF can be similarly obtained by means of the PDF transformation under a change of variables. The resulting PDF is given by The κ-lognormal PDF given by Equation (32) has heavier tails than the PDF in Equation (33) resulting from the Box-Cox transform of the Gaussian. This is confirmed by the parametric plots of the two parametric PDF families shown in Figure 8 for different values of λ = κ. Notice that for λ = 0 (left frame in Figure 8), and κ = 0 (right frame in Figure 8) the PDFs tend to the lognormal PDF. In order to compare the tails of the κ-lognormal with the tails of the lognormal distribution, let us define the PDF ratio R f (x; κ): First, we show that for κ > 0 and for x > 1 it holds that ln κ (x) > ln x. Let us define the function u(x) = ln κ (x) − ln x. (i) It holds that u(1) = ln κ (1) − ln 1 = 0. It suffices to show that u(x) is a monotonically increasing function for x > 1. (ii) The derivative of u(x) is du(x)/dx = (x κ + x −κ − 2)/(2x), and it is positive if x κ + x −κ > 2. By multiplying both sides with x κ (x κ > 0) and setting x κ = y, the inequality x κ + x −κ > 2 becomes equivalent to (y − 1) 2 > 0, which is true for any y ∈ R. Hence, it holds that lim x→∞ R f (x; κ) = 0, and thus the κ-lognormal has a lighter right tail than the lognormal distribution.
The proof of monotonicity of u(x) also holds for x < 1: by replacing x with x = 1/x where x > 1, it holds that x κ + x −κ = (x ) −κ + (x ) κ , and therefore du(x)/dx > 0; that is ln κ (x) > ln x for x < 1. Therefore, lim x→0 R f (x; κ) = 0, and thus the left tail of the κ-lognormal is also lighter than the lognormal's tail. Hence, the κ-lognormal concentrates more density in the middle of the distribution than the lognormal. The relative transfer of density from the tails to the body of the distribution is controlled by κ and becomes more pronounced as κ increases. This is confirmed by the parametric plots of R f (x; κ) shown in Figure 9 (bottom panel) for different values of κ.

Effective Permeability of Random Media
Single-phase, incompressible, steady-state flow in saturated random media is governed by the partial differential equation where s ∈ R d is the position vector, P(s) is the pressure, and ∇P(s) the pressure gradient, a · b denotes the inner product of the vectors a and b, and K(s) is the fluid permeability. The latter is assumed to be a scalar random field, i.e., a random function over the domain D ⊂ R d ; (in the case of anisotropic media K becomes a tensor). Equation (35) is the continuity equation which expresses the conservation of mass.
In the case of slow viscous flow through random media, the fluid velocity is given by Darcy's law, i.e., V(s) = −K(s)∇P(s)/µ where µ is the fluid viscosity [78], and thus the continuity equation becomes ∇ · V(s) = 0. The local variations of the velocity do not usually matter for the macroscopic flow behavior, i.e., for the average fluid velocity through a large domain. The macroscopic velocity is often determined in terms of an effective permeability; the latter connects the average pressure gradient to the average fluid velocity by means of E [V(s)] = −K e f f E [∇P(s)]. The ensemble average here is evaluated over the joint probability distribution of K(s), which represents the local variations (microstructure) of the medium. Similarly, effective measures can be defined for other properties (e.g., elasticity) of porous media by means of averages over the microstructural disorder [79][80][81].
The generalized mean K α with α = 1 − 2/d, also known as the power-law mean, is used to estimate the effective flow permeability (for single-phase, saturated flow), K e f f of random porous media with lognormal disorder and short-range correlations [16]. For a given power-law exponent α, the generalized mean is defined by means of the following expectation (assuming that it is well-defined mathematically): For α = −1, Equation (36) yields the harmonic mean, whereas for α = 0 the geometric mean, K G = exp[E (ln K)], is obtained as the limit lim α→0 K α . Furthermore, assuming that K = g(Y) where Y ∼ N m, σ 2 is a normally distributed random variable. In general, both Y and K can be random fields with spatial correlations. However, this does not affect the calculation of the generalized mean, which is a point (marginal) property. The generalized mean of K = g(Y) is now given by Thus, if K follows the lognormal distribution, namely, K = exp(Y), it holds that This equation, known as the Landau-Lifshitz-Matheron (LLM) ansatz, was first proposed for the dielectric permittivity of random dielectric mixtures [82]. In the case of electromagnetism, the continuity equation is embodied in Gauss's law; for zero free charge density the latter becomes ∇ · D(s) = 0, where D(s) = − (s) ∇φ(s) is the dielectric displacement field, φ(s) is the applied electric potential, and (s) is the dielectric permittivity of the medium. Notation differences aside, the mathematical form of the electrostatic equations is identical to those of the fluid flow problem; in both cases the continuity principle results in an elliptic partial differential equation with suitable boundary conditions. This reason underlies LLM's applicability to both fluid permeability and dielectric permittivity; as Richard Feynman wrote ( [83], Chap. 12.1): "there is a most remarkable coincidence: The equations for many different physical situations have exactly the same appearance". For d = 1 (e.g., for pipe flows), LLM yields the harmonic mean K H = 1/E [K −1 ]; in this case the flow is cut off if the permeability vanishes at a single point. For d = 2, K e f f = K 0 yields the geometric mean, K G = exp(E [ln K]), which coincides with the exact solution in two dimensions [84]. Finally, in d = 3 the expression K 1/3 follows from perturbative renormalization group analysis [16]. The LLM equation implies that the effect of disorder (as measured by the log-permeability variance σ 2 ) is reduced as the embedding dimension of the medium increases. The physical meaning of this dependence is that three-dimensional media include more permeable paths than one-or two-dimensional random media, thus enabling the bypass of flow bottlenecks.

Generalized Mean of the κ-Lognormal Distribution
This section focuses on the calculation of the generalized mean when K follows the κ-lognormal distribution. In this case, the respective ensemble average over the normal variable Y ∼ N m, σ 2 , defined by Equation (37), is given by where in deriving the last equality above we used the identity Note that for κ = 0 Equation (39) is equivalent to Equation (38) because the κ-exponential is reduced to the standard exponential. If, on the other hand, κ > 0, it holds that κ = 0. The expectation in Equation (39) can then be calculated by using the Taylor series expansion of exp κ/α (αY) around exp(αY); according to Equation (4), the expansion is given by and the polynomials ξ n (·) are defined in Equation (5). The power series in Equation (41) represents the correction of exp κ (αy) with respect to exp(αy). Based on Equation (7), it holds that exp κ (αy) = exp(αy) + O(κ 2 ). Notice that ξ n (κ ) − 1 = 0 for n = 0, 1, 2, whereas ξ n (κ ) − 1 < 0 for n = 3, 4. Hence, Hence, the expectation of the above is given by means of the following expression: The expectation E Y [Y n ] is given by means of the following expression by using the fluctuation Y = Y − m and Newton's binomial formula: The expectation over the fluctuations can be calculated by using the Wick-Isserlis theorem [52] E Y Y 2 = (2 )! σ 2 /(2 !), and E Y Y 2 +1 = 0 for ∈ N. Therefore the difference between the generalized mean of the κ-lognormal and the generalized mean of the lognormal (for the same value of α) is given by an infinite power series (correction factor), where δ j,2 is the Kronecker delta, as follows: Finally, returning to Equation (39) and using Equation (44) the generalized mean of the κ-lognormal distribution is given by The convergence of the power series in Equation (45) should be further investigated mathematically. To gain some insight into Equation (45), consider the case α = 1, which corresponds to the arithmetic (linear) mean. Then, the arithmetic mean of the κ-lognormal, i.e., K κ K α=1;κ is given by Hence, the arithmetic mean of the κ-lognormal is given by the standard arithmetic mean plus a correction factor which involves a double power series. Notice that when κ = 0, according to Equation (5) it holds ξ n (κ = 0) = 1 for all n ∈ N; hence, Equation (46) recovers the arithmetic mean when κ = 0 because the double power series vanishes.
On a more practical note, a numerical calculation of the generalized mean shows that the infinite series in K α;κ converges for −1 ≤ α ≤ 1 and 0 < κ < 1. Figure 10 shows parametric plots of the generalized mean obtained by a numerical evaluation of Equation (39). The expectations (for different α and κ) are calculated by using an ensemble of 10 4 random variates drawn from the standard normal distribution, i.e., Y ∼ N(0, 1).  The plots on the left of Figure 10 display the generalized mean as a function of κ for different values of the averaging exponent α. As evidenced in these plots, the difference between the harmonic mean (α = −1) and the arithmetic mean (α = 1) is reduced as κ increases. This behavior is due to the smaller tail weight of the κ-lognormal PDF for κ ↑. It is also observed that the geometric mean (α = 0) is independent of κ (for numerical reasons we use α = 0.001 instead of α = 0). This is also understood in light of Equation (39). Finally, the arithmetic mean (more generally, the generalized mean for α > 0) decays with increasing κ whereas the harmonic mean (more generally, the generalized mean for α < 0) increases. This behavior is due to the fact that the arithmetic mean reflects the diminishing right tail of the κ-lognormal as κ increases. On the other hand, the harmonic mean is more strongly influenced by the lower part of the distribution; according to Figure 9 the κ-lognormal has higher density in the left tail (except for values very close to zero) for higher values of κ.
The plots on the right of Figure 10 display the generalized mean as a function of α for different values of κ. All the curves exhibit an increase of the generalized mean with increasing α. This reflects the progression from the harmonic to the arithmetic mean according to the well-known ordering K H ≤ K G ≤ K. All the curves intersect at α = 0, marking the independence of the generalized geometric mean on κ. Finally, the slope of the curves is reduced with increasing κ due to the respective shrinking of the tails of the κ-lognormal.

Discussion
Section 3 introduced and investigated the properties of a nonlinear normalizing transform which is based on the κ-logarithm. The proposal generalizes the Box-Cox transform, and it can be used for normalizing skewed data before geostatistical methods are applied. An application to precipitation time series modeling was presented. Note that the nonlinear κ-logarithm transform could also be used in spatial models of precipitation in the framework of the censored latent Gaussian field approach [10]. At fine spatiotemporal scales, the correlations of dry/wet spells as well as storm autocorrelation patterns can be better captured by means of two-state models that use copulas to simulate the dependence structure [85]. We believe that the κ-Weibull and κ-lognormal distributions discussed herein will be useful in the framework of two-state models as well, e.g., for modeling the intensity of wet spells. Section 4 shows that the κ-Weibull model, which is a deformation of the classical Weibull with a power-law right tail, also respects the principle of weak-link scaling. Various extensions of the Weibull model have been proposed in the scientific literature, e.g., [86,87]. However, some of these models provide deformations of the Weibull expression that fail the weakest scaling equation (see Equation (28)) as pointed out by Zok [22]. Such modifications, although mathematically permissible, lack the physical justification of the classical Weibull model which is based on weakest-link scaling. Moreover, we established that the link statistical properties depend on the number of links in the system. This feature indicates a strongly interacting system; alternatively, it shows that the observed system is a part of a larger system.
The asymmetric κ-lognormal distribution introduced in Section 5 has lighter tails than the lognormal distribution. The deviation from the lognormal is controlled by the parameter κ. Smaller values of κ ≈ 0 imply small deviations, whereas larger values of κ ≈ 1 signify thinner tails than the lognormal. The κ-lognormal can be used as a model of fluid permeability for random porous media. We believe that the stochastic theory of single-phase, saturated fluid flow in κ-lognormal media can be derived, at least in the framework of perturbation analysis, and expressions for the effective permeability can be likewise obtained. An interesting question is how the parameter κ which controls tail behavior will impact flow properties.
We have reviewed the generalized mean for the lognormal distribution and its application in the estimation of the effective permeability of random media. We then studied the generalized mean for the κ-lognormal distribution. Note that our calculations do not prove that the effective permeability of random media with κ-lognormal disorder is given by the generalized mean. This intriguing hypothesis should be further explored in the framework of the stochastic theory of flow and transport [80].

Conclusions
We presented applications of the Kaniadakis κ-exponential and κ-logarithm functions in the modeling of mechanical strength and in earth science problems. In particular, we focused on κ deformations of classical distribution models such as the Weibull and the lognormal. The κ-Weibull distribution has a power-law tail which is useful for the modeling of mechanical strength, earthquake recurrence times, and properties of geological structures, among other applications. On the other hand, the κ-lognormal model has a tail lighter than the lognormal; this feature is of interest for skewed distributions which decline faster than the lognormal. The methodological applications of the Kaniadakis functions presented include the following: • The modified Box-Cox transform given by Equation (22). • Application of the modified Box-Cox transform to an autoregressive, intermittent model of precipitation as described in Section 3.3. • Connection between the κ-Weibull probability model with the theory of weakest-link scaling as shown in Section 4.2. • The study of the κ-lognormal distribution which is a generalization of the lognormal model with lighter tails. The PDF of this new model is given by Equation (32). • The calculation of the power-mean (generalized mean) of the κ-lognormal as shown in Section 5.2.
Further study of probability models based on the deformed exponential and logarithmic functions will lead to significant advances in different fields of earth science. The most obvious applications at this time include (i) modeling the mechanical strength of technological materials and geologic media, earthquake recurrence times, wind speed, and precipitation amounts, (ii) nonlinear transforms used for Gaussian anamorphosis in geostatistical and ensemble Kalman filtering applications [88], and (iii) the permeability of random porous media.