From Symmetry to Asymmetry on the Disc Manifold: Modeling of Marion Island Data

: The joint modeling of angular and linear observations is crucial as data of this nature are prevalent in multiple disciplines, for example the joint modeling of wind direction and another climatological variable such as wind speed or air temperature, the direction an animal moves and the distance moved, or wave direction and wave height. Hence, there is a need for developing ﬂexible distributions on the hyper-disc, which has support of the interior of the hyper-sphere, as it allows for modeling the combination of angular and linear observations. This paper addresses this need by developing ﬂexible distributions for the disc that have the ability to capture any inherent bimodality present in the data. A new class of bivariate distributions is proposed which has support on the unit disc in two dimensions that includes, as a special case, the existing Möbius distribution on the disc. This class is obtained by expressing the density function in a general form using a measurable function termed as generator. Special cases of this generator are considered to demonstrate the ﬂexibility. By applying a conformal mapping to the generator function a new Möbius distribution class emanates. This class of bivariate distributions on the disc is the ﬁrst to account for bimodality and skewness present in the data. The ﬂexible behavior of the proposed models in terms of bimodality and skewness is graphically demonstrated. Preliminary evidential analysis of the wind data observed at Marion Island reveals the absence of unimodality in the data. The ﬁt of the proposed models, which account for bimodality, to the Marion Island wind data were evaluated analytically and visually.


