Cylindrical Models Motivated through Extended Sine-Skewed Circular Distributions

Abstract

A class of cylindrical distributions, which include the Weibull-von Mises distribution as a special case, is considered. This distribution is obtained by combining the extended sine-skewed wrapped Cauchy distribution (marginal circular part) with the Weibull distribution (conditional linear part). This family of proposed distributions is shown to have simple normalizing constants, easy random number generation methods, explicit moment expressions, and identifiability in parameters. In particular, the marginal distribution of the circular random variable, and its conditional distribution given a linear random variable give relatively stronger skewness than those of existing cylindrical models. Some Monte Carlo simulations and real data analysis are performed to investigate the feasibility and tractability of the proposed models.

1. Introduction

Cylindrical data are a type of multivariate data that consist of a pair of variables, one representing an angular observation on the unit circle ${\mathbb{S}}^1$ and the other representing a real or non-negative observation on $[0,\infty )$. Because the variable representing an angle can be identified with a vector on the unit circle, these data are sometimes called linear–circular data. Such data can be found in various fields, such as environmental studies, biology, sports science, and others. Wind direction and SO₂ concentration have been explored by [1], and the directions and speed of ocean currents in the Adriatic Sea have been analyzed by [2] through hidden Markov models with cylindrical distributions. The Bayesian spatio-temporal modeling for wave heights and wave directions has been proposed by [3].

Some ways to model cylindrical data have been proposed by [4,5,6,7]. As an extension of the copula-based distribution of [7,8], we have proposed a cylindrical distribution called the Weibull sine-skewed von Mises (WeiSSVM) distribution, in which the marginal distribution of a circular random variable $\mathsf{\Theta }$ is the sine-skewed wrapped Cauchy distribution proposed by [9] and the conditional distribution of a linear random variable X given $\mathsf{\Theta }$ is a Weibull distribution. The WeiSSVM density has a simple normalizing constant and is amenable to calculations. This characteristic is referred to as tractability, which is one of the desirable properties of a class of flexible distributions in [10]. The WeiSSVM distribution also provides explicit expressions for the higher-order moments for X and $\mathsf{\Theta }$. This enables the representation of characteristics of the distribution, such as linear–circular correlation, with a high degree of interpretability in the parameter values.

On the other hand, when fitting the WeiSSVM distribution to observations, especially in cases where circular data exhibit significant asymmetry near the mode, the maximum likelihood estimator of the skewness parameter may approach the boundary of the parameter space. This phenomenon can be seen in Table 3 of [8]. This suggests that a distribution with stronger skewness is required for the circular random variable. In addition, the standard error cannot be appropriately estimated by the Hessian matrix of the negative log-likelihood function when the maximum likelihood estimate (MLE) takes values on the boundary of the parameter space, as the first-order conditions for maximizing the log-likelihood function are violated.

To overcome the shortcomings stated above, we consider a cylindrical distribution with the skewing function in the WeiSSVM distribution replaced by a Beta distribution function. We also show that the proposed distributions have simple normalizing constants, easy random number generation methods, explicit moment expressions, and identifiability in parameters. Through a Monte Carlo simulation study and two real data analyses, we observed that the estimated skewness parameter moves away from the boundary of the parameter space.

The rest of this paper is organized as follows. In Section 2, we introduce a class of cylindrical model whose marginal distribution of $\mathsf{\Theta }$ is the extended sine-skewed (ESS) wrapped Cauchy distribution, as proposed in [11]. In addition, the proposed class of cylindrical distribution has the Weibull distribution as the conditional distribution of the linear random variable X given $\mathsf{\Theta }$. The ESS-wrapped Cauchy distribution can exhibit a greater degree of skewness compared to the sine-skewed wrapped Cauchy distribution, making the proposed model more flexible than the existing cylindrical models. Section 2.1 presents some higher order moments of random variables X and $\mathsf{\Theta }$, Section 2.3 describes an easy random number generation method from the proposed distribution, and Section 2.4 presents an identifiability result for a family of the proposed distributions. In Section 3, the performance of the proposed models is studied through Monte Carlo simulations and, in Section 4 analysis of two real data sets is presented to demonstrate the performance of the proposed models. Section 5 presents concluding remarks. Some lemmas and proofs are given in Appendix A.

2. A New Model for Linear–Circular Data

The cylindrical distribution proposed by [8] features a simple normalizing constant, and its marginal distribution of the linear random variable X blue and the conditional distribution of the circular random variable $\mathsf{\Theta }$ become well known: the Weibull and sine-skewed wrapped Cauchy distributions. Hence, it is easy to calculate the log-likelihood function and to generate random numbers from this distribution. As can be seen from [12], the skewing function of the sine-skewed wrapped Cauchy distribution, which is the marginal distribution of the circular component $\mathsf{\Theta }$, is given by the distribution function of a uniform random variable. To extend the WeiSSVM distribution, we replace the skewing function with the distribution function of a density function on the interval $[-1,1]$ proposed in [13]:

$$\begin{array}{c}\hfill G_m\left(x\right)={\int }_{-1}^x{g}_m\left(t\right)dt,\end{array}$$

where

$$g_m\left(x\right)=\frac{\mathsf{\Gamma }\left(2\right(m+1\left)\right)}{2^{2m+1}\mathsf{\Gamma }{(m+1)}^2}{(1-x^2)}^m(-1\le x\le 1),$$

and $\mathsf{\Gamma }(·)$ represents the gamma function. $m\ge 0$ is the degree of the skewness parameter, which is assumed to be pre-specified for simplicity and is referred to as the “order”. We call this distribution the Weibull extended sine-skewed von Mises (WeiESSVM) distribution.

Here, we introduce a joint probability density of $(\mathsf{\Theta },X)$,

$$\begin{array}{cc}\hfill f_{\mathrm{WeiESSVM}}^{\left(m\right)}(\theta ,x)& =\frac{\alpha {\beta }^{\alpha }}{\pi \cosh \left(\kappa \right)}{G}_m\left(\lambda \sin (\theta -\mu )\right)x^{\alpha -1}\hfill \\ \hfill & \times \exp \left\{-{\left(\beta x\right)}^{\alpha }\left(1-\tanh \left(\kappa \right)\cos (\theta -\mu )\right)\right\},\hfill \end{array}$$

(1)

where $(x,\theta )\in [0,\infty )\times [-\pi ,\pi )$, $\tanh (·)$ denotes the hyperbolic-tangent function, $\alpha >0$, $\beta >0$, $-\pi \le \mu <\pi$, $\kappa >0$ and $-1\le \lambda \le 1$.

Although the proposed density (1) can be defined for any $m\ge 0$, it is tractable in the case where m is a nonnegative integer since the function $G_m(·)$ is expressed as

$$\begin{array}{cc}\hfill G_m\left(x\right)& =C_m\sum _{\mathcal{l}=0}^m\left(\genfrac{}{}{0pt}{}{m}{\mathcal{l}}\right)\frac{1}{2\mathcal{l}+1}{(-1)}^{\mathcal{l}}\left(x^{2\mathcal{l}+1}+1\right)\hfill \\ \hfill & =C_m\sum _{\mathcal{l}=0}^m\left(\genfrac{}{}{0pt}{}{m}{\mathcal{l}}\right)\frac{1}{2\mathcal{l}+1}{(-1)}^{\mathcal{l}}{x}^{2\mathcal{l}+1}+\frac{1}{2}\hfill \end{array}$$

