Estimation methods for the ratio of medians of three-parameter lognormal distributions containing zero values and their application to wind speed data from northern Thailand

The fiducial method

Hannig, Iyer & Patterson (2006) described the fiducial method based on the fiducial GPQ, a GPQ subclass introduced by Weerahandi (1993). Fiducial techniques have been used in constructing CIs in several research studies (Kharrati-Kopaei, Malekzadeh & Sadooghi-Alvandi, 2013; Li, Zhou & Tian, 2013; Maneerat, Niwitpong & Niwitpong, 2020). Let Y = (Y₁, Y₂, …, Y_n) be a random variable with probability density function f_Y(y; λ, θ), where λ and θ are vectors of the parameter of interest and the nuisance parameter, respectively. Moreover, let y = (y₁, y₂, …, y_n) be the observed values of Y, and assume that fiducial GPQ T(Y; y, λ, θ) is only a function of λ. This method is especially associated with the fiducial inference proposed by Fisher (1935). The fiducial GPQ CI depends on the fiducial GPQ defined by Hannig, Iyer & Patterson (2006) in Definition 1. Recently, Chankham, Niwitpong & Niwitpong (2022) recommended the fiducial GPQ-based CI for estimating the coefficient of variation of an inverse gaussian distribution when the sample size was small. Similarly, the performance of a fiducial GPQ CI in terms of expected length was the shortest when used to estimate the common coefficient of variation of delta-lognormal distributions by Yosboonruang, Niwitpong & Niwitpong (2022).

Definition 1.

A GPQ T(Y; y, λ, θ) for a parameter λ is called a fiducial generalized pivotal quantity (FGPQ) if it satisfies the following conditions: (FGPQ1) Given Y = y, the T(Y; y, λ, θ) distribution is free of all parameters. (FGPQ2) For every y ∈ R⁺, the observed pivotal T(Y; y, λ, θ) = λ.

If the T(Y; y, λ, θ) satisfies the conditions (FGPQ1) and (FGPQ2), then it is possible to construct the 100(1 − φ)% fiducial GPQ-based CI for λ is [T_λ(φ/2), T_λ(1 − φ/2)]; T_λ(φ) is the φ^th percentile of T(Y; y, λ, θ).

The NA method

According to probability theory, the concept behind this method is the assumption that the approximate distributions of all of the samples approach a normal distribution pattern if the sample size is sufficiently large. This idea is integrated with the central limit theorem, in which the distribution of a given sample mean is approximated as a normal pattern if the sample size is sufficiently large under the assumption that all of the samples are similar to each other regardless of the shape of the population distribution. Recently, Maneerat, Niwitpong & Nakjai (2022b) proposed an NA-based CI for the median of a TPLN distribution, which performed well for a large sample size. Thus, we also considered the NA method.

Given a set of observations, threshold γ_i can be estimated by using the Adam algorithm to find the MLEs of the mean and variance (${\hat{\mu }}_{X_i},{\hat{\sigma }}_{X_i}^2$). Here, the medians of TPLN models with zero observations can be log-transformed to become (24)${\displaystyle ln{\eta }_i=ln{\rho }_i+ln\left[{\gamma }_i+exp\left({\mu }_{X_i}\right)\right]}$

which can be approximated by using (${\hat{\gamma }}_i,{\hat{\mu }}_{X_i},{\hat{\rho }}_i$) to give $ln{\hat{\eta }}_i=ln{\hat{\rho }}_i+ln\left[{\hat{\gamma }}_i+exp\left({\hat{\mu }}_{X_i}\right)\right]$. Using the delta method, the variance of lnρ_i can be derived as (25)${\displaystyle V_{ln{\hat{\rho }}_i}=\frac{1-{\rho }_i}{n{\rho }_i}.}$

Likewise, McKean & Schrader (1984) estimated the variance of the median as a distribution-free estimate defined as follows: (26)${\displaystyle {\text{V}}_{{\hat{\eta }}_i}^{MS}={\left[\frac{X_{i\left(n_{i1}-c+1\right)}+X_{i\left(c\right)}}{2\left(1.96\right)}\right]}^2}$

