Close Form Design Quantiles Under Skewness and Kurtosis: A Hermite Approach to Structural Reliability

1. Introduction

In structural reliability, EN 1990 [1] adopts a simplified approach that uses only the mean and standard deviation to derive design values via the First Order Reliability Method (FORM). The Standard [1] introduces three basic probability distributions: the Gaussian, the lognormal, and the Gumbel probability density functions (pdfs). Structural resistance is most often modelled using the Gaussian [2,3] or lognormal [4,5] pdfs, while the Gumbel pdf is applied to represent extremes of dynamic response [6,7] or long-term, single-variable load actions [8,9,10,11]. In practice, EN 1990 does not prescribe a specific mapping of variables to probability distributions and provides closed-form design expressions only for these three two-parameter families, see Figure 1. The target reliability must therefore be demonstrated using a distribution and parameter values that reproduce the observed data. However, two-parameter pdfs may be inadequate for this purpose, see, e.g., [12,13,14].

Numerous studies have shown that the actual probability distribution of structural resistance often departs from simple Gaussian or lognormal models [15,16]. For example, recent material tests indicate that different failure modes or material properties may follow distinct distributions — tensile strength may fit a normal law, whereas compression or shear strength may follow Weibull or lognormal laws [15]. Strength data for high-performance fibres and polymers provide further evidence that, even when statistical tests “cannot reject normality”, the true strength histograms deviate from the ideal bell-shaped Gaussian curve [16]. In advanced composites, manufacturing defects can produce fat-tailed or skewed resistance distributions rather than the assumed thin-tailed form [17]. Similarly, accelerated ageing or damage can alter an initially normal resistance distribution into a Weibull-type distribution with heavier tails [18]. In such cases, structural strength variables exhibit skewness, heavy tails, or even multimodal behaviour, departing significantly from the classical normal or lognormal assumptions [17,18].

There are measurement cases in which the data distribution is only approximately Gaussian. However, when slightly larger deviations characterised by higher skewness or kurtosis occur, alternative distribution models should be considered [19,20]. Such cases have been observed both in experimental data [21,22] and in model outputs in engineering [23,24]. In both cases, the histogram approaches a bell curve, and the higher moments (skewness and kurtosis) approach zero or small values typical of a normal distribution. For such models, a suitable option is to employ a four-parameter Hermite probability density function (pdf) [25].

This article analyses the limit states of load-bearing structures using a Hermite model [25], which enables calibration to measured data exhibiting non-negligible skewness and kurtosis. The focus is placed on situations where the shape of the probability distribution materially influences the lower design quantiles of resistance and, consequently, the required level of reliability. Sobol sensitivity analysis is employed to quantify the contribution of each input variable and their interactions to the variance of the model output. In addition, it is used to evaluate how different choices of input probability distributions affect overall model sensitivity.

2. Hermite Models for Tail-Sensitive Reliability

Many studies report that the probability distributions of loads and resistances deviate from the Gaussian law and can exhibit non-zero skewness and kurtosis [26,27,28,29]. These higher-order moments can significantly influence failure probabilities and the reliability index [30,31]. Consequently, computed results depend on the assumed input distributions, as demonstrated for wind-loaded truss towers [32]. Therefore, more advanced probabilistic models are required.

The Hermite translation method serves as a tail-aware default when departures from normality are of interest. Its third-order mapping transforms a Gaussian variable into a non-Gaussian one, with shape parameters controlling skewness and kurtosis [25,33,34]. Evidence from wind engineering shows that the Hermite probability density function (pdf) improves tail representation for both mildly and strongly non-Gaussian pressure data on low-rise buildings [35]. The Hermite pdf has also been employed to obtain moment-based peak factors for pressure and displacement responses [36,37]. When full probability distributions are unavailable, reliability assessment and partial safety factor calibration can utilise moment methods based on the third to fifth moments [30,31,38]. Correlated inputs can be addressed using an L-moment-based normal transformation [39]. These approaches do not require full distributions and more faithfully capture non-Gaussian behaviour. As a general alternative, a unimodal four-moment density can represent a wide range of skewness and kurtosis and connects to both normal and lognormal limits [40].

Extensions of Hermite models target small to moderate deviations from normality. An adaptive Hermite model using probability-weighted moments estimates extreme value distributions for seismic reliability, while asymptotic analysis clarifies the monotonic domain and establishes connections to Gram–Charlier behaviour [41,42]. L-moment-based Hermite mappings enable efficient simulation of non-Gaussian random fields in structural [43] and geotechnical applications [44]. When higher-order moments are required but full distributions are not available, skewness and kurtosis can be estimated through a regression-based adaptation of the Cornish–Fisher expansion to improve shape recovery for asymmetric or heavy-tailed data [45]. In parallel, a flexible four-moment distribution employing a cubic normal transformation covers a broad region of the skewness–kurtosis plane and connects to the normal and lognormal limits [46]. For non-stationary, non-Gaussian responses, a unified Hermite mapping provides analytical mean upcrossing rates, and complementary studies quantify the transmission of kurtosis from excitation to response under both stationary and non-stationary loadings [47,48,49,50]. Hermite polynomial models have also been integrated into dynamic reliability formulations under non-Gaussian excitation [51].

Fatigue assessment is particularly sensitive to non-Gaussian features, as damage accumulation is driven by high-amplitude cycles in the distribution tails. Skewness and kurtosis can be embedded into damage estimators through Hermite-based mappings when full distributions are unavailable. For instance, Hermite transformations have been applied to fatigue damage estimation in broadband non-Gaussian processes using a conditional probability approach [52]. Time–frequency frameworks for non-stationary loading employ spectral methods for non-Gaussian processes and support fatigue evaluation [53]. Fast modal-domain estimators of response kurtosis enable efficient assessment of non-Gaussian effects in full finite element models used for random vibration fatigue [54]. Under simultaneous non-Gaussian processes, bias correction incorporating kurtosis improves fatigue damage estimates compared with linear summation [55]. For floating offshore wind turbine blades, non-Gaussian wind fields have been shown to reduce predicted fatigue life relative to Gaussian assumptions [56]. Fatigue in steel plate welded joints has also been evaluated using both non-Gaussian parametric and non-parametric formulations [57].

Direct applications of the Hermite pdf, or its explicit use for uncertainty modelling, span theoretical developments [33,34,41,42,47], simulation and random fields [43,44,58], wind loading and peak factor estimation [35,36,37,59], and reliability and fatigue analysis [25,51,52,56]. Several studies have employed quartic-order Hermite polynomial mappings to transform Gaussian variables into prescribed non-Gaussian distributions [31,44,58,59]. Although the quartic Hermite polynomial can improve the shape of the pdf, its monotonic domain is assessed locally and may be narrower over the practical working range. In contrast, the cubic model provides a transparent and easily verifiable monotonic boundary [44,59]. This article employs a cubic Hermite pdf model, the derivation of which is presented in the following chapter. The cubic model is suitable for formulating design quantiles and can surpass the accuracy of design values obtained using the classical FORM method.

3. Hermite Model: Map, Densities, and Properties

The Hermite translation model is introduced to represent moderate deviations from Gaussian behaviour while maintaining a compact four-parameter description defined by the mean, standard deviation, skewness, and kurtosis [25]. Throughout this chapter, a real-valued random variable X with finite moments is considered. The mean and variance are denoted by μ and σ². The central and standardised moments follow the conventions stated below.

$$\mu =E\left(X\right), {\sigma }^2=V\left(X\right), {\mu }_r^c=E\left({\left(X-\mu \right)}^r\right), r\ge 1.$$

(1)

$${\alpha }_3={\displaystyle \frac{{\mu }_3^c}{{\sigma }^3}}, {\alpha }_4={\displaystyle \frac{{\mu }_4^c}{{\sigma }^4}}.$$

(2)

For a Gaussian variable, one has α₄=3, departures from Gaussian kurtosis are therefore captured by α₄−3. Under these conventions, the standardisation in Equation (3) maps the observations to zero mean and unit variance, while the Hermite translation in Equation (4) introduces controlled adjustments to symmetry and tail weight through the calibrated coefficients defined in Equations (5)–(6) and Equation (14).

3.1. Polynomial Map and Parameterisation

A random variable X with mean μ, standard deviation σ, skewness α₃, and excess kurtosis α₄ ​is approximated by transforming a standard normal variable U∼N(0,1).

