Directional data analysis using the spherical Cauchy and the Poisson-kernel based distribution
Michail Tsagris


Federated Prototype-based Contrastive Learning for Privacy-Preserving Cross-domain Recommendation
    Li Wang, Quangui Zhang, Lei Sang, Qiang Wu, Senior Member, IEEE, and Min Xu^*, IEEE, Member


Li Wang, Qiang Wu, and Min Xu are with the School of Electrical and Data Engineering, University of Technology Sydney, Sydney 2000, Australia. Shoujin Wang is with the institute of data science, University of Technology Sydney, Sydney 2000, Australia.  Quangui Zhang is with the School of Artificial Intelligence, Chongqing University of Arts and Sciences, Chongqing 402160, China. *Corresponding author: Min Xu (e-mail: Min.Xu@uts.edu.au)
    September 9, 2024
=========================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================






Abstract


The spherical Cauchy distribution and the Poisson-kernel based distribution were both proposed in 2020, for the analysis of directional data. The paper explores both of them under various frameworks. Alternative parametrizations that offer numerical and estimation advantages, including a straightforward Newton-Raphson algorithm to estimate the parameters are suggested, which further facilitate a more straightforward formulation under the regression setting. A two-sample location test, based on the log-likelihood ratio test is suggested, completing with discriminant analysis. The two distributions are put to the test-bed for all aforementioned cases, through simulation studies and via real data examples comparing and illustrating their performance. 



Keywords: Directional data, maximum likelihood, regression, discriminant analysis




MSC: 62H11, 62H30 



§ INTRODUCTION

Directional data refers to multivariate data with a unit norm, and its sample space can be expressed as:


    𝕊^d={ x∈ℝ^d+1||| x||=1 },

where ||.|| denotes the Euclidean norm. Circular data, when d=1, lie on a circle, whereas spherical data, when d=2, lie on a sphere.
Circular data are met in various disciplines, such as political sciences <cit.>, criminology <cit.>, biology <cit.>, ecology <cit.> and astronomy <cit.> to name a few. Spherical data on the other hand are met in geology <cit.>, environmental sciences <cit.>, image analysis <cit.>, robotics <cit.> and space <cit.>.

Numerous spherical and hyper-spherical distributions have been proposed over time, with the von Mises-Fisher <cit.> and projected normal <cit.> distributions being among the earliest and most prevalent. The spherical Cauchy (SC) <cit.> and the Poisson-kernel based PKB) <cit.> distributions are two recently propositions. Despite these distributions assume rotational symmetry, which may restrict their applicability in certain scenarios, they have proved useful in many situations and for data on the sphere, they seem to perform well on some occasions <cit.>.

In this paper the SC <cit.> and the PKB <cit.> distributions are investigated with regards to five aspects. These are random vectors simulation, maximum likelihood estimation, hypothesis testing about locations, regression modelling and discriminant analysis. Most of the existing distributions cover these cases, with some drawbacks whatsoever, such lack of computational efficiency and lack of available techniques, for instance lack of proper regression models.  

Regarding simulation of random vectors, many distributions rely on rejection sampling, such as the von Mises-Fisher (vMF) <cit.> and Kent <cit.> distributions, while others, e.g. the projected normal, the elliptically symmetric angular Gaussian (ESAG) <cit.> and the spherical projected Cauchy <cit.> avoid this. As will be showed later, the SC is straightforward to simulate from, while the PKB requires rejection sampling which can hamper its use in cases where fast generation of random vectors is required. 

Efficient maximum likelihood estimation (MLE) is crucial for many reasons, not only for simulation study purposes but also when analysis of large scale data is involved. Estimate of mean direction of the VMF is available in closed form, and its concentration parameter can be estimated via a fast Newton-Raphson (NR) algorithm <cit.>. Estimation of the parameters of the Kent and the projected type of distributions is harder and thus one relies numerical optimizers which typically are slow, especially with increasing dimensions. As will be shown later, MLE with the SC relies on NR, whereas MLE with the PKB relies upon a numerical optimizer. Closely related to the MLE is the calculation of the normalizing constant, which in the cases under study exists in a closed form. 