where X_i(r) denotes the r^th order statistic in a random sample drawn from a TPLN model with zero observations of size n_i1, and $c=\left[\left(n_{i1}+1\right)/2\right]-1.96\left(\sqrt{n_{i1}/4}\right)$. Later, Hettmansperger & Sheather (1986) claimed that V^MS is a consistent estimator of the variance of the median. Similarly, the variance of $ln{\hat{\eta }}_i$ can be derived by applying the delta method as follows: (27)${\displaystyle {\text{V}}_{ln{\hat{\eta }}_i}\cong \frac{{\text{V}}_{{\hat{\eta }}_i}^{MS}}{{\left[{\text{E}}_{{\hat{\eta }}_i}\right]}^2}}$

where ${\text{E}}_{{\hat{\eta }}_i}={\hat{\gamma }}_i+exp\left[{\hat{\mu }}_{X_i}+{\hat{\sigma }}_{X_i}^2{\left(2n_{i1}\right)}^{-1}\right]$ is the expectation of ${\hat{\eta }}_i$. Estimated variance ${\widehat{\text{V}}}_{ln{\hat{\eta }}_i}$ is obtainable by replacing $\left({\hat{\gamma }}_i,{\hat{\mu }}_{X_i},{\hat{\sigma }}_{X_i}^2,{\hat{\rho }}_i\right)$ from the sample. Thus, by applying the central limit theorem , the random variable W_i can be defined as (28)${\displaystyle W_i=\frac{ln{\hat{\eta }}_i-ln{\eta }_i}{\sqrt{{\widehat{\text{V}}}_{ln{\hat{\eta }}_i}}}}$

which approaches a standard normal distribution as n → ∞. Subsequently, the 100(1 − φ)% NA-based CI for ω can be written as (29)${\displaystyle \left[l_{\omega }^{\text{N}},u_{\omega }^{\text{N}}\right]=exp\left[\left(ln{\hat{\eta }}_1-ln{\hat{\eta }}_2\right)\mp W_{1-\phi /2}\sqrt{{\widehat{\text{V}}}_{ln{\hat{\eta }}_1}+{\widehat{\text{V}}}_{ln{\hat{\eta }}_2}}\right]}$

where W_φ denotes the φ^th percentile of a standard normal distribution. Algorithm 4 was used to compute the CP of the NA (CP_N).

The Bayesian method

Bayesian methods are based on treating probability as beliefs rather than frequencies. Given unknown parameter θ, a prior distribution p(θ) represents the subjective belief as a subjective distribution formulated before the data are seen. The posterior distribution is obtained from a prior that is updated with the likelihood function (or sample information) by using Bayes’ rule (Casella & Berger, 2002). Importantly, the posterior distribution is considered to be a random quantity and can be used to make a statement about θ, For example, the point and interval estimates of θ can be computed by using its posterior. Equal-tailed Bayesian intervals based on the Jeffreys’ and uniform priors based on the posterior densities of the zero proportion and the variance have been shown to perform well in certain scenarios (Yosboonruang, Niwitpong & Niwitpong , 2022).

Here, the Bayesian CI for ω is formulated based on the Bayesian method. First, we define an informative prior for our objective assumption depending on the amount of information available in the data as (30)${\displaystyle P\left({\gamma }_{1i},{\mu }_{X_i},log{\sigma }_{X_i},{\rho }_i\right)=\text{constant}\ast {\left(1-{\rho }_i\right)}^{{\alpha }_i-1}{\rho }_i^{{\beta }_i-1}\Gamma \left({\alpha }_i+{\beta }_i\right){\left[\Gamma \left({\alpha }_i\right)\Gamma \left({\beta }_i\right)\right]}^{-1}.}$

The posterior densities of (γ_i, ρ_i) are obtained by obtaining Eq. (30) with the likelihood function (2) as