The standardised observation is

$$y={\displaystyle \frac{x-\mu }{\sigma }},$$

(3)

which has zero mean and unit variance. The Hermite translation model introduces a unit-variance surrogate defined as

$$Y=k\left[U+h_3\left(U^2-1\right)+h_4\left(U^3-3U\right)\right],$$

(4)

where H₂​(U)=U²−1 and H₃(U)=U³−3U are the second- and third-order Hermite polynomials, respectively. The third- and fourth-order base coefficients h₃ and h₄ are given by

$$h_3={\displaystyle \frac{{\alpha }_3}{6}}, h_4={\displaystyle \frac{{\alpha }_4-3}{24}},$$

(5)

where the constant k is a variance rescaling factor chosen such that V(Y)=1. Because the Hermite polynomials are orthogonal under the Gaussian inner product, E[U H_n(U)]=0, the centring is preserved, i.e. E(Y)=0. Equations (3)–(5) define the common basis of the model. The branch-specific calibrations ${\tilde{h}}_3$, ${\tilde{h}}_4$ and correction factors k_t are introduced below.

3.2. Density Derivations

3.2.1. Lower-Branch Density (α₄≤3)

For α₄≤3, the Winterstein calibration is applied:

$${\tilde{h}}_4=h_4-27h_4^2, {\tilde{h}}_3={\displaystyle \frac{h_3}{1+24{\tilde{h}}_4}}, k_t={\displaystyle \frac{1}{\sqrt{1+10{\tilde{h}}_3^2+42{\tilde{h}}_4^2}}}.$$

(6)

Dividing Equation (4) by k_t gives the working, dimensionless mapping:

$$y={\displaystyle \frac{Y}{k_t}}=U+{\tilde{h}}_3\left(U^2-1\right)+{\tilde{h}}_4\left(U^3-3U\right),$$

(7)

The density of Y follows from the change-of-variables formula:

$$f_Y\left(y\right)=\varphi \left(U\left(y\right)\right)\cdot \left|{\displaystyle \frac{dU}{dy}}\right|, \phi (u)={\frac{1}{\sqrt{2\pi }}}^{e\frac{-u^2}{2},}$$

(8)

and is obtained perturbatively for $\left|{\tilde{h}}_3\right|$, $\left|{\tilde{h}}_4\right|$≪ 1.

Applying a first-order Neumann inversion to Equation (7), and using the scalar rule (1+a)^-1 ≈ 1-a for ∣a∣≪1, the map U↦U+ε(U) can be inverted to first order as U≈y−ε(y). This results in:

$$U\left(y\right)\approx y-{\tilde{h}}_3\left(y^2-1\right)-{\tilde{h}}_4\left(y^3-3y\right).$$

(9)

Differentiating Equation (7) gives:

$${\displaystyle \frac{dy}{dU}}=1+2{\tilde{h}}_{3}U+3{\tilde{h}}_4\left(U^2-1\right) \Rightarrow \left|{\displaystyle \frac{dU}{dy}}\right|\approx 1-2{\tilde{h}}_{3}y-3{\tilde{h}}_4\left(y^2-1\right),$$

(10)

where U is replaced by y to first order. Substituting Equations (9)–(10) into Equation (8) then gives:

$$f_Y\left(y\right)\approx \varphi \left(U\left(y\right)\right)\cdot \left[1-2{\tilde{h}}_{3}y-3{\tilde{h}}_4\left(y^2-1\right)\right].$$

(11)

For a non-standardised X, it is convenient (and consistent with the calibration) to evaluate the model at:

$$y={\displaystyle \frac{x-\mu }{\sigma k_t}},$$

(12)

whence, by the one-dimensional change-of-variables rule,

$$f_X\left(x\right)={\displaystyle \frac{1}{\sigma k_t}}{f}_Y\left({\displaystyle \frac{x-\mu }{\sigma k_t}}\right).$$

(13)

Equations (9)–(13) define the lower-branch pdf used in practical applications. The approximation remains accurate while ∣α₃∣≪1 and the kurtosis surplus α₄−3 is non-positive or only mildly positive.

3.2.2. Upper-Branch Density (α₄>3)

For α₄>3 the calibrated moments are defined as:

$${\tilde{h}}_4={\displaystyle \frac{\sqrt{1+36h_4}-1}{18}}, {\tilde{h}}_3={\displaystyle \frac{h_3}{1+6{\tilde{h}}_4}}, k_t={\displaystyle \frac{1}{\sqrt{1+2{\tilde{h}}_3^2+6{\tilde{h}}_4^2}}}.$$

(14)

The working mapping remains as in Equation (7), with ${\tilde{h}}_3$,${\tilde{h}}_4$​and k_t obtained from Equation (14). In this branch, the inverse mapping U(y) can be obtained in the closed form. Let

$$B={\displaystyle \frac{1}{3{\tilde{h}}_4}}, A={\tilde{h}}_{3}B, Q=B-1-A^2, \xi ={\displaystyle \frac{3}{2}}B\left(A+y\right)-A^3.$$

(15)

and introduce the shift U=X−A. Substituting into Equation (7) yields the depressed cubic form:

$$X^3+3QX-2\xi =0.$$

(16)

where the shift U=X−A removes the quadratic term, yielding the depressed cubic with coefficients Q and ξ defined in Equation (15).

Define

$$R=\sqrt{{\xi }^2+Q^3},(\ge |\xi \left|\right).$$

(17)

By applying Cardano’s formula and selecting the principal real cube roots, the real solution for U(y) is obtained as:

$$U\left(y\right)=\sqrt[3]{R+\xi }-\sqrt[3]{R-\xi }-A.$$

(18)

This expression for U(y) satisfies Equation (16) (with X=U+A) and, consequently, the original mapping in Equation (7).

The Jacobian dU/dy is obtained by differentiating Equation (18). From Equations (15) and (17),

$${\displaystyle \frac{d\xi }{dy}}={\displaystyle \frac{3}{2}}B, {\displaystyle \frac{dR}{dy}}={\displaystyle \frac{\xi }{R}}{\displaystyle \frac{3}{2}}B.$$

(19)

and differentiating Equation (18) followed by the substitution of Equation (19) yields

$$\left|{\displaystyle \frac{dU}{dy}}\right|={\displaystyle \frac{B}{2}}\left[{\displaystyle \frac{\xi /R+1}{{\left(R+\xi \right)}^{2/3}}}+{\displaystyle \frac{1-\xi /R}{{\left(R-\xi \right)}^{2/3}}}\right].$$

(20)

Hence, from the change of variables relation in Equation (8), the probability density function is given by

$$f_Y\left(y\right)=\varphi \left(U\left(y\right)\right)\cdot {\displaystyle \frac{B}{2}}\left[{\displaystyle \frac{\xi /R+1}{{\left(R+\xi \right)}^{2/3}}}+{\displaystyle \frac{1-\xi /R}{{\left(R-\xi \right)}^{2/3}}}\right],$$

(21)

where U(y) is given by Equation (18) and A, B, Q, ξ, R are defined in Equations (15)–(17). For a non-standardised variable X,

$$y={\displaystyle \frac{x-\mu }{\sigma k_t}}, f_X(x)=\frac{1}{\sigma k_t}{f}_Y\left(\frac{x-\mu }{\sigma k_t}\right).$$

(22)

Equations (18)–(22) are exact; no terms are neglected. The use of principal real cube roots in Equation (18) is justified by the condition R±ξ≥0, which ensures continuity of U(y) across ξ=0.

3.3. Properties and Admissible Sets

This section summarises the conditions under which the Hermite pdf, derived in the previous sections, behaves well numerically and accurately reproduces the prescribed moments. The guidance is expressed in the native parameters α₃ and α₄ through the calibrated quantities ${\tilde{h}}_3$,${\tilde{h}}_4$​and k_t defined above.

3.3.1. Non-Negativity and Monotonicity

The lower-branch density in Equation (13) is a first-order change-of-variables approximation; neglected terms are of order O(${\tilde{h}}^2$). Moment matching is therefore only first-order accurate while $\left|{\tilde{h}}_3\right|$,$\left|{\tilde{h}}_4\right|$≪ 1. Using the calibrations of Equation (6), large values of ∣α₃∣ or ∣α₄-3∣ push the approximation outside its linear regime and bias the empirical moments by terms of order O(${\tilde{h}}^2$).