A plethora of hypothesis testing procedures for two or more mean directions has been proposed over the years. <cit.> performed a comparison of some methods for the case of two populations, where the computation of the p-value took place using both asymptotic theory and computer intensive techniques. The tests that did not allow for heterogeneity among the samples, did not perform accurately, i.e. they did not retain the nominal type I error, whereas tests that are independent of this assumption were more precise. The perk of the proposed log-likelihood ratio tests, based on the SC and the PKB distribution, do not assume equality of the concentration parameters. 

<cit.> related the mean direction of the projected normal, on the circle, to some covariates without assuming equality of the concentration parameter among the errors (the analogue of the homoscedasticity in linear regression). They allowed the concentration parameter to vary among the residuals. This concept was also used by <cit.> with the ESAG, Kent, and the vMF distribution, and by <cit.> for the projected Cauchy on the circle and the sphere. The same strategy is also proposed here with the two competing distributions. 

Finally, the case of discriminant analysis, a.k.a supervised learning, is another interesting topic in the field. <cit.> performed a comparison of maximum likelihood based classifiers and the k-NN algorithm and provided evidence that rotationally symmetric distributions sometimes perform equally well or better than distributions that are elliptically symmetric. As shown later, the SC and PKB perform equally well under the discrimination setting, however, discrimination based on the SC is much faster.

The SC and PKB distributions are presented in the next section, along with the five aforementioned aspects, for which new methods are proposed and discussed in detail. Section <ref> compares the performance of these two distributions, while Section <ref> illustrates their performance using real data. Finally, Section <ref> concludes the paper. 









§ THE SPHERICAL CAUCHY AND POISSON-KERNEL BASED DISTRIBUTIONS
 
A model that appears to be closely related to the classical vMF distribution is the SC distribution <cit.>, which can be seen as the generalisation of the wrapped Cauchy distribution <cit.> to the sphere (and hyper-sphere). 

The density of the SC on 𝕊^d is given by

    f( y)=C_d(1-ρ^2/1+ρ^2-2 y^⊤ m)^d,

where m∈𝕊^d is the location direction, that controls the mode of the density, ρ∈ [0, 1) plays the role of the concentration parameter and C_d=Γ[(d+1)/2]/2π^(d+1)/2 is the normalizing constant. For the most part of this paper though we will use an alternative parameterization that was used in <cit.> 

    f( y) =C_d(√(μ^2+1)- y^⊤μ)^-d=
    C_d(√(γ^2+1)-α)^-d,

where μ∈ℝ^d+1 is the unconstrained location parameter, α= y^⊤μ, γ=μ, m=μ/γ and ρ=(√(γ^2+1)-1)/γ (γ=2ρ/1-ρ^2). The benefit of this parameterization is that the maximisation with respect to the location parameter is unconstrained. 

A second, similar, distribution is the PKB distribution that was proposed by <cit.>, and can also be seen as the the generalisation of the wrapped Cauchy distribution, whose density is given by

    f( y)=C_d1-ρ^2/(1+ρ^2-2 y^⊤ m)^(d+1)/2,

where m, according to <cit.> is a vector orienting the center of the distribution. This sentence does not sound accurate and any reference to this parameter of the PKB, similarly to the SC, will be the location parameter. 

One may express this distribution in an alternative way, similar to the SC, as a function of the unconstrained location parameter μ  

    f( y)=C_d(√(γ^2+1)-α)^-(d+1)/2(1-(√(γ^2+1)-1)^2/γ^2)^-(d-1)/2.


Figure <ref> presents the contour plots for some location parameter and two values of the concentration parameter ρ for the SC distribution. The contour plots have the same shape, but the tails of the PKB decay slower than those of the SC distribution. The density value of the SC will be larger than the density value of the PKB, that is, if y^⊤μ<(1 + ρ^2)/2. 





 §.§ Simulation of directional vectors