(31)${\displaystyle f_{{\gamma }_i}=P\left({\gamma }_i|X,{\hat{\mu }}_{X_i},{\hat{\sigma }}_{X_i}^2\right)=\left[\text{constant}\ast {\left(1-{\rho }_i\right)}^{{\alpha }_i-1}{\rho }_i^{{\beta }_i-1}\Gamma \left({\alpha }_i+{\beta }_i\right){\left[\Gamma \left({\alpha }_i\right)\Gamma \left({\beta }_i\right)\right]}^{-1}\right]\ast \left[{\left(1-{\rho }_i\right)}^{n_{i0}}{\rho }_i^{n_{i1}}{\left(2\pi {\sigma }_{X_i}^2\right)}^{-n_{i1}/2}exp\left\{-\sum _{j:x_{ij}>0}ln\left(X_{ij}-{\gamma }_i\right)-\frac{1}{2{\sigma }_{X_i}^2}\sum _{j:x_{ij}>0}{\left[ln\left(X_{ij}-{\gamma }_i\right)-{\mu }_{X_i}\right]}^2\right\}\right]}$ (32)${\displaystyle \propto {\left({\hat{\sigma }}_{X_i}^2\right)}^{-n_{i1}/2}exp\left\{-\sum _{j:x_{ij}>0}ln\left(X_{ij}-{\gamma }_i\right)-\frac{1}{2{\hat{\sigma }}_{X_i}^2}\sum _{j:x_{ij}>0}{\left[ln\left(X_{ij}-{\gamma }_i\right)-{\hat{\mu }}_{X_i}\right]}^2\right\}}$ (33)${\displaystyle f_{{\rho }_i}=P\left({\rho }_i|X\right)=\left[\text{constant}\ast {\left(1-{\rho }_i\right)}^{n_{i0}+{\alpha }_i-1}{\rho }_i^{n_{i1}{\beta }_i-1}\Gamma \left({\alpha }_i+{\beta }_i\right){\left[\Gamma \left({\alpha }_i\right)\Gamma \left({\beta }_i\right)\right]}^{-1}\right]\ast \left[{\left(2\pi {\sigma }_{X_i}^2\right)}^{n_{i1}/2}exp\left\{-\sum _{j:x_{ij}>0}ln\left(X_{ij}-{\gamma }_i\right)-\frac{1}{2{\sigma }_{X_i}^2}\sum _{j:x_{ij}>0}{\left[ln\left(X_{ij}-{\gamma }_i\right)-{\mu }_{X_i}\right]}^2\right\}\right]}$ (34)${\displaystyle \propto {\left(1-{\rho }_i\right)}^{n_{i0}+{\alpha }_i-1}{\rho }_i^{n_{i1}+{\beta }_i-1}.}$

Next, we apply NA to the posterior distribution of κ_i = (μ_Xi, logσ_Xi); the logarithm of the posterior density is approximated by using a quadratic function of κ_i. The second derivatives of the log-posterior density are needed for constructing the approximation. From Eq. (3), the log-likelihood can be expressed as (35)${\displaystyle l\left({\kappa }_i|X\right)=n_{i0}ln\left(1-{\rho }_i\right)+n_{i1}ln{\rho }_i-\frac{n_{i1}}{2}ln\left(2\pi \right)+n_{i1}ln{\sigma }_{X_i}-\sum _{j:x_{ij}>0}ln\left(X_{ij}-{\gamma }_i\right)-\frac{1}{2{\sigma }_{X_i}^2}\left[\left(n_{i1}-1\right){\hat{\sigma }}_{X_i}^2+n_{i1}{\left({\hat{\mu }}_{X_i}-{\mu }_{X_i}\right)}^2\right]\propto \text{constant.}+n_{i1}ln{\sigma }_{X_i}-\frac{1}{2exp\left(2ln{\sigma }_{X_i}\right)}\left[\left(n_{i1}-1\right){\hat{\sigma }}_{X_i}^2+n_{i1}{\left({\hat{\mu }}_{X_i}-{\mu }_{X_i}\right)}^2.\right]}$

After that the first and second derivatives of (μ_Xi, logσ_Xi) respectively become