Positivity of the first-order form in Equation (13) holds as long as the bracket term remains non-negative over the range of interest. The bracket in Equation (11) is

$$g\left(y\right)=1-2{\tilde{h}}_{3}y-3{\tilde{h}}_4\left(y^2-1\right)=\left(1+3{\tilde{h}}_4\right)-2{\tilde{h}}_{3}y-3{\tilde{h}}_4{y}^2.$$

(23)

For platykurtic cases ${\tilde{h}}_4<0$, the corresponding parabola opens upward, and its minimum occurs at

$$y_{\min }=-{\displaystyle \frac{{\tilde{h}}_3}{3{\tilde{h}}_4}}, g_{\min }=1+3{\tilde{h}}_4+{\displaystyle \frac{{\tilde{h}}_3^2}{3{\tilde{h}}_4}}.$$

(24)

A sufficient condition for ∣dU/dy∣ ∝ g(y)>0 over the effective range is therefore

$${\tilde{h}}_3^2<-3{\tilde{h}}_4\left(1+3{\tilde{h}}_4\right).$$

(25)

When $\left|{\tilde{h}}_4\right|$≪ 1, Equation (25) reduces to a transparent rule of thumb,

$${\alpha }_3^2\lesssim 4.5(3-{\alpha }_4),$$

(26)

obtained by substituting ${\tilde{h}}_3$≃h₃=α₃/6 and ${\tilde{h}}_4$≃h₄=(α₄−3)/24.

Violation of Equations (25)–(26) identifies regions where g(y)<0. In exact theory, this would create negative lobes of f_Y in the far tails. In numerical implementations, such tails are often clipped to zero, which reduces the total mass ∫f_Y and shifts the empirical moments.

In practice, for a given pair α₃, α₄, the recommended procedure is to compute ${\tilde{h}}_3$ and ${\tilde{h}}_4$ and test Equation (25). If it is marginally violated, g(y)≥0 should also be verified over the effective range of ∣y∣. Failure of either test indicates that the first-order model may require post-normalisation (which perturbs moments) or the adoption of a different pdf.

The upper-branch density defined by Equations (21)–(22) is exact, employing the Cardano real root of the depressed cubic (16). The Jacobian factor (20) follows by direct differentiation, with no truncation error. Two numerical aspects merit attention.

First, monotonicity of the mapping U↦Y is guaranteed while the cubic has a single real branch. In terms of α₃ and α₄, this corresponds to the classic engineering criterion

$${\alpha }_3^2\lesssim 4.5({\alpha }_4-3),$$

(27)

beyond which the cubic becomes multi-valued and the Cardano transformation no longer defines a valid pdf.

Second, although Equation (20) remains finite for all admissible parameters (since R±ξ≥0), near R≈ξ≥0 the individual factors (R±ξ)^2/3 can become large. The balanced product in Equation (20),

$${\displaystyle \frac{B}{2}}\left[{\displaystyle \frac{\xi /R+1}{{\left(R+\xi \right)}^{2/3}}}+{\displaystyle \frac{1-\xi /R}{{\left(R-\xi \right)}^{2/3}}}\right],$$

(28)

avoids catastrophic cancellation and ensures numerical stability on implementation.

3.3.2. Realised Moments

This section presents the closed-form expressions for the realised first four moments of the Hermite pdf and explains how inputs may be specified either through the native parameters (μ, σ, α₃, α₄) or, on the upper branch, directly via the target moments (μ^*, σ^*,${\alpha }_3^{*}$,${\alpha }_4^{*}$).

The standardised surrogate variable Y defined in Equation (4) is employed together with the branch–specific calibrations of ${\tilde{h}}_3$,${\tilde{h}}_4$​and the scaling factor k_t given in Equations (6) (lower branch) and (14) (upper branch). The location–scale moments of a non-standardised quantity X=μ+σY then follow directly from those of Y.

On the upper branch (α₄>3), the variance rescaling defined in Equation (14) ensures exact centring and unit variance due to the orthogonality of the Hermite basis,

$$E\left(Y\right)=0, V\left(Y\right)=1.$$

(29)

Consequently,

$${\mu }^{*}=\mu , {\sigma }^{*}=\sigma .$$

(30)

Because Y is expressed as a polynomial in the Gaussian variable U, the realised skewness and kurtosis are obtained in closed form as functions of (${\tilde{h}}_3$,${\tilde{h}}_4$). These expressions are exact for the upper branch. Denoting ${\tilde{h}}_3$ and ${\tilde{h}}_4$ as defined in Equations (6) and (14), the realised skewness and kurtosis are given by

$${\alpha }_3^{*}={\displaystyle \frac{2{\tilde{h}}_3\left(4{\tilde{h}}_3^2+54{\tilde{h}}_4^2+18{\tilde{h}}_4+3\right)}{{\left(1+2{\tilde{h}}_3^2+6{\tilde{h}}_4^2\right)}^{3/2}}},$$

(31)

$${\alpha }_4^{*}={\displaystyle \frac{60{\tilde{h}}_3^4+2232{\tilde{h}}_3^2{\tilde{h}}_4^2+576{\tilde{h}}_3^2{\tilde{h}}_4+60{\tilde{h}}_3^2+3348{\tilde{h}}_4^4+1296{\tilde{h}}_4^3+252{\tilde{h}}_4^2+24{\tilde{h}}_4+3}{{\left(1+2{\tilde{h}}_3^2+6{\tilde{h}}_4^2\right)}^2}}.$$

(32)

The density on the upper branch is obtained by the exact Cardano inversion using the balanced Jacobian from Equation (28). When the target shape moments (${\alpha }_3^{*}$,${\alpha }_4^{*}$) are prescribed, the coefficients (${\tilde{h}}_3$,${\tilde{h}}_4$) can be recovered through a stable inversion of Equations (31) and (32). The forward relations given in Equation (14) provide an initial estimate that is first-order accurate in (α₃, α₄−3). Admissibility of the upper branch follows from the monotonicity criterion in Equation (27). It can be noted that similar explicit relations between higher-order moments and distribution parameters have recently been derived for the Weibull pdf and used to compute design quantiles in hydrological design [60].

On the lower branch, the Winterstein calibration in Equation (6), combined with the first-order Neumann inversion and the change of variables in Equation (13), yields an approximate density. In this regime, only first-order preservation of the first two moments is achieved,

$$E\left(Y\right)\approx 0, V\left(Y\right)\approx 1.$$

(33)

Consequently,

$${\mu }^{*}=\mu +O\left({\epsilon }^2\right), {\left({\sigma }^{*}\right)}^2={\sigma }^2\left[1+O\left({\epsilon }^2\right)\right], \epsilon =\max \left\{ \left|{\tilde{h}}_3\right|, \left|{\tilde{h}}_4\right|\right\}.$$

(34)

Using Equations (31) and (32) with the lower-branch values of (${\tilde{h}}_3$,${\tilde{h}}_4$) yields (${\alpha }_3^{*}$,${\alpha }_4^{*}$) to first-order; higher-order deviations increase with ∣α3∣ and ∣α4−3∣. Practical admissibility and non-negativity on the lower branch are governed by Equations (23)–(26), with the serviceability rule given by Equation (26).

Let θ=(μ, σ, α₃, α₄) denote the input parameters used to evaluate the Hermite pdf through Equation (4) and the branch-specific calibrations in Equation (6) or (14). Denote by

$$M\left(\theta \right)=\left({\mu }^{*}, {\sigma }^{*}, {\alpha }_3^{*}, {\alpha }_4^{*}\right),$$

(35)

the realised moments computed from the resulting density f_X(⋅∣θ). For the lower branch M(θ)−θ is of order O(${\tilde{h}}_3^2$,${\tilde{h}}_4^2$,${\tilde{h}}_3$${\tilde{h}}_4$)​. For the upper branch, the inversion is exact and the relations in Equations (31) and (32) express the realised skewness and kurtosis precisely as functions of (${\tilde{h}}_3$,${\tilde{h}}_4$).

Two short illustrations follow.

Lower branch: For the inputs θ=(0, 1, 0.05, 2.95), Equations (31) and (32) yield ${\alpha }_3^{*}$≈0.049996 and ${\alpha }_4^{*}$≈2.95334. Numerical integration of the pdf over the interval [−19, 19] yields ${\alpha }_3^{*}$≈0.0503 and ${\alpha }_4^{*}$≈2.9563, reflecting higher-order effects inherent to the first-order formulation.