In order to simulate random directional vectors y_i from the SC( y; m,ρ), <cit.> proposed a simple procedure that relies upon the uniform distribution and then perform some simple calculations. The fact that no rejection sampling is necessary, a strategy common with other distributions, is an appealing one. The two steps of the algorithm are described below.


  * Generate vectors u_i (i=…,n) from the uniform distribution in 𝕊^d.

  * Set y_i=( u_i + ρ m)(1-ρ^2∑_j=1^d+1 m_j^2)/∑_j=1^d+1( u_i + ρ m_j )^2 + ρ m.  


In order to simulate random directional vectors y_i from the PKB( y; m,ρ), <cit.> proposed a rejection sampling that is computationally more expensive than the simulation of the SC distribution. 


  Step 1 Set λ = 2 ρ / (1 + ρ^2).

  Step 2 Define  ω_d(λ, β) = 0.5 (d+1) log1 + √(1 - λ^2)/1 + √(1 - λ^2 / β) - 0.5 log(1-β) and find the β^* that minimizes ω_d(λ, β) in the interval (λ(2λ-1), 1).

  Step 4 Set β_1 = β^* / (1 - β^*) and β_2 =  - 1 + 1 / √((1 - β^*))

  Step 5 Simulate u ∼ U(0,1) and z=(z_1,…,z_d+1), where each z_i is simulated from a standard normal distribution, z_i ∼ N(0,1).

  Step 6 Set q=(μ^⊤ z + β_2μ^⊤ z)/
√( z^⊤ z+β_1(μ^⊤ z)^2).

  Step 7 If log(u)≤d+1/2[ - log(1 - λ * q) + log(1 - β^* q^2)-log2/1 + √(1 - λ^2/β^*)] 
 
set x←( z + β_2μ^⊤ zμ)/
√( z^⊤ z+β_1(μ^⊤ z)^2), otherwise return to Step 5.




 §.§ Maximum likelihood estimation
 
The log-likelihood for a sample of directional vectors y_i, i=1…,n of the SC distribution, using Eq. (<ref>) is given by 

    ℓ_SC = nlogC_d - d∑_i=1^nlog(√(γ^2+1)-α_i),

where C_d denotes the normalizing constant. To perform maximum likelihood estimation (MLE) the Newton-Raphson (NR) algorithm can be employed, to maximise ℓ (<ref>). The starting value for the NR algorithm is the sample mean vector for which the log-likelihood value is computed. At each successive step the estimated mean vector is updated via μ^t+1=μ^t+1 -  H^-1 J and the algorithm terminates when the difference between two successive log-likelihood values is smaller than some tolerance value ϵ (for instance ϵ=10^-6). 

A drawback of this method is that the concentration parameter ρ is embedded within the estimation of the mean vector μ. To deal with this problem an alternative optimization strategy is employed to disentangle the mean direction m from the concentration parameter. The method is a hybrid of the Brent algorithm and of the fixed points iteration algorithm. The relevant log-likelihood (excluding the normalizing constant) of the parameterization in Eq. (<ref>) is given by

    ℓ_SC = nlogC_d + ndlog(1 - ρ^2) - d∑_i=1^nlog(1+ ρ^2 - 2ρ y_i^⊤ m).

The steps of the hybrid algorithm are delineated below.


  Step 1 Start with an initial mean direction given by m̂=y̅/y̅, where y̅ denotes the sample mean vector.

  Step 2 Using this mean direction obtain ρ̂ that maximises the log-likelihood in Eq. (<ref>), using the Brent algorithm <cit.>.

  Step 3 For the estimated ρ̂ from Step 2, update the mean direction, using the fixed points iteration algorithm, by maximising the log-likelihood in Eq. (<ref>) under the constraint that the mean direction lies in 𝕊^d-1. The Lagrangian function takes the following form

    ℓ_SC = nlogC_d + ndlog(1 - ρ̂^2) - d∑_i=1^nlog(1+ ρ̂^2 - 2ρ̂ y_i^⊤ m) + λ( m^⊤ m - 1 ).