(2)

for $m=0,1,2,\dots$, where $C_m={\displaystyle \frac{\mathsf{\Gamma }\left\{2(m+1)\right\}}{2^{2m+1}\mathsf{\Gamma }{(m+1)}^2}}$.

For $m=0,1,\dots ,5$, $G_m\left(x\right)$ is explicitly defined by

$$\begin{array}{cc}\hfill G_0\left(x\right)& =\frac{x+1}{2},\hfill \\ \hfill G_1\left(x\right)& =\frac{1}{4}(-x^3+3x+2),\hfill \\ \hfill G_2\left(x\right)& =\frac{1}{16}(3x^5-10x^3+15x+8),\hfill \\ \hfill G_3\left(x\right)& =\frac{1}{32}(-5x^7+21x^5-35x^3+35x+16),\hfill \\ \hfill G_4\left(x\right)& =\frac{1}{256}(35x^9-180x^7+378x^5-420x^3+315x+128),\mathrm{and}\hfill \\ \hfill G_5\left(x\right)& =-\frac{1}{512}{(x+1)}^6\left(63x^5-378x^4+938x^3-1218x^2+843x-256\right).\hfill \end{array}$$

Notice that since the function $G_0\left(x\right)$ is the uniform distribution function on $[-1,1]$, the density (1) reduces to the distribution proposed by [8]. Figure 1 depicts plots of the functions $g_m\left(x\right)$ and $G_m\left(x\right)$ for different orders $m\in \{0,1,4\}$. As shown in Section 3.1 of [13], the greater the slope of the function $G_m$ near 0, the stronger the influence of the skewness parameter $\lambda$ on the marginal density of $\mathsf{\Theta }$.

Figure 2 shows contour plots of the density (1) with $\mu =0$, $\kappa =1$, $\alpha =2$ and $\beta =1$. The horizontal axis represents the variable x and the vertical axis represents the variable $\theta$. As Figure 2 illustrates, the shape of the density function $\mathsf{\Theta }$ of a model with $m=2$ is significantly skewed around $\mu =0$, in comparison to the model with $m=0$.

In what follows, we present marginal and conditional density functions of $\mathsf{\Theta }$ and X. The marginal density of $\mathsf{\Theta }$ is given by

$$\begin{array}{c}\hfill f\left(\theta \right)=\frac{1-{\tanh }^2(\kappa /2)}{\pi }\frac{G_m\left(\lambda \sin (\theta -\mu )\right)}{1+{\tanh }^2(\kappa /2)-2\tanh (\kappa /2)\cos (\theta -\mu )},\end{array}$$

(3)

indicating the ESS-Wrapped Cauchy distribution with a concentration parameter $\tanh (\kappa /2)$. The density (3) is derived in exactly the same manner as presented in [8]. The marginal density of X becomes

$$\begin{array}{cc}\hfill f\left(x\right)& =\frac{I_0\left(x^{\alpha }{\beta }^{\alpha }\tanh \left(\kappa \right)\right)}{\cosh \left(\kappa \right)}\alpha {\beta }^{\alpha }{x}^{\alpha -1}\exp \left\{-{\left(\beta x\right)}^{\alpha }\right\},\hfill \end{array}$$

(4)

where $I_0\left(\nu \right)={\left(2\pi \right)}^{-1}{\int }_0^{2\pi }\exp \left(\nu \cos \theta \right)d\theta$ is the modified Bessel function of the first kind and order zero. From Equation (2), $2G_m\left(\lambda \sin \left(\theta \right)\right)=1+\tau \left(\theta \right)$ where $\tau \left(\theta \right)$ is an odd function of $\theta$. Hence, the density (4) is derived in the same way as in [8].

Next, the conditional density of $\mathsf{\Theta }$ given $X=x$ becomes:

$$\begin{array}{cc}\hfill f\left(\theta \right|x)& =\frac{1}{\pi I_0\left(x^{\alpha }{\beta }^{\alpha }\tanh \left(\kappa \right)\right)}{G}_m\left(\lambda \sin (\theta -\mu )\right)\exp \left\{{\left(\beta x\right)}^{\alpha }\tanh \left(\kappa \right)\cos (\theta -\mu )\right\},\hfill \end{array}$$

(5)

which means that the conditional distribution follows an ESS-von Mises distribution of order m with concentration ${\left(\beta x\right)}^{\alpha }\tanh \left(\kappa \right)$. Finally, the conditional density of X given $\mathsf{\Theta }=\theta$ is expressed as

$$\begin{array}{cc}\hfill f\left(x\right|\theta )& =\alpha {\left[\beta {\left\{1-\tanh \left(\kappa \right)\cos (\theta -\mu )\right\}}^{1/\alpha }\right]}^{\alpha }{x}^{\alpha -1}\hfill \\ \hfill & \times \exp \left[-{\left\{\beta {(1-\tanh \left(\kappa \right)\cos (\theta -\mu ))}^{1/\alpha }x\right\}}^{\alpha }\right].\hfill \end{array}$$

(6)

This implies that the conditional distribution of X follows the Weibull distribution with $\alpha$ being a shape parameter and $1/\left\{\beta {(1-\tanh \left(\kappa \right)\cos (\theta -\mu ))}^{1/\alpha }\right\}$ being a scale parameter. It is noteworthy that the marginal density (4) of X and the conditional density (6) of X given $\mathsf{\Theta }=\theta$ are exactly the same as those of [8].

2.1. Moment Expressions

Throughout the paper, we use the notation $E_0[·]$ for the expectation under the density (1) with $\mu =0$. In this section, we present the expressions for moments $E_0\left[X^k\cos \left(s\mathsf{\Theta }\right)\right]$ and $E_0\left[X^k\sin \left(s\mathsf{\Theta }\right)\right]$ ($k\in \mathbb{N}$, $s\in \{0,1,2,\dots \}$) in explicit forms. To do that, we introduce the following symbol as an extension of the Legendre function,

$$\begin{array}{cc}\hfill A_{\nu }^s\left(x\right)& =\frac{1}{\pi }{\int }_0^{\pi }\cos \left(s\theta \right)\frac{1}{{(x+\sqrt{x^2-1}\cos \theta )}^{\nu +1}}d\theta .\hfill \end{array}$$

(7)

The function $A_{\nu }^s\left(x\right)$ is considered an extension of the Legendre function of the first kind. While the integral (7) is expressed by using the associated Legendre function, it is necessary to impose certain conditions on $\nu$ and s so that the integral is well-defined. Therefore, we use the symbol $A_{\nu }^s\left(x\right)$ instead of the associated Legendre function. In Appendix A.4, we present the conditions under which the integral (7) is expressed by the associated Legendre function.

First, we present $X^k$ times sth cosine moments for $s\in \{0,1,2,\dots \}$. It follows from Lemmas A1 and A2 in Appendix A that