Upper branch: For the inputs θ=(0, 1, −0.1, 3.1), numerical integration agrees with Equations (31) and (32) for the prescribed pair (${\tilde{h}}_3$,${\tilde{h}}_4$). When the forward relations of Equation (14) are used, the realised kurtosis ${\alpha }_4^{*}$≈3.113 slightly differs from the target due to the first-order character of that calibration. Solving Equations (31) and (32) for (${\tilde{h}}_3$,${\tilde{h}}_4$) enforces ${\alpha }_3^{*}$=-0.1 and ${\alpha }_4^{*}$=3.1 within numerical precision.

For reporting purposes, the moments obtained from Equations (31) and (32) should be regarded as the reference values, while those computed by numerical integration of the lower-branch pdf may serve as a first-order consistency check.

3.4. Implementation and Workflow

Equations (25)–(28) provide analytical criteria that define the admissible regions in α₃ and α₄, ensuring both positivity and monotonicity of the Hermite pdf. For practical applications, these conditions must be embedded within an implementable computational routine. The Pascal/Delphi function presented below encodes both branches of the density, using the Winterstein calibration (Equation (6)) for α₄≤3 and the exact Cardano inversion (Equations (18)–(22)) for α₄>3. The Jacobian factor is implemented in the balanced form given by Equation (28), which prevents numerical cancellation errors in the vicinity of R≈ξ.

The structure of the routine reflects the admissibility conditions. When α₄≤3, the Neumann inversion is truncated at first order, consistent with the approximation in Equation (13). When α₄>3, the exact cubic root is employed, and validity requires compliance with the monotonicity criterion of Equation (27). In both branches, the Gaussian kernel ϕ(U) is evaluated once, and rescaling by k_t ensures unit variance. A defensive clamp at the end of the routine enforces non-negativity of the pdf and prevents spurious negative values under extreme input conditions.

The Pascal (Delphi 6) implementation of the Hermite pdf is provided in Appendix A.

This implementation enables the direct numerical evaluation of the Hermite pdf without further algebraic manipulation. In practical applications, the computational workflow proceeds as follows:

• Compute ${\tilde{h}}_3$,${\tilde{h}}_4$ and k_t from α₃ and α₄.
• Select the branch: apply the Winterstein correction for α₄≤3, or the Cardano inversion for α₄>3.
• Verify the admissibility conditions, Equation (25) for the lower branch or Equation (27) for the upper branch.
• Evaluate the density using the function above.

If the admissibility conditions are violated, the Hermite pdf cannot be used directly. In such cases, an alternative parametric or non-parametric model, which is better suited to the empirical histogram, must be adopted.

4. Fitting and Validation on Materials Data

Recent experimental studies on reinforcing steel bars, welded joints and cementitious composites report non-Gaussian measurement profiles characterised by non-negligible skewness and excess kurtosis [61,62]. As a result, strength and geometric characteristics are often modelled using Weibull, lognormal or gamma distributions [63,64,65,66]. In fibre-reinforced polymer rods and in fibre-reinforced concrete exposed to high temperatures, Weibull fits effectively capture right-skewed tails and summarise the design-relevant variability [63,64]. For metakaolin masonry blocks and large regional datasets of green concrete, lognormal or gamma laws frequently provide the best or at least competitive fits [65,66]. In rolled steels and butt-welded seams, empirical skewness and kurtosis often deviate substantially from the normal benchmark, with several cases exhibiting kurtosis greater than three, motivating the explicit use of non-Gaussian modelling [26,61,62].

For many experimental datasets, the Hermite pdf can serve as an extension of the Gaussian model that accommodates small but non-zero skewness and variable excess kurtosis. For example, in weld toe geometry and reinforcement height distributions that remain close to symmetric yet exhibit mild excess kurtosis, a moment-matched Hermite fit can be applied instead of switching to lognormal or gamma families [61]. For reinforcing steel, where the sign of skewness varies across production lots, both positive and negative skewness can be represented within a single parametric family while preserving the observed variance and kurtosis [26]. For regional concrete datasets that are approximately normal, small departures from symmetry and tail weight can be captured using a parsimonious Hermite approximation that retains the Gaussian baseline for interpretability [66].

In the following section, Hermite pdfs are fitted to our experimentally obtained material characteristics of structural steel [67,68], which exhibit non-Gaussian skewness and kurtosis. The fits are validated using histograms and moment analysis.

4.1. S235 Yield Strength (Lower Branch)

In engineering practice, histograms of measured quantities may deviate, to varying degrees, from the Gaussian pdf. An illustrative example is provided by the yield strength f_y of structural steel grade S235 [67]. Based on n=562 valid observations of f_y, the sample statistics were estimated as mean of 297.3 MPa, standard deviation of 16.8 MPa, skewness of 0.325, and kurtosis of 2.5415 [67].

For a platykurtic input, the lower-branch relations given by Equations (6) to (13) apply. The rule in Equation (26) gives ${\alpha }_3^2$≪4.5(3−α₄), and the non-negativity check based on Equations (23)–(25) is satisfied across the effective range; hence, admissibility is confirmed.

When the pdf is evaluated using the input θ=(297.3 MPa, 16.8 MPa, 0.325, 2.541), the realised moments of the resulting Hermite pdf are M(θ)=(297.429 MPa, 16.207 MPa, 0.904, 3.821), where numerical integration was performed over the range [μ−19σ, μ+19σ]. The Kolmogorov–Smirnov (KS) test indicates a lack of fit with p=0.00148<0.05; see the red curve in Figure 2. This behaviour is consistent with the first-order character of Equations (9)–(13).

After a minor parameter refinement within the admissible region, θ=(297.5 MPa, 17.6 MPa, −0.071, 2.349) produces realised moments M(θ)=(297.337 MPa , 16.818 MPa, 0.328, 2.544), which closely approximate the sample moments. The KS test then reports p=0.2965>0.05; see the green curve in Figure 2. Both parametrisations satisfy Equation (26) and the positivity condition in Equation (25). The improved fit results from compensating for the first-order truncation error inherent in Equations (9)–(13).

This structure ensures that the first four moments are taken into account while maintaining conceptual transparency. If the admissibility conditions are violated, the recommended fallback is to revert to the Gaussian pdf or to adopt an alternative model.

4.2. S235 Ductility (Upper Branch)

Ductility measurements for S235 steel, based on n=562 observations, are analysed. The sample moments are μ=37.8%, σ=2.9%, α₃=−0.753 and α₄=5.18. Since α₄>3, the upper-branch formulation is applicable: calibration follows Equation (14), the exact cubic inversion is performed using Equations (18)–(22), and the Jacobian is evaluated in the balanced form given by Equation (28). The monotonicity requirement specified in Equation (27) is satisfied.

Using the sample moments directly as inputs θ=(37.8, 2.9, −0.753, 5.18) (red curve in Figure 3) yields the realised moments M(θ)=(37.8, 2.9, −0.751, 6.026) when evaluated either by numerical integration over the range [μ−19σ,μ+19σ] or by direct substitution into Equations (31) and (32). The discrepancy (${\alpha }_3^{*}$,${\alpha }_4^{*}$)≠(α₃,α₄) reflects the first-order nature of the polynomial translation in Equation (4). The Kolmogorov–Smirnov test gives p=0.0273<0.05; hence, the model is rejected.

A short refinement of the input parameters improves agreement between he realised and sample moments. With θ=(37.8, 2.9, −0.7605, 4.422) (green curve in Figure 3), the realised moments become M(θ)=(37.8, 2.9, −0.753, 5.180), which are very close to the sample moments. However, the Kolmogorov–Smirnov test still rejects the model at the 5% significance level. It can be noted that increasing the input kurtosis parameter to α₄=7.0 while keeping μ≈37.8, σ≈2.9 and α₃≈−0.762 yields p=0.121>0.05, i.e., non-rejection, but at the expense of a substantial deviation from the sample kurtosis.

A practical trade-off is implied: enforcing exact moment matching does not necessarily maximise agreement under the Kolmogorov–Smirnov criterion. Notably, in all configurations considered, the Hermite fits were not rejected by the Anderson–Darling goodness-of-fit test, indicating adequate tail behaviour despite the KS test’s sensitivity to central discrepancies for large n.