Equating the derivative of (<ref>), with respect to m, to zero, yields

    ∂ℓ_SC/∂ m= d∑_i=1^n2ρ y_i/1+ ρ̂^2 - 2ρ̂ y_i^⊤ m +2λ m =  0.

The updated mean direction is given by the unit vector parallel to ∑_i=1^nρ̂ y_i/1+ ρ̂^2 - 2ρ̂ y_i^⊤m̂.

  Step 4 Repeat Steps 2-3 until the log-likelihood in Eq. (<ref>) improves no more than some tolerance value ϵ.

The strategy employed in Step 3 was also employed by <cit.> in order to estimate the median direction and by <cit.> to obtain the eigenvectors of the Cauchy principal component analysis.  

As for the PKB distribution, the second representation (<ref>) eases the calculations, since the log-likelihood can be written as

    ℓ_PKB = nlog(C_d) -d+1/2∑_i=1^n log(√(γ^2+1)-α_i) - nd-1/2[log2 + log(√(γ^2+1)-1) - logγ^2].

The derivatives are very similar to those of the SC log-likelihood, with the exception of some extra terms. As for the hybrid algorithm, the mathematics are nearly the same in the SC case.



 §.§ Log-likelihood ratio test for equality of two location parameters

In order to test the equality of equal population location parameters, based on two samples[Evidently, the procedure can be generalised to the case of more than two groups.] that follow the SC distribution we will employ the log-likelihood ratio test, without assuming equality of the concentration parameters. Under the null hypothesis one must maximise the following log-likelihood (ignoring C_d) with respect to the common mean direction m_c and the two concentration parameters ρ_1 and ρ_2. 

    ℓ_0  = n_1dlog(1 - ρ_1^2) - d∑_i=1^n_1log(1+ ρ_1^2 - 2ρ_1  y_1i^⊤ m_c ) + n_2dlog(1 - ρ_2^2) - d∑_i=1^n_2log(1+ ρ_2^2 - 2ρ_2  y_2i^⊤ m_c ),

where n_j and y_ji refer to the j-th sample size and the i-th observation of sample j, respectively, for j=1,2. Under the alternative hypothesis, the two location parameters are not equal and hence the relevant log-likelihood to be maximised is given by

    ℓ_1  = n_1dlog(1 - ρ_1^2) - d∑_i=1^n_1log(1+ ρ_1^2 - 2ρ_1  y_1i^⊤ m_1 ) + n_2dlog(1 - ρ_2^2) - d∑_i=1^n_2log(1+ ρ_2^2 - 2ρ_2  y_2i^⊤ m_2 ),

where m_j refers to the location parameter of the j-th sample. In order to maximise ℓ_0 (<ref>) the hybrid maximisation approach is used, whereas for the maximisation of ℓ_1 (<ref>) this is accomplished via the NR algorithm, applied to each sample separately. If the null hypothesis is true, standard likelihood theory states that Λ = 2[ℓ_1(m̂_1, m̂_2, ρ̂_1,ρ̂_2)-ℓ_0(m̃_c, ρ̃_1,ρ̃_2)] χ^2_d, where m̂_1, m̂_2, ρ̂_1,ρ̂_2 denote the estimated parameters under H_1, while m̃_c, ρ̃_1,ρ̃_2 denote the estimated parameters under H_0.

The same strategy was adopted for the case of the PKB distribution as well, with the exception that the formulas are slightly different.



 §.§ Regression analysis

When a set of p covariates is present, we link the location parameter to the covariates by μ_i=  B^⊤ x_i, where B=(β_1^⊤,…,β_p^⊤) denotes the matrix of the regression coefficients and x_i denotes a row of the design matrix. The relevant log-likelihood, using (<ref>), becomes 

    ℓ_SC = log(C_d) - d∑_i=1^nlog(√(γ_i^2+1)-α_i),

where α_i= y_i^⊤μ_i and γ_i=μ_i. 

The advantage of the parametrization in Eq. (<ref>) is evident in the regression setting. Similarly to <cit.> and <cit.> the errors are anisotropic, because we do not assume a common concentration parameter. Each directional vector has its own concentration parameter γ_i that is linked to ρ_i as mentioned earlier. The matrix of the regression coefficients is again estimated via the NR algorithm. 