$$\begin{array}{cc}\hfill & E_0\left[X^k\cos \left(s\mathsf{\Theta }\right)\right]={\int }_{-\pi }^{\pi }{\int }_0^{\infty }{x}^k\cos \left(s\theta \right)f_{\mathrm{WeiESSVM}}(\theta ,x)dxd\theta \hfill \\ \hfill & ={\int }_{-\pi }^{\pi }{\int }_0^{\infty }\cos \left(s\theta \right)\frac{\alpha {\beta }^{\alpha }}{2\pi \cosh \left(\kappa \right)}2G_m\left(\lambda \sin \theta \right)x^{\alpha +k-1}\exp \left\{-{\left(\beta x\right)}^{\alpha }\left(1-\tanh \left(\kappa \right)\cos \theta \right)\right\}dxd\theta \hfill \\ \hfill & ={\int }_{-\pi }^{\pi }\cos \left(s\theta \right)\frac{\alpha {\beta }^{\alpha }}{2\pi \cosh \left(\kappa \right)}2G_m\left(\lambda \sin \theta \right){\displaystyle \frac{1}{\alpha {\beta }^{\alpha +k}}}\mathsf{\Gamma }\left({\displaystyle \frac{\alpha +k}{\alpha }}\right){\displaystyle \frac{1}{{\left(1-\tanh \left(\kappa \right)\cos \theta \right)}^{(\alpha +k)/\alpha }}}d\theta \hfill \\ \hfill & =\frac{\alpha {\beta }^{\alpha }}{2\pi \cosh \left(\kappa \right)}{\displaystyle \frac{1}{\alpha {\beta }^{\alpha +k}}}\mathsf{\Gamma }\left({\displaystyle \frac{k}{\alpha }}+1\right){\int }_{-\pi }^{\pi }\cos \left(s\theta \right)2G_m\left(\lambda \sin \theta \right){\displaystyle \frac{1}{{\left(1-\tanh \left(\kappa \right)\cos \theta \right)}^{(\alpha +k)/\alpha }}}d\theta \hfill \\ \hfill & =\frac{{\beta }^{\alpha }}{2\pi \cosh \left(\kappa \right)}{\displaystyle \frac{1}{{\beta }^{\alpha +k}}}\mathsf{\Gamma }\left({\displaystyle \frac{k}{\alpha }}+1\right){\int }_{-\pi }^{\pi }\cos \left(s\theta \right){\displaystyle \frac{1}{{\left(1-\tanh \left(\kappa \right)\cos \theta \right)}^{(\alpha +k)/\alpha }}}d\theta \hfill \\ \hfill & ={\displaystyle \frac{\cosh {\left(\kappa \right)}^{k/\alpha }}{{\beta }^k}}\mathsf{\Gamma }\left({\displaystyle \frac{k}{\alpha }}+1\right){(-1)}^s{A}_{k/\alpha }^s\left(\cosh \left(\kappa \right)\right).\hfill \end{array}$$

(8)

This expression is the same as that of [8] and independent of the order m. Next, we evaluate the expectation $E_0\left[X^k\sin \left(s\mathsf{\Theta }\right)\right]$.

$$\begin{array}{cc}\hfill & E_0\left[X^k\sin \left(s\mathsf{\Theta }\right)\right]={\int }_{-\pi }^{\pi }{\int }_0^{\infty }{x}^k\sin \left(s\theta \right)f_{\mathrm{WeiESSVM}}(\theta ,x)dxd\theta \hfill \\ \hfill & ={\int }_{-\pi }^{\pi }\sin \left(s\theta \right)\frac{\alpha {\beta }^{\alpha }}{2\pi \cosh \left(\kappa \right)}2G_m\left(\lambda \sin \theta \right){\int }_0^{\infty }{x}^{\alpha +k-1}\exp \left\{-{\left(\beta x\right)}^{\alpha }\left(1-\tanh \left(\kappa \right)\cos \theta \right)\right\}dxd\theta \hfill \\ \hfill & ={\int }_{-\pi }^{\pi }\sin \left(s\theta \right)\frac{\alpha {\beta }^{\alpha }}{2\pi \cosh \left(\kappa \right)}2G_m\left(\lambda \sin \theta \right)\hfill \\ \hfill & \times {\displaystyle \frac{1}{\alpha {\beta }^{\alpha +k}}}\mathsf{\Gamma }\left({\displaystyle \frac{\alpha +k}{\alpha }}\right){\displaystyle \frac{1}{{\left(1-\tanh \left(\kappa \right)\cos \theta \right)}^{(\alpha +k)/\alpha }}}d\theta \hfill \\ \hfill & =\frac{\alpha {\beta }^{\alpha }}{2\pi \cosh \left(\kappa \right)}{\displaystyle \frac{1}{\alpha {\beta }^{\alpha +k}}}\mathsf{\Gamma }\left({\displaystyle \frac{k}{\alpha }}+1\right){\int }_{-\pi }^{\pi }\sin \left(s\theta \right)2G_m\left(\lambda \sin \theta \right){\displaystyle \frac{1}{{\left(1-\tanh \left(\kappa \right)\cos \theta \right)}^{(\alpha +k)/\alpha }}}d\theta \hfill \\ \hfill & =\frac{\cosh {\left(\kappa \right)}^{k/\alpha }}{2\pi {\beta }^k}\mathsf{\Gamma }\left({\displaystyle \frac{k}{\alpha }}+1\right){\int }_{-\pi }^{\pi }\sin \left(s\theta \right)2G_m\left(\lambda \sin \theta \right){\displaystyle \frac{1}{{\left(\cosh \left(\kappa \right)-\sinh \left(\kappa \right)\cos \theta \right)}^{(\alpha +k)/\alpha }}}d\theta .\hfill \end{array}$$

(9)

Then, the moments for the cases $m=1$ and $m=2$ can be given by the following equations. The derivation details are provided in Appendix A.

The case of $m=1$: Using the result of Appendix A.1, we have

$$\begin{array}{cc}\hfill & E_0\left[X^k\sin \left(s\mathsf{\Theta }\right)\right]{|}_{m=1}\hfill \\ \hfill & ={\displaystyle \frac{\cosh {\left(\kappa \right)}^{k/\alpha }}{2{\beta }^k}}\mathsf{\Gamma }\left({\displaystyle \frac{k}{\alpha }}+1\right){(-1)}^{s+1}[-{\displaystyle \frac{{\lambda }^3}{8}}\{3A_{k/\alpha }^{|s-1|}\left(\cosh \left(\kappa \right)\right)-3A_{k/\alpha }^{s+1}\left(\cosh \left(\kappa \right)\right)\hfill \\ \hfill & -A_{k/\alpha }^{|s-3|}\left(\cosh \left(\kappa \right)\right)+A_{k/\alpha }^{s+3}\left(\cosh \left(\kappa \right)\right)\}+{\displaystyle \frac{3}{2}}\lambda \{A_{k/\alpha }^{|s-1|}\left(\cosh \left(\kappa \right)\right)-A_{k/\alpha }^{s+1}\left(\cosh \left(\kappa \right)\right)\left\}\right].\hfill \end{array}$$

(10)

The case of $m=2$: It follows from Appendix A.2 that