4.3. Extreme-Moment Inputs

This subsection examines inputs θ=(μ, σ, α₃, α₄) that lie outside the perturbative regime assumed in the derivation. In such cases, θ must be regarded as control parameters of the transformation rather than as guaranteed realised moments, and the resulting pdf f_X(⋅∣θ) must be validated a posteriori. Validation includes checks for non-negativity and normalisation, admissibility via Equation (25) for the lower branch and Equation (27) for the upper branch, and consistency of the realised moments M(θ)=(μ^*, σ^*, ${\alpha }_3^{*}$, ${\alpha }_4^{*}$) with the intended targets. Where discrepancies occur, a constrained numerical calibration of θ is required, followed by re-evaluation using Equations (13) for the lower branch or Equations (21)–(22) for the upper branch, with the balanced Jacobian defined in Equation (28).

For strongly platykurtic inputs (α₄<3), the first-order representation of the lower branch, f_Y(y)≈ϕ(U(y)) g(y) with $g\left(y\right)=1-2{\tilde{h}}_{3}y-3{\tilde{h}}_4\left(y^2-1\right)$ defined by Equations (11) and (23), governs admissibility and the shape of the pdf. When ${\tilde{h}}_4$<0, the factor g(y) forms a convex parabola whose minimum occurs at y_min and reaches g_min as given in Equation (24). Positivity requires g_min>0, corresponding to the bound in Equation (25). This bound explains two features observed in Figure 4. Firstly, very small target kurtosis depresses the centre of g(y) and may create a double-humped pdf even when the mapping remains monotonic. Secondly, large values of ∣α3∣ shift the minimum of g(y) and tighten the constraint imposed by Equation (25). Consequently, extremely low kurtosis cannot be combined with large skewness without violating non-negativity. In the mild-departure limit, the transparent rule of thumb expressed by Equation (26) is recovered.

For leptokurtic inputs (α₄>3), the upper branch employs the exact cubic inversion defined by Equations (18)–(22), and monotonicity reduces to the single inequality expressed in Equation (27). Because the right-hand side of Equation (27) increases with α₄−3, larger magnitudes of ∣α₃∣ become admissible as kurtosis increases. This behaviour is consistent with the broader family of attainable pdf shapes illustrated in Figure 5. The Jacobian remains finite when evaluated using the balanced product form in Equation (28). Nevertheless, even on this branch, the polynomial translation given in Equation (4) represents a cubic truncation. Consequently, the realised higher-order moments deviate from the input values once (α₃, α₄) move sufficiently far from (0, 3).

In practical modelling, the extreme-input regime can generate useful pdf shapes, including light-tailed or bimodal forms. The reliability of such fits depends on enforcing the admissibility conditions given by Equations (25)–(27), verifying the normalisation and non-negativity of f_X, and, where necessary, numerically calibrating θ so that M(θ) satisfies the desired moment constraints while the pdf remains admissible.

5. Reliability Application

In EN 1990, design values are obtained from the mean and standard deviation using the FORM [1]. When full probability distributions are specified, design quantiles are typically derived from Gaussian, lognormal or Gumbel laws [1]. A review of the literature reveals no published calculation of the 0.12% quantile directly under a Hermite probability density function. Existing applications of the Hermite model have focused primarily on peak factor estimation [37], dynamic reliability [51], and fatigue assessment [53], rather than on code-level design quantile evaluation. Motivated by these gaps, the following section derives closed-form design quantiles under a Hermite density, allowing skewness and kurtosis to be incorporated consistently into the estimation of design quantiles.

5.1. Design Quantiles

5.1.1. EN 1990 Quantile

Let the structural resistance R and the load effect F be modelled as independent random variables with means μ_R, μ_F and standard deviations σ_R, σ_F, see, e.g., [69]. The performance function is defined as

$$Z=R-F.$$

(36)

Under the Gaussian baseline, the reliability index

$$\beta ={\displaystyle \frac{{\mu }_{\mathrm{Z}}}{{\sigma }_{\mathrm{Z}}}}, {\mu }_{\mathrm{Z}}={\mu }_{\mathrm{R}}-{\mu }_{\mathrm{F}}, {\sigma }_{\mathrm{Z}}=\sqrt{{\sigma }_{\mathrm{R}}^2+{\sigma }_{\mathrm{F}}^2},$$

(37)

leads to the failure probability P_f=Φ(−β). Eurocode practice compares β with a target reliability index β_d (for the ultimate limit state: β_d≈3.8, corresponding to P_f≈7.2⋅10⁻⁵).

Within FORM, the contribution of R and F to σ_Z is described by sensitivity factors α_R, α_F>0 (for independent normal variables α_R=σ_R/σ_Z, α_F=σ_F/σ_Z, some standards adopt fixed values, e.g., α_R=0.8, α_F=0.7 [1]). The design values are then expressed as

$$R_{d,FORM}={\mu }_{\mathrm{R}}-{\alpha }_R{\beta }_d{\sigma }_R,$$

(38)

$$F_{d,FORM}={\mu }_{\mathrm{F}}+{\alpha }_F{\beta }_d{\sigma }_F.$$

(39)

The reliability requirement is verified through the inequality of the design quantiles:

$$F_{d,FORM}\le R_{d,FORM}.$$

(40)

The design resistance R_d,FORM represents the lower Gaussian quantile at probability level

$$p_d=\Phi \left(-{\alpha }_R{\beta }_d\right),$$

(41)

where, for the ultimate limit state, β_d=3.8 and α_R=0.8 are used for reliability class RC2 in [1]. This gives p_d≈Φ(-3.04)≈0.0012.

The design load action F_d,FORM corresponds to the higher Gaussian quantile at probability level

$$p_F=\Phi \left(+{\alpha }_F{\beta }_d\right),$$

(42)

where β_d=3.8 and α_F=0.7 are used for RC2 [1], giving p_d≈Φ(+2.66)≈0.996.

Symbols α₃, α₄ denote skewness and kurtosis used in the Hermite pdf, while α_R and α_F are FORM sensitivity coefficients (weight factors) and must not be confused.

5.1.2. Hermite Quantile (Closed Form)

Let U∼N(0,1) and define the Hermite translation as

$$T\left(U\right)=k_t\left[U+{\tilde{h}}_3\left(U^2-1\right)+{\tilde{h}}_4\left(U^3-3U\right)\right],$$

(43)

with branch-specific ${\tilde{h}}_3$,${\tilde{h}}_4$​ and k_t determined by (α_3,R, α_4,R) under the admissibility and monotonicity conditions stated earlier. The structural resistance is then

$$R={\mu }_R+{\sigma }_{R}T\left(U\right).$$

(44)

If T is strictly increasing (upper branch: exact, lower branch: valid in the small-departure regime), the cumulative density function (CDF) transforms as F_R(r)=Φ(U(r)). Hence, the p-quantile is obtained by inserting the Gaussian p-quantile z_p=Φ⁻¹(p) into the forward map:

$$q_p={\mu }_R+{\sigma }_R k_t\left[z_p+{\tilde{h}}_3\left(z_p^2-1\right)+{\tilde{h}}_4\left(z_p^3-3z_p\right)\right].,$$

(45)

where z_p=Φ⁻¹(p). For design purposes in the FORM sense, retain the target probability content

$$p_d=\Phi \left(-{\alpha }_R{\beta }_d\right), z_d=-{\alpha }_R{\beta }_d,$$

(46)

and define the Hermite-based design resistance as

$$R_{d,HERM}={\mu }_R+{\sigma }_R k_t\left[z_d+{\tilde{h}}_3\left(z_d^2-1\right)+{\tilde{h}}_4\left(z_d^3-3z_d\right)\right].$$

(47)

This expression reduces to the Gaussian design value R_d,FORM=μ_R−α_Rβ_dσ_R in the limit α_3,R→0, α_4,R→3, where ${\tilde{h}}_3$,${\tilde{h}}_4$→0 and k_t→1. Examples of design quantiles for normalised distributions with μ_R=1 σ_R=1 are shown in Figure 6.

When the load effect F is modelled as non-Gaussian, the same construction applies using parameters (μ_F, σ_F, α_3,F, α_4,F) and coefficients ${\tilde{h}}_3^{(F)}$,${\tilde{h}}_4^{(F)}$,$k_t^{(F)}$​. The design quantile corresponds to p_F=Φ(+α_Fβ_d) with z_F=+α_Fβ_d​, giving