The log-likelihood of the PKB regression model, using the parameterization from Eq. (<ref>) is written as follows

    ℓ_PKB=nlog(C_d) - nd-1/2log2 -d+1/2∑_i=1^n log(√(γ_i^2+1)-α_i) -d-1/2[∑_i=1^nlog(√(γ_i^2+1)-1) - ∑_i=1^nlogγ_i^2],

where α_i and γ_i are the same as in the case of the SC regression. 



 §.§ Discriminant analysis

Under the maximum likelihood discriminant analysis framework the rule is to allocate a new observation vector z∈𝕊^d in the group whose log-likelihood value has the highest value. In the case of the SC, with two groups, z is allocated to to group 1 iff

    log√(μ_2^2+1)- z^⊤μ_2/√(μ_1^2+1)- z^⊤μ_1 > 0

and to group 2 otherwise, where μ_1 and μ_2 denote the location parameter of the first and second group, respectively. Hence, at first the location parameter is estimated for each group separately and then the allocation rule is applied to the new observation z. The case of the PKB distribution the allocation rule is straightforward to write down analytically, and evidently it has a more complicated form

    log√(μ_2^2+1)- z^⊤μ_2/√(μ_1^2+1)- z^⊤μ_1 + log√(μ_2^2+1)-1/√(μ_1^2+1)-1 - logμ_2^2/μ_1^2 > 0.






§ SIMULATION STUDIES
 
Simulation studies were conducted to examine the computational efficiency of the two algorithms used for MLE of the SC distribution. Additionally, the SC distribution was compared to the PKB distribution using spherical and hyper-spherical data. The comparisons between the two distributions were carried out under the previously discussed topics: equality of two location parameters,  regression and discriminant analysis settings.



 §.§ Computational efficiency of the MLE algorithms

We compare the runtime of the two MLE algorithms, namely the hydrid and the NR described earlier, under a combination of various sample sizes and dimensionalities, when the data have been generetaed from the SC or the PKB distribution. The speed-up factor of NR compared to the hybrid algorithm appears in Table <ref>, where evidently the NR is to be preferred to its opponent, especially with increasing sample sizes. However, with increasing dimensions the NR seems to have a reduced dominance, something which does not come by surprise, since the NR requires the inverse of the Hessian matrix. However, we will note that with circular data (d=1), Kent and Tyler's algorithm <cit.> is faster than the NR algorithm[There is no reason to compare the NR algorithm to R's built in function optim(), and hence we skip all other comparisons.]. 

We also compared the time required to fit the SC distribution and the time required to fit the PKB distribution, when the data are generated either from the SC or the PKB distribution. Table <ref> shows that the MLE of the SC is faster than the MLE of the PKB distribution, for both cases. In this example, the data were generated from the PKB distribution, but the results are similar even if the data were generated from the SC distribution. 







 §.§ Hypothesis testing for two location parameters

According to the simulation studies of <cit.>, conducted for circular and spherical data, the heterogeneous approach <cit.>, that does not assume equality among the concentration parameters, was shown to be the optimal test in terms of size attainment. We implemented a smaller scale simulation study to estimate the type I error and the power of the SC log-likelihood ratio test, but we cannot compare it to the heterogeneous approach because the latter compares mean directions. Table <ref> presents the estimated type I error and the estimated power of the tests, for both tests when the data are generated from the SC or the PKB distribution, at various dimensions. 





 §.§ Regression analysis

For the regression analysis, we adhered to the methodology outlined by <cit.>. We utilized a single covariate to generate n values from a standard normal distribution, linking it to the response directional variable in a linear manner, μ_i= X_i B, for i=1,…,n. Subsequently, data were generated from either the SC or the PKB distribution, with the median and mean directions, respectively, being m_i=μ_i/μ_i. For the SC distribution, the concentration parameter was ρ_i=(√(μ_i^2 + 1) - 1 ) / μ_i, whereas for the PKB distribution, it was κ_i=μ_i <cit.>.