$$\begin{array}{cc}\hfill & E_0\left[X^k\sin \left(s\mathsf{\Theta }\right)\right]{|}_{m=2}=\frac{\cosh {\left(\kappa \right)}^{k/\alpha }}{{\beta }^k}\mathsf{\Gamma }\left({\displaystyle \frac{k}{\alpha }}+1\right){(-1)}^{s+1}\hfill \\ \hfill & \times \left[{\displaystyle \frac{3{\lambda }^5}{256}}\right\{\left(A_{k/\alpha }^{|s-5|}\left(\cosh \left(\kappa \right)\right)-A_{k/\alpha }^{s+5}\left(\cosh \left(\kappa \right)\right)\right)-5\left(A_{k/\alpha }^{|s-3|}\left(\cosh \left(\kappa \right)\right)-A_{k/\alpha }^{s+3}\left(\cosh \left(\kappa \right)\right)\right)\hfill \\ \hfill & +10\left(A_{k/\alpha }^{|s-1|}\left(\cosh \left(\kappa \right)\right)-A_{k/\alpha }^{s+1}\left(\cosh \left(\kappa \right)\right)\right)\}-{\displaystyle \frac{5{\lambda }^3}{32}}\{3\left(A_{k/\alpha }^{|s-1|}\left(\cosh \left(\kappa \right)\right)-A_{k/\alpha }^{s+1}\left(\cosh \left(\kappa \right)\right)\right)\hfill \\ \hfill & -\left(A_{k/\alpha }^{|s-3|}\left(\cosh \left(\kappa \right)\right)-A_{k/\alpha }^{s+3}\left(\cosh \left(\kappa \right)\right)\right)\}+{\displaystyle \frac{15}{16}}\lambda \left(A_{k/\alpha }^{|s-1|}\left(\cosh \left(\kappa \right)\right)-A_{k/\alpha }^{s+1}\left(\cosh \left(\kappa \right)\right)\right)].\hfill \end{array}$$

(11)

Since the marginal distribution of the circular part is the extended sine-skewed wrapped Cauchy density (3), we obtain the following simple moment expressions. The sth cosine moments of $\mathsf{\Theta }$ under the density (1) with $\mu =0$ are

$$\begin{array}{c}\hfill {\alpha }_s:=E_0\left\{\cos \left(s\mathsf{\Theta }\right)\right\}={\left\{\tanh (\kappa /2)\right\}}^{\left|s\right|}.\end{array}$$

(12)

Accordingly, ${\alpha }_s$ does not depend on the order m. In contrast, from Equation (2), the sth sine moments of $\mathsf{\Theta }$ under the density (1) with $\mu =0$ are given by

$$\begin{array}{c}\hfill {\beta }_s:=2C_m\sum _{\mathcal{l}=0}^m\left(\genfrac{}{}{0pt}{}{m}{\mathcal{l}}\right)\frac{{(-1)}^{\mathcal{l}}}{2\mathcal{l}+1}{\int }_{-\pi }^{\pi }{\left(\lambda \sin \theta \right)}^{2\mathcal{l}+1}\sin \left(s\theta \right)f_0(\theta ;\rho )d\theta ,\end{array}$$

(13)

where $C_m$ is defined in Equation (2), $f_0(\theta ;\rho )$ is the wrapped Cauchy density with $\mu =0$, $\rho =\tanh (\kappa /2)$ and $\left(\genfrac{}{}{0pt}{}{0}{0}\right)\equiv 1$ for convenience. For non-zero location $\mu$, the sth cosine and sine moments of $\mathsf{\Theta }$ are given by

$$\begin{array}{cc}\hfill {\alpha }_{s,\mu }& =E\left\{\cos \left(s\mathsf{\Theta }\right)\right\}=\cos \left(s\mu \right){\alpha }_s-\sin \left(s\mu \right){\beta }_s\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill {\beta }_{s,\mu }& =E\left\{\sin \left(s\mathsf{\Theta }\right)\right\}=\cos \left(s\mu \right){\beta }_s+\sin \left(s\mu \right){\alpha }_s,\hfill \end{array}$$

respectively. These equations are derived in [13], and thus details are omitted here.

2.2. Linear–Circular Correlation

In this subsection, we present some correlations and a linear–circular squared correlation of X and vector $(\cos \mathsf{\Theta },\sin \mathsf{\Theta })$. By Equations (8), (10), (12) and (13), the covariance of X and $\cos \mathsf{\Theta }$ with $\mu =0$, that of X and $\sin \mathsf{\Theta }$, and that of $\cos \mathsf{\Theta }$ and $\sin \mathsf{\Theta }$ are expressed as

$$\begin{array}{cc}\hfill {\mathrm{Cov}}_0(X,\cos \mathsf{\Theta })& =-\frac{\cosh {\left(\kappa \right)}^{1/\alpha }}{\beta }\mathsf{\Gamma }\left(\frac{1}{\alpha }+1\right)\left\{A_{1/\alpha }^1\left(\cosh \left(\kappa \right)\right)+A_{1/\alpha }^0\left(\cosh \left(\kappa \right)\right)\tanh (\kappa /2)\right\},\hfill \\ \hfill {\mathrm{Cov}}_0(X,\sin \mathsf{\Theta })& =\frac{\lambda \cosh {\left(\kappa \right)}^{1/\alpha }\mathsf{\Gamma }\left(1/\alpha +1\right)}{2\beta }[-\frac{3}{8}({\lambda }^2-4)A_{1/\alpha }^0\left(\cosh \left(\kappa \right)\right)\hfill \\ \hfill & +\frac{1}{2}({\lambda }^2-3)A_{1/\alpha }^2\left(\cosh \left(\kappa \right)\right)-\frac{{\lambda }^2}{8}{A}_{1/\alpha }^4\left(\cosh \left(\kappa \right)\right)\hfill \\ \hfill & -\frac{A_{1/\alpha }^0\left(\cosh \left(\kappa \right)\right)}{8{\cosh }^2(\kappa /2)}\left\{-3({\lambda }^2-4)+{\lambda }^2{\tanh }^2(\kappa /2)\right\}],\hfill \\ \hfill {\mathrm{Cov}}_0(\cos \mathsf{\Theta },\sin \mathsf{\Theta })& =\frac{\lambda \tanh (\kappa /2)}{32{\cosh }^2(\kappa /2)}\left\{3({\lambda }^2-4\lambda )+\frac{{\lambda }^2}{{\cosh }^2(\kappa /2)}\right\}.\hfill \end{array}$$

For simplicity, let

$$\begin{array}{cc}\hfill {\xi }_1(\alpha ,\kappa )& =A_{2/\alpha }^0\left(\cosh \left(\kappa \right)\right)\mathsf{\Gamma }\left(\frac{2}{\alpha }+1\right)-A_{1/\alpha }^0{\left(\cosh \left(\kappa \right)\right)}^2\mathsf{\Gamma }{\left(\frac{1}{\alpha }+1\right)}^2\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill {\xi }_2(\kappa ,\lambda )& =128{\cosh }^2(\kappa /2)-{\lambda }^3\tanh (\kappa /2)+3({\lambda }^3-4\lambda ).\hfill \end{array}$$

Using these results, we have the following correlations:

$$\begin{array}{cc}\hfill r_{XC}:& =\mathrm{Corr}(X,\cos \mathsf{\Theta })\hfill \\ \hfill & =\frac{\sqrt{2}\cosh (\kappa /2)\mathsf{\Gamma }(1/\alpha +1)\left\{-A_{1/\alpha }^1\left(\cosh \left(\kappa \right)\right)-A_{1/\alpha }^0\left(\cosh \left(\kappa \right)\right)\tanh (\kappa /2)\right\}}{\sqrt{{\xi }_1(\alpha ,\kappa )}},\hfill \\ \hfill r_{XS}:& =\mathrm{Corr}(X,\sin \mathsf{\Theta })\hfill \\ \hfill & =\frac{8\lambda {\cosh }^2(\kappa /2)\mathsf{\Gamma }\left(1/\alpha +1\right)}{\sqrt{{\xi }_1(\alpha ,\kappa ){\xi }_2(\kappa ,\lambda )}}[-\frac{3}{8}({\lambda }^2-4)A_{1/\alpha }^0\left(\cosh \left(\kappa \right)\right)\hfill \\ \hfill & +\frac{1}{2}({\lambda }^2-3)A_{1/\alpha }^2\left(\cosh \left(\kappa \right)\right)-\frac{{\lambda }^2}{8}{A}_{1/\alpha }^4\left(\cosh \left(\kappa \right)\right)\hfill \\ \hfill & -\frac{A_{1/\alpha }^0\left(\cosh \left(\kappa \right)\right)}{8{\cosh }^2(\kappa /2)}\left\{-3({\lambda }^2-4)+{\lambda }^2{\tanh }^2(\kappa /2)\right\}]\hfill \end{array}$$