$$F_{d,HERM}={\mu }_F+{\sigma }_F k_t^{(F)}\left[z_F+{\tilde{h}}_3^{(F)}\left(z_F^2-1\right)+{\tilde{h}}_4^{(F)}\left(z_F^3-3z_F\right)\right],$$

(48)

which reduces to the Gaussian F_d,FORM as α_3,F→0 and α_4,F →3. In the present study, only R_d,HERM is used in the examples.

5.2. Case Study: Pin-Jointed Tie System

A pin-connected node is supported by two axially loaded ties and subjected to a vertical gravity force F>0, see Figure 7.

Small-displacement statics and perfectly pinned joints are assumed; self-weight and bending effects are neglected. Static equilibrium of the node gives

$$\sum F_x=F_1\cos {\alpha }_1+F_2\cos {\alpha }_2=0\Rightarrow F_1=F_2,$$

(49)

$$\sum F_y=F_1\sin {\alpha }_1+F_2\sin {\alpha }_2-F=0\Rightarrow F={\displaystyle \frac{2}{\sqrt{5}}}\left(F_1+F_2\right)={\displaystyle \frac{4}{\sqrt{5}}}{F}_1,$$

(50)

where $\sin {\alpha }_1=\sin {\alpha }_2=2/\sqrt{5}$ and $\cos {\alpha }_1=-\cos {\alpha }_2=1/\sqrt{5}$. The axial force in each tie under given vertical load is F₁=F₂=$\sqrt{5}/4$F by Equation (50).

Let f_y1 and f_y2 denote the yield strengths of the two steel ties, and A₁ and A₂ their respective cross-sectional areas. The axial load-carrying capacities are

$$N_{R1}=f_{y1}{A}_1, N_{R2}=f_{y2}{A}_2,$$

(51)

where f_y is expressed in MPa and A in mm², giving N_R in N. The vertical load action that brings tie i to its ultimate limit state follows from Equation (50):

$$F_{Ri}={\displaystyle \frac{4}{\sqrt{5}}}{N}_{Ri} (i=1, 2),$$

(52)

The structural (system) resistance is governed by the weaker tie reaching yield first:

$$R=\min \left\{F_{R1}, F_{R2}\right\}={\displaystyle \frac{4}{\sqrt{5}}}\min \left\{f_{y1}{A}_1, f_{y2}{A}_2\right\}={\displaystyle \frac{4\pi }{\sqrt{5}}}\min \left\{f_{y1}{r}_1^2, f_{y2}{r}_2^2\right\},$$

(53)

where $A_i=\pi r_i^2$(i = 1, 2) denotes the cross-sectional area of the solid circular bar with radius r_i. Equation (53) is nonlinear, because the system resistance is governed by the weaker tie and the axial capacity of each tie depends on geometry and yield strength.

It can be noted, that comparable nonlinear input–output relations occur in nonlinear shell analysis [70,71], in the nonlinear reliability assessment of steel frames with Gaussian imperfections [31], and in nonlinear fatigue analysis under Gaussian loading [72]. Studies [31,70,72] show that nonlinear structural models with Gaussian input variables often yield non-Gaussian outputs due to the inherent nonlinearity of the input–output relationship. Skewness, kurtosis and other higher-order features of model output may therefore arise from the model itself, from the input distributions, or from a combination of both.

5.2.1. Inputs Random Quantities

Equation (53) is a function of four random quantities, where input vector is

$$\mathit{X}={\left(f_{y1}, r_1, f_{y2}, r_2\right)}^T,$$

(54)

and the model output is the structural resistance Y=R(X) obtained from Equation (53). In Equation (54), the four input variables are assumed to be mutually independent, and both members are assumed to possess identical statistical distributions of material and geometric properties.

The yield strength f_yi is treated as a random variable described by the histogram shown in Figure 2. The green Hermite pdf is adopted, as it reproduces the first four statistical moments of the histogram with numerical precision and passes the Kolmogorov–Smirnov goodness-of-fit test. The mean, standard deviation, skewness, and kurtosis of this Hermite pdf are consistent with the observed sample statistics, providing a mathematically admissible representation of yield strength while preserving the probabilistic characteristics of the measured data, see Table 1.

The radius r_i is treated as a random variable for hot-rolled steel bars with nominal value 5 mm. When direct measurements of bar radii are unavailable, the diameter can be taken as μ_r = 5 mm and the standard deviation σ_r is estimated assuming that 95 percent of realisations lie within the tolerance limits specified in EN 10060:2003, which gives σ_r = 0.1021 mm. Information on skewness and kurtosis is not provided in Eurocode design or tolerance standards. A Gaussian pdf is therefore adopted for r_i as a simplified model in the absence of higher-order statistics.

Nevertheless, statistical investigations of geometric characteristics in other rolled profiles, such as IPE sections [67,68], have reported skewness values between −0.5448 and 1.0545 and kurtosis values between 2.65 and 7.473, indicating that actual distributions may deviate from the Gaussian assumption. All input random quantities are summarised in Table 1.

In Table 1, the skewness and kurtosis of r_i are treated as uncertain parameters and their influence on the output statistics of R is examined. From an engineering perspective, configurations with a smaller mean and a larger standard deviation of R are considered less safe. Under fixed μ_r = 5 mm and σ_r = 0.1021 mm, such behaviour is promoted when the distribution of r_i exhibits negative skewness and leptokurtosis. The subsequent analyses therefore focus on deviations from the Gaussian baseline toward α_3,r<0, and α_4,r>3.

5.2.2. Sobol Sensitivity Analysis of Resistance

Sobol sensitivity analysis is used to assess how variations in the skewness and kurtosis of the radius r influence the variance of the model output R. In this case study, the pairs (α_3,r, α_4,r) = (0, 3), (−1, 4), (−2, 5), (−3, 6) are considered, see Table 1.

The first-order and total Sobol indices for input X_i are defined as

$$S_i={\displaystyle \frac{V\left(E\left(Y\left|X_i\right.\right)\right)}{V\left(Y\right)}}, S_{Ti}=1-{\displaystyle \frac{V\left(E\left(Y\left|X_{~i}\right.\right)\right)}{V\left(Y\right)}},$$

(55)

where E(·) denotes the mean value, V(·) the variance and X_~i represents the full set of input variables excluding X_i.

Expectations and variance V(Y) were evaluated by tensor-product midpoint quadrature on a uniform grid [lb_j, ub_j] with N = 400 nodes per dimension. For j=1,…,4 let dx_j = (ub_j−lb_j)/N and x_j(i_j) = lb_j + (i_j + 0.5) for i_j = 0,…,N−1. For independent inputs, the joint density is given by f₁(x₁)f₂(x₂)f₃(x₃)f₄(x₄). Then

$$E\left(Y\right)\approx {\displaystyle \sum _{i_1=0}^{N-1}{\displaystyle \sum _{i_2=0}^{N-1}{\displaystyle \sum _{i_3=0}^{N-1}{\displaystyle \sum _{i_4=0}^{N-1}Y\left(x_1\left(i_1\right),x_2\left(i_2\right),x_3\left(i_3\right),x_4\left(i_4\right)\right){\displaystyle \prod _{j=1}^4{f}_j}\left(x_j\left(i_j\right)\right)d{x}_{1}d{x}_{2}d{x}_{3}d{x}_4}}}},$$

(56)

$$E\left(Y^2\right)\approx {\displaystyle \sum _{i_1=0}^{N-1}{\displaystyle \sum _{i_2=0}^{N-1}{\displaystyle \sum _{i_3=0}^{N-1}{\displaystyle \sum _{i_4=0}^{N-1}Y{\left(x_1\left(i_1\right),x_2\left(i_2\right),x_3\left(i_3\right),x_4\left(i_{41}\right)\right)}^2{\displaystyle \prod _{j=1}^4{f}_j}\left(x_j\left(i_j\right)\right)d{x}_{1}d{x}_{2}d{x}_{3}d{x}_4}}}},$$

(57)

$$V\left(Y\right)\approx E\left(Y^2\right)-E{\left(Y\right)}^2.$$

(58)

For the first-order index S₁, the conditional mean on the grid of X₁ was evaluated using triple midpoint sums over the remaining coordinates:

$$E\left(Y\left|X_1=x_1\left(i_1\right)\right.\right)\approx {\displaystyle \sum _{i_2=0}^{N-1}{\displaystyle \sum _{i_3=0}^{N-1}{\displaystyle \sum _{i_4=0}^{N-1}Y\left(x_1\left(i_1\right),x_2\left(i_2\right),x_3\left(i_3\right),x_4\left(i_4\right)\right){\displaystyle \prod _{j=2}^3{f}_j}\left(x_j\left(i_j\right)\right)d{x}_{2}d{x}_{3}d{x}_4}}},$$

(59)

and the variance of the conditional mean and corresponding sensitivity index were obtained as

$$V\left(E\left(Y\left|X_1\right.\right)\right)\approx {\displaystyle \sum _{i_1=0}^{N-1}{\left(E\left(Y\left|X_1=x_1\left(i_1\right)\right.\right)\right)}^2}{f}_1\left(x_1\left(i_1\right)\right)d{x}_1-E{\left(Y\right)}^2$$

(60)

where E(Y)=E(E(Y∣X₁)). Kahan compensated summation was applied to both unconditional and conditional sums to reduce round-off error. The interval [lb_j, ub_j] was defined by lb_j =μ_i -(10-α₃) σ_i and ub_j =μ_i+(10+α₃) σ_i ensuring adequate coverage of the left tails arising from negative skewness (α₃<0). The remaining first-order indices and total indices were estimated analogously.

Numerical results are presented for the radius-moment pairs (α_3,r, α_4,r) = (0, 3), (−1, 4), (−2, 5), (−3, 6), as listed in Table 1. In Figure 8, the baseline case is depicted, where the radius r follows a Gaussian pdf with (α_3,r, α_4,r)=(0, 3). The Gaussian pdf of r is shown in the left panel. Its effect on the variance decomposition is shown in the two right-hand panels. The pie chart presents the first-order Sobol indices S_i, with the grey sector labelled “higher order” corresponding to 1−∑S_i. The bar chart displays the corresponding total indices S_Ti.

Figure 9, Figure 10 and Figure 11 present the non-Gaussian cases. The radius r is modelled by a Hermite pdf with μ_r = 5 mm and σ_r = 0.1021 mm, using the moment pairs (α_3,r, α_4,r) = (−1, 4), (−2, 5), (−3, 6), as listed in Table 1. Figure 8, Figure 9, Figure 10 and Figure 11 show, from the left to right, the adopted Hermite pdf of r, the first-order Sobol indices S_i in the central panel, and the total indices S_Ti in the right-hand panel.

The first-order Sobol indices in Figure 8, Figure 9, Figure 10 and Figure 11 exhibit a symmetric pattern that reflects the symmetry of the stochastic model, see Figure 7. In the Gaussian baseline case (α_3,r, α_4,r)=(0, 3), S_fy1=S_fy2≈0.25 and S_r1=S_r2≈0.14, with the remainder 1−∑_iS_i≈0.22 attributed to higher-order effects. For the Hermite cases (−1, 4), (−2, 5), (−3, 6) the indices of the radii increase monotonically to about 0.17, 0.19, and 0.21, while the indices of the yield strengths decrease to approximately 0.23, 0.22, and 0.21. The higher-order remainder decreases concurrently from about 0.20 to 0.18 and 0.17. Departures of the radius distribution toward negative skewness and leptokurtosis therefore increase the first-order contribution of r₁ and r₂ to V(Y), accompanied by a corresponding reduction in the higher-order remainder.

The total Sobol indices S_Ti in Figure 8, Figure 9, Figure 10 and Figure 11 quantify the overall contribution of each input to V(Y), including all interaction terms. In the Gaussian baseline case (α_3,r, α_4,r)=(0, 3) the yield strengths dominate, while the radii have comparatively smaller total effects. As the radius distribution departs progressively from the Gaussian shape toward (−1, 4),(−2, 5), and (−3, 6), the total effects of r₁ and r₂ increase, whereas the total effects of f_y1 and f_y2 decrease. In the extreme case (−3, 6) the total effect of each radius approaches that of the corresponding yield strength indicating a comparable overall influence on V(Y). The differences S_Ti−S_i show that interaction effects remain more pronounced for the yield strengths, whereas the radii gain importance primarily due to their first-order contributions.

5.2.3. Output Statistics and Design Quantiles

For the resistance R, it is essential to identify cases in which the random output exhibits higher standard deviation, as these cases are potentially more critical from a design perspective.

The results of the sensitivity analysis indicate that the standard deviation of R increases when deviations from the Gaussian distribution are accompanied by increasing kurtosis and decreasing skewness of the radius r. This variation in skewness and kurtosis of the input radius r has relatively little influence on the mean value of resistance R but a considerable effect on its standard deviation. Consequently, an increase in kurtosis combined with a decrease in skewness of the radius r leads to a higher coefficient of variation of resistance R, and therefore to a lower design resistance R_d,FORM, which is evaluated according to Equation (38), where R_d,FORM corresponds to the target reliability index β_d=3.8, which in turn represents a failure probability of P_f = 7.2⋅10⁻⁵ [1].

The input random variables are defined in Table 1. The statistical analysis of R was performed by numerical integration, following the same procedure used for the evaluation of the first two statistical moments in the previous chapter. The yield strength is assumed to follow the green Hermite pdf introduced in Section 3.1., see Figure 2. The results of the statistical analysis of the output resistance R are presented in Table 2.

Table 2 reports the realised moments of R and the corresponding design resistance R_d,FORM for the Gaussian baseline and for Hermite pdfs with moment pairs (α_3,r, α_4,r)=(−1,4),(−2,5),(−3,6). The mean value of R changes only marginally (by approximately +0.07%), whereas the standard deviation increases monotonically by about 9%. The output skewness shifts from positive to negative, and the output kurtosis increases markedly, indicating the development of heavier tails. As a result of this increase in the standard deviation of R, the design value R_d,FORM decreases from 33.3048 to 32.7103 (a reduction of approximately 1.8%).

In Equation (38), the design value R_d,FORM corresponds to the 0.12% quantile of a Gaussian pdf, as defined in EN 1990 [1]. According to [1], a Gaussian model is admissible for both loads and structural resistance. However, large skewness and kurtosis and deviations from the Gaussian pdf may make the FORM concept insufficient to ensure reliable design, see, e.g., [73].

A more accurate 0.12% quantile can be obtained by using all four statistical moments listed in the last row of Table 2. For example, for the input parameters θ=(40.17, 2.4539, −0.7, 6.93), the 0.12% quantile can be obtained by approximating R using the Hermite pdf function as 27.4154 kN, which is significantly different from the Gaussian design value R_d,FORM=32.7103 kN obtained from Equation (38).

Table 3 compares the design resistances R_d,FORM obtained from Hermite pdfs and from the FORM expression in Equation (38). The first column lists the parameters (μ, σ, α₃, α₄) of the Hermite models for R, determined numerically to achieve the four moments of R given in Table 2. The second column provides the design resistance R_d,HERM-INT evaluated as the 0.12% quantile computed by numerical integration of the Hermite pdf. The third column lists the design value R_d,HERM calculated using Equation (47). It is evident that the values R_d,HERM-INT and R_d,HERM ​​are in excellent agreement.

The fourth column reports the design value R_d,FORM computed by the FORM method, which depends solely on the mean μ and standard deviation σ. Comparison of the design values ​​from the Hermite distribution R_d,HERM with the corresponding R_d,FORM value (last column) shows a significant influence of skewness and kurtosis on the design quantile. The discrepancy increases as α₃ becomes more negative and α₄ increases.

In Table 3, several modelling specifics deserve mention. In Table 2, skewness of 0.2161 and kurtosis of 3.0073 do not satisfy the admissibility condition in Equation (27). Therefore, the lower branch was used. This approach reproduces the target moments with sufficient accuracy, as shown in the first row of Table 3. In this case, a minor calibration was performed using the parameters θ=(40.1443, 2.2517, 0.2051, 2.898) of Hermite pdf, which yields realised moments M(θ)=(40.1431, 2.2494, 0.2160, 3.0074), practically identical to the moments listed in Table 2. For the other leptokurtic cases (α₄>3), the first two Hermite parameters match the first two moments in Table 2 exactly, while the third and fourth moments were slightly calibrated. With these settings, all Hermite pdf functions in Table 3 reproduce the moments from Table 2 with an accuracy of two decimal places.