This procedure was iterated for various combinations of sample sizes and dimensionalities, with regression coefficients being estimated using both the SC and PKB regression models. The entire process was repeated 1,000 times, and the average fit of the two models, as measured by the quantity ∑_i=1^n y_i^⊤ŷ_i/n is presented in Table <ref> The fit takes values from 0 up to 1, where higher values indicate better fit. The estimated fits of the regression models are nearly nearly identical.





 §.§ Discriminant analysis

Following <cit.> we also simulated data from the SC and the PKB distributions, assuming two groups whose mean directions differ by a specific angle and performed a 10-fold cross-validation protocol to estimate the percentage of correct classification[A key note difference from <cit.> is that we moved on to higher dimensions.]. This process was repeated 1,000 times and the average percentages are presented in Table <ref>. Evidently, regardless of the true data generation mechanism, the classification capabilities of either distribution are nearly the same. 







§ REAL DATA ANALYSIS
 



 §.§ Hypothesis testing for two location parameters

The Ordovician dataset <cit.> consists of two groups of 50 measurements, on the sphere, each from L_0^1 axes (intersections between cleavage and bedding planes of F, folds) in Ordovician turbidites, collected in the same sub-domain. Figure <ref> presents the data on the sphere with colours indicating the two groups and Table <ref> presents the MLE for the SC and PKB distributions. Evidently, there are small differences between the two models. Both the SC and PKBD log-likelihood ratio tests though provided high p-values, 0.733 and 0.856, respectively, and so were their relevant bootstrap based p-values, 0.825 and 0.914 for the SC and PKB, respectively. 








 §.§ Regression analysis

Data regarding crop productivity in the Greek NUTS II region of Thessaly during the 2017-218 cropping year were supplied by the Greek Ministry of Agriculture, also known as farm accountancy data network (FADN) data. The data refer to a sample of 487 farms and initially they consisted of 20 crops, but after aggregation they were narrowed down to 10 crops[A larger version of this dataset was used in <cit.>.]. For each of the 487 farms the cultivated area and the production in each of the 10 crops is known. However, the goal of the paper is to relate the composition of the production (simplicial response, Y) to the composition of the cultivated area (simplicial predictor, X) and for this reason were scaled to sum to unity[The raw data cannot be distributed due to disclosure restrictions.]. The square root was then applied to the compositional data (both the response and the predictor variables) so that they are mapped onto the 10-dimensional sphere.  

Including the constant terms, there are 110 regression parameters to be estimated. The SC based regression model required 0.11 seconds to complete, whereas the PKB regression model face numerical instability problems with the Hessian matrix and therefore we had to rely upon a numerical optimizer (the function optim() in R), which required more than 5 minutes to complete. This does not come by surprise as the such optimizers are not designed to work with high dimensional problems.

With regards to the fitting performance, the quantity ∑_i=1^n y_i^⊤ŷ_i/n was equal to 0.958 for the SC regression model and 0.955 for the PKB regression model. Due to the vast computational time required by the PKB regression model we did not perform the 10-fold cross-validation procedure. 



 §.§ Discriminant analysis

For this task we will consider the Wireless Indoor Localization data set, that is publicly available in the https://archive.ics.uci.edu/dataset/422/wireless+indoor+localizationUCI Machine Learning Repository’s website. The data were collected in indoor space by observing signal strengths of seven WiFi signals visible on a smartphone. The data consist of 2,000 measurements on 7 variables that report the measurements of the WiFi signal strength received from 7 Wi-Fi routers in an office location in Pittsburgh (USA). The grouping variable is one of the four rooms with 500 observations from each room. The WiFi signal strength is measured in dBm, decibel milliwatts, which is expressed as a negative value ranging from -100 to 0. In order to apply the discriminant analysis we first normalized the data (projected them onto the hyper-sphere).

We applied the 10-fold cross-validation process, repeated 50 times to quantify the variance caused by the different splits, computing the percentage of correct classification at each repetition. 
Figure <ref> presents the results. The average and median percentage of correct classification were equal to 0.9792 and 0.9790, respectively, for the SC distribution and equal to 0.9775 and 0.9775, respectively, for the PKB distribution.   