(14)

and

$$\begin{array}{cc}\hfill r_{CS}:& =\mathrm{Corr}(\cos \mathsf{\Theta },\sin \mathsf{\Theta })\hfill \\ \hfill & =\frac{\lambda \sinh (\kappa /2)\left\{3({\lambda }^2-4)+{\lambda }^2/{\cosh }^2(\kappa /2)\right\}}{\sqrt{2}\sqrt{{\xi }_2(\kappa ,\lambda )}}.\hfill \end{array}$$

(15)

Then, the linear–circular squared correlation of X with $(\cos \mathsf{\Theta },\sin \mathsf{\Theta })$ that is defined in [14,15] can be obtained by

$$\begin{array}{cc}\hfill R_{X\mathsf{\Theta }}^2& =\frac{r_{XC}^2+r_{XS}^2-2r_{CS}{r}_{XC}{r}_{XS}}{1-r_{CS}^2}.\hfill \end{array}$$

(16)

It is seen from these expressions that $R_{X\mathsf{\Theta }}^2$ does not depend on the parameter $\beta$ and the sign of $\lambda$. In addition, if $\lambda =0$ then $R_{X\mathsf{\Theta }}^2=r_{XC}^2$ holds. It should be noted that when $\lambda =0$, the linear–circular squared correlation coefficients for any order m are identical.

Similar to the cylindrical distribution proposed by [8], the linear–circular squared correlation $R_{X\mathsf{\Theta }}^2$ depends only on $\alpha$, $\kappa$, and the absolute value $\left|\lambda \right|$. The behavior of $R_{X\mathsf{\Theta }}^2$ for $\alpha$ and $\kappa$ is shown by the contour plots in Figure 3. As depicted in Figure 3, when $\lambda$ and $\kappa$ are fixed, $R_{X\mathsf{\Theta }}^2$ increases as $\alpha$ increases. Concerning the effects of $\kappa$ on $R_{X\mathsf{\Theta }}^2$, the shape of the functions of $\kappa$ is humped (rising then falling) when $\alpha$ and $\lambda$ are fixed. These characteristics resemble those of the distribution of [8]. The linear–circular squared correlations for models with $m=2,3,4$, and beyond can be obtained by using numerical or Monte Carlo integration methods.

2.3. Random Number Generation

Combining the results from Section 2 with the random number generation method proposed by [8], which is based on the decomposition $f(x,\theta )=f\left(x\right|\theta \left)f\right(\theta )$, we introduce a simple random number generation algorithm. This involves initially generating $\mathsf{\Theta }\sim f\left(\theta \right)$ and then $X|\mathsf{\Theta }\sim f(x\left|\theta \right)$. The algorithm proceeds as follows.

Step 1

Generate a random variable ${\mathsf{\Theta }}_1$ from the wrapped Cauchy distribution with location $\mu$ and concentration $\tanh (\kappa /2)$ and generate independently U from the uniform distribution $U[0,1]$.

Step 2

Set $\mathsf{\Theta }$ as