In Table 3, as skewness becomes more negative and kurtosis increases, the deviation between the design values R_d,FORM and R_d,HERM increases substantially. This confirms that the FORM approach, which neglects higher-order moments, systematically overestimates the design resistance when the output distribution of R departs from the Gaussian pdf due to negative skewness and higher kurtosis.

The value of R_d,HERM is a more accurate estimate of the design resistance compared to R_d,FORM. For example, the negative skewness and higher kurtosis of the resistance R place more probability mass near smaller values of R. With μ_R and σ_R held constant, these changes reduce the left-hand quantiles of R and thicken the lower tail of the distribution. For independent load effects F, this results in an increased failure probability P_f=P(F>R) relative to the Gaussian model used in FORM with the same μ_R and σ_R. In such cases, the reliability of the design is more accurately ensured by the value of R_d,HERM, since R_d,HERM < R_d,FORM. If the resulting P_f exceeds the EN 1990 target value of 7.2⋅10⁻⁵ for the ultimate limit state, corresponding to the reliability index β_d=3.8, the design no longer satisfies the prescribed reliability level.

Under the opposite regime, where α_3,R>0 and α_4,R<3, the failure probability decreases. With μ_R and σ_R again held constant, the left-hand quantiles of R increase and the lower tail becomes thinner. For independent load effects F, the failure probability P_f=P(F>R) is therefore smaller than for the moment-matched Gaussian model used in FORM, and the 0.12% quantile of R is higher. In this situation, the FORM-based design is conservative.

The results in Table 2 were verified by Monte Carlo simulation using n=10⁶ independent samples of (f_y1, r₁, f_y2, r₂), with inputs as listed in Table 1. For each sample, the resistance R was computed, and the output moments were estimated. The comparison between Table 2 and Table 4 shows agreement to three significant digits for both ${\mu }_R^{*}$ and ${\sigma }_R^{*}$. The estimates ∣${\alpha }_{3,R}^{*}$∣ and ${\alpha }_{4,R}^{*}$ in Table 4 are smaller, which is consistent with finite-sample estimation in cases where the distribution of R exhibits a heavy lower tail.

Figure 12 shows the design value R_d,MC. It is defined as the empirical 0.12% quantile and corresponds to the 1200th member in the ascending order of the Monte Carlo sample of R for n=10⁶.

The estimate R_d,MC=29.4 kN is higher than R_d,HERM-INT≈R_d,HERM≈27.5 kN in Table 3, which is consistent with the lower sample kurtosis ${\alpha }_{4,R}^{*}$=4.64 in Table 4, compared with ${\alpha }_{4,R}^{*}$=7.85 in Table 2.

6. Discussion

This work provides a Hermite-based framework suitable for reliability analysis and the limit states method. A closed-form expression for the design quantile under a Hermite pdf was derived for both resistance and load effects. The formulation uses the Gaussian design quantile as an input and produces the Hermite design value through a calibrated cubic translation. Agreement with numerical integration is within two per mille in all cases examined, validating the analytical construction. The closed-form expressions are given in Equations (47) and (48). Admissibility and implementation rules for the Hermite pdf were also established. The monotonicity bound, the positivity check for the lower branch, and the balanced Jacobian for the upper branch form a practical checklist for engineering use.

The case study demonstrates the significant influence of higher-order moments in the limit state method. The input quantities were the yield strength f_y and the radius r of a circular cross-section, while the output quantity was the structural resistance R. The mean resistance changed negligibly across all scenarios, whereas the standard deviation increased by approximately nine percent as the input radius moved toward negative skewness and higher kurtosis. For the most leptokurtic setting, the Hermite design resistance was 27.4957 kN, compared with the FORM value of 32.7103 kN based on the same mean and standard deviation. This result confirms that neglecting skewness and kurtosis can lead to a non-conservative design when lower tails become heavier.

The Sobol sensitivity analysis showed that the variance of R is strongly influenced by the shape of the distribution of r, even when the first two moments are fixed. As kurtosis increased and skewness became more negative, the first-order sensitivity indices of r grew, while the higher-order remainder decreased. The total Sobol indices quantified the overall contribution of each input r and f_y to the variance of R. In the most skewed and heavy-tailed case of r, the total indices for r and f_y were nearly equal. This finding highlights that the distributional characteristics of r can significantly thicken the lower tail of R and thereby reduce the structural reliability.

Monte Carlo simulation with one million samples confirmed the mean values and standard deviations but yielded smaller magnitudes of skewness and kurtosis, consistent with the slow convergence of tail quantiles and higher-order moments. The 0.12% quantile obtained from Monte Carlo (29.4 kN) was higher than the Hermite design value (27.5 kN), owing to the smaller kurtosis realised in the simulated sample compared with the calibrated analytical model. The FORM design value remained higher than both quantile estimates. These results demonstrate that standard Monte Carlo simulation may be inadequate when higher-order moments significantly influence reliability, as the estimation of tail quantiles and higher-order moments converges slowly. Importance sampling or tail-focused variance-reduction techniques are therefore recommended for future validation of very small quantiles.

It can be concluded that a FORM design value based solely on μ and σ may be insufficient whenever R or F depart from the Gaussian pdf, whether that departure arises from non-Gaussian inputs or from model mappings of Gaussian inputs. In such cases, a calibrated Hermite pdf should be considered, accompanied by admissibility checks and, where necessary, a brief numerical refinement of the input parameters to reproduce the target moments accurately. If R or F depart from the Gaussian pdf through skewness and kurtosis and exhibit a lower or upper tail, the Hermite quantile defined by Equations (47) and (48) should be reported alongside the FORM value.

In engineering applications, the results of the reliability analysis depend on how well the input pdfs reproduce skewness and kurtosis. For the two input random material characteristics studied (yield strength and ductility), the results show that the cubic Hermite polynomial model can represent the shape of the histogram with high accuracy. When the four input moments of the Hermite pdf are interpreted as density parameters, numerical adjustment of these parameters can bring the realised four moments of the Hermite pdf into close agreement with the moments of the approximated histogram. For small skewness and kurtosis, the prescribed parameters of the Hermite pdf are close to the realised moments of the Hermite pdf. Although the derivation treats the first four moments as parameters of the Hermite pdf, the assumption of small skewness and kurtosis need not hold in technical applications. On the lower branch, the realised first four moments of the Hermite pdf may differ slightly from the prescribed targets. For the lower branch (platykurtosis), fitting of the Hermite pdf to an empirical histogram was demonstrated for yield strength. The example showed that numerical tuning of the Hermite parameters can yield an excellent approximation, as confirmed by goodness-of-fit tests. On the upper branch (leptokurtosis), the first two moments are matched exactly, while small deviations may appear in the third and fourth moments. For this branch, fitting of the Hermite pdf to an empirical histogram was presented for steel ductility.

7. Conclusion

A Hermite-based framework for reliability assessment within the limit states method was presented. A closed-form expression for the design quantile under a Hermite density was derived and shown to agree closely with results obtained by numerical integration. Admissibility and implementation conditions were defined, providing a practical checklist for engineering use.

The closed-form Hermite design quantiles are backward compatible. In the Gaussian limit α₃→0 and α₄→3, they reduce to the FORM design values [1], allowing the method to be incorporated without disrupting existing design practice. Under extreme input conditions, the Hermite model can produce useful light-tailed or even bimodal shapes. However, these cases require explicit verification of normalisation and positivity before use.

The case study demonstrated that departures from Gaussian behaviour due to skewness and kurtosis can significantly reduce the design resistance when the lower tail becomes pronounced. Variance-based sensitivity analysis showed that the distributional shape of the model inputs affects the variance of the structural response even when the first two moments are fixed. The convergence of estimators for tail quantiles and higher-order moments was found to be slow under basic Monte Carlo simulation. For very small quantiles, tail-focused importance sampling or subset simulation is recommended.

In engineering practice, the Hermite design quantile should be reported together with the FORM design value based on the mean and the standard deviation. The resulting failure probability should then be verified against the target reliability specified in the relevant design code. The applicability of the Hermite-based approach remains subject to the stated admissibility checks. Further development is expected to focus on dependent input variables, joint non-Gaussian models for resistance and load effects, and direct calibration to target quantiles using efficient variance-reduction techniques.