§ CONCLUSIONS
 
We investigated two recently proposed distributions, namely the spherical Cauchy and the Poisson-kernel based distribution. Specifically, we suggested an alternative, hybrid, method to estimate the parameters of both distributions and in particularly for the SC distribution we suggested a new parameterization that allows for application of the Newton-Raphson algorithm. We also re-parameterised the density function of the PKB distribution, in the same manner as with the SC, but it did not prove useful for the maximum likelihood estimation of its parameters. However, the new parameterizations facilitated the implementation of regression modelling escaping the (hyper-)spherical constraint and allowing for anisotropic errors as in the case of the ESAG distribution <cit.>. The benefit of the re-parameterized SC distribution is that Newton-Raphson was again implemented for the regression setting, thus yielding computationally efficient estimation of the regression parameters. The aforementioned hybrid method for estimation of the parameters of either distribution enabled the hypothesis testing of location directions between two populations, based on the log-likelihood ratio test, without assuming equal concentration parameters. Finally, we explored the maximum likelihood discriminant analysis using either distribution.    

The simulation studies and the real data examples, showcased the performance of each distribution, and provided evidence that both distributions perform similarly. The fact that the SC distribution is easier to simulate values from, is faster when it comes to estimating its parameters, with and without covariates, renders it a better choice, for practitioners and researchers.  

Regarding future work we can mention the following. Extension of the hypothesis testing for more than two location parameters is straightforward, the only difficulty is the running time, since as the number of groups increase, so does the computational cost. Secondly, rejection sampling for the PKB using the SC as an envelope function did not work very satisfactorily. We matched the parameters of the SC distribution to those of the PKB distribution and estimated the bound between the ratio of the two distributions. However, the bound increases with increasing dimensionality, plus the accuracy of this method seems less than the rejection sampling already proposed by <cit.>. Model based clustering though, using either distribution is something we are currently working on. 



§ APPENDIX








 §.§ Difference in the log densities of the SC and PKB distributions

The difference of the log-densities, of PKB and SC, can be written as

    logf_SC-logf_PKB   =   d-1/2[ -log(√(γ^2+1)-α) + log(1-(√(γ^2+1)-1)^2/γ^2)] 
       =   d-1/2[log1+ρ^2-2 y^⊤ m/1-ρ^2+log(1-ρ^2)] 
       =   d-1/2log(1+ρ^2-2 y^⊤ m).




 §.§ Derivatives of Eq. (<ref>)


    J   =   ∂ℓ_SC/∂μ = - d∑_i=1^nμ/√(γ^2+1)- y_i√(γ^2+1)-α_i
     H   =   ∂^2ℓ_SC/∂μ∂μ^⊤ = - d∑_i=1^n I_d+1√(γ^2+1)-μμ^⊤/√(γ^2+1)/γ^2+1(√(γ^2+1)-α_i)-(μ/√(γ^2+1)- y_i)(μ/√(γ^2+1)- y_i)^⊤(√(γ^2+1)-α_i)^2.




 §.§ Derivatives of Eq. (<ref>)


    J = ∂ℓ_PKB/∂μ   =    - d+1/2∑_i=1^nμ/√(γ^2+1)- y_i√(γ^2+1)-α_i - nd-1/2( μ/√(γ^2+1)√(γ^2+1)-1 -2μ/γ^2)
    
     H = ∂^2ℓ_PKB/∂μ∂μ^⊤   =    - d+1/2∑_i=1^n I_d+1√(γ^2+1)-μμ^⊤/√(γ^2+1)/γ^2+1(√(γ^2+1)-α_i)-(μ/√(γ^2+1)- y_i)(μ/√(γ^2+1)- y_i)^⊤(√(γ^2+1)-α_i)^2
           - nd-1/2[ I_d+1√(γ^2+1)-μμ^⊤/√(γ^2+1)/γ^2+1(√(γ^2+1)-1)-(μ/√(γ^2+1))(μ/√(γ^2+1))^⊤(√(γ^2+1)-1)^2 - 2 I_d+1γ^2-4μμ^⊤/γ^4].




 §.§ Derivatives of Eq. (<ref>)


    ∂ℓ_SC/∂β_k   =    - d∑_i=1^nμ_ik x_i/√(γ_i^2+1)- y_ik x_i√(γ_i^2+1)-α_i
    ∂^2ℓ_SC/∂β_k∂β_l^⊤   =   {[ - d∑_i=1^n x_i x_i^⊤√(γ_i^2+1)- x_iμ_ikμ_ik x_i^⊤/√(γ_i^2+1)/γ_i^2+1(√(γ_i^2+1)-α_i)-( x_iμ_ik/√(γ_i^2+1)- y_ik x_i)( x_iμ_ik/√(γ_i^2+1)- y_ik x_i)^⊤/(√(γ_i^2+1)-α_i)^2,                                                                                                                                                                    if k=l;                     - d∑_i=1^n- x_iμ_ikμ_il x_i^⊤/√(γ_i^2+1)/γ_i^2+1(√(γ_i^2+1)-α_i)-( x_iμ_ik/√(γ_i^2+1)- y_ik x_i)( x_iμ_il/√(γ_i^2+1)- y_il x_i)^⊤/(√(γ_i^2+1)-α_i)^2,                                                                                                                                                                  if k ≠ l;                                                                                                                                                                           ]}




 §.§ Derivatives of Eq. (<ref>)