$$\begin{array}{c}\hfill \mathsf{\Theta }=\left\{\begin{array}{c}{\mathsf{\Theta }}_1\mathrm{if}U<G_m\left(\lambda \sin ({\mathsf{\Theta }}_1-\mu )\right)\hfill \\ -{\mathsf{\Theta }}_1\mathrm{if}U\ge G_m\left(\lambda \sin ({\mathsf{\Theta }}_1-\mu )\right).\hfill \end{array}\right.\end{array}$$

Then, $\mathsf{\Theta }$ follows the sine-skewed wrapped Cauchy distribution.

Step 3

Generate X from the Weibull distribution (6) with a shape parameter $\alpha$ and a scale parameter $1/\left\{\beta {(1-\tanh \left(\kappa \right)\cos (\theta -\mu ))}^{1/\alpha }\right\}$.

2.4. Identifiability

Identifiability is an essential condition for giving a unique interpretation to estimators and for deriving the consistency of estimators. If parameter estimation is performed using a non-identifiable family of distributions, two different estimates will represent the same model, resulting in difficulty in interpreting the estimated parameters. If we consider a family of distributions (1) that includes the order m as a parameter, the family becomes non-identifiable with respect to m when $\lambda =0$. For this reason, it is necessary to discuss the identifiability of the subfamily ${\mathcal{F}}_{\mathrm{cyl}}^{\left(m\right)}:=\left\{f_{\mathrm{WeiESSVM}}^{\left(m\right)}(\theta ,x;\mathit{\eta })\right|\mathit{\eta }\in \mathit{H}\}$ for which the order m is fixed. Here, $\mathit{H}=\left\{\right(\mu ,\kappa ,\lambda ,\alpha ,\beta \left)\right|-\pi \le \mu <\pi ,\kappa >0,-1\le \lambda \le 1,\alpha >0,\mathrm{and}\beta >0\}$ denotes the parameter space. Then the following holds.

Proposition 1.

For each $m\in \{0,1,2,\dots \}$, the family ${\mathcal{F}}_{\mathit{cyl}}^{\left(m\right)}$ is identifiable.

The proof is given in Appendix A.

3. Statistical Inference

Assume that observations $({\theta }_1,x_1),\dots ,({\theta }_n,x_n)$ are sampled from the density (1). To estimate the parameter vector $\mathit{\eta }$, we employ the maximum likelihood estimator ${\hat{\mathit{\eta }}}_n\in {\mathrm{argmax}}_{\mathit{\eta }\in \mathit{H}}{\mathcal{l}}_m\left(\mathit{\eta }\right)$ based on the log-likelihood function

$$\begin{array}{cc}\hfill {\mathcal{l}}_m\left(\mathit{\eta }\right)& =(\alpha -1)\sum _{i=1}^n\log x_i-{\beta }^{\alpha }\sum _{i=1}^n{x}_i^{\alpha }\left\{1-\tanh \left(\kappa \right)\cos ({\theta }_i-\mu )\right\}\hfill \\ \hfill & +\sum _{i=1}^n\log G_m\left(\lambda \sin ({\theta }_i-\mu )\right)+n\left\{\alpha \log \beta +\log \alpha -\log \left(\pi \cosh \right(\kappa \left)\right)\right\}.\hfill \end{array}$$

After all parameters are converted to unconstrained variables using variable conversions, such as $\tilde{\kappa }=\log \kappa$, optimization is carried out by the BFGS quasi-Newton method in R’s ‘optim’ function. Because the density (1) has the hyperparameter m, we need to specify it. Here, we consider two methods for selecting order m. The first is the maximum log-likelihood approach, which is given by $m={\mathrm{argmax}}_{m=0,1,...,M}{\mathcal{l}}_m\left({\hat{\mathit{\eta }}}_n\right)$, where $M>0$ is a prespecified upper bound. The second is to find the order m that minimizes the Takeuchi information criterion (TIC), which is proposed by [16]:

$$\begin{array}{c}\hfill \mathrm{TIC}\left(m\right)=(-2){\mathcal{l}}_m\left({\hat{\mathit{\eta }}}_n\right)+2tr\left({\hat{J}}^{-1}\hat{K}\right),\end{array}$$

(17)

where

$$\begin{array}{c}\hfill \hat{J}=-\frac{1}{n}\frac{{\partial }^2}{\partial \mathit{\eta }\partial {\mathit{\eta }}^T}{\mathcal{l}}_m\left({\hat{\mathit{\eta }}}_n\right)\end{array}$$

and

$$\begin{array}{c}\hfill \hat{K}=\frac{1}{n}\sum _{i=1}^n{\displaystyle \frac{\partial }{\partial \mathit{\eta }}}\log f_{\mathrm{WeiESSVM}}^{\left(m\right)}({\theta }_i,x_i;{\hat{\mathit{\eta }}}_n){\displaystyle \frac{\partial }{\partial {\mathit{\eta }}^T}}\log f_{\mathrm{WeiESSVM}}^{\left(m\right)}({\theta }_i,x_i;{\hat{\mathit{\eta }}}_n).\end{array}$$

In the right-hand side of Equation (17), the first term $(-2){\mathcal{l}}_m\left({\hat{\mathit{\eta }}}_n\right)$ measures the goodness of fit of a model, and the second term $2tr({\hat{J}}^{-1}\hat{K})$ measures how much the model used for estimation deviates from the true model.

Monte Carlo Simulation Study

In this section, we demonstrate the model fitting performance of the WeiESSvM distribution (1). To investigate the finite sample performances of the MLE, we generate random samples from WeiESSVM distributions with two different orders $m=2$ and $m=3$. The data-generating processes (DGP) considered here are as follows. In the first DGP, hereafter referred to as DGP 1, the parameter values are chosen as $\alpha =2,\beta =0.5,\mu =\pi /6\approx 0.5236,\kappa =1,\lambda =-0.8$ with $m=2$ and for the second DGP (DGP 2), we set $\alpha =2,\beta =1,\mu =0,\kappa =1.5,\lambda =0.5$ with $m=3$. The random samples are generated according to the methods provided in Section 2.3. The sample size is selected as $n=100$ and $n=300$, and the Monte Carlo simulations were repeated 1000 times. We computed average bias and root mean squared error (RMSE). Since the true order $m_{\mathrm{true}}$ is unknown, we fit possible candidates for $m_{\mathrm{fitted}}$ that varies from 0 to 4, and the unknown order m is estimated by both maximized likelihood values and TIC, which are defined by

$$\begin{array}{c}\hfill {\hat{m}}^{\left(\mathrm{MLL}\right)}=\underset{m\in \{0,1,\dots ,4\}}{\mathrm{argmax}}{\mathcal{l}}_m\left({\hat{\mathit{\eta }}}_n\right)\mathrm{and}{\hat{m}}^{\left(\mathrm{TIC}\right)}=\underset{m\in \{0,1,\dots ,4\}}{\mathrm{argmin}}\mathrm{TIC}\left(m\right).\end{array}$$

Table 1. Simulation results of the mean bias and RMSE of the MLE for DGP 1. The true order is $m=2$. The RMSE values are reported in parentheses, and the best results among models are highlighted in boldface. True parameter values are $\alpha =2$, $\beta =0.5$, $\mu =\pi /6$, $\kappa =1$, $\lambda =-0.8$.

$\mathit{n}=100$	$\mathit{\alpha }$	$\mathit{\beta }$	$\mathit{\mu }$	$\mathit{\kappa }$	$\mathit{\lambda }$
$m=0$	0.0667	0.0235	$-0.0873$	0.0749	$-0.1895$
	(0.1721)	(0.0475)	(0.1157)	(1.0784)	(0.1915)
$m=1$	0.0486	0.0076	$-0.0202$	0.0254	$-0.1213$
	(0.1648)	(0.0404)	(0.0864)	(1.0603)	(0.1518)
$m=2$	0.0427	0.0031	–0.0026	0.0107	–0.0243
	(0.1634)	(0.0395)	(0.0895)	(1.0546)	(0.1266)
$m=3$	0.0405	0.0015	0.0033	0.0055	0.0584
	(0.1628)	(0.0393)	(0.0922)	(1.0523)	(0.1459)
$m=4$	0.0396	0.0010	0.0057	0.0034	0.1242
	(0.1627)	(0.0393)	(0.0935)	(1.0515)	(0.1818)
$n=300$
$m=0$	0.0373	0.0203	$-0.0806$	0.0652	$-0.1932$
	(0.0948)	(0.0318)	(0.0912)	(1.0409)	(0.1937)
$m=1$	0.0197	0.0046	$-0.0141$	0.0168	$-0.1230$
	(0.0888)	(0.0243)	(0.0499)	(1.0234)	(0.1378)
$m=2$	0.0154	0.0013	0.0003	0.0058	–0.0105
	(0.0877)	(0.0240)	(0.0508)	(1.0190)	(0.0774)
$m=3$	0.0142	0.0002	0.0038	0.0026	0.0793
	(0.0876)	(0.0239)	(0.0515)	(1.0179)	(0.1086)
$m=4$	0.0131	$-0.0003$	0.0060	0.0007	0.1465
	(0.0874)	(0.0239)	(0.0520)	(1.0168)	(0.1620)

and

Table 2. Simulation results of the mean bias and RMSE of the MLE for DGP 2. The true order is $m=3$. The RMSE values are reported in parentheses, and the best results among models are highlighted in boldface. True parameter values are $\alpha =2$, $\beta =1$, $\mu =0$, $\kappa =1.5$, $\lambda =0.5$.

$\mathit{n}=100$	$\mathit{\alpha }$	$\mathit{\beta }$	$\mathit{\mu }$	$\mathit{\kappa }$	$\mathit{\lambda }$
$m=0$	0.0459	0.0111	0.0173	0.0282	0.3881
	(0.1490)	(0.0874)	(0.0553)	(0.5640)	(0.4032)
$m=1$	0.0422	0.0033	0.0031	0.0170	0.2042
	(0.1483)	(0.0870)	(0.0580)	(0.5605)	(0.2475)
$m=2$	0.0412	0.0015	–0.0000	0.0142	0.0897
	(0.1481)	(0.0871)	(0.0595)	(0.5596)	(0.1587)
$m=3$	0.0407	0.0007	$-0.0013$	0.0131	0.0159
	(0.1480)	(0.0872)	(0.0601)	(0.5591)	(0.1214)
$m=4$	0.0404	0.0003	$-0.0021$	0.0124	$-0.0358$
	(0.1480)	(0.0872)	(0.0605)	(0.5589)	(0.1177)
$n=300$
$m=0$	0.0126	0.0106	0.0162	0.0156	0.3921
	(0.0799)	(0.0492)	(0.0338)	(0.5186)	(0.3981)
$m=1$	0.0096	0.0044	0.0044	0.0067	0.1930
	(0.0795)	(0.0482)	(0.0325)	(0.5157)	(0.2090)
$m=2$	0.0088	0.0028	0.0016	0.0045	0.0777
	(0.0795)	(0.0481)	(0.0328)	(0.5149)	(0.1066)
$m=3$	0.0084	0.0021	0.0004	0.0035	0.0045
	(0.0794)	(0.0481)	(0.0331)	(0.5145)	(0.0664)
$m=4$	0.0082	0.0018	–0.0003	0.0029	$-0.0467$
	(0.0794)	(0.0481)	(0.0333)	(0.5143)	(0.0767)

presented the results of the MLEs for DGP 1 and DGP 2, respectively. From these tables, we observed the following facts. First, when a fitted value of $m_{\mathrm{fitted}}$ is equal to the true value of $m_{\mathrm{true}}$, that is $m_{\mathrm{fitted}}=m_{\mathrm{true}}=2$ for DGP 1, and $m_{\mathrm{fitted}}=m_{\mathrm{true}}=3$ for DGP 2, it is observed that as n increases, the bias and the RMSE decrease, which indicates the consistency of the MLE. In contrast, even for the misspecified case with $m_{\mathrm{fitted}}\ne m_{\mathrm{true}}$, we can confirm that the bias and the RMSE of the parameters of the Weibull distribution, such as $\alpha$, $\beta$, and $\kappa$, decrease as n becomes large. Hence, it can be confirmed that the consistency of the parameters of the Weibull distributions holds, particularly when the order of fitted models is $m_{\mathrm{fitted}}=5$, the efficiency of the estimators is convincingly demonstrated. These facts can be explained by the structures of the Fisher information matrix of the log-likelihood function, as the function $\log G_m(·)$ is independent of parameters $\alpha$, $\beta$, and $\kappa$. For the case with underestimated $m_{\mathrm{fitted}}$, that is $m_{\mathrm{fitted}}<m_{\mathrm{true}}$ for DGP 1, the estimated $\hat{\lambda }$ tends to lie on the boundary of its parameter space, yielding negative bias in estimation, and vice versa for $m_{\mathrm{fitted}}>m_{\mathrm{true}}$ in DGP 2. Similarly, in the case with $m_{\mathrm{fitted}}>m_{\mathrm{true}}$ for DGP 1, the bias of the estimated parameter $\hat{\lambda }$ became positive, and vice versa for $m_{\mathrm{fitted}}<m_{\mathrm{true}}$ in DGP 2.

Figure 4 and Figure 5 show boxplots of the estimated skewness parameter $\hat{\lambda }$ for changing values of $m_{\mathrm{fitted}}$ for DGP1 and 2, respectively. It can be observed that for DGP 1, at least 3/4 and 1/2 of the $\hat{\lambda }$ lies on the boundary $-1$ when $m_{\mathrm{fitted}}=0$ for sample size $n=100$ and $n=300$, respectively. As for DGP2, at least 1/4 of the $\hat{\lambda }$ lies on the boundary 1 when $m_{\mathrm{fitted}}=0$ for $n=100$. This indicates that the existing sine-skewed Weibull-von Mises distribution corresponding to the case with $m=0$, could not sufficiently capture the skewness of the angular variables, and our extended sine-skew modeling overcomes this drawback. The unbiasedness property is observed with correctly specified case $m_{\mathrm{true}}=m_{\mathrm{fitted}}=2$, which indicates that the specifying unknown order m is crucial for estimating $\lambda$.

Table 3. Simulation results of the ratio of selected models by MLL and TIC for DGP 1. The true order is $m=2$.

	$\mathit{m}=0$	$\mathit{m}=1$	$\mathit{m}=2$	$\mathit{m}=3$	$\mathit{m}=4$
${\hat{m}}^{\left(\mathrm{MLL}\right)}$
$n=100$	0.177	0.324	0.122	0.068	0.309
$n=300$	0.072	0.392	0.139	0.058	0.339
${\hat{m}}^{\left(\mathrm{TIC}\right)}$
$n=100$	0.170	0.457	0.085	0.035	0.253
$n=300$	0.047	0.442	0.086	0.055	0.370

and

Table 4. Simulation results of the ratio of selected models by MLL and TIC for DGP 2. The true order is $m=3$.

	$\mathit{m}=0$	$\mathit{m}=1$	$\mathit{m}=2$	$\mathit{m}=3$	$\mathit{m}=4$
${\hat{m}}^{\left(\mathrm{MLL}\right)}$
$n=100$	0.426	0.103	0.033	0.011	0.427
$n=300$	0.340	0.125	0.045	0.014	0.476
${\hat{m}}^{\left(\mathrm{TIC}\right)}$
$n=100$	0.330	0.188	0.012	0.028	0.442
$n=300$	0.320	0.110	0.040	0.025	0.505

summarize the results of the selected orders by using ${\hat{m}}^{\left(\mathrm{MLL}\right)}$ and ${\hat{m}}^{\left(\mathrm{TIC}\right)}$. According to Table 3, both the maximized likelihood method and TIC select $m=1$ model for both sample sizes of 100 and 300. The ratio of selecting orders $m=1$ and $m=4$ increases as n increases, whereas the ratio of selecting true model $m=2$ remains low and stable for different model selection methods and sample sizes. The same findings are observed from Table 4, as the ratio of selecting orders $m=4$ increases whereas the ratio of the selecting true model $m=3$ remains low.

4. Real Data Analysis

4.1. Periwinkle Data Set

We consider the data set of size 31 with direction and distance blue periwinkles have moved after re-positioning them at a specific point on the ground. The data set is given in Table 1 of [17] and available in R packages NPCirc and circular.

Table 5. Estimated parameters with maximized log-likelihood function and TIC for periwinkle data. The maximum MLL and minimum TIC values are in boldface.

	$\hat{\mathit{\alpha }}$	$\hat{\mathit{\beta }}$	$\hat{\mathit{\mu }}$	$\hat{\mathit{\kappa }}$	$\hat{\mathit{\lambda }}$	MLL	TIC
$m=0$	2.011	0.050	1.242	1.682	1.000	$-168.466$	350.383
$m=1$	2.012	0.047	1.192	1.637	1.000	$-168.009$	344.615
$m=2$	2.003	0.046	1.157	1.596	1.000	$-167.871$	345.010
$m=3$	1.991	0.045	1.129	1.558	1.000	−167.838	345.251
$m=4$	1.988	0.044	1.121	1.547	0.936	$-167.845$	345.385

shows the maximum likelihood estimators, the log-likelihood functions evaluated at the maximum likelihood estimators, and the TIC values. The linear–circular association measure defined by (16) was $R_{X\mathsf{\Theta }}^2(\hat{\alpha },\hat{\kappa },\hat{\lambda })=0.263$ for $m=1$, whereas the sample estimate was ${\hat{R}}_{X\mathsf{\Theta }}^2=0.294$, which was calculated by substituting $r_{XC}^2$, $r_{XS}^2$, and $r_{CS}^2$ with their sample moments ${\hat{r}}_{XC}^2$, ${\hat{r}}_{XS}^2$, and ${\hat{r}}_{CS}^2$. This subtle difference between model and sample association measures comes from the drawbacks of the proposed models with fitted skewness parameters at the boundary $\hat{\lambda }=1$. The estimated linear–circular association measure $R_{X\mathsf{\Theta }}^2(\hat{\alpha },\hat{\kappa },\hat{\lambda })$ for models $m=2,3,4$ becomes $R_{X\mathsf{\Theta }}^2(\hat{\alpha },\hat{\kappa },\hat{\lambda };m=2)=0.292$, $R_{X\mathsf{\Theta }}^2(\hat{\alpha },\hat{\kappa },\hat{\lambda };m=3)=0.296$, $R_{X\mathsf{\Theta }}^2(\hat{\alpha },\hat{\kappa },\hat{\lambda };m=4)=0.297$, respectively. It can be seen that for orders m from 0 to 3, the estimates for the skewness parameter $\lambda$ appear on the boundary of the parameter space. However, as m increases, the estimates of $\lambda$ are away from the boundary, and become interior points of the parameter space.

The maximum log-likelihood method selects the model with $m=3$, whereas the TIC minimization method selects the model with $m=1$. If we want to conduct a hypothetical test on the asymmetry of the distribution (i.e., the null hypothesis is $\lambda =0$) using the asymptotic method, the model with $m=4$, in which the estimate of $\lambda$ is less than 1, is theoretically preferable.

Figure 6 shows a scatter plot of the data and a contour plot of the fitted WeiESSVM density of order $m=3$, which has the maximum log-likelihood (MLL). The vertical axis indicates the angle in radians, and the horizontal axis indicates the length.

Now we compare it with cylindrical models of [4,6,8]. The cylindrical density of [6] is given by

$$\begin{array}{c}\hfill f_{KS}(\theta ,x)\propto \exp \left[-\frac{{\{x-\mu \left(\theta \right)\}}^2}{2{\sigma }^2}+{\kappa }_1\cos (\theta -{\mu }_1)+{\kappa }_2\cos \left\{2(\theta -{\mu }_2)\right\}\right],\end{array}$$

where $-\pi \le \theta <\pi$, $x>0$, $\sigma >0$, ${\kappa }_1>0$, ${\kappa }_2>0$, $-\pi \le {\mu }_1<\pi$, $-\pi /2\le {\mu }_2<\pi /2$, $\mu \left(\theta \right)={\mu }^{\prime }+\lambda \cos (\theta -\nu )$, ${\mu }^{\prime }\in \mathbb{R}$, $\lambda >0$, and $-\pi \le \nu <\pi$.

It is reported by [8] that the MLLs of the fitted densities of [4,6] are $-176.88$ and $-168.46$, respectively. As shown in Table 5, all MLLs of the WeiESSVM model with orders from 1 to 4 are greater than these two values. Because the numbers of parameters in the models of [4,6] are 6 and 8, respectively, and greater than those of the WeiESSVM distribution of any order, the Akaike information criterion and Bayesian information criterion of the WeiESSVM model of any order are better than those of the two models. Table 5 demonstrates that all WeiESSVM models of order from 1 to 6 enhance the MLL value compared to the WeiSSVM distribution proposed in [8], which corresponds to the WeiESSVM model with $m=0$.

4.2. Occurrence Date of Typhoon and the Reciprocal of Minimum Pressure Data Set

We fit the bivariate data of the occurrence date of the typhoon and the reciprocal of minimum pressure for the second example. The dataset is sourced from the website “Digital Typhoon” [18] at the National Institute of Informatics. The data we used here are from 2014 to 2022 and the sample size is $n=231$. The occurrence date of the typhoon is transformed into $[0,2\pi )$ such that the number of days from 1 January to 31 December corresponds to 1 through 365 into radians by ${\theta }_i={\mathrm{days}}_i\times 2\pi /365$. Since the typhoon season is generally from May to October, we set 1 July at 0 radian so that the mode of the angular data should be located at the origin. The minimum pressure of typhoon data is transformed into reciprocals so that the smaller pressures correspond to larger values. This transformation is made by ${\tilde{x}}_i=1/\log \left({\mathrm{pressure}}_i\right)$ and $x_t=1000({\tilde{x}}_i-\min \left({\tilde{x}}_i\right)+c)$, where c is a small constant to avoid zero values; here, we set $c=0.0001$. The resulting bivariate cylindrical data are plotted in Figure 7 with contour plots with estimated best-fitting cylindrical density function whose parameters were summarized in

Table 6. Estimated parameters with maximized log-likelihood function and TIC for typhoon data. The maximum MLL and minimum TIC values are in boldface.

	$\hat{\mathit{\alpha }}$	$\hat{\mathit{\beta }}$	$\hat{\mathit{\mu }}$	$\hat{\mathit{\kappa }}$	$\hat{\mathit{\lambda }}$	MLL	TIC
$m=0$	1.3729	1.0213	2.3993	0.2353	$-0.8749$	$-554.45$	$1119.95$
$m=1$	1.3733	1.0082	2.4753	0.1930	$-0.6820$	$-554.11$	$\mathbf{1119}.\mathbf{29}$
$m=2$	1.3718	1.0053	2.4965	0.1798	$-0.5710$	$-554.06$	$1119.40$
$m=3$	1.3717	1.0037	2.5086	0.1736	$-0.5001$	$-554.04$	$1119.49$
$m=4$	1.3705	1.0030	2.5153	0.1691	$-0.4500$	$-\mathbf{554}.\mathbf{03}$	$1119.54$

. The estimated skewness order is $\hat{m}=1$ by TIC, whereas the MLL method selects $\hat{m}=4$, and the clear negative skewness is observed with moderately larger values of $\hat{\lambda }=-0.682$ with $\hat{m}=1$. Similar to the previous example, the Mardia–Sutton and Kato–Shimizu densities were estimated, which yielded the MLL as $-591.02$ and $-590.43$, respectively, indicating that our proposed model outperforms existing cylindrical distributions. The linear–circular association measure defined by (16) was $R_{X\mathsf{\Theta }}^2(\hat{\alpha },\hat{\kappa },\hat{\lambda })=0.017$ for $m=1$, whereas the sample estimate was ${\hat{R}}_{X\mathsf{\Theta }}^2=0.015$. This indicates that the model and sample association measures are quite close to each other. For models $m=0,2,3,4$, $R_{X\mathsf{\Theta }}^2(\hat{\alpha },\hat{\kappa },\hat{\lambda };m=0)=0.025$, $R_{X\mathsf{\Theta }}^2(\hat{\alpha },\hat{\kappa },\hat{\lambda };m=2)=0.015$, $R_{X\mathsf{\Theta }}^2(\hat{\alpha },\hat{\kappa },\hat{\lambda };m=3)=0.014$, and $R_{X\mathsf{\Theta }}^2(\hat{\alpha },\hat{\kappa },\hat{\lambda };m=4)=0.013$, respectively. This indicates that as m increases, the linear–circular association of the cylindrical data calculated by the estimated model parameters decreases. From this fact, the most adequate model, namely the model with $\hat{m}=1$, can explain the linear–circular association of the data, and this model is selected by a TIC.

5. Concluding Remarks

In this paper, we have proposed a tractable and flexible model for data on the cylinder. As an advantage, this model achieves a better fit by selecting an appropriate order m, even when dealing with a cylindrical dataset in which the circular part exhibits significant skewness. Consequently, if we incorporate this distribution as a component in finite mixture models or hidden Markov models, they can effectively capture local asymmetries in the observed data. Although [19] employs the hidden Markov model with the cylinder distribution proposed by [8] as a component, replacing it with the WeiESSVM distribution could enhance data fitting.

In the analysis of real data in Section 4, the Hessian matrix of minus the log-likelihood function of this model did not degenerate numerically. On the other hand, it is not clear at this point whether the Fisher information matrix, which is the population version of the Hessian matrix, is regular or not. Such a study would be well worthwhile since the regularity of the Fisher information matrix is an important factor in the asymptotic normality of the maximum likelihood estimator.