Introduction
The analytical description of wind climate is usually confined to the description of wind speed; however, the joint description of the directional and linear wind characteristics is also essential. This joint modeling of wind speed and direction is important for climatology and a variety of ocean engineering and ecological applications. Further emphasis on the necessity of bivariate modeling of wind speed and wind direction was raised by [1]. One approach to capture both the linear and directional variables is by means of a disc. There are very few distributions available with support on the unit disc. To the best of the author's knowledge, the only literature that investigates distributions on the hyper-disc are those by Jones [2,3] and Uesu, et al. [4]. In [2], the bivariate beta distribution was proposed for the disc which is a symmetric distribution. In [3], the Möbius distribution on the disc was proposed, which allowed for skewness. The multivariate extension of the Möbius distribution on the disc to the hyper-disc was then proposed by Uesu, et al. [4]. A shortcoming in the existing distributions on a disc is the ability to account for bimodality present in the data. There exists flexible distributions for cylindrical models as proposed in [5][6][7], but not the disc. The need for jointly modeling a linear and angular variable on a disc manifold stems from the cyclic nature of the angular variable which cannot be modeled the same as a linear variable. Many distributions have been proposed for modeling angular variables on a circle and other manifolds. For a full summary of these models, refer to the works of Mardia, et al. [8], Ley, et al. [9] and Pewsey, et al. [10].
In this paper, we propose a general framework for distributions which has support on the unit disc and includes the above-mentioned Möbius distribution on the disc [3] as a special case. This proposed class of the distributions accounts for both bimodality and skewness. We define this distribution by using a generator function and applying a conformal mapping, Möbius transformation, to obtain a new Möbius distribution class. This new class allows for bimodality when the generator function exhibits bimodal behavior with asymmetry introduced by the Möbius transformation.
An advantage of the proposed model is that it extends the shape characteristics of the model on the unit disc proposed by Jones [3] to account for more flexibility. This approach to modeling wind data differs from those in [1,[11][12][13][14][15], by means of the correlation structure and method for obtaining the joint distribution. The aforementioned models have been built by means of: (i) the Johnson-Wehrly (JW) [7] method; and (ii) a copula approach where the univariate models of best fit are chosen to obtain a joint bivariate distribution for the wind speed and wind direction. The proposed model in this paper has an embedded correlation structure present to jointly account for the wind speed and wind direction. It was pointed out in [16] that the joint modeling of wind speed and wind direction is more significant than a univariate analysis of the wind speed. The univariate distribution of the wind speed or wind direction is not considered since Soukissian, et al. [1] pointed out that the best univariate distribution for wind speed does not provide the best bivariate fit for the wind speed and wind direction. There has been much attention of wind power assessment worldwide, however, there has been no analysis of the wind speed and wind direction analysis at Marion Island. The proposed models was used to analyse the wind data for Marion Island. These models were compared to the only known model on the disc, which is a special case of the proposed class. The main contributions of this paper can be summarized as follows: • A new class of distributions on the disc for the joint modeling of a linear and angular variable is introduced. This class includes as a special case the Möbius distribution on the disc. The class can assume a parametric or semi-parametric model based on the choice of the practitioner.

•
In the semi-parametric model, a general class is introduced. For this model, we propose a computational algorithm which simplifies to a univariate kernel density function in order to obtain the parameter estimates.

•
The proposed class of distributions includes an embedded correlation structure present in the model to account for any relation between the linear and angular variables. This is addresses a shortcoming in the joint modeling of linear and angular variables. The embedded correlation structure reduces any assumption of the univariate models for the linear and angular variables.

•
The need for flexible models on a disc is addressed by introducing a class of distributions that can account for bimodality and skewness present in the data. An advantage of the proposed model is that it extends the shape characteristics of the model on the unit disc proposed by Jones [3] to account for more flexibility.

•
The modeling of the wind speed and wind direction data at Marion Island was jointly analyzed for the first time in this study.
This new class of distributions is defined in Section 3 and special cases of the generator function are highlighted in Section 4. In Section 5, an illustrative example of fitting the different models emanating from the class, to Marion Island wind data is presented. The numerical results along with the evaluation of the bivariate fits are also presented. Subsequently, a background discussion of the Marion Island weather data ia provided in Section 2. Section 6 concludes with a summary and some final remarks.

Background
The Prince Edward Island group (PEI) consist of two small volcanic islands: (1) Marion Island; and (2) Prince Edward Island. These islands are located in the sub-Antarctic Indian Ocean at latitude 46 • 49 -46 • 59 South and longitude 37 • 35 -37 • 55 East (see Figure 1a). Both island territories belong to South Africa, which is situated 1900 km to the northwest of the island group [17]. The recently declared PEIs Marine Protected Area (MPA) is one of the largest of its kind in the world and serves as a research area for South Africa as well as the world.
With an area of about 290 km 2 , a circumference of 72 km and the highest peak at 1300 m, Marion Island is the larger of the Prince Edward Islands. With its geographic isolation, Marion Island's climate is highly oceanic in nature, coupled with the influence of passing frontal weather systems. This creates conditions of high temperature variability and where gale force winds can blow during most weeks of the year. In 1948, South Africa established a fully operational weather station on Marion Island (see Figure 1b) with continuous weather data records available since the early 1950s [18]. Due to its duration and the scarcity of weather data from the sub-Antarctic region, these weather records provide valuable and unique opportunities for research.
Data from the Marion Island describe climate variability collected by a weather station that is located near the coast on the eastern side of the island [19]. Wind differences between the eastern (predominantly leeward winds) and western (windward) parts of the island, both within the larger scale prevalence of westerly winds, have long been speculated about [19][20][21]. It has been suggested that winds at the western side of the island, where no consistent data collection exists, are stronger and less variable, in comparison to winds at the eastern weather station. Despite these orographically induced differences, wind recordings of the weather should capture the larger scale variability of climate at Marion Island. The data used in this study were obtained from Marion Island and were recorded daily at 08:00, 14:00 and 20:00 SAST (relates to the main synoptic hours). The daily averages were obtained and the wind speed and wind direction investigated for nine different years. These nine years were considered in two-year intervals from 2001 to 2017. Figure 2 provides visual displays of the data to the reader with some background to the behavior of the data. A scatter plot of the wind direction and wind speed is given in Figure 3. In Table 1, the values of the main circular statistics [8] of the wind direction for the nine years are given. Specifically, the mean resultant length (r), mean direction (θ ), circular variance (V Θ ), circular standard deviation (ν Θ ), circular dispersion (δ), circular skewness (ŝ) and circular kurtosis (k) are shown.

Construction Methodology
In this section, we define the general Möbius class and consider an algorithm to estimate the generator function and the parameters using a semi-parametric approach.

Möbius Transformation
The Möbius transformation is defined on the extended complex plane denoted byC [22]. Hence, we rewrite the cartesian coordinates in polar coordinates. Considering the polar coordinate form of the unit disc, we write z = (x, y) and w = (u, v) as complex numbers with w = re iθ , 0 ≤ r ≤ 1, − π < θ ≤ π. The Möbius transformation is defined by [23] as where c = ae iµ , 0 ≤ a < 1 andc denotes the complex conjugate. The only conformal mapping from the unit disc to itself is the Möbius transformation [23]. Using this Möbius transformation, we can develop probability distributions [3,24,25].

Cartesian Coordinate Configuration of the General Class
Consider the general bivariate density function, where 0 ≤ x 2 + y 2 ≤ 1 and g(·) is a continuous function termed as the "generator". The normalizing constant C can be obtained from the polar form of Equation (1) as when the g(·) function is infinitely differentiable at the point zero. Numerical computation of the integral in Equation (2) is given by where E[·] denotes the expected value of a random variable. Writing the Taylor's series expansion of g (u) in Equation (2) around zero gives which allows the computation of the normalizing constant by derivatives of g(·) around zero.

A New Möbius Distribution Class on the Disc
By applying the Möbius transformation in Section (3.1), we obtain the class of general Möbius distributions on the disc as, where 0 ≤ a < 1, −π < µ ≤ π and 0 ≤ r ≤ 1, −π < θ ≤ π. We denote this class as (r, θ) ∼ GM (a, µ, g). Trivial calculations reveal for a strictly increasing generator function, These bounds are achieved when θ = µ(+π). For the limiting case a → 0, the joint distribution tends to Crg r 2 /π. We can rewrite the density function in Equation (4) in the Cartesian space by performing the transformation: u = r cos θ and v = r sin θ. Then, the resulting density function is where 0 ≤ u 2 + v 2 ≤ 1, 0 ≤ a < 1 and −π < µ ≤ π.

Maximum Likelihood Estimation
In this section, we derive the maximum likelihood (ML) estimates of the parameters.
, where all parameters are unknown. The likelihood function has the form Computational stages for estimating g(·) and the parameters are given by Algorithm 1.

Algorithm 1 ML estimation algorithm
Step 1. Set the initial values for the unknown parametersâ = 0 andμ = −π. These values are chosen as the lower bound as the algorithm performs an ascending sequential search through all possible parameter values.
Step 4. Repeat Steps 2 and 3, using an optimization technique to achieve convergence for the parameters a and µ. The final estimated values are those that maximize the likelihood function.
In Algorithm 1, the bandwidth h can be obtained using existing methods since the values in the argument of g(·) function belong to the Cartesian space. Hence, in our numerical studies, we applied the commonly used cross validation method to optimize h. Furthermore, we considered the standard Gaussian model for the kernel function, K(·). In Algorithm 1, the computation time can be improved by considering a grid approach in Step 4 over the domain of the two parameters.

Special Cases
In this section, we consider three special cases of the generator function, g(·). The first generator function, Pearson Type II distribution [26], results in the well known Möbius distribution on the disc proposed by Jones [3]. The second generator function, Kummer beta distribution [27], results in the Kummer beta Möbius distribution and the third generator function, beta type III distribution [28], results in the beta type III Möbius distribution. The Pearson Type II distribution was used for comparsion purposes and to show that the proposed new class of distributions include the well known distribution on the disc [3]. The other special cases are considered due to their flexibility and the capacity of the distribution to handle bimodality [27].

Möbius Distribution on the Disc
By defining the generator function as we obtain the bivariate spherically symmetric beta (Pearson Type II) distribution [26], with density function, where C = γ π > 0, 0 ≤ x 2 + y 2 ≤ 1. The density function and marginal distribution graphs in Figures 4 and 5 illustrate the symmetrical and unimodal behavior of the generator function in Equation (8).  After applying the Möbius transformation, this choice of generator function in Equation (8) simplifies to the well known Möbius distribution on the disc defined by Jones [3].
The parameters can be interpreted as γ controls the concentration, a controls the skewness (or the asymmetry for the length from the center of the disc) and µ controls the orientation of the distribution (see Figure 6).

Maximum Likelihood Estimation
The log-likelihood of the polar form density function, Equation (9), is obtained for maximum likelihood estimation as:

Kummer Beta Möbius Distribution on the Disc
Consider the following generator function where γ, β, λ > 0. By substituting the generator function in Equation (1), the resulting distribution is a bivariate Kummer beta distribution [27] with density function where γ, β, λ > 0, 0 ≤ x 2 + y 2 ≤ 1 and is the normalizing constant where 1 F 1 (·) is the confluent hypergeometric function of the first kind [29]. The density function and marginal distribution graphs in Figures 7 and 8 illustrate the symmetrical and bimodal behavior of the generator function.  After applying the inverse Möbius transformation in Section 3.1, the Kummer beta Möbius distribution is obtained from Equation (4) with the density function in terms of the polar coordinates as where 0 ≤ a < 1, −π < µ ≤ π, 0 ≤ r ≤ 1, −π < θ ≤ π and C given by Equation (13).
The parameters can be interpreted as γ controls the concentration, while β contributes to the concentration by adding bimodality. As the value of β increases, the model introduces another mode in the form of an antimode, as seen in the plots of the generator function in Figure 7. The parameter λ controls the steepness of the concentration, a controls the skewness (or the asymmetry for the length from the center of the disc) and µ controls the orientation of the distribution (see Figure 9).  (c) Density function and contour plots (16) for a = 0.2, µ = 0, γ = 2, β = 2 and λ = 0, 1, 5 (from left to right).

Maximum Likelihood Estimation
The log-likelihood of the polar form density function, Equation (15), is obtained for maximum likelihood estimation as:

Beta Type III Möbius Distribution on the Disc
Consider the following generator function where γ, β > 0. By substituting the generator in Equation (1), the resulting distribution is a bivariate Beta type III distribution [28] with density function where γ, β > 0, 0 ≤ x 2 + y 2 ≤ 1 and is the normalizing constant. The density function and marginal distribution graphs in Figures 10 and 11 illustrate the symmetrical and bimodal behavior of the generator function.  After applying the inverse Möbius transformation in Section (3.1), the Beta type III Möbius distribution is obtained from Equation (4) with the density function in terms of the polar coordinates as where 0 ≤ a < 1, −π < µ ≤ π, 0 ≤ r ≤ 1, −π < θ ≤ π and C as given in Equation (18).
The parameters can be interpreted as γ controls the steepness of the concentration, while β contributes to the concentration by adding bimodality. As the value of β increases (β > 1), the model introduces another mode in the form of an antimode, as seen in the generator function plots in Figure 10. The parameters a controls the skewness (or the asymmetry for the length from the center of the disc) and µ controls the orientation of the distribution (see Figure 12).

Maximum Likelihood Estimation
The log-likelihood of the polar form density function, Equation (19), is obtained for maximum likelihood estimation as:

Remark 2.
The marginal distributions for our proposed models are intractable. Obtaining the marginal distribution via numerical integration is computationally expensive.

Illustrative Example: Marion Island Data
For evaluation of the models, in addition to the AIC and BIC, we considered the performance measure proposed by Mathisen, et al. [30], which is used to evaluate the fits of bivariate models. This method used for linear-circular models in [1] is the normalized deviation. This method is based on the normalized deviation, d ij between the observed number of data points, N ij . The mathematical expression for the normalized deviation, d ij , is given by, where σ (e) ij is the normalizing factor and the expected standard deviation for the number of data points falling in the cell. Values of normalized deviation close to zero indicate a good fit of the estimated models. The normalized deviation identify the areas in the (r, θ) plane where the model either underperforms or overperforms. The normalized deviation, d ij , is presented in Figure 13, for the model of best fit for each year according to the performance measures in Table 2. In Table 3, the range and median of the d ij for the models of best fit for each year is given. The values if d ij close to zero indicate a good fit of the model. Positive and negative values indicate underestimation and overestimation, respectively. In Table 4, a summary of the normalized deviations for the model of best  fit for each year, based on Table 2, is provided. Figure 13 illustrates the normalized deviation plots for the model of best fit for each year.
In Table 2, the estimated values of the parameters of the models (Equations (4), (9), (15) and (19)) and their performances are summarized. The expressions for the maximum likelihood estimation for the parametric models are given in Appendix A.
The following inferences can be drawn from Table 2: 1.
The general Möbius distribution (Equation (4)) performs fairly well most of the time and has the advantage of no distributional assumptions. However, it is worth noting that the selection of the bandwidth plays an important role in the directional analysis.

2.
The Möbius distribution on the disc (Equation (9) The Beta type III Möbius distribution in Equation (19) outperforms when there is bimodal behavior inherent in the data. This can be seen for the years 2001, 2009 and 2011. It is also worth noting that the Beta type III Möbius distribution has the ability to model unimodality as well as bimodality. Hence, the performance of this model fairs well against the Möbius distribution on the disc.

4.
It is interesting to note that the parameters introduced by the Möbius transformation, a and µ, have the same estimated values for the three parametric models. This highlights the influence of the Möbius transformation to the generator models.  For the year 2017, consider the estimated density function plot of the Möbius distribution on the disc in Equation (9) in Figure 14. In Table 4, the range and median of the normalized deviation for the proposed parametric models are compared and the plots of these normalized deviation are illustrated in Figure 15.  Table 4. Normalized deviation ranges for the parametric models for 2017.

Model Range Median
Möbius distribution on disc (9)

Conclusions
In this paper, we have developed a new general Möbius distribution class on the disc which has support on the unit disc that includes the well-known Möbius distribution on the disc proposed by [3] as a special case. This new class can assume either a semi-parametric or parametric form depending on the choice of the practitioner. This class of distributions is suitable when there exist bimodal patterns hidden in the data structure. Using the Marion Island climate data, model evaluation showed the performance of the different proposed models and the flexibility of this new class of distributions in capturing the bimodality and asymmetry present in the data. The Möbius transformation introduces asymmetry to the model in terms of off-centeredness as well as an orientation for the asymmetry. The ability of this new class to capture the bimodal behavior is a result of the generator function. One of the reviewers commented on that the low kurtosis values in Table 1 indicate a multimodal distribution. This needs further exploration. This work can be extended for modeling wind power density in wind energy analysis. To this end, it is sufficient to multiply the wind power density with the general Möbius density given by Equation (4) (for further details, refer to [12]).

Acknowledgments:
The authors would like to thank the anonymous reviewers for their insightful comments which led to an improvement of this paper. The authors would also like to thank Elsa de Jager from the National Climate Center, South African Weather Service for her assistance with the data description.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A Appendix A.1. Möbius Distribution on This Disc
For the Möbius distribution on the disc, assuming that a = 0, the maximum likelihood estimators can be obtained by estimating the equations below derived from Equation (12),

Appendix A.3. Beta Type III Möbius Distribution
For the beta type III Möbius distribution, assuming that a = 0, the maximum likelihood estimators can be obtained by estimating the equations below derived from Equation (22), where Ψ(·) is the digamma function.