The vector of the first derivative is given by

    ∂ℓ_PKB/∂β_k =
    - d+1/2∑_i=1^nμ_ik x_i/√(γ_i^2+1)- y_ik x_i√(γ_i^2+1)-α_i - d-1/2∑_i=1^nμ_ik x_i/√(γ_i^2+1)√(γ_i^2+1)-1 + d-1/2∑_i=1^n2μ_ik x_i/γ_i^2.

The Jacobian matrix of the second derivatives comprises of

    ∂^2ℓ_SC/∂β_k∂β_k^⊤   =   
    - d+1/2∑_i=1^n x_i x_i^⊤√(γ_i^2+1)- x_iμ_ikμ_ik x_i^⊤/√(γ_i^2+1)/γ_i^2+1(√(γ_i^2+1)-α_i)-( x_iμ_ik/√(γ_i^2+1)- y_ik x_i)( x_iμ_ik/√(γ_i^2+1)- y_ik x_i)^⊤/(√(γ_i^2+1)-α_i)^2
           - d-1/2∑_i=1^n x_i x_i^⊤√(γ_i^2+1)- x_iμ_ikμ_ik x_i^⊤/√(γ_i^2+1)/γ_i^2+1(√(γ_i^2+1)-1)-( x_iμ_ik/√(γ_i^2+1))( x_iμ_ik/√(γ_i^2+1))^⊤/(√(γ_i^2+1)-1)^2
           + d-1/2∑_i=1^n2 x_i x_i^⊤γ_i^2-4μ_ik x_iμ_ik x_i^⊤/γ_i^4, 
    ∂^2ℓ_SC/∂β_k∂β_l^⊤   =   
    - d+1/2∑_i=1^n- x_iμ_ikμ_il x_i^⊤/√(γ_i^2+1)/γ_i^2+1(√(γ_i^2+1)-α_i)-( x_iμ_ik/√(γ_i^2+1)- y_ik x_i)( x_iμ_il/√(γ_i^2+1)- y_il x_i)^⊤/(√(γ_i^2+1)-α_i)^2
           - d-1/2∑_i=1^n- x_iμ_ikμ_il x_i^⊤/√(γ_i^2+1)/γ_i^2+1(√(γ_i^2+1)-1)-( x_iμ_ik/√(γ_i^2+1))( x_iμ_il/√(γ_i^2+1))^⊤/(√(γ_i^2+1)-1)^2
           + d-1/2∑_i=1^n-4μ_ik x_iμ_il x_i^⊤/γ_i^4.



apalike