(36)${\displaystyle \left(\frac{\partial l\left({\kappa }_i|X\right)}{\partial {\mu }_{X_i}},\frac{{\partial }^{2}l\left({\kappa }_i|X\right)}{\partial {\mu }_{X_i}^2}\right)=\left(\frac{n_{i1}\left({\hat{\mu }}_{X_i}-{\mu }_{X_i}\right)}{{\sigma }_{X_i}^2},-\frac{n_{i1}}{{\sigma }_{X_i}^2}\right)}$ (37)${\displaystyle \left(\frac{\partial l\left({\kappa }_i|X\right)}{\partial ln{\sigma }_{X_i}},\frac{{\partial }^{2}l\left({\kappa }_i|X\right)}{\partial {\left[ln{\sigma }_{X_i}\right]}^2}\right)=\left(-n_{i1}+{\sigma }_{X_i}^{-2}\left[\left(n_{i1}-1\right){\hat{\sigma }}_{X_i}^2+n_{i1}{\left({\hat{\mu }}_{X_i}-{\mu }_{X_i}\right)}^2\right]\right.,\left.-2{\sigma }_{X_i}^{-1}\left[\left(n_{i1}-1\right){\hat{\sigma }}_{X_i}^2+n_{i1}{\left({\hat{\mu }}_{X_i}-{\mu }_{X_i}\right)}^2\right]\right).}$

The point estimates of (μ_Xi, logσ_Xi) are derived after setting their first derivatives to zero. Meanwhile, the variances of their estimates are respectively obtained by using an inverse Fisher information matrix as follows:

(38)${\displaystyle \left(\begin{array}{c}\hfill {\hat{\mu }}_{X_i}\hfill \\ \hfill ln{\hat{\sigma }}_{X_i}\hfill \end{array}\right)=\left(\begin{array}{c}\hfill {\hat{\mu }}_{X_i}\hfill \\ \hfill ln\sqrt{\left(n_{i1}-1\right)n_{i1}^{-1}{\hat{\sigma }}_{X_i}^2}\hfill \end{array}\right)}$ (39)${\displaystyle \left(\begin{array}{c}\hfill {\text{V}}_{{\hat{\mu }}_{X_i}}\hfill \\ \hfill V_{ln{\hat{\sigma }}_{X_i}}\hfill \end{array}\right)=\left(\begin{array}{c}\hfill {\hat{\sigma }}_{X_i}^2/n_{i1}\hfill \\ \hfill 2n_{i1}\hfill \end{array}\right).}$

After using the Jacobian to back-transform from $ln{\hat{\sigma }}_{X_i}$ to ${\sigma }_{X_i}^2$, the posterior densities of (${\mu }_{X_i},{\sigma }_{X_i}^2$) are respectively approximated as normal distribution as follows:

(40)${\displaystyle f_{{\mu }_{X_i}}=p\left({\mu }_{X_i}|X,{\hat{\sigma }}_{X_i}^2,{\hat{\gamma }}_i\right)=\frac{1}{{\left(2\pi {\sigma }_{X_i}^2/n_{i1}\right)}^{1/2}}exp\left\{-\frac{n_{i1}}{2{\sigma }_{X_i}^2}{\left[ln\left(X_{ij}-{\hat{\gamma }}_i\right)-{\hat{\mu }}_{X_i}\right]}^2\right\}}$ (41)${\displaystyle f_{{\sigma }_{X_i}^2}=p\left({\sigma }_{X_i}^2|X\right)={\left(\frac{n_{i1}+2}{4\pi {\tilde{\sigma }}_{X_i}^4}\right)}^{1/2}exp\left\{-\frac{n_{i1}+2}{4{\tilde{\sigma }}_{X_i}^4}{\left[ln\left(X_{ij}-{\hat{\gamma }}_i\right)-{\tilde{\sigma }}_{X_i}^2\right]}^2\right\}}$

where ${\tilde{\sigma }}_{X_i}^2=n_{i1}{\hat{\sigma }}_{X_i}^2/\left(n_{i1}+2\right)$. Therefore, the posterior density of ω becomes (42)${\displaystyle f_{\omega }=\frac{f_{{\rho }_1}\left[f_{{\gamma }_1}+exp\left(f_{{\mu }_{X_1}}\right)\right]}{f_{{\rho }_2}\left[f_{{\gamma }_2}+exp\left(f_{{\mu }_{X_2}}\right)\right]}}$

Finally, the 100(1 − φ)% Bayesian-based CI for ω is (43)${\displaystyle \left[l_{\omega }^B,u_{\omega }^B\right]=\left[f_{\omega }\left(\phi /2\right),f_{\omega }\left(1-\phi /2\right)\right]}$

where f_ω(φ) denotes the φ^th percentile of f_ω. The CP of Bayesian CI (CP_B) can be computed by using Algorithm 5.