A flexible wrapped Lindley-type distribution for angular data modelling

A FL E X I B L E W R A P P E D L I N D L E Y - T Y P E D I S T R I B U T I O N F O R A N G U L A R D A T A M O D E L L I N G Johan Ferr eira School of Statistics and Actuarial Science, Univ ersity of the W itwatersrand, South Africa johan.ferreira@wits.ac.za Delene van W yk-de Ridder Department of Statistical Sciences, Univ ersity of Cape T o wn, South Africa Department of Statistics, Univ ersity of Pretoria, South Africa delene.vanwyk@uct.ac.za Janet van Niekerk Department of Statistics, Univ ersity of Pretoria, South Africa janet.vanniekerk@up.ac.za March 18, 2026 A B S T R AC T Flexible distrib utions for modelling angular data have recei ved considerable attention in recent years, with ongoing work extending e xisting circular models to provide greater fle xibility in capturing div erse angular beha viours. In this paper , we introduce and study the w3PL distrib ution, a circular model obtained by extending the wrapped Lindley distrib ution by incorporating two additional shape parameters. The proposed generalisation increases flexibility in modelling concentration and sk ewness while preserving analytical tractability and encompassing e xisting circular models as special cases. Closed-form expressions for the probability density function, cumulati ve distrib ution function, and trigonometric moments are deri ved, allo wing key distrib utional properties to be studied analytically . The distributional modality is characterised, and the nature of in v ariance is inv estigated for the newly proposed circular model. Parameter estimation is de veloped within a re gularised maximum likelihood framew ork, and a simulation study demonstrates reliable parameter recovery and stable finite-sample performance. Applications to angular datasets from geology , marine biology , and finance illustrate the model’ s practical significance and show improv ed fit relativ e to existing circular alternati ves. K eywords directional, SDG goals, regularised maximum likelihood estimation, trigonometric moments, ske wness and kurtosis 1 Introduction Angular data naturally arise in many scientific and applied conte xts where the v ariable of interest is inherently directional or periodic. Unlike linear data, where data are measured simply as positive ( x > 0 ) or unbounded ( −∞ < x < ∞ ), circular data wraps around (i.e., 0 ° and 360 ° represent the same direction), rendering traditional linear statistical models inadequate or misleading. When data are measured in a "non-linear" w ay , the accompanying probabilistic inference should correctly account for this beha viour - in this case, where data are measured as angles. Consequently , there has been a resurgence in recent years in the study of the potential probability structures underlying data distributed on the circle. Circular distrib utions are essential for modelling biological, geological and financial data when the variables in question exhibit such inherent periodicity , such as time of day , wind direction, animal migration orientation, or seasonal cycles. Moreo ver , inspired by the Sustainable Dev elopment Goals (SDGs), it is imperative that research focuses on A P R E P R I N T - M A R C H 1 8 , 2 0 2 6 addressing and understanding the goals in meaningful ways. Proper use of circular statistics for circularly distributed data, especially in the en vironmental sciences, enables more accurate analyses and predictions, which are critical for addressing SDGs such as Goal 13: Climate Action, Goal 14: Life Below W ater , and Goal 15: Life On Land. For example, understanding wind and ocean current directions using circular models supports climate modelling and marine conserv ation planning. This precision supports informed decision-making for sustainability and resilience in the face of en vironmental change. T o illustrate the practical relev ance of circular models, this paper considers three datasets from distinct application areas. As one of the real data examples, we consider the dataset on trav el directions of n = 31 blue periwinkles, as in vestigated by [ 1 ]. In this case, it would be rele vant to marine biologists working in this field to establish the extent to which periwinkles exhibit directional preferences in their movement patterns [ 2 ]. This is meaningful as it can help determine whether periwinkles tend to move in particular directions in response to en vironmental cues such as tidal cycles, moisture gradients or food resources. Naturally , it is understandable that a standard "linear" PDF would not suffice to model the intricacies of data such as this, which emphasises the continued importance of wrapped models for proper statistical inference. Where circular models are formulated in radians (when θ ∈ [0 , 2 π ) ), data measured in degrees should be con verted to radians by multiplying the data by π 180 prior to model fitting. Many authors have described and implemented the methodology of obtaining (and subsequently in vestig ating) a distribution on the unit circle via wrapping; recent examples in the literature include [ 3 , 4 , 5 , 6 , 7 , 8 , 9 , 10 , 11 ]. In particular , [ 12 ] completed a re view of recent discrete and continuous models for the circle via wrapping, and two important monographs within the circular realm are those by [ 13 ] and [ 14 ]. For a linear continuous random v ariable X with (continuous) infinite support with probability density function (PDF) f ( x ; α ) and corresponding parameter vector α , the corresponding circular PDF of the random variable θ = X mo d 2 π such that θ ∈ [0 , 2 π ) is giv en by g ( θ ; α ) = ∞ X m = −∞ f ( θ + 2 mπ ; α ) . If X has (continuous) positiv e support, then the circular PDF is giv en by g ( θ ; α ) = ∞ X m =0 f ( θ + 2 mπ ; α ) . (1) One particular example of a linear distribution that has been wrapped (and recei ved subsequent attention, see [ 15 ]) is the Lindley distrib ution (see [16]), with (linear) PDF f ( x ; κ ) = κ 2 1 + κ (1 + x ) exp ( − κx ) , x > 0 , κ > 0 (2) (where α ≡ κ ) and subsequent circular PDF g ( θ ; κ ) = κ 2 1 + κ exp ( − κθ )  1 + κ 1 − exp ( − 2 π κ ) + 2 π exp ( − 2 π κ ) (1 − exp ( − 2 π κ )) 2  , θ ∈ [0 , 2 π ) , κ > 0 . (3) The cumulativ e distribution function (CDF) for this distrib ution is gi ven by G ( θ ; κ ) = 1 1 − e xp ( − 2 π κ )  1 − exp ( − κθ ) − κθ exp ( − κθ ) κ + 1  − 2 π κ κ + 1 (1 − exp ( − κθ )) exp ( − 2 π κ ) (1 − exp ( − 2 π κ )) 2 . Although it is generally preferred that a distrib ution hav e fe wer parameters than more, the single parameter α in this model may fall short of capturing the true intricacies of the underlying data. Here, the perceptiv e practitioner can foresee that κ is the single parameter that has to accurately describe the (circular) mean, variance, ske wness, and kurtosis (see Appendix A.2 for these expressions) of sample data for a good fit. In the spirit of relieving this single parameter of this burden, it is of both theoretical and inferential v alue to in vestigate models which i) generalises the model in (3) where the generalised model nests this wrapped Lindley as a special case, and ii) subsequently study such a generalisation to establish the behaviour of additional parameters in a model fitting context – in order to verify such a model’ s validity and feasibility . [ 17 ] introduced and studied a three-parameter Lindle y distribution (referred to as the 3PL distribution); for a random variable X follo wing this distrib ution, the PDF is giv en by f ( x ; κ, α, β ) = κ 2 ακ + β ( α + β x ) exp ( − κx ) , x > 0 , β , κ > 0 , ακ + β > 0 . (4) 2 A P R E P R I N T - M A R C H 1 8 , 2 0 2 6 with the CDF giv en by F ( x ; κ, α, β ) = 1 −  1 + κβ x ακ + β  exp ( − κx ) , x > 0 . (5) It follo ws immediately that, in (4) , when α = β = 1 , the usual Lindley distrib ution (2) is obtained. When α = 1 , the two-parameter Lindley distribution of [ 18 ] is obtained. Considering this distribution in a circular context would be impactful, as it nests the usual Lindley distrib ution as a special case and provides meaningful parameter enrichment, offering additional fle xibility that may better reflect the true underlying data characteristics. In this paper , we deriv e and present a novel circular distribution, the wrapped 3PL (w3PL) distribution, for data distributed on the circle. The distrib ution and its key properties and characteristics are presented in Section 2, where the PDF , CDF , trigonometric moments, modality and property of in v ariance are established and in vestigated. In Section 3, aspects of estimation and sampling are discussed, and Section 4 presents numerical studies, including a simulation study and model fitting to three real datasets. Section 5 concludes the paper . 2 A wrapped 3PL distribution In this section, we deriv e and in vestigate a circular companion to the linear PDF in (4) via wrapping (see (1) ). The proofs for these results are outlined in Appendix A.1. Theorem 2.1. Suppose X follows the 3PL distribution with the PDF (4) with parameters ( κ, α, β ) . Then, the corr esponding cir cular w3PL variable Θ has the following PDF: g ( θ ; κ, α , β ) = κ 2 ακ + β exp( − κθ ) 1 − exp( − 2 π κ )  α + β θ + 2 π β exp ( − 2 κπ ) 1 − exp ( − 2 κπ )  (6) wher e θ ∈ [0 , 2 π ) and κ, β > 0 and ακ + β > 0 . Wher e rele vant, this distribution will be denoted by w3PL. W e note her e that, as in the linear case, the "usual" wrapped Lindle y model with PDF (3) is obtained fr om (6) when α = β = 1 (see [15]). Theorem 2.2. Suppose X follows the 3PL distribution with the PDF (4) with parameters ( κ, α, β ) . Then, the CDF of the corr esponding cir cular w3PL variable Θ is given by: G ( θ ; κ, α, β ) =  1 − exp( − κθ ) − κβ θ ακ + β exp( − κθ )  1 1 − e xp( − 2 πκ ) +  2 π κβ ακ + β { 1 − e xp( − κθ ) }  exp( − 2 π κ ) (1 − exp( − 2 π κ )) 2 . (7) As befor e, the "usual" wr apped Lindley model with CDF (5) is obtained fr om (7) when α = β = 1 . Figure 1 represents the PDF (6) and CDF (7) for different parameter configurations for κ , α , and β , respecti vely . 2.1 T rigonometric moments W e present a Lemma that is useful for establishing the trigonometric moments of (6). Lemma 2.1. Suppose a random variable X follows the 3PL distribution with the PDF (4) with parameter s ( κ, α, β ) . Then, the char acteristic function of X is given by ϕ X ( t ) = κ 2 ακ + β  α ( κ − it ) + β ( κ − it ) 2  . (8) wher e i = √ − 1 . Using (8) , the trigonometric moments of order p of the distribution in (6) can be established via the relation ϕ ( p ) = ρ p exp ( iµ p ) = ϕ X ( t ) . This means that the characteristic function of the distrib ution in (6) is giv en by ϕ ( p ) = κ 2 ακ + β  α ( κ − ip ) + β ( κ − ip ) 2  . Let z = a − ib . The polar representation of a complex number gi ves z = r exp( iθ ) 3 A P R E P R I N T - M A R C H 1 8 , 2 0 2 6 (a) w3PL (6) PDF plots (b) w3PL (7) CDF plots Figure 1: w3PL PDF and CDF plots for dif ferent parameter configurations where r = | z | = √ a 2 + b 2 and θ = arg( z ) , as expressed in [19]. For z = a − ib we have θ = arg( a − ib ) = − arctan  b a  , so that a − ib = r exp  − i arctan  b a  . Using the exponential form of De Moi vre’ s formula, ( r exp( iθ )) m = r m exp( imθ ) , which follows from Euler’ s identity exp( iθ ) = cos θ + i sin θ , we obtain ( a − ib ) − m = r − m exp  im arctan  b a  . (9) Using (9), the characteristic function can be written as ϕ ( p ) = κ 2 ακ + β (( ακ + β ) 2 + α 2 p 2 ) 1 2 κ 2 + p 2 exp  i  2 arctan  p κ  − arctan  αp ακ + β  which leads to the identification of µ p = 2 arctan  p κ  − arctan  αp ακ + β  , (10) ρ p = κ 2 ακ + β (( ακ + β ) 2 + α 2 p 2 ) 1 2 κ 2 + p 2 . (11) Subsequently , the noncentral trigonometric moments of the distribution in (6) can be obtained. Since the definition of trigonometric moments is giv en by ϕ p = a p + ib p where a p = ρ p cos ( µ p ) and b p = ρ p sin ( µ p ) , the noncentral trigonometric moments follow as a p = κ 2 ακ + β (( ακ + β ) 2 + α 2 p 2 ) 1 2 κ 2 + p 2 cos  2 arctan  p κ  − arctan  αp ακ + β  , (12) b p = κ 2 ακ + β (( ακ + β ) 2 + α 2 p 2 ) 1 2 κ 2 + p 2 sin  2 arctan  p κ  − arctan  αp ακ + β  . (13) 4 A P R E P R I N T - M A R C H 1 8 , 2 0 2 6 In particular, the central trigonometric moments are given by ¯ a p = ρ p cos ( µ p − pµ 1 ) and ¯ b p = ρ p sin ( µ p − pµ 1 ) leaving: ¯ a p = κ 2 ακ + β (( ακ + β ) 2 + α 2 p 2 ) 1 2 κ 2 + p 2 × cos  2 arctan  p κ  − arctan  αp ακ + β  − p  2 arctan  1 κ  − arctan  α ακ + β  , ¯ b p = κ 2 ακ + β (( ακ + β ) 2 + α 2 p 2 ) 1 2 κ 2 + p 2 × sin  2 arctan  p κ  − arctan  αp ακ + β  − p  2 arctan  1 κ  − arctan  α ακ + β  . As important summary measures in the circular context, the follo wing quantities are of immediate relev ance (see [ 13 ]): • Mean direction (using (10) where p = 1 ): µ = µ 1 = 2 arctan  1 κ  − arctan  α ακ + β  ; (14) • Circular variance (using (11) where p = 1 ): V 0 = 1 − ρ 1 = 1 − ρ = 1 − κ 2 ακ + β (( ακ + β ) 2 + α 2 ) 1 2 κ 2 + 1 ; (15) • Circular deviation: σ 0 = p − 2 log (1 − V 0 ) = v u u t − 2 log κ 2 ακ + β (( ακ + β ) 2 + α 2 ) 1 2 κ 2 + 1 ! ; • Coefficient of ske wness: η 1 = ¯ b 2 V 3 2 0 ; (16) and • Coefficient of kurtosis: η 2 = ¯ a 2 − (1 − V 0 ) 4 V 2 0 . (17) T able 1 illustrates representativ e parameter configurations chosen to highlight distinct behavioural regimes of the w3PL distribution. V alues of the (circular) mean (14) , variance (15) , ske wness (16) , and kurtosis (17) are calculated for arbitrary v alues of the parameters κ, α , and β . The shaded line illustrates, for e xample, agreement with the results of [ 15 ], sho wing the model proposed in that paper as a special case of the one presented in this paper . K ey observ ations are summarised in the dif ferent groups. In Group A the concentration parameter is fixed at a small v alue κ = 0 . 25 . In this regime, the distrib ution approaches a near -uniform baseline, reflected by lar ge values of V 0 and shape measures η 1 and η 2 close to zero. Although variations in α and β influence the distrib ution, their impact is relativ ely mild compared to the dominant effect of a small κ . Group B exhibits the cases where the concentration parameter is fixed at κ = 1 , while the shape parameters α and β are v aried. Changes in these parameters lead to visible shifts in the mean direction, together with increasing departures from symmetry . As the parameter configurations become more extreme, larger magnitudes of ske wness and kurtosis are observed, illustrating ho w the model accommodates increasingly asymmetric and heavy-tailed circular beha viour . Lastly , in Group C, the parameters α and β are held constant while κ increases across v alues surrounding a critical transition point. As κ increases, the mean direction mo ves to ward zero, concentration increases (smaller V 0 ), and the magnitude of the higher-order shape measures gro ws substantially . This behaviour highlights the control that κ exerts ov er modal structure and distributional sharpness. T o further illustrate these relationships, the dependence of the summary measures on the concentration parameter is examined graphically . The summary measures sho wn in T able 1 can be observed in Figure 2, where the k ey summary measures with respect to the concentration parameter κ is illustrated for selected parameter configurations. Figure 2a displays the mean direction µ as a function of κ . Across all parameter configurations, the mean direction decreases 5 A P R E P R I N T - M A R C H 1 8 , 2 0 2 6 T able 1: Representative summary measures of the w3PL distrib ution illustrating low-concentration, shape-variation, and concentration-transition regimes. Group κ α β µ V 0 η 1 η 2 A 0.25 1 3 2.353 0.938 -0.014 -0.015 0.25 1 2 2.233 0.936 -0.017 -0.016 0.25 2 1 1.724 0.902 -0.050 -0.011 0.25 3 1 1.609 0.883 -0.066 -0.003 B 1.00 1 1 1.107 0.441 -0.683 0.526 1.00 1 4 1.373 0.490 -0.497 0.267 1.00 4 1 0.896 0.360 -0.973 1.124 1.00 1 50 1.551 0.500 -0.453 0.230 1.00 50 1 0.795 0.300 -1.228 1.680 C 1.60 6 3 0.673 0.204 -1.815 2.657 2.00 6 3 0.547 0.138 -2.441 4.006 2.40 6 3 0.458 0.099 -2.982 5.236 monotonically and gradually approaches zero as κ increases, confirming the concentration parameter’ s dominant influence on the distrib ution’ s location. Differences between curves reflect the influence of α and β , which primarily shifts the le vel of the mean direction while preserving the ov erall trend. Figure 2b sho ws the circular variance V 0 . For all configurations, V 0 decreases rapidly with increasing κ , indicating increasing concentration around the modal direction. The curves show similar decay patterns, indicating that concentration is largely go verned by κ , while the shape parameters introduce moderate adjustments to dispersion. The behaviour of circular ske wness η 1 is illustrated in Figure 2c. Here it can be seen that as κ increases, the magnitude of ske wness gro ws steadily , moving aw ay from zero and indicating increasing asymmetry in the distribution. P arameter configurations associated with stronger shape effects produce more pronounced departures from symmetry , consistent with the patterns observed in Group B of T able 1. Finally , Figure 2d presents circular kurtosis η 2 . Here, we observe an increase in kurtosis as κ increases, indicating that higher concentration is associated with heavier tails and a sharper modal structure. The difference between the curves highlights how v ariations in α and β magnify higher-order shape characteristics, e ven when the o verall concentration trend is controlled by κ . While the circular mean, v ariance, ske wness and kurtosis provide v aluable insight into the ov erall location, dispersion and shape of the distribution, they do not directly describe the structural form of the PDF itself. A natural next step is therefore to in vestigate the distrib ution’ s modality and determine the conditions under which the w3PL e xhibits a unique mode. 2.2 Modality of the w3PL In this section, the modality of the w3PL distribution with PDF in (6) is formally in vestigated, and conditions go verning the existence and location of the mode are established. Theorem 2.3. Suppose Θ follows the w3PL distribution with the PDF (6) with parameters ( κ, α, β ) . Then ther e exists one unique mode at the point θ mode = 1 κ − α β − 2 π exp( − 2 π κ ) 1 − exp( − 2 π κ ) . (18) From (18), the critical point of the PDF in (6) is giv en by θ mode = S ( κ ) − α β , S ( κ ) := 1 κ − 2 π exp(2 π κ ) − 1 . An interior mode, therefore, requires 0 < θ mode = S ( κ ) − α β < 2 π , which is equiv alent to 0 < S ( κ ) − α/β < 2 π , and particularly requires S ( κ ) > α/β . T o in vestigate the validity of this condition, we in vestigate the limiting beha viour of S ( κ ) as κ → 0 + and κ → ∞ . Consider first lim κ → 0 + S ( κ ) = lim κ → 0 +  1 κ − 2 π exp (2 π κ ) − 1  . Using the T aylor expansion exp (2 π κ ) = 1 + 2 π κ + 2 π 2 κ 2 + O ( κ 3 ) , 6 A P R E P R I N T - M A R C H 1 8 , 2 0 2 6 (a) (b) (c) (d) Figure 2: Summary measures of the w3PL distrib ution as functions of κ for selected ( α, β ) configurations we obtain exp (2 π κ ) − 1 = 2 π κ  1 + π κ + O ( κ 2 )  . Hence, 2 π exp (2 π κ ) − 1 = 1 κ 1 1 + π κ + O ( κ 2 ) . Applying the T aylor expansion of (1 + x ) − 1 about x = 0 , (1 + x ) − 1 = 1 − x + O ( x 2 ) , | x | < 1 , yields 1 1 + π κ + O ( κ 2 ) = 1 − π κ + O ( κ 2 ) . Therefore, lim κ → 0 + S ( κ ) = lim κ → 0 +  1 κ − 1 κ  1 − π κ + O ( κ 2 )   = lim κ → 0 +  π + O ( κ )  = π . 7 A P R E P R I N T - M A R C H 1 8 , 2 0 2 6 Next, as κ → ∞ , we note that exp (2 π κ ) → ∞ , and hence lim κ →∞ S ( κ ) = lim κ →∞  1 κ − 2 π exp (2 π κ ) − 1  = 0 . Consequently , S ( κ ) is a positi ve, strictly decreasing function that maps (0 , ∞ ) onto (0 , π ) . Let r := α β , 0 < r < π. Then there exists a unique cutof f v alue κ r > 0 satisfying S ( κ r ) = r . For 0 < κ < κ r , the stationary point θ mode = S ( κ ) − r lies in (0 , 2 π ) and corresponds to an interior mode. When κ ≥ κ r , the stationary point leaves the interior of the support and the maximum of the PDF occurs at the boundary θ = 0 . Accordingly , the mode of the distribution on [0 , 2 π ) may be written in the piecewise form Mod ( θ ) =    1 κ − 2 π exp(2 π κ ) − 1 − α β , 0 < κ < κ r , 0 , κ ≥ κ r or r > π . The following remark clarifies the e xistence and uniqueness of the cutof f κ r . Remark 2.1. Rewriting the condition of the interior mode yields 0 < S ( κ ) − α β < 2 π ⇐ ⇒ β  S ( κ ) − 2 π  < α < β S ( κ ) . Since α > 0 , the admissible range r educes to max  0 , β  S ( κ ) − 2 π  < α < β S ( κ ) . F or β > 0 , the upper bound is meaningful pr ovided S ( κ ) > 0 . Indeed, S ( κ ) > 0 ⇐ ⇒ 1 κ > 2 π exp (2 π κ ) − 1 ⇐ ⇒ exp (2 π κ ) > 1 + 2 π κ, which holds for all κ > 0 . Consequently , for any r = α/β ∈ (0 , π ) ther e exists a unique cutof f κ r > 0 satisfying S ( κ r ) = r , and the interior mode occur s pr ecisely for 0 < κ < κ r . In Figures 3 and 4 we observe (6) for a fixed α = 6 and β = 4 for two dif ferent values of κ relativ e to the cutoff κ r = 0 . 6102 , which satisfies S ( κ r ) = 1 . 5 . In Figure 3 we observe the PDF for κ = 0 . 2441 < κ r , so that the interior-mode regime is illustrated where the critical point and maximum, θ mode = 0 . 8681 , lies within [0 , 2 π ) . In contrast, Figure 4 shows the PDF for κ = 0 . 7322 > κ r so that the boundary-mode regime is illustrated where the critical point leav es the interior of the support and the PDF is maximised at θ mode = 0 . While the preceding analysis establishes the modality of the w3PL distrib ution, an additional property of interest in circular models is their behaviour under transformations of the reference system. It is particularly desirable that a circular distribution maintains a coherent probabilistic structure under rotations on the circle. 2.3 The property of in variance In this section, the in v ariance property of the w3PL distribution is e xamined using the framew ork dev eloped by [20]. Theorem 2.4. Let T = [ a, b ) with b − a = 2 π , and let Θ ∈ T be a cir cular random variable with PDF g Θ ( · ; ψ ) , wher e ψ ∈ Ψ is a parameter vector . Define a transformed variable Θ ∗ = δ (Θ + ϵ ) , δ ∈ {− 1 , 1 } , ϵ ∈ T , and let g Θ ∗ ( · ; ψ ∗ ) be the PDF of Θ ∗ , with ψ ∗ ∈ Ψ ∗ . Then g Θ is invariant under changes in orientation (ICO) and in variant under chang es in initial dir ection (ICID) if and only if g Θ ∗ ( θ ∗ ; ψ ∗ ) = g Θ ( θ ∗ ; ψ ) for almost every θ ∗ ∈ T , with Ψ ∗ ≡ Ψ , in the sense that the char acteristic functions of Θ and Θ ∗ must belong to the same functional family . 8 A P R E P R I N T - M A R C H 1 8 , 2 0 2 6 Figure 3: Interior -mode regime when 0 < κ < κ r Figure 4: Boundary-mode re gime when κ ≥ κ r W e will no w sho w that the w3PL distrib ution is not in variant, b ut can be transformed to meet the requirements. The characteristic function of θ ∗ for the w3PL distribution is pro vided by ϕ θ ∗ ( p ) = E (exp( ipθ ∗ )) = E (exp( ip [ δ ( θ + ϵ )])) = E (exp( iδ pθ + iδ pϵ )) = exp( iδ pϵ ) E (exp( iδ pθ )) = exp( iδ pϵ ) ϕ θ ( δ p ) It was pre viously stated that ϕ θ ( p ) = a p + ib p so that ϕ θ ( δ p ) = a δ p + ib δ p . Note that a p and b p were defined in (12) and (13), respectiv ely , and that a δ p = κ 2 ακ + β (( ακ + β ) 2 + α 2 p 2 ) 1 2 κ 2 + p 2 cos  2 arctan  δ p κ  − arctan  αδ p ακ + β  b δ p = κ 2 ακ + β (( ακ + β ) 2 + α 2 p 2 ) 1 2 κ 2 + p 2 sin  2 arctan  δ p κ  − arctan  αδ p ακ + β  . Using Euler’ s formula (see [19]) we obtain exp( iδ pϵ ) = cos( δ pϵ ) + i sin( δ pϵ ) , so that ϕ θ ∗ ( p ) = (cos( δ pϵ ) + i sin( δ pϵ ))( a δ p + ib δ p ) = a δ p cos( δ pϵ ) − b δ p sin( δ pϵ ) + i ( a δ p sin( δ pϵ ) + b δ p cos( δ pϵ )) which will be equal to ϕ θ ( p ) if and only if δ = 1 and ϵ = 0 . Within these specific requirements, we conclude that g Θ ( · ; ψ ) is not ICO and ICID. The PDF g Θ ∗ ( · ; ψ ∗ ) , with Θ ∗ = δ (Θ + ϵ ) , δ ∈ {− 1 , 1 } , ϵ ∈ T , and ψ ∗ = ( ψ , δ, ϵ ) ∈ Ψ ∗ is, ho wev er , ICO and ICID, according to [ 20 ]. W e can rewrite Θ ∗ = δ (Θ + ϵ ) in terms of Θ so that Θ = δ − 1 Θ ∗ − ϵ . Therefore, with this transformation, 9 A P R E P R I N T - M A R C H 1 8 , 2 0 2 6 g Θ ∗ ( θ ∗ ; δ, ϵ, κ, α, β ) =     dθ dθ ∗     g Θ ( δ − 1 θ ∗ − ϵ ; κ, α, β ) =     d ( δ − 1 θ ∗ − ϵ )) dθ ∗     g Θ ( δ − 1 θ ∗ − ϵ ; κ, α, β ) = | δ − 1 | g Θ ( δ − 1 θ ∗ − ϵ ; κ, α, β ) = 1 · g Θ ( δ θ ∗ − ϵ ; κ, α, β ) since δ ∈ {− 1 , 1 } = κ 2 ακ + β exp( − κ ( δ θ ∗ − ϵ ))) 1 − e xp( − 2 πκ )  α + β ( δ θ ∗ − ϵ ) + 2 π β exp ( − 2 κπ ) 1 − exp ( − 2 κπ )  , so that the in v ariant family can be expressed as { g Θ ( δ θ − ϵ ; κ, α, β ) : δ ∈ {− 1 , 1 } , ϵ ∈ [0 , 2 π ) } . Having established ke y structural properties of the w3PL distribution, attention no w turns to statistical inference, where procedures for parameter estimation and sampling from the model are dev eloped. 3 Estimation and sampling In this section, we outline the estimation frame work for the w3PL distribution and describe the sampling and simulation design used to assess parameter recovery . W e first present the maximum likelihood estimator and its regularised extension in Section 3.1, and thereafter detail the sampling scheme employed to generate data for ev aluating finite- sample performance in Section 3.2. 3.1 Estimation using maximum likelihood Suppose that ( θ 1 , ..., θ n ) is a random sample from the w3PL distribution with PDF in (6) . The lik elihood function is then giv en by L ( κ, α, β | θ 1 , ..., θ n ) = n Y i =1 g ( θ i ; κ, α, β ) = n Y i =1 κ 2 ακ + β exp ( − κθ i )  α + β θ i 1 − exp ( − 2 κπ ) + 2 β π exp ( − 2 κπ ) (1 − e xp ( − 2 κπ )) 2  = κ 2 n ( ακ + β ) n exp − κ n X i =1 θ i ! n Y i =1  α + β θ i 1 − exp ( − 2 κπ ) + 2 β π exp ( − 2 κπ ) (1 − exp ( − 2 κπ )) 2  from where the log-likelihood (LL) follo ws as LL ( κ, α, β | θ 1 , ..., θ n ) = 2 n log ( κ ) − n log( ακ + β ) − κ n X i =1 θ i + n X i =1 log  α + β θ i 1 − exp ( − 2 κπ ) + 2 β π exp ( − 2 κπ ) (1 − exp ( − 2 κπ )) 2  where log denotes the natural logarithm. Closed-form solutions to the score equations are not tractable; therefore, maximum likelihood estimates (MLEs) are obtained via numerical optimisation. The parameter space is restricted to κ > 0 , α > 0 , and β > 0 . In practice, optimisation is performed in log-parameters to enforce positivity constraints. Computation is carried out in R using the optim function with the L-BFGS-B algorithm. T o ensure inference is reliable and ef ficient, we employ a re gularised maximum likelihood approach to estimate the parameters of the w3PL distrib ution. Regularisation plays an important role in likelihood-based inference for directional and circular models, particularly in settings where the likelihood may be irregular , flat or poorly behav ed. As noted by [ 21 ], regularisation, or a penalised likelihood approach, is commonly used in likelihood-based estimation to cir cumvent unboundedness in mixture likelihoods and to possibly reduce o verfitting. It can impro ve small-sample generalisation and stabilise numerical optimisation, thereby yielding better-beha ved estimators. 10 A P R E P R I N T - M A R C H 1 8 , 2 0 2 6 In mixtures of v on Mises-Fisher distrib utions, [ 22 ] sho ws that their likelihood function is unbounded from abov e so that a global MLE fails to exist. T o overcome this, a penalty term is introduced so that the resulting penalised lik elihood function is bounded and the existence of the penalised MLE is guaranteed. Similarly , in testing homogeneity for mixtures of von-Mises distributions, [ 23 ] adds a penalty term to the LL function to obtain a tractable limiting distribution. In multi variate circular modelling, [ 24 ] highlights that penalising model comple xity helps av oid ov erfitting when the number of samples is relativ ely small compared to the number of model parameters, and proposes Frobenius-norm regularisation for the multiv ariate von Mises distribution. These results illustrate that penalised likelihood methods are well established in directional statistics for controlling model complexity and ensuring stable estimation. In line with this literature, we no w introduce an L2-regularised maximum likelihood estimator (L2-RMLE) for the w3PL distribution and in vestigate its beha viour through simulation and parameter recovery . W e propose the L2-RMLE, which optimises the LL in a numerically stable manner while adding a quadratic shrinkage term to mitigate instability in finite samples for inferring from a w3PL distribution. The regularised objecti v e augments the negati ve LL with a L2-RMLE in log-space, shrinking parameters toward chosen tar gets ( α 0 , β 0 ) , defined as follows: ∆ = − LL ( κ, α, β | θ 1 , ..., θ n ) + λ α (log α − log α 0 ) 2 + λ β (log β − log β 0 ) 2 . Since the penalty is quadratic in the log-parameters, the estimator remains an L2-RMLE, but applied in the transformed parameter space. Consequently , multiplicati ve de viations, rather than additiv e deviations, are shrunk, which is natural for positiv e scale-like parameters. The penalty weights stabilise the optimisation and constrain the estimator away from regions where the lik elihood is nearly flat or the data pro vide weak identification of ( α, β ) combinations. By pulling estimates to ward sensible targets, the penalty reduces v ariance and often lowers the MSE, e ven though it introduces bias. This approach trades a small, controlled bias for substantial v ariance reduction. In this approach, the penalty weights are scaled by the sample size so that shrinkage weakens as n increases, ensuring that the estimator asymptotically approaches the classic MLE as n → ∞ , under identifiability . By adopting a vanishing re gularisation schedule λ i = c n for i = α, β , we obtain stabilisation in small samples while ensuring that the penalty influence diminishes as n increases. The constant c is selected to achiev e stable con vergence with minimal shrinkage. It is important to note that L2-RMLE is sensiti ve to the choice of penalty weights. If λ i is too large or does not decay with an increase in the sample size n , the estimator can be dominated by the penalty and track the ’tar get’ rather than the data. In ex ecuting the L2-RMLE, a multistart approach is employed, providing a practical globalisation strategy that reduces dependence on initialisation and helps avoid con vergence to flat ridges of the likelihood surface, an observed feature of this model. The multistart approach is based on K random initial points drawn from broad parameter ranges and a bounded quasi-Newton optimiser (L-BFGS-B) run from each initial point. The solution with the lo west objectiv e value is retained. 3.2 Sampling and simulation design T o ev aluate estimator performance under controlled conditions, samples are generated from the w3PL distribution using a numerical implementation of the in verse transform method. Since repeated root-solving of the theoretical CDF in (7) prov ed computationally unstable, inv ersion is performed numerically by constructing a monotone grid-based approximation to the in verse CDF . For fix ed ( κ, α, β ) the analytic CDF G ( θ ) is ev aluated on a dense grid θ 1 , ..., θ m spanning [0 , 2 π ) . Since finite-precision ev aluation can introduce non-monotonicities or v alues slightly outside [0 , 1] , the grid CDF values are clamped to [0 , 1] , and monotonicity is then enforced by applying a cumulativ e maximum operator , ensuring a nondecreasing CDF . The grid is then thinned by remo ving near-duplicate CDF le vels so that in version is stable. The resulting strictly increasing pairs { ( F i , θ i ) } define a numerical approximation to the in verse CDF via linear interpolation. Sampling reduces to dra wing U ∼ Unif (0 , 1) and returning θ = ˆ F − 1 ( U ) . T o examine parameter recovery , samples of sizes n = 50 , 100 , 250 and 500 are generated under fixed parameter configurations. For each of these configurations, we obtain the corresponding L2-regularised maximum likelihood estimates of κ , α , and β . The procedure is repeated N = 200 times for each sample size. Let λ ∗ denote the true parameter value. Performance is ev aluated using empirical bias and mean squared error (MSE), defined as bias ( ˆ λ ) = 1 N N X i =1 ( ˆ λ − λ ∗ ) MSE ( ˆ λ ) = 1 N N X i =1 ( ˆ λ − λ ∗ ) 2 These metrics quantify systematic deviation from the truth and o verall estimation accurac y , respectiv ely . 11 A P R E P R I N T - M A R C H 1 8 , 2 0 2 6 4 Numerical examples In this section, we e v aluate the finite-sample performance of the w3PL distrib ution through a simulation e xperiment, followed by applications to real circular datasets. The simulation study in Section 4.1 focuses on parameter reco very under both interior- and boundary-mode re gimes to assess estimation accuracy and stability across distinct re gions of the parameter space. In addition, we examine how estimation beha viour changes with increasing sample size, providing insight into con ver gence trends. Three real data applications are presented in Section 4.2 to illustrate the model’ s practical applicability . 4.1 Simulation studies W e ev aluate the performance of the w3PL distribution in a simulation setting. T ables 2 to 5 inv estigate different parameter configurations, using grid-in version with L2-RMLE, to illustrate both interior- and boundary-mode re gimes for the model. In all cases, initial v alues were generated within a multistart frame work. F or the L2-RMLE, the shrinkage targets in log-space were set to the logarithm of the initial values, i.e., log( α 0 ) and log( β 0 ) . This choice centres the penalty locally around plausible parameter magnitudes while a voiding shrinkage toward arbitrary fix ed constants. Based on preliminary sensiti vity analyses, a constant c = 0 . 5 was selected across all configurations to achiev e stable con ver gence with minimal shrinkage. Each table is based on N = 200 Monte Carlo replicates, with sample sizes ranging from n = 50 to 500 . W e consider four representative parameter configurations that deliberately span both interior- and boundary-mode regimes. By including both regimes, we ensure that the simulation study accurately reflects the model’ s geometric structure. The parameter configurations considered are: • ( κ, α, β ) = (1 , 2 , 5) , an interior-mode regime where κ < κ r = 2 . 5 (T able 2) • ( κ, α, β ) = (4 , 1 , 8) , an interior-mode regime where κ < κ r = 8 (T able 3) • ( κ, α, β ) = (2 , 4 , 1) , a boundary-mode regime where r > π (T able 4) • ( κ, α, β ) = (8 , 1 , 6) , a boundary-mode regime where κ > κ r = 6 (T able 5) T able 2: P arameter recov ery for the interior-mode configuration, ( κ, α, β ) = (1 , 2 , 5) . n Parameter Estimate SD Bias MSE 50 κ 0.942 0.206 -0.058 0.045 α 2.572 1.364 0.572 2.151 β 6.721 7.714 1.721 11.279 100 κ 0.926 0.162 -0.074 0.031 α 2.566 1.594 0.566 2.811 β 5.417 2.627 0.417 6.938 250 κ 0.885 0.090 -0.115 0.021 α 2.675 0.765 0.675 1.030 β 4.238 1.187 -0.762 1.962 500 κ 0.900 0.063 -0.100 0.014 α 2.489 0.518 0.489 0.502 β 4.391 0.859 -0.609 1.095 Across the four configurations, the simulation results rev eal a distinction between interior- and boundary-mode regimes. In the interior-mode settings, captured in T ables 2 and 3, parameter recovery is stable and impro ves consistently with increasing sample size. For both configurations, bias and v ariability decrease as n increases, with the concentration parameter κ generally well estimated e ven in moderate samples. The parameters α and β exhibit more dispersion in small samples, with variance contracting steadily , leading to reductions in MSE and con ver gence tow ard the true values. The boundary-mode configurations, captured in T ables 4 and 5, display greater instability , especially for α and β . While κ remains reasonably well recovered across sample sizes, estimation of especially α shows ele v ated variability and slower MSE reduction, reflecting the flatter , ridge-like geometry of the likelihood when the mode lies on or near the boundary of the support. Despite these observations, increasing the sample size improves stability and reduces error . Overall, the results indicate that the w3PL distribution is reliably estimable in interior-mode re gimes and practically workable in boundary-mode settings, despite greater finite-sample v ariability . 12 A P R E P R I N T - M A R C H 1 8 , 2 0 2 6 T able 3: P arameter recov ery for the interior-mode configuration, ( κ, α, β ) = (4 , 1 , 8) . n Parameter Estimate SD Bias MSE 50 κ 4.143 0.651 0.143 0.435 α 1.087 0.454 0.087 0.210 β 9.532 5.670 1.532 13.847 100 κ 4.115 0.456 0.115 0.217 α 0.987 0.392 -0.013 0.151 β 9.870 3.781 1.871 7.511 250 κ 3.971 0.258 -0.029 0.066 α 1.074 0.219 0.074 0.052 β 8.153 1.704 0.153 2.868 500 κ 3.997 0.183 -0.003 0.033 α 1.034 0.156 0.034 0.025 β 8.308 1.252 0.308 1.631 T able 4: P arameter recov ery for the boundary-mode configuration, ( κ, α , β ) = (2 , 4 , 1) . n Parameter Estimate SD Bias MSE 50 κ 2.302 0.507 0.302 0.343 α 2.133 1.209 -1.867 4.916 β 1.589 1.207 0.589 1.775 100 κ 2.272 0.440 0.272 0.264 α 2.836 1.714 -1.164 4.232 β 1.257 0.934 0.257 0.921 250 κ 2.053 0.318 0.052 0.102 α 3.950 2.249 -0.050 4.960 β 0.845 0.585 -0.155 0.360 500 κ 2.064 0.232 0.064 0.057 α 4.059 2.816 0.058 2.775 β 0.831 0.495 -0.169 0.269 T able 5: P arameter recov ery for the boundary-mode configuration, ( κ, α , β ) = (8 , 1 , 6) . n Parameter Estimate SD Bias MSE 50 κ 8.092 1.514 0.092 2.254 α 1.077 0.679 0.077 0.458 β 7.195 5.222 1.195 18.150 100 κ 8.237 1.418 0.237 2.028 α 1.213 0.905 0.213 0.848 β 6.426 3.506 0.426 12.226 250 κ 7.802 0.943 -0.198 0.911 α 1.157 0.609 0.157 0.388 β 5.430 1.925 -0.570 3.954 500 κ 7.990 0.634 -0.010 0.394 α 1.076 0.811 0.076 0.651 β 5.629 1.313 -0.371 1.827 4.2 Real data In this section, we inv estigate the behaviour of the w3PL distribution by applying it to three real-world datasets, namely the cross-bed measurements (Section 4.2.1), tra vel directions of blue periwinkles (Section 4.2.2), and credit-card transaction timestamps (Section 4.2.3). Model performance is ranked via the Akaike information criterion (AIC; [ 25 ]), AIC = 2 k − 2 LL ( ˆ κ, ˆ α, ˆ β ) , the corrected Akaike information criteria (AICc; [26]) for small samples, AICc = AIC + 2 k ( k + 1) n − k − 1 , 13 A P R E P R I N T - M A R C H 1 8 , 2 0 2 6 and the Bayesian information criterion (BIC; [27]), BIC = log( n ) k − 2 LL ( ˆ κ, ˆ α, ˆ β ) , where k is the number of parameters and n is the number of observations, ˆ κ, ˆ α , and ˆ β are the ML estimates of κ, α , and β , respecti vely , and where LL ( ˆ κ, ˆ α, ˆ β ) is the maximized LL v alue for the model. In this formulation of both the AIC and BIC cases, the model with the smallest value is considered the preferred model. Where relev ant, this is indicated in the results with the smallest value indicated in bold . T o illustrate the viability of the model in (6) as a meaningful alternati ve within the circular frame work, we benchmark it to similar models as those considered in [ 15 ] and [ 5 ]. Specifically , we compare it via the AIC and BIC with the wrapped Lindley (wL) model (3), the wrapped XLindle y (wXL) model of [5] with PDF g ( θ ; κ ) (19) = κ 2 exp ( − θ κ ) ( κ + 1) 2 (1 − exp ( − 2 π κ )) 2 [(1 − exp ( − 2 π κ ))( κ + θ + 2) + 2 π exp ( − 2 π κ )] , θ ∈ [0 , 2 π ) , κ > 0 , and the wrapped exponential (wExp) model of [28] with PDF g ( θ ; κ ) = κ exp ( − κθ ) 1 − exp ( − 2 π κ ) , θ ∈ [0 , 2 π ) , κ > 0 . (20) The notation of the parameter κ is kept uniform throughout these models for notational simplicity; ho wev er , we note that, in a model such as (20) , the role of κ dif fers by design from that in other considered models, where it serves as the "original" Lindley parameter . The datasets used in this paper were initially introduced and made av ailable by [ 13 ] and are easily accessible in R via the circular package. 4.2.1 Cross-bed measurements Cross-stratification is a primary sedimentary structure used to infer palaeoflo w directions. The geometry of trough or planar cross-beds is linked to the migration of dunes under a pre v ailing flo w , so that collections of measured orientations yield a directional sample that can be analysed to reconstruct flow re gimes and sediment transport patterns, as described by [ 29 ]. The dataset under consideration contains n = 104 observations of cross-bed orientations from the Himalayan molasse in Pakistan, specifically the Chaudan Zam large bedforms, and is av ailable as cross.beds2 in the R package, NPCirc , reported in [13]. (a) Empirical histogram and fitted PDFs (b) Rose diagram and fitted PDFs Figure 5: Model fits to the cross-bed measurement data W e fit the proposed w3PL distribution alongside the considered contenders, with results summarised in T able 6. Based on both AIC and BIC, the w3PL provides the strongest empirical fit for these data. The impro vement in information criteria is substantial; relativ e to the next-best contender (wXL), the AIC difference e xceeds 10 units, providing strong evidence in fav our of the additional flexibility introduced by the w3PL ’ s two shape parameters. The fitted model 14 A P R E P R I N T - M A R C H 1 8 , 2 0 2 6 T able 6: Results of MLE to the cross-bed data according to the AIC and BIC. Parameter/Measure w3PL (6) wL (3) wXL (19) wExp (20) κ 0.637 0.200 0.361 0.174 α 0.833 - - - β 119.992 - - - LL -177.094 -236.521 -184.259 -186.125 AIC 360.189 475.041 370.517 374.250 BIC 368.122 477.686 373.162 376.894 comparisons are shown in Figure 5. The empirical histogram and rose diagram suggest a unimodal pattern with moderate dispersion and mild asymmetry . amongst fitted models, the w3PL PDF most closely follows the observed peak location and overall pattern, while the competing models either under -represent the concentration around the dominant direction or produce less accurate tail beha viour . This is consistent with the superior AIC and BIC values reported in T able 6. 4.2.2 T rav el directions of blue periwinkles [ 30 ] describes periwinkles as ab undant intertidal grazers that influence community structure by feeding on algae and other organisms on rock y shores and marsh grasses. Their grazing can alter algal ab undance and distribution, af fecting habitat structure and other species in the food web . Movement studies hav e revealed that distance and direction are often correlated with habitat structure, sensory cues, and ecological pressures, making this type of data biologically significant, as explained by [ 2 ]. The dataset under consideration contains n = 31 blue periwinkles, for which the direction of displacement was recorded after individuals were transferred do wnshore from the heights at which they normally liv e, as described in [1]. The data are a vailable in the periwinkles dataset of the R package, NPCirc . (a) Empirical histogram and fitted PDFs (b) Rose diagram and fitted PDFs Figure 6: Model fits to the tra vel directions of blue periwinkles T able 7: Results of MLE to the periwinkle data according to the AIC and BIC. Parameter/Measure w3PL (6) wL (3) wXL (19) wExp (20) κ 1.155 1.165 0.708 0.508 α 0.944 - - - β 105.974 - - - LL -37.988 -51.837 -45.409 -46.462 AICc 82.865 105.811 92.955 95.062 BIC 86.278 107.108 94.251 96.358 W e fit the proposed w3PL, alongside the competing wrapped distributions, with results summarised in T able 7. Based on both the AICc and BIC, the w3PL pro vides the strongest empirical fit amongst the considered models. The improv ement 15 A P R E P R I N T - M A R C H 1 8 , 2 0 2 6 in AICc relati ve to the next-best competitor is substantial, indicating strong support for the additional fle xibility pro vided by the two-parameter shape structure of the w3PL. Gi ven the relati v ely small sample size, the use of AICc is appropriate, and the agreement between AICc and BIC strengthens the robustness of the selection. The fitted parameters indicate a moderate concentration, along with asymmetry in the data. Figure 6 highlights the superiority of the w3PL fit, where the empirical histogram and corresponding rose diagram suggest a unimodal but ske wed pattern. The w3PL most closely tracks the observed peak location and asymmetric decay from the mode, while the competing models either ov er-concentrate around the peak or fail to capture the tail beha viour in the sample. 4.2.3 Credit-card transaction timestamp data T ransaction timestamp data reflect regular patterns in consumer beha viour and system usage, with time-of-day having shown to be an informativ e predictor in credit card fraud detection and transaction modelling, as highlighted by [ 31 ] and [ 32 ]. Because the 24-hour clock is inherently periodic, timestamps reduced to hour-of-day form a circular v ariable, where 23:00 and 00:00 are adjacent. The dataset considered here consists of n = 349 transactions with timestamps occurring between 1 January and 29 July 2020, obtained via the PyCircular package in Python . Transactions were filtered by transaction type (types 1 and 2) and con verted to circular observ ations representing time of day; i.e., transaction timestamps were parsed to a date-time format and binned to the hour of day h ∈ { 0 , ..., 23 } . Hourly frequencies were computed, and each hour was encoded as an equally spaced circular angle via θ ( h ) = 2 π h/ 24 , equiv alent to 15 h degrees. An observation-le vel circular sample w as formed by repeating each hour-angle according to its hourly count. The resulting sample represents transaction intensity over the 24-hour c ycle and may exhibit asymmetric concentration due to peak activity periods. (a) Empirical histogram and fitted PDFs (b) Rose diagram and fitted PDFs Figure 7: Model fits to the credit-card transaction timestamp data T able 8: Results of MLE to the credit-card transaction timestamp data according to the AIC and BIC. Parameter/Measure w3PL (6) wL (3) wXL (19) wExp (20) κ 0.842 0.990 0.592 0.392 α 5.536 - - - β 18.063 - - - LL -487.804 -608.870 -498.228 -503.193 AIC 981.609 1219.741 998.457 1008.385 BIC 992.828 1223.481 1002.197 1012.125 W e fit our proposed model, the w3PL, alongside the considered contenders, and summarise the results in T able 8. Based on both AIC and BIC, the w3PL pro vides the strongest empirical fit. Notably , the AIC dif ference between the w3PL and the next-best competitor e xceeds 15 units, providing strong e vidence in fa v our of the additional flexibility introduced by the model’ s two-parameter shape structure. The BIC results lead to the same conclusion, reinforcing the robustness of the selection. The estimated w3PL parameters indicate a moderate concentration and a pronounced asymmetric component, gov erned by α and β , within the interior -mode regime. The findings in T able 8 are reflected in Figure 7, 16 A P R E P R I N T - M A R C H 1 8 , 2 0 2 6 where the empirical histogram and corresponding circular rose diagram indicate a unimodal distribution of transaction times with moderate concentration and noticeable sk ewness. Amongst the fitted models, the w3PL PDF aligns most closely with the observed peak location and overall pattern, capturing the asymmetric shape more effecti vely than the competing wrapped distributions. 5 Conclusion In this paper , we deri ved, presented, and studied a pre viously unconsidered flexible wrapped distrib ution, namely the w3PL distribution. An immediate generalisation of the wrapped Lindley model of [ 15 ], our proposition uses the 3PL distribution of [ 33 ] as the linear basis in a wrapping context and nests the wrapped Lindley model as a special case. In particular , the PDF , CDF , and trigonometric moments were deriv ed in closed form, yielding expressions for calculating key circular characteristics, including the circular mean, v ariance, ske wness, and kurtosis. The w3PL distribution (6) achiev es increased flexibility through the additional tw o shape parameters α and β , in addition to the parameter κ , which the model inherited from the classical Lindley and wrapped Lindley distributions. The inclusion of the parameters α and β allows additional parameters to gov ern the higher-order moments of the distribution, thereby providing greater control over skewness and kurtosis while retaining the structural simplicity of the wrapped Lindley framework. At the same time, the model remains tractable and computationally feasible. Simulation from the w3PL distribution is computationally feasible, and likelihood-based parameter estimation can be implemented in a stable and practical manner within a regularised framew ork. Considering these advantages, it is worth noting that the proposed distribution closely aligns with the desirability guidelines for flexible distrib utions outlined by [34]. The theoretical properties were juxtaposed with those of the wrapped Lindley model, sho wcasing the advantage of the additional parameters α and β . In particular: • α and β giv e access to a broader range of ske wness and kurtosis for the practitioner , enhancing descriptiv e flexibility; • the model is reliable in parameter recovery across dif ferent sample sizes and regime configurations, including both interior- and boundary-mode settings; and • The model is conv enient to implement in practice for fitting real circular datasets. In the simulation study , we observ ed that estimation accuracy impro ves with increasing sample size, with empirical bias and MSE generally decreasing as n grows. The L2-RMLE approach stabilised estimation in regions where the likelihood surface exhibits flat ridges, while maintaining asymptotic consistency . Overall, the simulation results confirm that the w3PL distribution can numerically recov er the true parameter values with satisfactory accurac y across a range of configurations. Finally , across three real data fitting scenarios, the w3PL distribution proved a viable contender in the circular modelling space. Across the applications, the w3PL sho wed superior performance, as measured by common model performance metrics: AIC, AICc, and BIC. The improvements were particularly e vident in datasets exhibiting asymmetric concentration, where the additional shape flexibility of the w3PL translated into tangible g ains in descriptiv e accuracy . The theoretical tractability , simulation performance, and empirical competiti veness of the w3PL distrib ution suggest that it provides a useful and fle xible addition to practitioners’ toolkits of wrapped distrib utions for modelling asymmetric circular data. Acknowledgments Ferreira has been partially supported by the National Research Foundation (NRF) of South Africa (SA), grants RA201125576565 (nr 145681) and CSRP250220299317, as well as the DSI-NRF Centre of Excellence in Mathematical and Statistical Sciences (CoE-MaSS), South Africa. V an W yk-de Ridder acknowledges the support of the Faculty of Science at the Univ ersity of Cape T own, as well as the Centre for Statistics in Ecology , the Environment and Conserv ation, based at the University of Cape T o wn. The opinions expressed and conclusions arri ved at are those of the authors and are not necessarily to be attributed to the NRF . 17 A P R E P R I N T - M A R C H 1 8 , 2 0 2 6 Data availability statement All datasets considered in this paper are freely av ailable, particularly av ailable in the R package circular or NPCirc (which can be found in hardcopy in [ 13 ]) or in Python . The names of the datasets in this package are reported in the Real Data Application section. Disclosure statement The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article. A uthor contrib utions All authors contributed to the study conception and design. Material preparation and data analysis were performed by Ferreira, van W yk-de Ridder and v an Niekerk. The first draft of the manuscript was written by Ferreira, and all authors commented and contributed to earlier v ersions. All authors read and appro ved the final manuscript. A ppendix A Proofs and other expressions A.1 Characteristics of the w3PL distribution (6) Pr oof of Theor em 2.1. By using (1) and (4), we see that g ( θ ; κ, α , β ) = ∞ X m =0 f ( θ + 2 mπ ; κ, α, β ) = ∞ X m =0 κ 2 ακ + β ( α + β ( θ + 2 mπ )) exp ( − κ ( θ + 2 mπ )) = κ 2 ακ + β exp ( − κθ )  ( α + β θ ) 1 − exp ( − 2 π κ ) + 2 β π exp ( − 2 κπ ) (1 − exp − 2 κπ ) 2  which leav es the final result. Pr oof of Theor em 2.2. By using the definition of a circular CDF (see [35]) and (5), we see that G ( θ ; κ, α, β ) = ∞ X m =0 F ( θ + 2 π m ; κ, α, β ) − F (2 π m ; κ, α , β ) = ∞ X m =0  1 −  1 + κβ ( θ + 2 π m ) ακ + β  exp ( − κ ( θ + 2 π m )) −  1 −  1 + κβ (2 π m ) ακ + β  exp ( − κ (2 π m ))  = 1 1 − exp( − 2 π κ )  1 − exp( − κθ ) − κβ θ ακ + β exp( − κθ )  + exp( − 2 π κ ) (1 − exp( − 2 π κ )) 2  2 π κβ ακ + β { 1 − exp( − κθ ) }  which leav es the final result. Pr oof of Lemma 2.1. By using the definition of the characteristic function and (4), we see that ϕ X ( t ) = E (exp ( itX )) = Z ∞ 0 exp ( itx ) κ 2 ακ + β ( α + β x ) exp ( − κx ) dx = κ 2 ακ + β  Z ∞ 0 α exp ( − x ( κ − it )) dx + Z ∞ 0 β x e xp ( − x ( κ − it )) dx  = κ 2 ακ + β  α ( κ − it ) + β ( κ − it ) 2  . which leav es the final result. 18 A P R E P R I N T - M A R C H 1 8 , 2 0 2 6 Pr oof of Theor em 2.3. By using (6) and the chain rule, we determine the deriv ati ve and set it equal to 0 : d dθ g ( θ ; κ, α , β ) = d dθ  κ 2 ακ + β exp( − κθ ) 1 − e xp( − 2 πκ )  α + β θ + 2 π β exp ( − 2 κπ ) 1 − exp ( − 2 κπ )  = d dθ  ακ 2 ακ + β exp( − κθ ) 1 − e xp( − 2 πκ ) + β θ κ 2 ακ + β exp( − κθ ) 1 − exp( − 2 π κ ) + κ 2 ακ + β 2 π β exp( − 2 π κ ) 1 − exp( − 2 π κ ) exp( − κθ ) 1 − exp( − 2 π κ )  = κ 2 ( ακ + β )(1 − exp( − 2 π κ )) exp( − κθ )  − ακ + β − β κθ − 2 π κβ exp( − 2 π κ ) 1 − e xp( − 2 πκ )  which implies exp ( − κθ ) | {z } =0  − ακ + β − β κθ − 2 π κβ exp( − 2 π κ ) 1 − exp( − 2 π κ )  | {z } =0 = 0 . Solving for the second bracket lea ves the final result. T o confirm that the critical point θ mode in (18) is indeed a maximum, let C = κ 2 ( ακ + β )(1 − exp( − 2 πκ )) and A = 2 π β exp( − 2 π κ ) 1 − exp( − 2 π κ ) so that d dθ g ( θ ; κ, α , β ) = C exp( − κθ ) ( β − κ ( α + β θ + A )) . From this expression, we note that when we set β − κ ( α + β θ mode + A ) = 0 , we end up with κ ( α + β θ mode + A ) = β . (21) Using the chain rule, we find d 2 g ( θ ) d 2 θ = d dθ ( C exp( − κθ ) ( β − κ ( α + β θ + A ))) = − κC exp( − κθ ) ( β − κ ( α + β θ + A )) − κβ C exp( − κθ ) = κC exp( − κθ )( κ ( α + β θ + A ) − 2 β ) If we ev aluate d 2 g ( θ ) d 2 θ at θ mode (via (21)), we hav e κ ( α + β θ mode + A ) − 2 β = β − 2 β = − β , so that d 2 g ( θ mode ) d 2 θ mode = C exp( − κθ mode ) κ ( − β ) = − κβ C exp( − κθ mode ) < 0 ∀ κ, α, β > 0 , confirming that the mode is a strict local maximum. A.2 Expressions of (circular) mean, variance, skewness, and kurtosis for the wrapped Lindley distribution (wL) (3) [15] • Mean direction: µ = 2 arctan  1 κ  − arctan  1 κ + 1  ; • Circular variance: V 0 = 1 − κ 2 κ + 1 (( κ + 1) 2 + 1) 1 2 κ 2 + 1 ; • Circular deviation: σ 0 = p − 2 log (1 − V 0 ) = v u u t − 2 log κ 2 κ + 1 (( κ + 1) 2 + 1) 1 2 κ 2 + 1 ! ; 19 A P R E P R I N T - M A R C H 1 8 , 2 0 2 6 • Coefficient of ske wness: η 1 = ¯ b 2 V 3 2 0 where ¯ b 2 = = κ 2 κ + 1 (( κ + 1) 2 + 4) 1 2 κ 2 + 4 × sin  2 arctan  2 κ  − arctan  2 κ + 1  − 2  2 arctan  1 κ  − arctan  1 κ + 1  ; and • Coefficient of kurtosis: η 2 = ¯ a 2 − (1 − V 0 ) 4 V 2 0 where ¯ a 2 = κ 2 κ + 1 (( κ + 1) 2 + 4) 1 2 κ 2 + 4 × cos  2 arctan  2 κ  − arctan  2 κ + 1  − 2  2 arctan  1 κ  − arctan  1 κ + 1  . 20 A P R E P R I N T - M A R C H 1 8 , 2 0 2 6 References [1] Nicholas I Fisher and Alan J Lee. Regression models for an angular response. Biometrics , pages 665–677, 1992. [2] Peter S Petraitis. Occurrence of random and directional mo vements in the periwinkle, littorina littorea (l.). Journal of Experimental Marine Biology and Ecology , 59(2-3):207–217, 1982. [3] W illiam Bell and Saralees Nadarajah. The wrapped Rama distribution. Scientific Reports , 14(1):31936, 2024. [4] Abdullah Y ilmaz and Cenker Biçer . A ne w wrapped exponential distribution. Mathematical Sciences , 12:285–293, 2018. [5] E Zinhom, MM Nassar, SS Radwan, and A Elmasry . The wrapped XLindley distribution. En vir onmental and Ecological Statistics , 30(4):669–686, 2023. [6] B Sango Otieno and Christine M Anderson-Cook. Measures of preferred direction for en vironmental and ecological circular data. En vir onmental and Ecological Statistics , 13:311–324, 2006. [7] Lu Chen, V ijay P Singh, Shenglian Guo, Bin Fang, and Pan Liu. A new method for identification of flood seasons using directional statistics. Hydr ological Sciences Journal , 58(1):28–40, 2013. [8] Sahana Bhattacharjee, Inzamul Ahmed, and Kishore Kumar Das. Wrapped two-parameter Lindle y distribution for modelling circular data. Thailand Statistician , 19(1):81–94, 2021. [9] Said Benlakhdar , Mohammed Rziza, and Rachid Oulad Haj Thami. In-depth analysis of von Mises distribution models: Understanding theory , applications, and future directions. Statistics, Optimization & Information Computing , 12(4):1210–1230, 2024. [10] Sakhr AS Alotaibi, Nora Muda, Ahmad MH Al-Khazaleh, and Shawkat Alkhazaleh. The wrapping Epanechnikov exponential distrib ution: A nov el flexible model for asymmetric circular data. PLoS One , 20(9):e0331906, 2025. [11] Ahmed M. Gemeay and Ahmed M. T . Abd El-Bar . On wrapping a new v ersion of Lindley distribution for circular data analysis. En vir onmental and Ecological Statistics , 2026. [12] W illiam Bell and Saralees Nadarajah. A review of wrapped distributions for circular data. Mathematics , 12(16):2440, 2024. [13] Nicholas I. Fisher . Statistical Analysis of Cir cular Data . Cambridge Uni versity Press, Cambridge, 1995. [14] Arthur Pewse y , Markus Neuhäuser, and Graeme D. Ruxton. Cir cular Statistics in R . Oxford University Press, Oxford, 2013. [15] Savitri Joshi and KK Jose. Wrapped Lindley distribution. Communications in Statistics-Theory and Methods , 47(5):1013–1021, 2018. [16] Mohamed E Ghitany , Barbra Atieh, and Saralees Nadarajah. Lindley distrib ution and its application. Mathematics and Computers in Simulation , 78(4):493–506, 2008. [17] Rama Shanker , Kamlesh Kumar Shukla, Ra vi Shanker , and T ekie Asehun Leonida. A three-parameter Lindley distribution. American Journal of Mathematics and Statistics , 7(1):15–26, 2017. [18] Rama Shanker , Shambhu Sharma, and Ravi Shanker . A two-parameter Lindley distrib ution for modeling waiting and surviv al times data. Applied Mathematics , 4(2):363–368, 2013. [19] Lars V . Ahlfors. Complex Analysis . McGraw-Hill, New Y ork, 3 edition, 1979. [20] Gianluca Mastrantonio, Giov anna Jona Lasinio, Antonello Maruotti, and Gianfranco Calise. Inv ariance properties and statistical inference for circular data. Statistica Sinica , 29(1):67–80, 2019. [21] Peter J Bickel, Bo Li, Ale xandre B Tsybakov , Sara A v an de Geer , Bin Y u, T eófilo V aldés, Carlos Ri vero, Jianqing Fan, and Aad V an Der V aart. Regularization in statistics. T est , 15(2):271–344, 2006. [22] T in Lok James Ng. Penalized maximum likelihood estimator for mixture of von Mises–Fisher distributions. Metrika , 86(2):181–203, 2023. [23] Y uejiao Fu, Jiahua Chen, and Pengfei Li. Modified likelihood ratio test for homogeneity in a mixture of von Mises distributions. Journal of Statistical Planning and Infer ence , 138(3):667–681, 2008. [24] Luis Rodriguez-Lujan, Pedro Larranaga, and Concha Bielza. Frobenius norm regularization for the multiv ariate von Mises distrib ution. International Journal of Intelligent Systems , 32(2):153–176, 2017. [25] Hirotugu Akaike. A new look at the statistical model identification. IEEE Transactions on Automatic Contr ol , 19(6):716–723, 1974. [26] Clif ford M Hurvich and Chih-Ling Tsai. Regression and time series model selection in small samples. Biometrika , 76(2):297–307, 1989. 21 A P R E P R I N T - M A R C H 1 8 , 2 0 2 6 [27] Gideon Schwarz. Estimating the dimension of a model. The Annals of Statistics , 6(2):461–464, 1978. [28] Sreeniv asa Rao Jammalamadaka and T omasz J Kozubo wski. New f amilies of wrapped distributions for modeling ske w circular data. Communications in Statistics-Theory and Methods , 33(9):2059–2074, 2004. [29] Jasper Maars, Gijs v an Dijk, Mark J Dekkers, F Ja vier Hernández-Molina, Federico Andreetto, Francisco J Rodríguez-T ov ar , and W out Krijgsman. New palaeocurrent analysis approach from two-dimensional trough cross-strata using photographs and anisotropy of magnetic susceptibility . The Depositional Recor d , 10(2):298–311, 2024. [30] AJ Underwood and MG Chapman. Experimental analyses of the influences of topography of the substratum on mov ements and density of an intertidal snail, littorina unifasciata. J ournal of Experimental Marine Biology and Ecology , 134(3):175–196, 1989. [31] Richard J Bolton and David J Hand. Statistical fraud detection: A revie w . Statistical Science , 17(3):235–255, 2002. [32] Alejandro Correa Bahnsen, Aleksandar Stojanovic, Djamila Aouada, and Björn Ottersten. Improving credit card fraud detection with calibrated probabilities. In Pr oceedings of the 2014 SIAM International Conference on Data Mining , pages 677–685. SIAM, 2014. [33] Mohammed S Eliwa, Emrah Altun, Mohamed El-Dawoody , and Mahmoud El-Morshedy . A new three-parameter discrete distribution with associated IN AR(1) process and applications. IEEE access , 8:91150–91162, 2020. [34] Christophe Ley , Sla ¯ dana Babi ´ c, and Domien Craens. Flexible models for complex data with applications. Annual Revie w of Statistics and Its Application , 8(1):369–391, 2021. [35] S. Rao Jammalamadaka and Ambar Sengupta. T opics in Circular Statistics . Series on Multi v ariate Analysis. W orld Scientific, Singapore, 2001. 22