arXiv CSopen access

ANGLE: Angular Neural Generative Learning via Engression

ANGLE: Angular Neural Generative Learning via Engression Rajdeep Pathak1,2, * Archi Roy3,*

Tanujit Chakraborty1,2,†

1

LPSM, Sorbonne Université, Paris, France SAFIR, Sorbonne University Abu Dhabi, UAE 3 Indian Institute of Management, Kozhikode, India

arXiv:2607.12833v1 [stat.ML] 14 Jul 2026

2

July 15, 2026

Abstract Circular data, representing angles or directions, are frequently encountered in computer vision, biology, geology, and meteorology. Traditional regression targets the conditional mean, which is often geometrically misleading for circular responses under multimodal, skewed, or asymmetric data structures. To address these limitations, a lightweight deep generative framework, namely ANGLE, is introduced for non-parametric distributional regression on the circle. The full conditional distribution of an angular response, given Euclidean and circular covariates, is learned through a generative map optimized via a generalized circular energy score (GCES) loss. Desirable theoretical properties including the strict propriety of the loss and the rotational equivariance of the estimators are established. Furthermore, both pre- and post-additive noise models are accommodated. A unified toolbox is provided for advancing previously underexplored challenges in circular statistics: extrapolation, sufficient dimension reduction, and conditional distribution equality testing. The framework’s efficacy is demonstrated through extensive simulations and real-world applications. Specifically, the proposal is utilized for object pose estimation from imagery and wind direction prediction, which are integral to surveillance, autonomous vehicles, and energy systems, respectively. Superior predictive performance and robust uncertainty quantification of the proposed method in these tasks are revealed. Software implementation is made publicly available1 via our Python package anglepy2 .

Keywords: Circular data, distributional regression, Engression, extrapolation, pose estimation * Joint first authors. †

Corresponding author: [email protected] https://github.com/PyCoder913/anglepy 2 Documentation: https://anglepy.readthedocs.io 1

1

1

Introduction

Circular data refers to observations that take the form of angles or directions on the unit circle S1 , and they arise naturally in many scientific disciplines. Common examples of circular data include times of day or days of a calendar year, which wrap naturally onto a 24-hour or 12-month clock, and directions or orientations, which are measured as angles on the unit circle. Several interesting applications of circular data are presented in Gill & Hangartner (2010), where the authors model the timing of domestic terrorism events as points on the calendar circle; and in Villarini (2016), where the day of annual peak river flow is modeled as an angle to capture seasonal flooding patterns. Beyond these examples, circular statistical methods have also been extensively employed in domains like biological rhythm discovery (Gorczyca & Sefas 2025, Downs & Mardia 2002) and movement analysis (Ranalli & Maruotti 2020). In most of these applications, the primary statistical objective is to estimate the conditional mean of the angular response, given a set of explanatory variables. However, there has been increasing recognition that a single summary statistic such as the mean often fails to capture the full structure of the conditional relationship, and that the entire conditional distribution is the more relevant inferential target. This need is particularly acute for circular responses where the conditional distribution can be multimodal or asymmetric in ways that render point summaries uninformative and geometrically misleading (Di Marzio, Fensore, Panzera & Taylor 2016). As a motivating example, consider pose estimation, a well-established task in computer vision. Recovering object orientation from imagery or videos is critical for head pose (Prokudin et al. 2018) and gaze direction (Nonaka et al. 2022) detection, autonomous steering (Rasib et al. 2021), projectile tracking, robotic manipulation, augmented reality, and in-orbit operations (see Liu et al. (2024) for a comprehensive review). However, occlusion, poor illumination, low resolution, or symmetries frequently result in multiple plausible poses. Estimating the full conditional distribution

2

aeroplane

aeroplane

bicycle

boat

boat

0°

0°

0°

0°

0°

bus

car

chair

motorbike

motorbike

0°

0°

0°

0°

0°

Ground Truth

ANGLE Prediction

CLR Prediction

95% Prediction Interval

Figure 1: Object pose estimation with uncertainty quantification with the proposed ANGLE. A detailed discussion is provided in Sec. 5.1.

captures this uncertainty; as shown in Fig. 1, ambiguous inputs (e.g., the left-facing motorbike) necessitate wider predictive intervals to bound the true pose, which point estimators like circular linear regression (CLR) fails to provide. An analogous challenge arises in meteorology, where atmospheric dynamics, terrain-induced channeling, land-sea interactions, and competing weather systems create complex, multimodal wind direction distributions (Rad et al. 2022). Because single point predictions obscure these underlying regimes, estimating the full distribution is vital for generating the probabilistic predictions (predictive intervals) required by safety-critical operations like aircraft takeoff and landing (Khattak et al. 2023), wildfire spread prediction (Quill et al. 2019), and maritime navigation (Ning et al. 2025). Motivated by the inadequacy of point summaries in many such applications involving angular responses, we develop ANGLE - Angular Neural Generative Learning via Engression, a deep generative distributional regression framework for circular data. Given covariates X ∈ Rd and an angular response Y ∈ [0, 2π), the method learns the conditional distribution of Y |X through a generative map g(X, ε), where ε ∼ U(0, 1), estimated by minimizing a generalized circular energy score (GCES) loss. While our work draws conceptual inspiration from Shen & Meinshausen

3

(2025), the extension to circular responses is far from straightforward. Notably, standard deep networks cannot be used directly to fit bounded and periodic angular responses. Moreover, the Euclidean energy score fails to accommodate the periodic wrapping effect (0 ≡ 2π), leading to erroneous predictive penalization. Unlike the Euclidean setting, where the energy score is known to be strictly proper, the circular setting lacks a canonical notion of distance, making distributional learning substantially more challenging. To address this issue, we introduce a family of GCES losses based on geodesic and chordal distances on S1 and characterize conditions under which they remain strictly proper. The resulting framework is nonparametric, bandwidth-free, and rotationally equivariant, ensuring that both distributional estimates and derived functionals respect the geometry of the circle. Beyond conditional distribution estimation, the learned generative object naturally enables a range of downstream statistical tasks. We develop methodologies for point prediction with uncertainty quantification - specifically, the estimation of circular mean, median, mode, density, dispersion, and quantiles; sufficient dimension reduction; extrapolation; and testing equality of conditional distributions, accompanied by desirable theoretical guarantees. The practical utility of the proposed framework is demonstrated through applications to object pose estimation from images using the PASCAL3D+ benchmark (Xiang et al. 2014) and wind direction prediction using datasets from India and Germany (Cohen et al. 2025). In summary, this work provides a unified lightweight generative framework for distributional learning on circular data. We summarize our key contributions below. (a) ANGLE framework and theoretical foundations. We propose ANGLE, a deep generative distributional regression method for circular data, establishing population-level guarantee and rotational equivariance for the derived estimators (Theorem 1; Lemma 1). A family of GCES losses is introduced and conditions for their strict propriety are characterized. (b) Comprehensive downstream applications. We develop a unified toolkit for angular responses

4

with theoretical guarantees across all tasks, enabling point prediction with uncertainty quantification (Lemma 1), sufficient dimension reduction (Theorems 2-4), extrapolation (Theorem 5), and conditional distribution equality testing (Algorithm 1). (c) Architectural flexibility and computational efficiency. The framework adapts to both preadditive (covariate noise) and post-additive (response noise) models, remaining significantly more computationally lightweight than existing classical, Bayesian, and modern deep learning alternatives available in circular statistics. (d) Applications to image and tabular data. We demonstrate the method’s efficacy on image data (object pose estimation, SDR, and distribution equality testing) and tabular data (wind direction prediction and extrapolation), alongside extensive simulation experiments. The proposed ANGLE methodology provides model-intrinsic uncertainty quantification for circular responses, the efficacy of which is systematically benchmarked against state-of-the-art alternatives using the circular continuous ranked probability score (CRPS).

1.1

Background

To estimate the conditional distribution of an Euclidean target variable given covariates, Engression was recently introduced as a deep distributional regression framework (Shen & Meinshausen 2025). Engression optimizes a strictly proper scoring rule within an additive noise model, providing theoretical extrapolation guarantees under mild monotonicity assumptions. Recent works have utilized this methodology for problems with temporal (Kraft et al. 2026) and spatiotemporal (Pathak & Chakraborty 2026) data. While it was introduced for Euclidean responses with the main narrative around extrapolation, the original framework did not explore other downstream statistical tasks such as SDR and conditional distribution equality testing. We bridge this gap for the circular data domain by introducing a comprehensive and unified toolkit designed for distributional regression

5

and subsequent downstream analysis of angular data. In this paper, we develop the ANGLE methodology and evaluate its performance across various data modalities, including both imagery and tabular datasets.

1.2

Organization and Notations

The remainder of the paper is organized as follows. Following the notation setup, Sec. 2 outlines the ANGLE methodology and model architecture. Sec. 3 develops downstream statistical tasks: conditional functional estimation (Sec. 3.1), sufficient dimension reduction (Sec. 3.2), extrapolation (Sec. 3.3), and conditional distribution equality testing (Sec. 3.4). Empirical validation via simulations is presented in Sec. 4, followed by applications to object pose and wind direction prediction in Sec. 5. Sec. 7 concludes with future directions. Proofs, software demonstrations, implementation details, and supporting materials are deferred to the Appendix. d

Throughout, = denotes equality in distribution and ι =

√ −1. For a real-valued random variable X,

let ψX (r) = E[exp(ιrX)], r ∈ Z, denote its r-th circular moment; similarly, for a function f , let ψf (r) = E[exp(ιrf (ε))]. The k-th component of a vector u is denoted by u(k) . For a measurable set A, U(A) denotes the uniform distribution on A. For a measure µ, Iµ denotes its support, f |A the restriction of a function f to A, and µ∗ ν the push-forward of a measure ν under µ, defined by µ∗ ν(A) = ν(µ−1 (A)). We write Ptr for the training-data distribution and X ⊂ Rd for the training support. Throughout, S1 is parameterized by [0, 2π), with circular addition and subtraction defined by a ⊕ b = (a + b) mod 2π and a ⊖ b = (a − b) mod 2π. For a random variable X, σ(X) denotes the generated sigma-field. We write atan2(s, c) = arg(c + ιs), where arg : C \ {0} → (−π, π] is the principal argument, C \ {0} denotes the entire complex plane excluding the origin, and (s, c) ̸= (0, 0). ∥ · ∥ denotes the Euclidean norm unless stated otherwise. The notation X ⊥ ⊥Y denotes that the random variable X is independent of Y .

6

2

Methodology

We consider a generative approach to modeling a univariate angular response Y ∈ [0, 2π) given the Euclidean covariates X ∈ Rd . Rather than directly specifying a parametric family for the conditional distribution PY |X , we characterize it implicitly through a generative process as follows. We posit a function g belonging to a prescribed class M, and a real-valued noise variable ε independent of X d

such that Y |X = g(X, ε). For each g ∈ M and x ∈ X ⊂ Rd , we write Pg (·|x) as the distribution induced by g(x, ε), where the randomness arises solely through ε. This generative formulation allows for two desirable features. First, it does not require explicit density specification on the circle, and second, it enables simulation-based approximation of the conditional distribution and sample-based estimation of circular statistics. We further illustrate this in detail in Section 3. Our proposed framework is inspired by the pre-ANM and post-ANM formulations of engression (Shen & Meinshausen 2025), where latent noise is injected either before or after the non-linear link, and the resulting response is mapped onto the circle. We consider the following classes of circular generative models:

Mpost =

(1)

n 





(2)

n 





Mpost =

o

g β ⊤ x + η mod 2π | g : R → R, g ∈ G, β ∈ B , o

g β ⊤ x + η mod 2π | g : R → [0, 2π), g ∈ G, β ∈ B ,

n 



n 



o

⊤ M(1) pre = g β x + η mod 2π | g : R → R, g ∈ G, β ∈ B ,

o

⊤ M(2) pre = g β x + η | g : R → [0, 2π), g ∈ G, β ∈ B .

(2.1)

(1)

Note that Mpost and M(1) pre follow a ‘wrapped modeling’ approach, where the link function g is unbounded, but the final output is modulated by 2π. In contrast, g is bounded in [0, 2π) for models d

(2)

in Mpost and M(2) pre . In each class, the noise η enters through a transformation η = h(ε) for some (1)

(2)

measurable h ∈ H, where we take ε ∼ U(0, 1) to be independent of X. The classes Mpost and Mpost

7





(2) introduce noise after the signal g β ⊤ x is formed (post-ANM), while in M(1) pre and Mpre the noise

is directly introduced in the covariate level (pre-ANM). The parameter space B and the function classes G, H are chosen depending on the statistical objective. For standard prediction tasks, it n

o

suffices to take B = β ∈ Rd : ∥β∥ < ∞ , with G and H classes of smooth functions. However, as we demonstrate in Section 3.3, stronger structural assumptions on the parameter and function spaces are required when the goal is extrapolation. Although the response variable is circular, we model the latent noise variable η as real-valued rather than circular. This is primarily motivated by many real-world phenomena, where the observed angle arises as the circular projection of an underlying Euclidean system subject to ordinary real-valued perturbations. While considering circular noise is also possible in the same framework, we leave the theoretical details for a future work. For simplicity and ease of presentation, we restrict our attention to single-index additive noise structures in the theoretical derivations. The results extend straightforwardly to the multiple-index setting, where β ∈ Rd is replaced by B ∈ Rd×b , with suitable modifications to the identifiability conditions imposed whenever required. An important element of any distributional regression method is a quantitative metric to assess a distributional fit. We introduce the generalized circular energy score (GCES) loss, which, for any direction θ, is given as ◦ 1 ES(P, θ) = E[F (d(θ, Θ))] − E[F(d(Θ, Θ∗ ))], 2

(2.2)

where Θ, Θ∗ are independent copies of circular random variables coming from the distribution P, d(θ1 , θ2 ) = arccos ⟨θ1 , θ2 ⟩ = min{|θ1 − θ2 |, 2π − |θ1 − θ2 |} is the geodesic metric on S1 , and F represents a continuous function on [0, π] such that k(x, y) = F(d(x, y)) is an isotropic kernel on S1 . As we illustrate later in Theorem 1, for the population level guarantee of the proposed methodology, we require (2.2) to be a strictly proper scoring rule. Because we frame our objective

8

as a loss function, we adopt the negatively oriented convention: a scoring rule is strictly proper if its expected penalty is minimized if and only if the predicted probability distribution exactly matches the true underlying distribution. According to Theorem 1.1 and 4.5 of Steinwart & Ziegel (2021), this holds if and only if −F belongs to the class of functions that induce a strictly positive definite kernel on S1 . Examples of such functions are given in Gneiting (2013, Table 1). Note that a simpler version of the circular energy score loss was considered in Grimit et al. (2006), taking F to be the identity function, where the authors showed the metric to be proper. We need the above reformulation (2.2) since strict propriety is crucial for this study. As illustrated in Lemma A.1.2 (Appendix A), it is possible to extend the construction beyond the geodesic metric. Specifically, (2.2) continues to define a strictly proper energy score loss when F is the identity function and d is chosen to be the chordal distance. The population version of the proposed method can then be written as the minimizer of the GCES loss (2.2), i.e., 1 ge ∈ arg min E F (d (Y, g(X, ε))) − F (d (g(X, ε), g(X, ε∗ ))) 2 g∈M 





,

(2.3)

where (X, Y ) ∼ Ptr (x, y) and ε, ε∗ ⊥ ⊥ X are independent copies drawn from the U(0, 1) distribution. The first step now is to establish a population level guarantee for the proposed method as a general distributional approach. That is, we need to show that when the model is correctly specified, (2.3) learns the true conditional distribution for all x ∈ X .

Theorem 1 (Population level guarantee). Assume that there exists a g ∈ M such that g(x, ε) ∼ Ptr (y|x) for all x ∈ X , and that the generalized circular energy score loss is strictly proper. Let ge be any population minimizer in (2.3) that minimizes the GCES loss. Then ge(x, ε) ∼ Ptr (y|x) PX -almost everywhere. Having established the population level guarantee for our methodology, we now move on to finite sample estimation. For each observation in an i.i.d. sample {(Xi , Yi )|i = 1, 2, . . . , n}, we sample 9

n

o

(j)

an ε from the U(0, 1) distribution m times: ϵi |i = 1, 2, . . . , n, j = 1, 2, . . . , m . Based on the finite samples, we define the empirical version of circular engression as: ( −1

b g ∈ arg min n g∈M

n X i=1

"

m

1 X F m

 

d Yi , g



(j) Xi , ϵi



m

X 1 − 2m(m − 1)

j=1

m X

F

   d g

(j) Xi , ϵi



,g



(k) Xi , ϵi



#) .

j=1 k=1,k̸=j

(2.4)

For each fixed g, the objective (2.4) is an unbiased estimator of the population risk (2.3). In particular, the second term in the right hand side of (2.4) is written as a U-statistic over pairs j ̸= k to obtain the unbiased estimate of the self-dispersion term. We parameterize the model g using a broad class of neural networks, as detailed in the subsequent section.

2.1

Model Architecture

The ANGLE framework seamlessly accommodates single-index, multiple-index, and fully connected architectures, allowing the model complexity to be tailored to the task at hand. The framework also supports both pre- and post-ANMs, as well as bounded circular outputs and unconstrained real-valued outputs projected onto [0, 2π), thereby providing a unified implementation of all model classes in (2.1). The architecture utilizes a multilayer perceptron (MLP) backbone. The network’s hidden representation incorporates fully connected transformations, batch normalization, and residual connections, augmented by stochastic noise injection (sampled from either U(0, 1) or N (0, 1)) to enable conditional distribution estimation via the engression methodology (Shen & Meinshausen 2025). To obtain circular outputs, the terminal layer uses one of the three transformations: (I) A bivariate output (u, v) ∈ R2 can be mapped through θb = atan2(v, u) (mod 2π) with a small numerical stabilization to avoid the singular point (u, v) = (0, 0); (II) An unbounded scalar output u ∈ R can be reduced modulo 2π before evaluating the loss and at inference time;

10

(III) An unconstrained scalar output can be mapped to [0, 2π) using a scaled sigmoid transformation. (1)

The modulo transformation implements the wrapped-output variants (Mpre,post ), whereas the atan2 (2)

and scaled-sigmoid transformations implement the bounded circular output variants (Mpre,post ). The distinction between Mpre and Mpost is determined by whether the stochastic noise is injected before or after the nonlinear map. Notably, the ANGLE architecture also accommodates circular covariates by replacing each angular covariate with its trigonometric embeddings before the initial layer. Finally, the network parameters are optimized by minimizing the empirical GCES loss defined in (2.4), using either geodesic or chordal distance as the underlying circular discrepancy. For geodesic distance, the function F is chosen as F(t) = −k(t) with k selected from Gneiting (2013, Table 1). The different variants of the proposed model are compared in an ablation study in Section 6.

3

Statistical Inference with ANGLE

Having established the generative model classes and ANGLE framework for circular responses, we now turn to the theoretical and practical utilities of the proposed framework for statistical inference. A key advantage of the proposed approach is that, once the generative map g(x, ε) has been estimated, the full conditional distribution PY |X=x becomes accessible through forward simulation. Specifically, for any covariate value x, we can draw i.i.d. ε1 , . . . , εm ∼ U(0, 1) and generate approximate conditional samples from the fitted model: Yb1 = g(x, ε1 ), . . . , Ybm = g(x, εm ). These simulated samples can then be used to estimate a wide range of conditional circular summaries without requiring an explicit closed-form expression for the conditional density. In the remainder of this section, we investigate the various inferential tasks enabled by the proposed framework.

11

3.1

Point Prediction of Conditional Functionals

The proposed generative formulation enables straightforward estimation of a broad class of conditional functionals for circular responses, such as the mean and median directions, measures of circular dispersion, conditional quantiles, cumulative distribution functions, and densities. Even under standard parametric circular models, many of these quantities do not admit simple closed-form expressions and are therefore most naturally estimated by sampling from the fitted generative map. For completeness, Appendix E reviews existing approaches to estimating these functionals under classical circular models. In this section, we show how each of these quantities can be estimated directly from the fitted generative map in (2.3). A fundamental requirement for statistical inference on the circle is rotational equivariance, which ensures that inference is invariant to the arbitrary choice of reference direction. Lemma 1 establishes that, under correct model specification, the conditional functionals induced by ANGLE satisfy this property. Specifically, if the response space is rotated by an angle α, then the estimated mean direction, median direction, mode, and density rotate by the same angle α, while the circular dispersion remains unchanged.

Lemma 1 (Rotational equivariance). Let (X, Y ) ∼ Ptr (x, y), ge be the population ANGLE object defined in (2.3), and (X, Y α ) denote the rotated training problem where Y α = Y ⊕ α for any α ∈ R with the corresponding training distribution Pαtr (y, x). Denote by geα the population ANGLE object fitted on (X, Y α ). Assume that the underlying model class M is correctly specified for both problems, i.e., there exist g, gα ∈ M such that g(x, ε) ∼ Ptr (·|x), gα (x, ε) ∼ Pαtr (·|x) for all x ∈ X. (A) Assume that the mean resultant length |Eε [exp(ιge(x, ε))]| > 0 PX -almost everywhere.

12

Define the conditional mean direction of Y |X = x as

µeg (x) = atan2 (Eε [sin ge(x, ε)] , Eε [cos ge(x, ε)]) ,

(3.1)

where atan2 is the two-argument arctangent operator. Then µegα (x) = µeg (x) ⊕ α PX almost everywhere. Hence, the conditional mean direction under the rotated problem is equal to the original conditional mean direction rotated by the same angle. (B) Define the circular median as meg (x) = arg minθ∈[0,2π) {Eε [d(θ, ge(x, ε))]}. Then megα (x) = meg (x) ⊕ α PX -almost everywhere. For notational simplicity, we may assume uniqueness of the median; however, the rotational equivariance property remains valid for each of the medians in the non-unique case. (C) Define the circular dispersion measure (Demni et al. 2026) as σeg (x)

=

Medε [d(ge(x, ε), meg (x))], where Medε denotes the median with respect to the random variable ε. Then σegα (x) = σeg (x) PX -almost everywhere, i.e., the conditional circular dispersion is rotationally invariant. (D) Assume that the conditional density of Y |x is square-integrable on [0, 2π). Following the Fourier inversion theorem for circular variables, we define the truncated conditional density representation of the population engressor as Peg (θ|x) = (2π)−1

PR

r=−R ψe g (r|x) exp(−ιrθ) where ψe g (r|x)

= Eε [exp(ιrge(x, ε))],

and R > 0 is a specified integer-valued truncation level. Then Pegα (θ|x) = Peg (θ ⊖ α|x) for all θ ∈ [0, 2π), PX -almost everywhere. (E) Under the same assumption stated in part (D), define the conditional mode n

and antimode respectively as ῡeg (x) = arg maxθ∈[0,2π) Peg (θ|x) n

o

o

and υeg (x) =

arg minθ∈[0,2π) Peg (θ|x) . Then ῡegα (x) = ῡeg (x) ⊕ α and υegα (x) = υeg (x) ⊕ α. For

13

notational simplicity, we may assume uniqueness of the mode and antimode; however, the rotational equivariance remains valid for each mode in the non-unique case.

Proof of Lemma 1 is provided in Appendix A. Next, following Di Marzio, Panzera & Taylor (2016), we define the conditional circular cumulative distribution function (CDF) estimator as

Feg (θ|x) =

Z θ 0

Peg (u|x)du,

θ ∈ [0, 2π),

where the conditional density Peg (·|x) is defined in part (D) of Lemma 1. This definition depends on a choice of reference direction, which is unavoidable for CDF-type summaries on the circle. Following the same line of work, we define the conditional circular quantile at level 100(1 − α)%, for α ∈ (0, 1), as

qα (x) = arg min Eε [ρα (ge(x, ε), q)] ,

(3.2)

q∈[0,2π)

where the circular check loss function is given by

ρα (θ1 , θ2 ) =

     αδ(θ1 , θ2 ),

if δ(θ1 , θ2 ) ≥ 0,

    −(1 − α)δ(θ1 , θ2 ),

if δ(θ1 , θ2 ) < 0

,

(3.3)

where δ is the signed shortest angular distance on the circle, wrapped to the interval [−π, π). Further properties of the circular CDF and quantile estimators, including rotational equivariance, can be established following similar arguments as in Di Marzio, Panzera & Taylor (2016). Finite sample consistency and convergence properties of these estimators require a separate asymptotic analysis and are beyond the scope of the present work. Notably, the rotational equivariance results established in Lemma 1 rely on the assumption of correct model specification, which is an idealized population-

14

level condition. In practice, however, the expressive power of neural network generators such as ANGLE allows the approximation error to be made arbitrarily small by increasing model capacity. Nevertheless, model misspecification may still arise due to finite network capacity, optimization error, or the choice of latent noise representation.

3.2

Sufficient Dimension Reduction

In many circular regression problems, the covariates are high dimensional, making sufficient dimension reduction (SDR) an important objective. Unlike unsupervised methods such as PCA, which ignore the response and may discard information relevant for prediction (Artemiou & Li 2009), SDR seeks low-dimensional representations of the covariates that preserve all information about the response. Although SDR has recently been studied for general non-Euclidean responses (see Huang et al. (2026) among others), relatively little work has explicitly developed the SDR methodology tailored to circular-valued responses. In this section, we show that ANGLE provides a flexible and principled framework to this end. Given the response-covariate pair (Y, X) ∈ [0, 2π) × X , the task of SDR is to find a sub-sigma field GY |X ⊆ σ(X), such that

Y ⊥ ⊥ X|GY |X ,

(3.4)

i.e., Y and X are independent conditioned on the information contained in GY |X . Note that the choice of GY |X in (3.4) is not unique. The countable intersection of all such choices, termed as the central sigma field, is denoted as G∗Y |X . The existence result for G∗Y |X (Theorem 2) is a direct analogue of the corresponding SDR result for general responses in Lee et al. (2013, Theorem 1); since circular responses take values in a standard measurable space, the proof carries over to the present setting under the same regularity condition. We include the proof in Appendix A for 15

completeness. n

o

Theorem 2. Assume that the family of probability measures PX|Y (·|y)|y ∈ S1 exists and is dominated by a σ-finite measure on X . Then there exists a unique minimal sufficient sub-sigma field G∗Y |X of σ(X) such that Y ⊥ ⊥ X|GY |X , and if G0 is another sub-sigma field satisfying (3.4), then G∗Y |X ⊆ G0 , upto PX -null sets. Typically, GY |X is abstract and difficult to estimate directly. However, in the specific model structures considered in (2.1), we exploit the single-index form Y = f (β ⊤ X, η) with η ⊥ ⊥ X, to convert this infinite dimensional problem of estimating G∗Y |X to a finite dimensional problem of estimation of the true parameter vector β. Under suitable identifiability conditions (see Assumptions 1-4 and Appendix A.2), this parameter β can be uniquely recovered and the corresponding central sigma field G∗Y |X can be determined uniquely up to PX -null sets. We state this formally in Theorem 3 and Theorem 4. 



Theorem 3. Let the true data generating process (DGP) be Y |x ∼ f0 β0⊤ x, ε , where f0 ∈ M takes the structural form f0 (t, ε) = g(t) + h(ε) or f0 (t, ε) = g(t + h(ε)), with β0 being the true DGP parameter. Assume that the ANGLE objective (2.3) is well-defined at the true parameter values, n

o

and that the map t 7→ f0 (t, ·) is injective on β0⊤ x | x ∈ X (i.e., the conditional law is injective with respect to the linear predictor). Then the central sigma field is given by G∗Y |X = σ(β0⊤ X), upto PX -null sets.

Using Theorem 3, we now propose the following methodology for performing SDR with ANGLE.

16

Let f0 (x, ε) = g0 (β0⊤ x, h0 (ε)) be the true DGP. We recover these parameters through ANGLE:

e h) e ∈ (ge, β,



arg min g∈G, h∈H, β∈B

   1  E F d(Y, g(β ⊤ X, h(ε))) − F d(g(β ⊤ X, h(ε)), g(β ⊤ X, h(ε∗ ))) 2 



(3.5)

Theorem 4. Let X be absolutely continuous on Rd and the assumptions in Theorem 3 are true. Let βe be the index parameter obtained by minimizing the objective (3.5). Then βe = cβ0 for some c ̸= 0, and consequently G∗Y |X = σ(βe⊤ X). Proofs for Theorems 2-4 are provided in Appendix A. The recovery of the central sigma field proceeds similarly for both the pre- and post-ANMs in our framework, since both are estimated through the same energy score optimization. However, outside this common estimation strategy, recovering the central sigma-field in pre-additive models is generally more challenging; see Chen (2021) for related developments in the Euclidean setting. Note that the single-index mode can be extended by replacing β ∈ Rd with B ∈ Rd×b , b < d, to obtain a reduced b-dimensional subspace through SDR. Following the dimension reduction step, one can readily fit a stage-two ANGLE on the reduced covariates for downstream tasks, or perform standard post-SDR inferential tasks (Kim et al. 2026).

3.3

Extrapolation

A key advantage of engression over classical distributional regression techniques is its ability to extrapolate, i.e., to produce reliable predictions at covariate values lying outside the support of the training data. In the Euclidean case, extrapolation is facilitated by the fact that the signal g(β ⊤ x) can continue to evolve continuously outside the observed covariate range (Shen & Meinshausen 2025), while the additive noise structure remains preserved. On S1 , however, the response necessarily

17

.

wraps around the circle, so that a continuously increasing signal does not move progressively farther away from the observed responses in any meaningful geometric sense, but instead repeatedly cycles through the same angular values. To make this notion precise, we quantify the degree to which two conditional circular distributions that are indistinguishable on the training support may nevertheless differ at nearby extrapolation points, formalized through the following definition.

Definition 1 (Distributional extrapolability). For a δ > 0, define the distributional extrapolation uncertainty as

UP (δ) =

sup {x∗ : inf{|x∗ −x|:x∈X }≤δ}}

      

sup

P,P∗ ∈P D(P,P∗ )=0 ∀x∈X

   

{D (P(y|x∗ ), P∗ (y|x∗ ))} ,   

where D is a suitable measure of distance between probability distributions and P = {P(y|x)} is a class of conditional distributions. P is said to be distributionally extrapolable upto δ0 if ∃ 0 < δ < δ0 such that UP (δ) = 0.

As identifiability of the model classes serve as a basis for extrapolability (Shen & Meinshausen 2025), we introduce some conditions to study identifiability of (2.1). Specifically, we take G = {g|g satisfies Assumption 2} , and H = {h|h satisfies Assumption 3} . Assumption 1. The training support X is a compact subset of Rd . The covariate X admits a density that is bounded and strictly positive on an open convex set X ′ ⊆ X .

Assumption 2.

(A) The function g is non-constant and smooth.

(B) The support of g in the model classes defined in (2.1) is compact and convex with a non-empty interior. (2) (C) For model classes M(1) pre and Mpre with pre-additive noise, g is strictly monotone.

18

Assumption 3. h : [0, 1] → R is a strictly monotonically increasing, bounded differentiable function with h

  1 2

= 0 and |h(1) − h(0)| < 2π.

Assumption 4. The set of coefficient vectors is B = {β ∈ Rd | β(1) > 0, ∥β∥ = 1}.

The assumptions and identifiability of (2.1) are discussed in Appendix A.2 in detail. We demonstrate in Theorem 5 that all four model classes in (2.1) are extrapolable.

Theorem 5. Let the domain of the functions in G be I ⊂ R, which is an open connected interval. Suppose that Assumptions 1-4 hold and all g ∈ G are real analytic on I. (A) Let Tpost = {β ⊤ x | x ∈ X , β ∈ B}. Assume Tpost ⊂ I. Then there exists a δ0 > 0 such that for all 0 < δ ≤ δ0 , UM(1) (δ) = UM(2) (δ) = 0. post

post

(B) Let Tpre = {β ⊤ x + h(ε) | x ∈ X , ε ∈ [0, 1], β ∈ B, h ∈ H}. Assume Tpre ⊂ I. Then there exists a δ0 > 0 such that for all 0 < δ ≤ δ0 , UM(1) (δ) = UM(2) (δ) = 0. pre

pre

The proof of Theorem 5 is given in Appendix A. Under the identifiability assumptions, the real analytic restriction makes the model classes (locally) distributionally extrapolable. It is useful to contrast our extrapolability result with the engression framework of Shen & Meinshausen (2025, Theorem 1) for classical regression setup, where functional extrapolability of post-ANMs is necessary and sufficient for distributional extrapolability. They operate under a twice-differentiable (C 2 ) assumption which does not inherently guarantee functional extrapolability, and require unbounded noise. In our circular framework, the additive noise is intrinsically bounded because of the bounded transformation h. To achieve extrapolability, we instead constrain our model class G to real analytic functions. By the identity theorem for analytic functions, agreement on the training support dictates global agreement across the connected domain I. Consequently, the identifiability assumptions and real analytic constraint inherently renders the function class functionally extrapolable on the

19

extended domain I, naturally satisfying the criteria for distributional extrapolability. Another distinction concerns the covariate space itself. While Shen & Meinshausen (2025) allow the training support X to be potentially unbounded, we assume that X is compact. Compactness is entirely consistent with practice, since observed covariates are often finite valued and necessarily lie within a bounded region. We assume that the feature space I extends beyond the observed training data as an open connected but bounded domain. This is a natural assumption in essentially all real applications, since covariates are always constrained by physical, geometric, or measurement limitations and therefore cannot vary arbitrarily.

3.4

Testing Equality of Conditional Distributions

In this section, we demonstrate how the ANGLE estimator can be applied for testing of equality of conditional distributions. Consider two independent samples S1 = {(Y1i , X1i )|i = 1, 2, . . . , n1 } (1)

and S2 = {(Y2i , X2i )|i = 1, 2, . . . , n2 } generated from the true conditional distributions PY |X and (2)

PY |X , respectively. We test

(1)

(2)

(1)

(2)

H0 : PY |X = PY |X vs. H1 : PY |X ̸= PY |X .

(3.6)

Tests of this form are widely used in various domains, especially in the context of covariate shift (Qiu et al. 2024) and causal inference (Fan et al. 2024). We take inspiration from the GCA-CDET test of Zheng et al. (2025), adapted with minimal changes to incorporate circular responses. The pseudoscope of the test is given in Algorithm 1. The test reduces (3.6) to an unconditional two(1)

sample problem using the idea that under H0 , synthetic responses generated from PY |X at covariate values from S2 should be indistinguishable from the true responses in S2 . The sample split in Step 2 of Algorithm 1 ensures the independence required for the angular two-sample test ATS, such as Watson’s test (Fornasin et al. 2025) or Kuiper’s test (Anichini et al. 2023). Consistency

20

of the procedure follows from Zheng et al. (2025) under standard smoothness and sample size conditions, provided the test ATS is itself consistent. Notably, several conditional independence testing procedures originally developed for Euclidean data can, in principle, be extended to the circular setting; see Ren et al. (2025), Gao et al. (2026), Zhang et al. (2026) among others. Algorithm 1 Conditional distribution equality test 2 1 , S2 = {(Y2i , X2i )}ni=1 , level α, circular two-sample test ATS. Require: S1 = {(Y1i , X1i )}ni=1 (1) e1 , such that g e1 (x, ε) ∼ PY |X=x 1: Fit the ANGLE network on S1 to obtain the generative model g under correct specification. n2 /2 2: Randomly partition S2 , independently of S1 , into two equal halves S21 = {(Y21,i , X21,i )}i=1 n2 /2 and S22 = {(Y22,i , X22,i )}i=1 . 3: for i = 1, . . . , n2 /2 do 4: Draw noise εi ∼ Pε independently of S1 and S2 . 5: Generate synthetic circular response Ye21,i = ge1 (X21,i , εi ). 6: end for n2 /2 7: Form the synthetic dataset Se21 = {(Ye21,i , X21,i )}i=1 . Under correct specification and suffi(1) (2) (2) (2) ciently large sample sizes, Se21 ∼ PY |X ⊗ PX and S22 ∼ PY |X ⊗ PX . Consequently, under e H0 , the datasets  S21 andS22 are approximately identically distributed. 8: return ATS Se21 , S22 , α and reject H0 if the p-value < α.

4

Simulation Experiments

We conduct a comprehensive simulation study across varied settings to empirically assess the predictive accuracy and extrapolative properties of ANGLE. The following sections detail the data generation procedures, model specification, and empirical results for standard predictive tasks specifically, the estimation of the circular conditional mean and distribution. The simulation study for extrapolation is subsequently presented in Section 4.3.

4.1

Circular Data Generation

We generate synthetic circular data using four approaches derived from two base mapping mechanisms: stereographic projection (arctangent link) and the Projected Normal (PN) distribution 21

(Bhuyan et al. 2026), with linear and nonlinear latent structures. Euclidean and circular covariates are sampled from standard normal and von Mises distributions, respectively, and perturbed with preadditive uniform or Gaussian noise. Each mechanism is evaluated under four configurations varying covariate dimensions and structural parameters, yielding 16 simulation settings in total. Each setting contains 2,000 training and 200 test observations. All generation settings are summarized in Table A.4 (Appendix F).

4.1.1

Approach 1 and 2: Arctangent Link

We define the linear and circular covariates as XL ∼ N (0, 1)n×pL and ΘC ∼ VM(0, 1)n×pC , respectively, where n denotes the number of data points. To simulate pre-additive measurement uncertainty, we introduce independent noise matrices EL and EC with elements drawn from f = X +E and Θ f = (Θ +E ) U(−σ, σ) or N (0, σ 2 ). The perturbed covariates are given by X L L L C C C

(mod 2π). To bypass boundary discontinuities, the noisy circular features are trigonometrically f = [X f , cos(Θ f ), sin(Θ f )] ∈ Rn×(pL +2pC ) . Given embedded, yielding the unified design matrix X L C C

independent coefficient vectors β1 , β2 ∈ RpL +2pC , we compute the base linear combinations zk = f Xβ k for k ∈ {1, 2}. To govern functional complexity, a scalar parameter a introduces coupled

non-linear cross-terms to form the latent signals:

ξ1 = z1 + a tanh(z2 ),

ξ2 = z2 + a tanh(z1 )

We set a = 0 for a strictly linear latent mapping (Approach 1) and a = 4 to induce strong non-linear interactions (Approach 2). Treating ξ1 and ξ2 as orthogonal Cartesian coordinates, we project these unbounded signals onto the circular support [0, 2π)n via a full-quadrant arctangent transformation:

θ = atan2(ξ2 , ξ1 ) mod 2π.

22

This methodology establishes a robust data-generating process that natively accommodates heterogeneous predictors under specified noise and varying structural non-linearities.

4.1.2

Approach 3 and 4: Projected Normal Distribution

f ∈ Rn×p (where p = p + 2p ), our second dataBuilding upon the unified design matrix X L C

generating mechanism derives from the PN distribution (Bhuyan et al. 2026). Given independent f coefficient vectors β1 , β2 ∈ Rp , we compute base linear combinations zk = Xβ k for k ∈ {1, 2}.

To introduce structural non-linearity, the latent Euclidean mean vectors are coupled via a scalar parameter a: µ1 = z1 + a tanh(z2 ),

µ2 = z2 + a tanh(z1 ).

Again, we set a = 0 for a strictly linear latent mapping (Approach 3) and a = 2 to induce non-linear cross-interactions (Approach 4). To capture the intrinsic structural uncertainty of the isotropic PN distribution, bivariate Cartesian targets are sampled as Yk ∼ N (µk , In ), where k ∈ {1, 2} and In denotes the n × n identity matrix. Finally, these unbounded coordinates are radially projected onto the standard circular domain:

θ = atan2(Y2 , Y1 ) mod 2π.

This procedure implicitly models the conditional response as θi |x̃i ∼ PN 2 (µi , I2 ).

4.1.3

Data Generation Settings and Parameters

Details such as number of covariates and the corresponding parameters governing these simulations are given in Table A.4 (Appendix F). For both the arctangent and PN formulations, the true regression coefficient vectors βk expand to dimension pL + 2pC to accommodate the trigonometric embedding. Specifically, the first pL elements govern the linear covariates, while the remaining 23

2pC elements serve as paired weights for the cosine and sine components of the circular features. For high-dimensional configurations (Settings 1.4, 2.4, 3.4, and 4.4) where explicit specification is intractable, the 125-dimensional coefficient vectors are drawn independently from U(−10, 10). Fig. 2 illustrates the first setting (1.1) with linear arctangent link and 2 Euclidean covariates; analogous plots for the remaining settings are omitted for brevity. Scatter Plot with Cyclic Colormap

Training Data Simulated from Setting 1.1 Quiver Plot (Arrows indicate direction)

2

4 3

1 0 1

1 0 1

3

1

0

1

Feature 1 (Linear)

2

3

0

2

1 0

1

2

3

2

1

0

1

Feature 1 (Linear)

2

3

0

3

4

r)

2

3

a ne

3

/2

/2 0

(Li

2

2 3 /2

e2 tur

2

3

3 /2

2

Fea

/2

Feature 2 (Linear)

3 /2

2

Circular Response ( )

Feature 2 (Linear)

3

Combined 3D Scatter and Quiver

2

Circular Response ( )

4

1 2 1 0 (Linear) 3 2 eature 1

3

F

Figure 2: 2000 training points generated from setting 1.1; Scatter plot (left), quiver plot with arrows depicting the direction of the observations in 2D (middle), and the same in 3D (right).

4.2

Model Specification and Simulation Results

All simulated datasets are fitted using a 3-layer (pre-ANM) ANGLE network with hidden dimension 100 and noise dimension 64. To evaluate the effect of geometry on distributional learning, we train separate models using geodesic and chordal distance-based energy score losses. Training is performed for 500 epochs with learning rate 0.05 and no batch normalization. Performance is assessed using mean absolute angular deviation (MAAD), circular mean directional error (CMDE or cosine loss), and circular continuous ranked probability score (CRPS; see Appendix C for definitions). Since ANGLE is generative, inference is done based on 100 samples per test point. Further, the point metrics are computed from the circular mean, while CRPS uses the full predictive distribution. We compare against the baseline models listed in Table 1, whose architectural and implementation details are deferred to Appendix B. We exclude Bayesian implementations due to their substantially higher computational cost; however, Bayesian von-Mises Quasi Processes

24

(vMQP; Cohen et al. (2025)) is included for comparison during real data analysis. Table 1: Comparison of baseline models for circular data. Columns indicate native support for linear/circular covariates, input/output ranges, intrinsic uncertainty quantification (UQ), computational efficiency, and evaluation tasks (simulations, object pose estimation (OPE), and wind direction prediction (WDP)). Model

Linear

Circular

Input/output

covariates

covariates

range

UQ

Lightweight

Simulations

OPE

WDP

Circular Linear Regression (CLR; Presnell et al. (1998))

✓

✗

(−π, π]

✗

✓

✓

✗

✓

Nonparametric Kernel CLR (Taylor 2012)

✓

✓

(−π, π]

✗

✓

✓

✗

✓

SCBR-NLL (Hassanzadeh 2021)

✓

✗

[0, 2π)

✗

✓

✓

✗

✗

SCBR-MCPE (Hassanzadeh 2021)

✓

✗

[0, 2π)

✗

✓

✓

✓

✓

Lifted EP-MGvM (Wu Navarro 2019)

✓

✗

(−π, π]

✗

✓

✓

✓

✓

MCR2 (Jha & Biswas 2017)

✗

✓

[0, 2π)

✗

✗

✓

✓

✓

vMQP (Cohen et al. 2025)

✓

✗

(−π, π]

✓

✗

✗

✗

✓

Parametric Mixture Models (Prokudin et al. 2018)

-

-

[0, 2π)

✓

✗

✗

✓

✗

ANGLE (Proposed)

✓

✓

[0, 2π)

✓

✓

✓

✓

✓

(A) Donut Plot

(B) Spokeplot N (0 rad)

Clockwise Deviation

Anticlockwise Deviation

True (Observed)

Predicted

Figure 3: (A) Donut plot and (B) Spokeplot of ANGLE-Chordal predictions on simulated data, setting 1.1. In (A), points near the outer circumference (radius = 2) indicate good fit; points inside the unit disc indicate poor fit. Filled and empty circles denote clockwise and anticlockwise deviations, respectively.

Table A.5 (Appendix F) reports the mean and standard deviation of the evaluation metrics over 200 independent trials. Across most settings, the energy score based on the chordal distance achieves the best performance. Predictions by ANGLE on the simulated test set of setting 1.1 is presented through (A) donut plot (Jha & Biswas 2017) and (B) spokeplot in Fig. 3. For fairness, circular covariates were represented by their sine and cosine embeddings for baseline methods that do not natively support directional predictors. Overall, ANGLE consistently outperforms all 25

ANGLE-Chordal

3.983

ANGLE-Geodesic

4.529

Circular Linear Regression

6.854

Kernel CLR

30.806

Lifted EP-MGvM

33.954

SCB-NLL

37.268

SCB-MCPE

37.501

MCR2

42.724 0

5

10

15

20

25

30

35

40

Mean CRPS across all simulation settings (lower is better)

Figure 4: Comparison of model performance on simulated datasets using mean CRPS across all 16 simulation settings. competitors (see Fig. 4) in terms of CRPS across all 16 simulation settings. While circular linear regression (CLR) remains competitive under linear data-generating mechanisms (Approaches 1 and 3), ANGLE substantially outperforms it in the nonlinear settings (Approaches 2 and 4), highlighting the advantage of neural networks in capturing complex nonlinear relationships.

4.3

Simulations for Extrapolation

Theorem 5 establishes the extrapolability of the model classes (2.1) under some identifiability conditions and the assumption that the generative function g is real analytic. To achieve extrapolation, (1)

(2)

we qualitatively evaluate the implementations of Mpre, post , and Mpre, post with scaled sigmoid output. (1)

For Mpre, post , by defining g : R → R as an unbounded scalar output and applying the modulo operator post-hoc during inference, the network is permitted to grow linearly outside the training support, which maps to a continuous ‘wrapping’ behavior on the circular manifold. Theoretically, M(1) pre is not fully identifiable due to the 2kπ equivalence class, but is extrapolable under the real analytic assumption (Theorem 5). We demonstrate extrapolability empirically using simulated pre-ANM data with one predictor in Fig. 5. 5000 training covariates X ∈ R are sampled uniformly from a bounded support and perturbed 26

Figure 5: Simulation results for extrapolation with ANGLE: unbounded output with 2π modulation (first two columns) and output bounded in [0, 2π) with scaled sigmoid (last two columns). As a baseline, the kernel CLR is considered. with pre-additive Gaussian noise, η ∼ N (0, 1). These noisy inputs are passed through two base functions - linearly asymptotic g = Softplus, and sub-linear g = Log, and projected onto the circular manifold to generate the angular response Y . The true function g is chosen to be strictly monotonic, i.e., not wrapping around the circle. To establish a ground-truth extrapolation baseline, the exact conditional circular mean is numerically approximated over an extended evaluation grid via Monte Carlo integration utilizing 10,000 noise samples per grid point. The training and 27

evaluation supports for the Softplus base function are defined as [0, 2] and [−2, 5], respectively; for the logarithmic configuration, these boundaries are restricted to [0, 1.5] and [−1, 3]. We visualize the extrapolation fits by training ANGLE over the simulated datasets under both pre-ANM (panels (A)-(H)) and post-ANM (panels (I)-(P)) formulations. Within these plots, the training samples are denoted by gray scatter points bounded within the training support. To approximate the model’s predicted conditional distribution, we generate 100 samples from the trained ANGLE at each evaluation coordinate, illustrated by the blue scatter clouds. The theoretical ground-truth conditional mean and ANGLE’s predicted circular mean are superimposed as solid red and dashed navy lines, respectively. The dash-dotted indigo line represents the fit of kernel CLR, which we consider as a baseline. It is visually evident that although the baseline fit (indigo) closely approximates the true conditional mean (red) within the training support, it asymptotes outside the support and fails to extrapolate. (1)

Conversely, the proposed ANGLE with the Mpost,pre -like architecture smoothly preserves the structural trend of the data, providing a better approximation of the true conditional mean, both inside and beyond the training support (first two columns in Fig. 5). Finally, extrapolation capabilities (2)

of the Mpost,pre -like architecture with scaled sigmoid is presented in the last two columns in Fig. 5, where the predicted circular mean (dashed navy line) closely aligns with the true conditional mean (red line).

5

Real Data Applications

To empirically validate the practical utility of our proposed framework for prediction with uncertainty quantification, sufficient dimension reduction, extrapolation, and testing equality of conditional distributions, we conduct comprehensive evaluations across both image and tabular data modalities. For prediction, we address two distinct tasks: object pose estimation from images and location-based

28

wind direction prediction. For the former, we utilize the PASCAL3D+ benchmark dataset (Xiang et al. 2014) comprising 36,697 images across 12 object categories with annotated ground-truth orientations. The proposed architecture for this task combines a fine-tuned visual encoder with a class-specific ANGLE head for targeted object pose estimation. For wind direction prediction, we evaluate the model using two distinct meteorological datasets: the German wind direction dataset previously analyzed by Cohen et al. (2025) and district-level mean wind directions in India for December 2025. The Indian wind dataset is also used to showcase extrapolability of ANGLE. The proposed models are implemented in fully dense mode with atan2 in the final layer. Point predictions are obtained by computing the circular mean of samples simulated from the estimated ANGLE object. To ensure the statistical robustness of our findings, all experimental setups are executed over 50 independent trials, with performance metrics reported as the empirical mean and standard deviation.

5.1

Object Pose Estimation

Pose estimation is a fundamental task in computer vision, which deals with predicting the direction towards which an object in a given image or video is pointing. It is an important building block in systems that aim to understand complex scenes and arises in many applications such as estimation of head pose (Prokudin et al. 2018), gaze direction (Nonaka et al. 2022), object class orientation, and so on. We evaluate our approach for pose detection of objects on the PASCAL3D+ dataset (Xiang et al. 2014), a rigorous benchmark for pose estimation comprising 12 object categories (see Table 2) sourced from PASCAL VOC 2012 and ImageNet (Deng et al. 2009). The dataset provides 3D annotations including angles for azimuth, elevation, and distance. We model the azimuth, which denotes the horizontal rotation angle defining the object’s orientation relative to a 0◦ reference. Consistent with prior literature (Prokudin et al. 2018), we utilize the ImageNet and PASCAL training images for training, ImageNet validation images for validation, and PASCAL 29

validation images for testing (data split statistics are provided in Table 2). Notably, we analyze the raw images without applying any prior preprocessing and data augmentation. Table 2: Number of images per class in train, validation, and test splits of the PASCAL3D+ dataset. aeroplane

bicycle

boat

bottle

bus

car

chair

diningtable

motorbike

sofa

train

tvmonitor

Train

1467

1094

1843

1515

865

4043

1996

1533

1019

Validation

979

658

1262

761

536

2780

514

1153

613

Test

484

380

491

733

320

1173

1449

374

Total class-wise

2930

2132

3596

3009

1721

7996

3959

3060

Total

1150

989

1041

18555

707

646

623

11232

376

387

329

414

6910

2008

2244

1964

2078

36697

Since the proposed ANGLE expects data in the form of arrays or tensors, feature extraction is performed using state-of-the-art visual encoders: Inception-v3 (Szegedy et al. 2016) and ConvNeXt (Liu et al. 2022), initialized with pre-trained ImageNet weights from torchvision3 . As classificationbased pre-training encourages rotation invariance, we fine-tune these models using PASCAL3D+ azimuth annotations to extract rotation-sensitive representations. Hence, the object pose estimation task is divided into two stages: (1) fine-tuning pre-trained visual encoder to perform as feature extractor and (2) fitting an ANGLE model on the extracted features. For the first stage, we replace the standard classification blocks of the image models with a custom multi-class biternion regression head. To circumvent the periodic discontinuity of angular targets at 0◦ /360◦ , we encode the groundtruth azimuth angle θ as a continuous 2D directional vector, or biternion, v = [cos θ, sin θ]⊤ . The modified final linear layer maps the extracted image features to a tensor of size C × 2, where C = 12 is the total number of target object classes. A feature-wise L2 -normalization is then applied to the outputs to ensure the predicted 2D vectors lie on the unit circle. During fine-tuning, we employ a conditional routing strategy: the network extracts the predicted biternion v̂c corresponding strictly to the object class c. The entire network is fine-tuned end-to-end by minimizing the cosine distance loss between the predicted and target biternions: L = 1 − v̂c · v. Optimization is performed for 20 epochs using Adam with a learning rate of 1 × 10−4 and a weight decay of 0.01. During training, 3

https://github.com/pytorch/vision

30

83%

75%

17%

51.

68.08%

48.08%

69.

33%

00%

49.50%

92%

64.

42.91

40 32.97

Model

37.96

36.04 31.76

32.53 33.52

33.90

30 24.4324.34

20

10

66.

50%

65.

Mean Circular CRPS (Lower is better)

69.

(B) Mean Circular CRPS over 12 Classes 44.38

53.

58% 70.

53.25%

(A) Mean Accuracy /6 over 12 Classes

23.06

Mixture-Inception Mixture-DenseNet Mixture-MobileNet Inceptionv3-SCBR Inceptionv3-LEP-mGvM Inceptionv3-MCR2 ConvNeXt-SCBR ConvNeXt-LEP-mGvM ConvNeXt-MCR2 Inception-v3-SDR-ANGLE Inception-v3-ANGLE ConvNeXt-ANGLE

0

Figure 6: Mean (A) Accuracyπ/6 and (B) CRPS obtained by the proposed and baseline models across all 12 classes on the PASCAL3D+ test set.

we utilize the validation images for early stopping. A single global encoder processes images across all 12 classes, generating 2048-dimensional (Inception-v3) or 768-dimensional (ConvNeXt) feature vectors. In the second stage, these features are subsequently passed to the ANGLE framework. While feature extraction is global, we train independent, class-specific models to perform targeted pose estimation. We compare our approach against the parametric mixture models of Prokudin et al. (2018) and standard baselines reported in Table 1. CLR and its kernel variant are excluded due to their inability to scale to high-dimensional regimes where features outnumber samples. Evaluation metrics include MAAD, CMDE, circular CRPS, Accuracyπ/6 , and median error (MedErr; see Appendix C for the mathematical formulations). Quantitative class-wise comparisons are presented in Table A.7 (Appendix F), and the test-set predictions with 95% intervals are visualized in Fig. 1. The prediction intervals are constructed by extracting the 0.025 and 0.975 circular quantiles using (3.2) from 100 predictive samples from the estimated conditional distribution of azimuth angles. Fig. 6 illustrates the mean scores obtained by the models on Accuracyπ/6 and CRPS across 12 classes. The proposed model with ConvNeXt as the visual encoder (ConvNeXt-ANGLE) achieves the best mean scores across all 12 classes for all metrics, followed by Inception-v3-ANGLE.

31

(A) Metric: CRPS

(B) Metric: MedErr

(C) Metric: MAAD

Inception−v3−SCBR − 11.33

Inception−v3−SCBR − 10.67

Inception−v3−SCBR − 9.75

ConvNeXt−SCBR − 10.33

ConvNeXt−SCBR − 9.58

Mixture−MobileNet − 9.58

ConvNeXt−LEP−mGvM − 8.67

ConvNeXt−LEP−mGvM − 8.33

Mixture−Inception − 8.67

Mixture−MobileNet − 7.79

ConvNeXt−MCR2 − 7.83

ConvNeXt−SCBR − 8.42

ConvNeXt−MCR2 − 7.75

Inception−v3−MCR2 − 7.50

Mixture−DenseNet − 8.42

Inception−v3−MCR2 − 7.38

Mixture−MobileNet − 7.00

ConvNeXt−MCR2 − 6.08

Inception−v3−LEP−mGvM − 6.75

Inception−v3−LEP−mGvM − 6.92

Inception−v3−MCR2 − 5.92

Mixture−DenseNet − 5.83

Mixture−Inception − 6.58

ConvNeXt−LEP−mGvM − 5.83

Mixture−Inception − 5.25

Mixture−DenseNet − 6.08

Inception−v3−SDR−ANGLE − 5.00

Inception−v3−SDR−ANGLE − 2.92

Inception−v3−SDR−ANGLE − 3.29

Inception−v3−LEP−mGvM − 4.83

Inception−v3−ANGLE − 2.33

Inception−v3−ANGLE − 2.21

Inception−v3−ANGLE − 3.33

ConvNeXt−ANGLE − 1.67

ConvNeXt−ANGLE − 2.00

ConvNeXt−ANGLE − 2.17

0

2

4

6

8

10

12

Mean Ranks of Models (All 12 Classes)

0

2

4

6

8

10

12

Mean Ranks of Models (All 12 Classes)

0

2

4

6

8

10

12

Mean Ranks of Models (All 12 Classes)

Figure 7: MCB test results to statistically evaluate model performances in object pose estimation using (A) CRPS, (B) MedErr, and (C) MAAD metrics. The y-axis in each plot denotes the models and their average ranks (lower the better) corresponding to that metric.

To statistically test the performance of the proposed models against baselines on this task, we conduct the non-parametric Multiple Comparisons with the Best (MCB) test (Koning et al. 2005). The results presented in Fig. 7 indicate that ANGLE is statistically superior to the baselines and achieve the lowest average ranks with respect to CRPS, MedErr, and MAAD metrics.

5.2

Wind Direction Prediction

The modeling of wind direction is a well-established problem in circular statistics. As a second application, we leverage spatial coordinates to predict wind direction at unobserved locations across two meteorological datasets. Predicting wind direction at unobserved locations is useful in several applications, including pollutant dispersion tracking, extreme weather forecasting, and the optimization of renewable wind energy infrastructure. The German dataset, sourced from the Deutscher Wetterdienst (DWD)4 , utilizes a single cross-sectional observation from 260 stations on the final day of a calm weather period (Cohen et al. 2025). Observations for this data are in the range (−π, π], which is converted to the range of [0, 2π) for modeling with ANGLE. However, modeling solely using spatial coordinates is rather simplistic and the procedure can be made more 4

https://www.dwd.de

32

robust by incorporating external covariates that affect wind direction, such as temperature, pressure, and humidity. Another dataset considered in this study is the Indian wind data5 , which comprises the mean wind direction across 737 districts for December 2025. The observations lie in [0, 2π), with the northern direction indicating 0 radians and angles increasing clockwise. The distribution of wind directions for both datasets is presented as rose plots in Fig. A.11 (Appendix F). The German and Indian datasets were randomly partitioned into 208/52 and 552/185 train/test splits, respectively (Fig. A.12 in Appendix F). ANGLE model comprising 2 layers (128 hidden units, 64 noise dimension) was trained for 500 epochs on the German data, and a 3-layer architecture (100 hidden units, 64 noise dimension) was trained for 800 epochs on the Indian data. We benchmark predictive performance against the baselines listed in Table 1. The proposed approach attains superior CRPS scores across both datasets (Fig. 8), signifying better modeling of the target conditional distribution. Quantitative results are summarized in Table A.6 (Appendix F). Notably, the Indian wind dataset presents a more challenging modeling task than the German dataset due to its larger volume and wider observational variance. As a result, although baseline models yield comparable results on the German data, the proposed ANGLE demonstrates a significant performance advantage on the Indian data. Furthermore, our proposal proves to be a remarkably computationally efficient conditional distribution estimator; while Bayesian vMQP requires 21,677 seconds for its sampling procedure on the German dataset, our models complete training and inference in just 1.96-2.2 seconds (close to CLR), representing a 99.99% decrease in computational overhead. 5

https://iced.niti.gov.in/energy/fuel-sources/wind/speed

33

(A) Circular CRPS: Indian Wind Data

(B) Circular CRPS: German Wind Data

80

20

70 18

60 50

16

40 14

30 20

al

sic de AN

GL

E-C

ho rd

QP GL

Model

AN

M

R2

E-G

eo

vM

MC

BR SC

P-m Gv

AN

GL

LE

LR

R

el C

CL

al

sic de

E-C

ho rd

QP GL

Model

AN

M

R2

E-G

eo

vM

MC

BR

LR

P-m Gv

LE

SC

el C

Ke rn

CL

R

12

Ke rn

Circular CRPS (Lower is better)

90

Figure 8: CRPS obtained by the models on (A) Indian and (B) German wind datasets over 50 independent runs.

5.3

ANGLE for Sufficient Dimension Reduction

We demonstrate the practical utility of sufficient dimension reduction (SDR) using ANGLE for the PASCAL3D+ dataset. In object pose estimation, feature extractors such as Inception-v3 generate high-dimensional representations (here, d = 2048). To mitigate redundancy and isolate the predictive subspace that is ‘sufficient’ to explain the angular responses, we estimate a d × b index parameter βb (b ≪ d) through ANGLE, following Theorem 4. Setting b = 5, we project the original features to obtain the reduced 5-dimensional covariates X ′ = βb⊤ X. A stage-two model, termed Inception-v3-SDR-ANGLE, is subsequently trained on X ′ . Despite a 99.7% reduction in feature dimensionality, this approach yields highly competitive predictive performance and consistently outperforms all baselines (see Table A.7 and Fig. 6), achieving the third lowest average rank in the MCB test (Fig. 7) based on MedErr and CRPS. The estimated 5-dimensional central subspace sufficiently captures the structural information required for azimuth prediction. Beyond computational efficiency, SDR naturally facilitates visual diagnostics and exploratory analysis.

34

5.4

ANGLE for Extrapolation

To evaluate the extrapolation capabilities of the ANGLE framework, we implement a targeted spatial block-splitting strategy on the Indian wind dataset. Specifically, the training set is restricted to observations whose latitude and longitude are both below their respective 85th percentiles, with the remaining data allocated to the test set. This configuration establishes an out-of-distribution (OOD) spatial block, confining the training data to a bounded south-western region of the coordinate space and forcing the model to extrapolate toward the northern and eastern peripheries. Hence, the northern and north-eastern states are isolated as the held-out test set (Fig. A.13 in Appendix F) with 517 training and 220 test observations. Table 3: Model performance on Indian wind dataset under the extrapolation regime. The best results are highlighted. MedErr, CRPS, and MAAD are reported in degrees. Metric

vMQP

CLR

Kernel CLR

SCBR

LEP-MGvM

MCR2

ANGLE-Geodesic

ANGLE-Chordal

Accuracy

0.20 ± 0.02

0.36 ± 0.00

0.22 ± 0.00

0.26 ± 0.01

0.22 ± 0.00

0.21 ± 0.03

0.30 ± 0.04

0.31 ± 0.02

MedErr

89.14 ± 1.45

52.96 ± 0.41

78.72 ± 0.00

60.75 ± 6.02

73.00 ± 0.00

95.22 ± 7.90

51.27 ± 9.56

45.74 ± 1.55

CRPS

57.84 ± 0.34

59.16 ± 0.25

78.62 ± 0.00

65.15 ± 3.68

74.43 ± 0.00

91.96 ± 3.61

38.84 ± 4.58

37.81 ± 0.67

MAAD

88.71 ± 0.61

59.16 ± 0.25

78.62 ± 0.00

65.15 ± 3.68

74.43 ± 0.00

91.96 ± 3.61

58.25 ± 7.18

55.90 ± 0.87

CMDE

0.99 ± 0.01

0.60 ± 0.00

0.85 ± 0.00

0.65 ± 0.06

0.80 ± 0.00

1.03 ± 0.05

0.56 ± 0.11

0.52 ± 0.01

The predictive performance under this extrapolation paradigm is summarized in Table 3. Notably, some conditions required for the theoretical guarantees of extrapolability established in Section 3.3, such as the monotonicity assumption, are violated in this empirical setting; broad-scale meteorological phenomena over the Indian subcontinent dictate complex circulating wind patterns that do not vary monotonically with latitudinal or longitudinal progression. Although the scores substantially degrade as compared to the prediction (interpolation) regime (Table A.6), the results demonstrate that ANGLE substantially outperforms the baseline models on this OOD task.

35

(B) Trigonometric MANOVA p-Values

(C) MANOVA: Conditional Equality

boat 2.0e-17 5.0e-01

bottle 1.0e-100 1.0e-100 1.0e-100

bottle 1.0e-100 1.0e-100 2.4e-98

10 16

bottle

bus 6.3e-08 3.1e-18 2.6e-22 9.6e-93

bus 3.1e-05 2.0e-07 3.2e-15 4.7e-34

10 28

bus

car 2.6e-16 1.6e-27 1.5e-43 1.0e-100 5.3e-02

car 4.8e-18 4.0e-23 7.3e-52 3.7e-33 9.2e-02

10 40

car

chair 5.0e-07 1.8e-25 1.1e-33 1.0e-100 1.5e-06 8.6e-33

chair 6.6e-07 6.0e-30 2.1e-39 1.8e-55 2.0e-04 7.7e-31

10 52

chair

diningtable 3.3e-36 1.0e-60 2.4e-25 1.0e-100 3.5e-13 7.2e-21 1.3e-25

diningtable 1.2e-41 9.5e-70 4.6e-37 5.6e-57 2.0e-11 8.7e-04 3.0e-37

motorbike 1.9e-31 2.6e-05 4.9e-02 1.0e-100 9.5e-12 3.4e-39 9.7e-17 3.8e-51

motorbike 3.9e-27 2.7e-04 1.1e-01 1.0e-100 4.9e-05 1.2e-24 7.4e-17 7.4e-57

sofa train

sofa

train

motorbike

diningtable

tvmonitor car

10 100

motorbike

chair

sofa

10 88

train

motorbike

diningtable

car

chair

bus

boat

bottle

bicycle

aeroplane

sofa

train

motorbike

diningtable

car

tvmonitor 8.5e-18 4.5e-89 1.7e-95 2.2e-25 3.9e-16 5.5e-47 3.3e-04 9.9e-01 3.2e-77 3.0e-04 8.1e-07

chair

tvmonitor 3.7e-23 2.4e-75 2.2e-77 1.6e-60 2.2e-22 1.2e-37 3.2e-11 1.4e-07 1.6e-74 9.8e-13 4.0e-11 bus

train 2.1e-06 2.8e-18 3.8e-44 1.2e-22 9.5e-18 2.6e-11 4.8e-09 1.9e-01 4.0e-18 5.3e-02

boat

train 7.5e-13 6.7e-35 4.8e-54 5.4e-61 3.1e-20 4.9e-13 5.9e-26 1.5e-10 1.7e-33 3.3e-05

bottle

sofa 1.3e-40 3.8e-61 3.1e-64 5.8e-37 1.1e-17 3.6e-19 6.0e-01 1.6e-04 2.6e-71

bicycle

sofa 4.6e-36 8.5e-54 2.7e-50 2.1e-73 1.4e-20 6.4e-45 3.1e-21 5.1e-22 1.0e-57

diningtable

bus

10 76

boat

10 64

boat

bottle

p-value (Log Scale)

boat 3.9e-17 3.0e-02

aeroplane

Equal (p 0.05) Different (p < 0.05)

bicycle 10 4

bicycle

bicycle 3.0e-27

aeroplane

(A) Kuiper's Test p-Values bicycle 7.1e-34

Figure 9: Testing pairwise equality of conditional distributions of 12 objects in PASCAL3D+ training set using ANGLE.

5.5

ANGLE for Testing Equality of Conditional Distributions

To evaluate conditional distribution equality using ANGLE, we tested whether azimuth distributions given Inception-v3 features are preserved across the 12 PASCAL3D+ categories. For each class pair, we tested the null hypothesis (H0 ) of identical conditional pose distributions. ANGLE was fitted on the larger class (S1 ). The smaller class (S2 ) was bisected: one half generated synthetic azimuths via the fitted model, while the other provided the observed azimuths. These samples were compared using Kuiper’s test and Trigonometric MANOVA. The latter serves as our primary testing criterion due to its superior power and Type-I error control for circular data (Landler et al. 2021). Fig. 9 summarizes the exhaustive pairwise tests, with the symmetric upper triangle and diagonal masked. Panels (A) and (B) visualize p-values for the Kuiper and MANOVA tests, respectively; lighter shades denote severe generative mismatches (p ≪ 0.05). Panel (C) binarizes the MANOVA results at α = 0.05: green cells indicate statistically indistinguishable distributions (failing to reject H0 ), while red cells denote significant divergence. Results indicate that structurally analogous objects (e.g., car/bus, sofa/chair) share equivalent conditional distributions. However, the prevalence of rejections highlights a high degree of class-specificity, confirming distinct feature-to-orientation mappings across most pairs. This divergence might be due to varying geometric symmetries, photographic framing biases (e.g., frontal cars versus lateral aeroplanes), or fundamental differences 36

in the underlying neural feature representations.

6

Ablation Study

Implementation of the proposed ANGLE framework affords significant architectural flexibility. Within the GCES loss (2.2), the distance metric can be specified as either chordal or geodesic; for the latter, the function F may be the negative of any of those outlined in Gneiting (2013, Table 1) to ensure strict propriety. Architecturally, the network supports both single-, multiple-index, and fully dense modes, alongside three output head formulations: unbounded (modulated by 2π during inference), scaled sigmoid, or a biternion representation utilizing atan2. To compare the contribution of these components in model performance, we conducted an ablation study on the Indian wind dataset utilizing a distinct train-test split from Sec. 5.2. By varying the distance metric, network mode, kernel, and output head, we compared 60 different model variants. Full results are detailed in Table A.8 (Appendix F), and the top 10 configurations based on CRPS are highlighted in Fig. 10. The results demonstrate that the fully dense mode consistently surpasses the single-index mode, capitalizing on the superior expressive power of fully connected layers. Moreover, models employing an unbounded output modulated by 2π exhibit high instability, characterized by large standard deviations. We observe that the dense variant combining the chordal distance with a scaled sigmoidal output head achieves the highest overall performance, while the C 2 -Wendland kernel proves most effective among models utilizing the geodesic distance.

7

Discussion

In this paper, we proposed a deep generative framework called ANGLE for distributional learning on circular data, that combines the flexibility of deep learning with the inferential capabilities of distributional modeling. The framework is nonparametric, computationally efficient, produces 37

Ablation Study on Indian Wind Data: Top 10 Model Variants based on CRPS Chordal-Dense-Sigmoid

16.948 ± 0.80 18.180 ± 1.02

Geodesic-Dense-C 2-Wendland-Sigmoid Geodesic-Dense-C 4-Wendland-Sigmoid

18.460 ± 1.30

Chordal-Dense-atan2

18.578 ± 0.95

Geodesic-Dense-Sine Power-Sigmoid

18.786 ± 1.20

Geodesic-Dense-C 2-Wendland-atan2

20.464 ± 1.22

Geodesic-Dense-C 4-Wendland-atan2

20.599 ± 1.15

Geodesic-Dense-Multiquadric-Sigmoid

20.698 ± 1.24

Chordal-Dense-Modulated

20.884 ± 4.08

Geodesic-Dense-Sine Power-atan2

21.604 ± 1.03 0

5

10

15

CRPS (Lower is Better)

20

25

Figure 10: Top 10 ANGLE model variants based on CRPS in the ablation study. The mean and standard deviation of the CRPS scores are reported over 50 independent training and evaluation runs. rotationally equivariant estimators, and offers extrapolation along with uncertainty quantification while naturally accommodating both Euclidean and circular covariates. Through simulations and real data applications to object pose estimation and wind direction prediction, we demonstrated its ability to recover complex conditional distributions and provide model-intrinsic uncertainty quantification. Beyond prediction, ANGLE yields a unified toolbox for supervised learning with angular responses, including sufficient dimension reduction and testing equality of conditional distributions. Several directions remain open for future research. Methodologically, extending the method for probabilistic forecasting of circular time series (Harvey et al. 2024) and space-time data (Wang & Gelfand 2014, Pathak & Chakraborty 2026) and to more general non-Euclidean response spaces, such as spheres, tori, and manifolds, would considerably broaden its scope. On the applied side, performance in computer vision tasks can be improved through stronger feature extractors and data augmentation. Extending the framework to real-time video analysis would enable applications such as dynamic gaze estimation (Nonaka et al. 2022) and trajectory tracking, further broadening the scope of ANGLE.

38

Appendix Table of Contents A Mathematical Proofs

40

A.1 Auxiliary Lemmas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 A.2 Identifiability Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 A.3 Proof of Theorem 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 A.4 Proof of Lemma 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 A.5 Proof of Theorem 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 A.6 Proof of Theorem 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 A.7 Proof of Theorem 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 A.8 Proof of Theorem 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 B Baseline Models and Implementation

56

C Evaluation Metrics

58

D Software Demonstration

59

D.1 Architectural Configurations and Sufficient Dimension Reduction . . . . . . . . . 60 D.2 Output Heads . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 D.3 Distance Metrics in GCES Loss . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 D.4 Noise Injection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 D.5 Usage Example: Training and Inference . . . . . . . . . . . . . . . . . . . . . . . 61 E Literature Review on Estimation of Circular Conditional Functionals

63

E.1 Conditional Mean Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 E.2 Conditional Mode and Concentration Functions . . . . . . . . . . . . . . . . . . . 64 E.3 CDF, Quantiles, and Density Functions . . . . . . . . . . . . . . . . . . . . . . . . 64 F Tables and Figures

65

39

A

Mathematical Proofs

A.1

Auxiliary Lemmas

Lemma A.1.1. Let h : [0, 1] → R be a strictly increasing continuous function and ε ∼ U(0, 1). Then the α-quantile of h(ε) is h(α). Proof. Note that as h is strictly increasing and continuous, it is injective. For any ε ∼ U(0, 1), Fh(ε) (t) = P(h(ε) ≤ t) = P(ε ≤ h−1 (t)) = h−1 (t), where FX denotes the CDF of the random variable X. So the quantile function of h(ε) is given as n

o

n

o

−1 (α) = inf t | Fh(ε) (t) ≥ α = inf t|h−1 (t) ≥ α = h(α), Fh(ε)

since h is strictly increasing. Lemma A.1.2 (Strict propriety of the chordal circular energy score). Let P denote the class of all Borel probability measures on [0, 2π). Let u : [0, 2π) → R2 be the embedding u(θ) = (cos θ, sin θ)⊤ , and let dc (θ1 , θ2 ) = ∥u(θ1 ) − u(θ2 )∥ be the chordal distance. For P ∈ P and θ ∈ [0, 2π), define the chordal circular energy score loss by taking F to be the identity function and d to be dc in (2.2). ◦

Then ES is a strictly proper scoring rule relative to P, i.e., for any Q ∈ P,  ◦

 ◦





EΘ0 ∼Q ES(Q, Θ0 ) ≤ EΘ0 ∼Q ES(P, Θ0 ) , with equality if and only if P = Q. Proof. Let P, Q ∈ P be arbitrary. The map u : [0, 2π) → R2 is continuous and strictly injective on [0, 2π), hence Borel measurable. Define the pushforward measures u∗ P and u∗ Q on S1 . Since u is an injective Borel measurable map between standard Borel spaces, the pushforward operator is injective on Borel probability measures, so that u∗ P = u∗ Q if and only if P = Q. Since the support of u∗ P is contained in S1 , every X ∼ u∗ P satisfies ∥X∥ = ∥u(Θ)∥ = 1 almost surely. So, ∥X − Y ∥ ≤ 2 almost surely for any Y ∼ u∗ Q, and hence E(∥X − Y ∥) < ∞. Substituting dc (θ1 , θ2 ) = ∥u(θ1 ) − u(θ2 )∥ into (2.2) under Q and applying the change of variables X = u(Θ) ∼ u∗ P, X ′ = u(Θ′ ) ∼ u∗ P, and Y = u(Θ0 ) ∼ u∗ Q, we obtain ◦ 1 EΘ0 ∼Q ES(P, Θ0 ) = EY ∼u∗ Q EX∼u∗ P [∥X − Y ∥] − EX,X ′ ∼u∗ P [∥X − X ′ ∥] 2 = EY ∼u∗ Q [ES(u∗ P, Y )] ,





where ES(·, ·) denotes the standard multivariate energy score loss on R2 defined in Shen & Meinshausen (2025). Now due to Józiak (2015, Theorem 4.2), we know that the Euclidean norm is conditionally strictly negative definite on R2 , and hence the same holds for S1 . Then, using Waghmare & Ziegel (2025, Theorem 4), the associated energy score loss must be strictly proper.

40

Therefore we can write,  ◦



EY ∼u∗ Q [ES(u∗ Q, Y )] ≤ EY ∼u∗ Q [ES(u∗ P, Y )] = EΘ0 ∼Q ES(P, Θ0 ) , with equality if and only if u∗ P = u∗ Q. Since u∗ P = u∗ Q if and only if P = Q, equality holds if and only if P = Q. This completes the proof. Lemma A.1.3. Let X be a compactn subset of Rd , and let h : [0, 1] → o R be a continuous function. d ⊤ For any vector β ∈ R , the set S = β x + h(ε) | x ∈ X , ε ∈ [0, 1] is a compact subset of R. Proof. Let Z = X × [0, 1] denote the Cartesian product of the domain spaces. Because X is compact by assumption and the closed interval [0, 1] is compact, their product Z is a compact subset of Rd+1 by Tychonoff’s theorem. Now, define the mapping g : Z → R as g(x, ε) = β ⊤ x+h(ε). The function g can be decomposed into the sum of two functions: the linear projection (x, ε) 7→ β ⊤ x, which is universally continuous, and the scalar function (x, ε) 7→ h(ε), which is continuous by the hypothesis on h. As the sum of two continuous functions, g is continuous on Z. By definition, the set S is precisely the image of Z under g, meaning S = g(Z). Since the continuous image of a compact topological space is compact, it follows that S is a compact subset of R.

A.2

Identifiability Analysis

As identifiability of the model classes serve as a basis for extrapolability (Shen & Meinshausen 2025), we first investigate the identifiability of (2.1). For this, we impose some assumptions on the covariate X, coefficient vectors β ∈ B, and the function classes G and H. The assumptions are stated in Sec. 3.3, where we take, G = {g|g satisfies Assumption 2} ,

H = {h|h satisfies Assumption 3} .

The assumptions establish the necessary regularity conditions to guarantee model identifiability and extrapolability. Assumption 1 imposes a compact, well-behaved support for the covariates with a strictly positive density. This is a standard regularity condition in nonparametric estimation that naturally aligns with the bounded nature of real-world measurements, ensuring stable estimation across the domain. Assumption 2 ensures that the data-generating function g is smooth and, crucially for pre-ANMs, strictly monotone. Following Shen & Meinshausen (2025), requiring monotonicity provides a structural middle ground. It avoids the highly restrictive assumption of strict linearity (which trivially guarantees identifiability) while imposing sufficient shape constraints to track the generative mechanism outside the training support. Assumption 3 characterizes the noise distribution, mapping a uniform base measure via the strictly increasing function h. The condition h(1/2) = 0 anchors the noise with a zero median. More importantly, for directional data, the strict bound |h(1) − h(0)| < 2π ensures the noise support cannot span the S1 manifold entirely. Without this constraint, the noise would fully wrap around the circle, obscuring the structural signal and rendering the model fundamentally unidentifiable. Finally, Assumption 4 specifies the parameter space for the index vector β. Fixing the Euclidean norm (∥β∥ = 1) and the sign of the first component (β(1) > 0) resolves the inherent scale and sign ambiguities standard to single-index 41

models. In practice, this structural constraint can be maintained during optimization via manifold parameterization or projected gradient techniques. (2)

Through the subsequent results, we show that the classes Mpost and M(2) pre are fully identifiable, i.e., two models in these classes that induce the same conditional distribution PY |x for all x ∈ X must necessarily share the same index β, the same regression function g on the training support (1) almost surely, and the same noise function h almost surely on [0, 1]. In contrast, the classes Mpost and M(1) pre exhibit what we term as the 2kπ-identifiability trap. Since the response is observed only modulo 2π, shifting the regression function g by any integer multiple of 2π leaves the conditional distribution of Y |X unchanged. Consequently, parameter configurations differing only by such shifts are observationally indistinguishable, and g is identifiable only up to an additive constant of the form 2kπ, where k ∈ Z. This lack of full identifiability arises because g is allowed to take values in all of R rather than being constrained to [0, 2π). The restriction g : R → [0, 2π) imposed (2) in Mpost and M(2) pre restores full identifiability. Intuitively, once g is constrained to take values within a single period of the circle, the usual ambiguity caused by adding multiples of 2π disappears. Thus, the range restriction selects a unique representative among all wrapped equivalent functions, thereby restoring identifiability. Theorem 6 (Identifiability). (2)

(A) (Identifiability of Mpost ). Let β1 , β2 ∈ B, g1 , g2 ∈ G (g1 , g2 : R → [0, 2π)), and h1 , h2 ∈ H such that 





d



g1 β1⊤ x + h1 (ε) mod 2π =









g2 β2⊤ x + h2 (ε) mod 2π

(A.1)

for all x ∈ X , ε ∈ [0, 1]. Then β1 = β2 =: β, g1 (β ⊤ x) = g2 (β ⊤ x) almost surely in X , and h1 (ε) = h2 (ε) almost surely in [0, 1]. (1)

Let β1 , β2 ∈ B, g1 , g2 ∈ G (g1 , g2 : R → R), and

(B) (2kπ-identifiability of Mpost ). h1 , h2 ∈ H such that 





d



g1 β1⊤ x + h1 (ε) mod 2π =









g2 β2⊤ x + h2 (ε) mod 2π

(A.2)

for all x ∈ X , ε ∈ [0, 1]. Then g1 (β1⊤ x) = g2 (β2⊤ x) + 2kπ for some k ∈ Z almost surely in X , and h1 (ε) = h2 (ε) almost surely in [0, 1]. Hence g is not uniquely identified, and the model class falls into a 2kπ-identifiability trap. (C) (2kπ-identifiability of M(1) pre ). Assume g1 , g2 ∈ G (g1 , g2 : R → R), h1 , h2 ∈ H, and β1 , β2 ∈ B, with the map ε 7→ g(β ⊤ x + h(ε)) supported on an interval of length strictly less than 2π. Suppose 

d







g1 β1⊤ x + h1 (ε) mod 2π = g2 β2⊤ x + h2 (ε) mod 2π

(A.3)

for all x ∈ X , ε ∈ [0, 1]. Then β1 = β2 , h1 = h2 , and g1 (β1⊤ x + h1 (ε)) = g2 (β2⊤ x + h2 (ε)) + 2kπ for some k ∈ Z and all ε ∈ [0, 1], so that M(1) pre falls into a 2kπ-identifiability trap. 42





(D) (Identifiability of M(2) pre ). Assume g1 , g2 ∈ G g1 , g2 : R → [0, 2π) , h1 , h2 ∈ H, and β1 , β2 ∈ B. Suppose 

g1 β1⊤ x + h1 (ε)



d



= g2 β2⊤ x + h2 (ε)



(A.4)

for all x ∈ X , ε ∈ [0, 1]. Then β1 = β2 = β, g1 (β ⊤ x + h1 (ε)) = g2 (β ⊤ x + h2 (ε)), and h1 (ε) = h2 (ε) for all ε ∈ [0, 1]. Proof. The proof is as follows. Proof of (A). For a fixed x ∈ X , (A.1) implies the equality of all trigonometric moments on the unit circle (Bell & Nadarajah 2024), which means for all r ∈ Z: h

n





h

n



oi

E exp ιr g1 (β1⊤ x) + h1 (ε) (mod 2π) =⇒ E exp ιr g1 (β1⊤ x) + h1 (ε) n

o

h

oi

h

n





= E exp ιr g2 (β2⊤ x) + h2 (ε) (mod 2π) n



= E exp ιr g2 (β2⊤ x) + h2 (ε)

n

oi

oi

(A.5)

o

=⇒ exp ιrg1 (β1⊤ x) ψh1 (r) = exp ιrg2 (β2⊤ x) ψh2 (r). Due to Assumption 3, we must have ψh (1) ̸= 0. Taking r = 1 in (A.5) allows us to write: o

n 

exp ι g1 (β1⊤ x) − g2 (β2⊤ x) 

Let δ0 = arg

ψh2 (1) ψh1 (1)



=

ψh2 (1) . ψh1 (1)

(A.6)

∈ (−π, π] be the principal argument of the constant complex fraction. Since

the right-hand side of (A.6) is independent of x, evaluating the complex exponential implies that the arguments must be equal up to an integer multiple of 2π: g1 (β1⊤ x) − g2 (β2⊤ x) = δ0 + 2πκ(x),

(A.7)

where κ(x) ∈ Z is some integer-valued function of x. Substituting this relationship back into (A.5) for an arbitrary integer r ∈ Z yields: exp {ιr(δ0 + 2πκ(x))}ψh1 (r) = ψh2 (r). Because both r and κ(x) are integers, the term exp(ιr2πκ(x)) is identically 1 for all x ∈ X . Hence, exp (ιrδ0 )ψh1 (r) = ψh2 (r),

∀r ∈ Z.

(A.8)

By the uniqueness of trigonometric moments, (A.8) implies that the wrapped distributions are globally identical: d (h2 (ε)) mod 2π = (h1 (ε) + δ0 ) mod 2π. By Assumption 3(A), h is strictly increasing, so its image is the interval [h(0), h(1)]. By Assumption 3(B), |h(1) − h(0)| < 2π, which ensures this image interval has a length strictly less than 2π. Consequently, the mapping to the circle is injective and contains no overlapping probability mass.

43

Therefore, the unwrapped distributions must be equal up to a single global integer shift k ∈ Z: d

h2 (ε) + 2kπ = h1 (ε) + δ0 .

(A.9)

Because ε ∼ U(0, 1) and both h1 , h2 are strictly monotonically increasing, the median of the distribution occurs exactly at ε = 0.5. From Assumption 3(A), h1 (0.5) = h2 (0.5) = 0. Taking the median of both sides of (A.9) yields 2kπ = δ0 . Because δ0 ∈ (−π, π] is a principal argument, the only integer k that satisfies this equality is k = 0, which consequently forces δ0 = 0. Substituting this into (A.9) reveals that h1 (ε) and h2 (ε) are identical in distribution. Because h1 and h2 are strictly monotonically increasing functions mapping from the same uniform random variable ε, they uniquely define the quantile functions of their distributions (see Lemma A.1.1). Identical distributions possess identical quantile functions, forcing pointwise equality: h1 (ε) = h2 (ε)

almost surely in [0, 1].

Finally, we substitute δ0 = 0 back into the spatial equation (A.7) to obtain: g1 (β1⊤ x) − g2 (β2⊤ x) = 2πκ(x). (2)

By the definition of the model class Mpost , both g1 and g2 strictly map to the codomain [0, 2π). Therefore, their absolute difference is strictly bounded: |g1 (β1⊤ x) − g2 (β2⊤ x)| < 2π. This forces |2πκ(x)| < 2π. Since κ(x) is an integer, the only valid solution is κ(x) = 0 for all x ∈ X . Hence: g1 (β1⊤ x) = g2 (β2⊤ x).

(A.10)

Due to Assumption 2(A), Assumption 2(B), and the normalization constraints in Assumption 4, (A.10) implies β1 = β2 and g1 = g2 almost surely on X by Theorem 1 of Lin & Kulasekera (2007). Proof of (B). For a fixed x ∈ X , (A.2) implies the equality of all trigonometric moments on the unit circle (Bell & Nadarajah 2024). Because the additive noise structure is identical to that of (2) model class Mpost , we can apply the exact same deductions utilized in the proof of Theorem 6(A). By isolating the first trigonometric moments and unwrapping the bounded noise distributions, we immediately recover both the pointwise equality of the noise functions, h1 (ε) = h2 (ε) almost surely on [0, 1], and the global spatial phase relationship: g1 (β1⊤ x) = g2 (β2⊤ x) + 2kπ, (2)

for some k ∈ Z. (1)

Now, unlike model class Mpost , the functions g1 , g2 ∈ Mpost map to the unbounded codomain R. Because their outputs are not restricted to [0, 2π), their absolute difference |g1 (β1⊤ x) − g2 (β2⊤ x)| is not geometrically constrained to be strictly less than 2π. Hence, k = 0 is not forced, and the functions are only identifiable up to an integer multiple of 2π, establishing that the model class falls into a 2kπ-identifiability trap.

44

Proof of (C). Define a push-forward probability measure on the target space [0, 2π) as: µj (A) = U

n





a ∈ [0, 1] | gj βj⊤ x + hj (a) mod 2π ∈ A

o

,

for any Borel set A ∈ B([0, 2π)).

The equality in distribution (A.3) implies that the two random variables induce the identical probability measure on the target space [0, 2π). Therefore, for every measurable set A ∈ B([0, 2π)), h





i

h





i

P g1 β1⊤ x + h1 (ε) mod 2π ∈ A = P g2 β2⊤ x + h2 (ε) mod 2π ∈ A .

(A.11)

Because ε ∼ U(0, 1), we have P(ε ∈ B) = U(B) for any measurable B ⊆ [0, 1]. Thus, 







P gj βj⊤ x + hj (ε) mod 2π ∈ A = µj (A). Substituting this into (A.11) yields µ1 (A) = µ2 (A) for all Borel sets A, meaning µ1 ≡ µ2 as measures on [0, 2π). Next, define an unwrapped push-forward measure on R as: νj (B) = U

n





a ∈ [0, 1] | gj βj⊤ x + hj (a) ∈ B

o

,

for any Borel set B ∈ B(R).

(A.12)

Let ω : R → [0, 2π) denote the Borel measurable wrapping map ω(t) = t mod 2π. By assumption, the mapping ε 7→ gj (βj⊤ x + hj (ε)) is supported on an interval strictly shorter than 2π. Let Iνj ⊂ R denote this interval. The measure νj is therefore entirely supported on Iνj . Because the length of Iνj is strictly less than 2π, the restriction of the wrapping map to this interval, ω|Iνj , is strictly injective. For any measurable A ∈ B([0, 2π)), µj (A) = U

n





o

a ∈ [0, 1] | gj βj⊤ x + hj (a) ∈ ω −1 (A)

= νj (ω −1 (A)) = ω∗ νj (A), which establishes that µj is the push-forward of νj under ω. Conversely, since ω is injective on Iνj , we can write for any measurable B ∈ B(R): νj (B) = νj (B ∩ Iνj ) (since νj is supported on Iνj ) 





= νj ω −1 ω(B ∩ Iνj ) 

(by injectivity of ω on Iνj )



= µj ω(B ∩ Iνj ) (since µj = ω∗ νj ).

(A.13)

Because µ1 ≡ µ2 , and both ν1 and ν2 are supported on intervals strictly shorter than 2π, their supports can only differ by an exact integer multiple of 2π. That is, there exists a unique k ∈ Z such that Iν1 = Iν2 + 2πk. Consequently, the wrapped images ω(Iν1 ) and ω(Iν2 ) perfectly coincide on [0, 2π). Due to the 2π-periodicity of ω, for any measurable B: ω(B ∩ Iν1 ) = ω (B ∩ (Iν2 + 2πk)) = ω ((B − 2πk) ∩ Iν2 ) .

45

Therefore, utilizing (A.13), we obtain: ν1 (B) = µ1 (ω(B ∩ Iν1 )) = µ2 (ω(B ∩ Iν1 )) = µ2 (ω((B − 2πk) ∩ Iν2 )) = ν2 (B − 2πk).

(A.14)

Equation (A.14) establishes that ν1 and ν2 are identical measures up to a spatial shift of 2πk. Returning to the definition of νj in (A.12), this implies that for some k ∈ Z: 

d







g1 β1⊤ x + h1 (ε) = g2 β2⊤ x + h2 (ε) + 2πk. Fixing x ∈ X , the composition ε 7→ g(β ⊤ x + h(ε)) is strictly monotonically increasing in ε due to Assumption 2(C) and Assumption 3(A). Because strictly increasing functions of a uniform variable uniquely define the quantile functions of their distributions, equality in distribution necessitates pointwise equality almost surely: g1 (β1⊤ x + h1 (ε)) = g2 (β2⊤ x + h2 (ε)) + 2πk,

∀ε ∈ [0, 1].

(A.15)

By taking the gradient of both sides of this equation with respect to x, we obtain the vector equation: 











g1′ β1⊤ x + h1 (ε) β1 = g2′ β2⊤ x + h2 (ε) β2 .

(A.16)

Taking the Euclidean ℓ2 -norm of both sides yields: 



g1′ β1⊤ x + h1 (ε) ∥β1 ∥2 = g2′ β2⊤ x + h2 (ε) ∥β2 ∥2 . By Assumption 4, the coefficient vectors are normalized such that ∥β1 ∥2 = ∥β2 ∥2 = 1, reducing the equality to:     g1′ β1⊤ x + h1 (ε) = g2′ β2⊤ x + h2 (ε) . By Assumption 2(C), g1 and g2 are strictly monotone, meaning non-zero  their derivatives  are strictly   ′ ⊤ ′ ⊤ and maintain a constant sign. Consequently, we must have g1 β1 x + h1 (ε) = ±g2 β2 x + h2 (ε) . Substituting this into (A.16), 







±g2′ β2⊤ x + h2 (ε) β1 = g2′ β2⊤ x + h2 (ε) β2 . Because the derivative g2′ (·) is strictly non-zero, we can divide both sides by this scalar term, isolating the coefficient vectors to obtain ±β1 = β2 . Assumption 4 further dictates that the first coordinate of any β ∈ B is strictly positive (β(1) > 0). This strictly rules out the negative root β1 = −β2 , forcing the conclusion that β1 = β2 = β. Having established β1 = β2 = β, we return to the intermediate equality (A.15). By Assumption 3(A), the noise functions are anchored at the median such that h1 (0.5) = h2 (0.5) = 0.

46

Evaluating (A.15) at ε = 0.5 yields: 







g1 β ⊤ x = g2 β ⊤ x + 2kπ. Because this relationship holds for all x ∈ X , it describes a global constant shift between the two functions. We can substitute this functional relationship g2 (t) + 2kπ = g1 (t) back into (A.15) for an arbitrary ε: 







g1 β ⊤ x + h1 (ε) = g1 β ⊤ x + h2 (ε) . Because g1 is strictly monotone (Assumption 2(C)), it is strictly injective. The equality of the evaluated outputs therefore guarantees the equality of their inputs: β ⊤ x + h1 (ε) = β ⊤ x + h2 (ε) =⇒ h1 (ε) = h2 (ε). This confirms that h1 = h2 for all ε ∈ [0, 1]. We have thus established β1 = β2 , h1 = h2 , and g1 (·) = g2 (·) + 2kπ, proving that M(1) pre falls into the 2kπ-identifiability trap. Proof of (D). Fix x ∈ X . Due to Assumption 2(C) and Assumption 3(A), the composition ε 7→ g(β ⊤ x + h(ε)) is strictly monotone in ε. By Lemma A.1.1, the a-th quantile of h(ε) is h(a), so that 







gj βj⊤ x + Qa (hj (ε)) = gj βj⊤ x + hj (a) ,

(A.17)

where Qa (·) denotes the a-level quantile function. Now (A.4) implies that all the quantiles of the random variables on both sides should be equal, and hence using (A.17), we must have 







g1 β1⊤ x + h1 (a) = g2 β2⊤ x + h2 (a)

(A.18)

for all x ∈ X , a ∈ [0, 1]. By Assumption 3(A), h1 (0.5) = h2 (0.5) = 0. Evaluating (A.17) at a = 0.5 yields: 







g1 β1⊤ x = g2 β2⊤ x . Because this holds for all x ∈ X , under Assumption 1, Assumption 2(B), and the coefficient constraints in Assumption 4, we can invoke Theorem 1 of Lin & Kulasekera (2007) to conclude that β1 = β2 = β

and

g1 = g2 .

Having established the equality of the index parameters and the spatial functions, we substitute β1 = β2 = β and g2 = g1 back into the pointwise equality (A.18) for an arbitrary ε ∈ [0, 1] to get 







g1 β ⊤ x + h1 (ε) = g1 β ⊤ x + h2 (ε) . Because g1 is strictly monotone (Assumption 2(C)), it is injective. The equality of the evaluated

47

outputs therefore guarantees the equality of their inputs: β ⊤ x + h1 (ε) = β ⊤ x + h2 (ε) =⇒ h1 (ε) = h2 (ε). This confirms that h1 (ε) = h2 (ε) for all ε ∈ [0, 1]. We have thus established β1 = β2 , g1 = g2 , and h1 = h2 , completing the identifiability proof for M(2) pre .

A.3

Proof of Theorem 1

Proof. Given any model g ∈ M with g ∼ Ptr (y|x), assume for the sake of contradiction that there exists a subset X1 ⊆ X with P(X ∈ X1 ) > 0, such that Peg (y|x) ̸= Ptr (y|x) for all x ∈ X1 . Now, since (2.2) is a strictly proper scoring rule, we must have for all x ∈ X1 ,  ◦

 ◦ 



EY ∼Ptr (y|x) ES (Ptr (y|x), Y ) < EY ∼Ptr (y|x) ES Peg (y|x), Y



,

∀x ∈ X1 .

(A.19)

,

∀x ∈ X1 .

(A.20)

Since we have g ∼ Ptr (y|x), (A.19) can be written as  ◦

 ◦ 



EY ∼Ptr (y|x) ES (Pg (y|x), Y ) < EY ∼Ptr (y|x) ES Peg (y|x), Y



Taking expectation on both sides of (A.20) with respect to Ptr (since X1 has non-zero measure) gives  ◦

 ◦ 



EPtr ES (Pg (y|X), Y ) < EPtr ES Peg (y|X), Y



,

which is a contradiction since ge is the minimizer in (2.3). So, we must have P(X ∈ X1 ) = 0, which means that ge ∼ Ptr (y|x). This completes the proof.

A.4

Proof of Lemma 1

Proof. By the population level guarantee (Theorem 1) and our correct specification assumption, the minimizers satisfy ge(x, ε) ∼ Ptr (·|x) and geα (x, ε) ∼ Pαtr (·|x) almost everywhere. By definition, Pαtr (·|x) is the distribution of Y α = Y ⊕ α, where Y ∼ Ptr (·|x). This means that Pαtr is the push-forward measure of Ptr under the map Tα : θ 7→ θ ⊕ α, i.e., Pαtr = (Tα )∗ Ptr . Using the change of variables for push-forward measures, we can write EΘ∼(Tα )∗ Ptr [ϕ(Θ)] = EΘ∼Ptr [ϕ(Tα (Θ))] = EΘ∼Ptr [ϕ(Θ ⊕ α)],

(A.21)

for any Borel measurable function ϕ. Note that the map Tα is piecewise continuous on [0, 2π), making it measurable with respect to the Borel σ-algebra on [0, 2π). Thus, the push-forward measure is well-defined for any Borel measurable ϕ.

48

Proof of (A). Since geα ∼ Pαtr , we can write Eε [sin geα (x, ε)] = EΘ∼Pαtr (·|x) [sin Θ],

Eε [cos geα (x, ε)] = EΘ∼Pαtr (·|x) [cos Θ].

(A.22)

Applying (A.21) with ϕ(θ) = sin θ and ϕ(θ) = cos θ to (A.22), and noting that the 2π-periodicity of the trigonometric functions allows us to drop the modulo operator (ϕ(Θ ⊕ α) = ϕ(Θ + α)), we obtain EΘ∼Pαtr (·|x) [sin Θ] = EΘ∼Ptr (·|x) [sin(Θ + α)],

EΘ∼Pαtr (·|x) [cos Θ] = EΘ∼Ptr (·|x) [cos(Θ + α)]. (A.23)

Using the sum identities sin(θ + α) = sin θ cos α + cos θ sin α and cos(θ + α) = cos θ cos α − sin θ sin α, and denoting s̄ = EΘ∼Ptr (·|x) [sin Θ] and c̄ = EΘ∼Ptr (·|x) [cos Θ], (A.23) expands to EΘ∼Pαtr (·|x) [cos Θ] = c̄ cos α − s̄ sin α = c̄α . (A.24)

EΘ∼Pαtr (·|x) [sin Θ] = s̄ cos α + c̄ sin α = s̄α ,

Using (A.24), we can write c̄α + ιs̄α = (c̄ cos α − s̄ sin α) + ι(s̄ cos α + c̄ sin α) = (c̄ + ιs̄)(cos α + ι sin α) = (c̄ + ιs̄) exp(ια), and hence, e µe α (x) = atan2(s̄α , c̄α ) = atan2(s̄, c̄) ⊕ α = µ(x) ⊕ α,

which completes the proof. Proof of (B). For the sake of simplicity, assume that the circular median is unique. However, the proof holds even if the set of minimizers contains multiple values. Note that we can write for any θ1 , θ2 ∈ [0, 2π), d(θ1 , θ2 ) = D(θ1 ⊖ θ2 ),

(A.25)  

where D(u) = π − |π − u| for geodesic metric and D(u) = 2 sin u2 for chordal distance metric (both metrics are rotationally invariant). Applying (A.21) with ϕ(θ) = d(θ, ψ) for any fixed ψ ∈ [0, 2π), we have Eε [d(geα (x, ε), ψ)] = EΘ∼Pαtr (·|x) [d(Θ, ψ)] = EΘ∼Ptr (·|x) [d(Θ ⊕ α, ψ)] = EΘ∼Ptr (·|x) [d(Θ, ψ ⊖ α)] (using d(θ1 , θ2 ) = D(θ1 ⊖ θ2 )) = Eε [d(ge(x, ε), ψ ⊖ α)] (since ge(x, ε) ∼ Ptr (·|x)).

49

Taking arg min over [0, 2π) on both sides, we get arg min {Eε [d(geα (x, ε), ψ)]} = arg min {Eε [d(ge(x, ε), ψ ⊖ α)]} ψ∈[0,2π)

ψ∈[0,2π)

= arg min {Eε [d(ge(x, ε), ψ)]} ⊕ α, ψ∈[0,2π)

or equivalently, megα (x) = meg (x) ⊕ α PX -almost everywhere. This completes the proof. Proof of (C). We have σegα (x) = Medε [d(geα (x, ε), megα (x))] h

i

= MedΘ∼Pαtr (·|x) d(Θ, megα (x)) (since geα ∼ Pαtr (·|x)) h

i

= MedΘ∼Ptr (·|x) d(Θ ⊕ α, megα (x)) (using (A.21) with ϕ(θ) = d(θ, megα (x))) h

i

= MedΘ∼Ptr (·|x) d(Θ ⊕ α, meg (x) ⊕ α) (using part (B))

(A.26)

Due to (A.25), we must have d(Θ⊕α, meg (x)⊕α) = D((Θ⊕α)⊖(meg (x)⊕α)). Using the properties of modulo 2π arithmetic over [0, 2π), the α terms cancel, yielding (Θ ⊕ α) ⊖ (meg (x) ⊕ α) = Θ ⊖ meg (x). Therefore, d(Θ ⊕ α, meg (x) ⊕ α) = D(Θ ⊖ meg (x)) = d(Θ, meg (x)). Substituting this back into our median expression, we obtain: h

i

σegα (x) = MedΘ∼Ptr (·|x) d(Θ, meg (x)) h

i

= Medε d(ge(x, ε), meg (x)) = σeg (x), which holds PX -almost everywhere. This completes the proof. Proof of (D). By definition of the expected trigonometric moments, we can write: ψeg (r|x) = Eε [exp(ιrge(x, ε))] = EΘ∼Ptr (·|x) [exp(ιrΘ)] , ψegα (r|x) = Eε [exp(ιrgeα (x, ε))] = EΘ∼Pαtr (·|x) [exp(ιrΘ)] . Applying the push-forward change of variables (A.21) with ϕ(θ) = exp(ιrθ), we have: ψegα (r|x) = EΘ∼Pαtr (·|x) [exp(ιrΘ)] = EΘ∼Ptr (·|x) [exp(ιr(Θ ⊕ α))]. Because the complex exponential function exp(ιr·) is 2π-periodic for any integer r ∈ Z, the modulo addition is equivalent to standard addition in the exponent. Thus, exp(ιr(Θ ⊕ α)) = exp(ιr(Θ + α)) = exp(ιrα) exp(ιrΘ). Factoring out the constant, we have ψegα (r|x) = exp(ιrα)EΘ∼Ptr (·|x) [exp(ιrΘ)] = exp(ιrα)ψeg (r|x).

50

Substituting this into the truncated density equation, we obtain: Pegα (θ|x) = (2π)−1 −1

= (2π)

R X r=−R R X r=−R

= (2π)−1

R X r=−R

ψegα (r|x) exp(−ιrθ) ψeg (r|x) exp(ιrα) exp(−ιrθ) ψeg (r|x) exp(−ιr(θ − α)).

Since the exponential function evaluates evaluating the phase modulo 2π, we can replace (θ − α) with (θ ⊖ α), resulting in Peg (θ ⊖ α|x) in the RHS. This completes the proof. Proof of (E). For the sake of simplicity, assume the conditional mode is unique. Using part (D), we have: n

o

n

o

ῡegα (x) = arg max Pegα (θ|x) = arg max Peg (θ ⊖ α|x) . θ∈[0,2π)

θ∈[0,2π)

Let θ′ = θ ⊖ α. Because the mapping θ 7→ θ ⊖ α is a bijection on [0, 2π), evaluating the maximum over θ is equivalent to evaluating the maximum over θ′ and applying the inverse transformation θ = θ′ ⊕ α. Therefore: n

o

n

′

! o

arg max Peg (θ ⊖ α|x) = arg max Peg (θ |x)

⊕α

θ′ ∈[0,2π)

θ∈[0,2π)

= ῡeg (x) ⊕ α. This establishes that ῡegα (x) = ῡeg (x) ⊕ α. The proof for the conditional antimode, υegα (x) = υeg (x) ⊕ α, follows identically by replacing arg max with arg min. This completes the proof.

A.5

Proof of Theorem 2 ◦

◦

Proof. Let Πy = PX|Y (·|y) and P = {Πy |y ∈ S1 }. Since P is dominated by a σ-finite measure, ◦

◦

◦

it contains a countable subset Q = {Qk |k = 1, 2, . . .} such that P ≡ Q (we use ≡ to denote that one measure dominates the other). Let {ck }k=1,2,... be a sequence of positive numbers adding up ◦

◦

◦

◦

◦

to 1 and let Q0 = ∞ k=1 ck Qk . Then Q0 is a probability measure on X such that {Q0 } ≡ Q ≡ P. dΠy Let πy = ◦ and let G be a sub-sigma field of σ(X). We claim that Y ⊥ ⊥ X|G if and only if πy is P

dQ0

◦

essentially measurable with respect to G for all y ∈ S1 modulo Q0 . To show the necessary part, let

51

B ∈ σ(X). Then, h

E◦

Q0

i

πy (X)1{X∈B} = EΠy [1{X∈B} ] h



i

= EΠy EΠy 1{X∈B} |G h

= E◦



i

πy (X)EΠy 1{X∈B} |G

Q0

.

(A.27)

Now, since Y ⊥ ⊥ X|G, Πy (B|G) is same for all y ∈ S1 . Hence Πy (B|G) = Qk (B|G) for all k, ◦

implying that Πy (B|G) = Q0 (B|G). So, we can write the RHS of (A.27) as 

E◦

Q0

πy (X)E ◦



Q0





1{X∈B} |G

= E◦

Q0



1{X∈B} E ◦ (πy (X)|G) . Q0

Using this in (A.27), we can write E◦

h

Q0



i

πy (X)1{X∈B} = E ◦



1{X∈B} E ◦ (πy (X)|G) ,

Q0

(A.28)

Q0

◦

which implies πy (X) = E ◦ [πy (X)|G]-a.s. Q0 , since (A.28) holds for all B ∈ σ(X). Next, we Q0

prove the sufficiency part. Assume that πy is essentially measurable with respect to G for all y ∈ S1 ◦

modulo Q0 , which is equivalent to E ◦ [πy (X)|G] = πy (X).

(A.29)

Q0

Then for any A ∈ G, we can write 

Eπy 1{X∈A} E ◦

Q0





1{X∈B} |G



= E◦

1{X∈A} πy (X)E ◦

Q0



Q0



1{X∈B} |G





= E ◦ 1{X∈A} 1{X∈B} E ◦ (πy (X)|G) Q0

Q0

h

i

= E ◦ 1{X∈A} 1{X∈B} πy (X) (due to (A.29)) Q0

= Eπy [1{X∈A} 1{X∈B} ] = Πy (X ∈ A ∩ B). Therefore, we can conclude that E ◦

Q0

h

i

◦

1{X∈B} |G = Q0 (B|G) is the conditional probability

Πy (B|G), which means that Πy (B|G) does not depend on y, hence Y ⊥ ⊥ X|G. Now let G∗ be the intersection of all SDR sigma fields G. Then G∗ is itself a sigma field. Moreover, since πy is measurable with respect to all SDR sigma fields for all y ∈ S1 , it should also be measurable with respect to G∗ for all y ∈ S1 . Hence, G∗ is itself an SDR sigma field, which implies that it is also the smallest SDR sigma field. If G0 is another smallest SDR sigma field then we know that G∗ ⊆ G0 and G0 ⊆ G∗ . Therefore, G∗ is unique.

52

A.6

Proof of Theorem 3

Proof. Since ε is independent of X, under correct specification, PY |X depends on X only through β0⊤ X. That is, for any Borel measurable set B ⊂ [0, 2π), P(Y ∈ B|X = x) = P(Y ∈ B|β0⊤ X = β0⊤ x), which means that we should have Y ⊥ ⊥ X|σ(β0⊤ X). This means that σ(β0⊤ X) is always an SDR sigma field and therefore, we must have G∗Y |X ⊆ σ(β0⊤ X). Next, we need to show that σ(β0⊤ X) ⊆ G∗Y |X . Suppose for the sake of contradiction that σ(β0⊤ X) ⊈ G∗Y |X . We can write G∗Y |X = σ(R(X)) for some measurable function R such that σ(β0⊤ X) ⊈ σ(R(X)). Now since σ(R(X)) is an SDR sigma field, the conditional distribution PY |X depends on X only through R(X), i.e., PY |X=x = PY |R(X)=R(x) ,

PX -a.s..

Therefore, we can write R(x1 ) = R(x2 ) =⇒ PY |R(X)=R(x1 ) = PY |R(X)=R(x2 ) =⇒ PY |X=x1 = PY |X=x2 =⇒ β0⊤ x1 = β0⊤ x2 , where the last step holds because the map t 7→ PY |β0⊤ X=t is injective (note that injectivity of f0 guarantees the injectivity of the induced measure PY |· ). The above chain means that β0⊤ X is a function of R(X) almost surely. Hence, using Doob-Dynkin’s lemma, β0⊤ X is measurable with respect to σ(R(X)), i.e., σ(β0⊤ X) ⊆ σ(R(X)), which is a contradiction. Thus, we must have σ(β0⊤ X) = G∗Y |X , which completes the proof.

A.7

Proof of Theorem 4

Proof. By Theorem 1, under correct specification, the population minimizer satisfies   ⊤ e e e g β x, h(ε) ∼ PY |x for all x ∈ X . Because the conditional distribution of Y given X = x depends on x strictly through the linear combination βe⊤ x, it follows that Y ⊥ ⊥ X | βe⊤ X, meaning that σ(βe⊤ X) is an SDR sigma field. From Theorem 3, the unique central subspace satisfies G∗Y |X = σ(β0⊤ X). By the minimality of the central subspace, any other sub-sigma field that 







satisfies the SDR condition must contain it. Therefore, we must have σ β0⊤ X ⊆ σ βe⊤ X . By the Doob-Dynkin factorization lemma, this inclusion  implies the existence of a Borel measurable function ψ : R → R such that β0⊤ X = ψ βe⊤ X almost surely. Suppose for the sake of contradiction that β0 and βe are linearly independent vectors in Rd . Define the linear transformation h i e β ⊤ , which by assumption has full row rank 2. Because X is absolutely continuous matrix B = β, 0 with respect to the Lebesgue measure on Rd , its projection onto this full-rank two-dimensional  ⊤ subspace, given by the random vector Z = βe⊤ X, β0⊤ X , must be absolutely continuous with respect to the two-dimensional Lebesgue measure λ2 on R2 . Consequently, for any Borel set A ⊂ R2 with λ2 (A) = 0, we must have P(Z ∈ A) = 0. Define Ψ = {(u, v) ∈ R2 | v = ψ(u)} as the graph of the function ψ. Because ψ is Borel measurable, Fubini’s Theorem guarantees that its graph has a 2D Lebesgue measure of  zero (λ2 (Ψ) = 0). Utilizing the absolute continuity of Z, this implies P



βe⊤ X, β0⊤ X

⊤

∈ Ψ = 0. However, our previous deduction states that 53





β0⊤ X = ψ βe⊤ X almost surely, meaning the random vector resides on the graph with probability one: P



 ⊤ ⊤ ⊤ e β X, β X ∈ Ψ = 1. This is a contradiction. We conclude that β

0 and β cannot

0

e

be linearly independent. Since neither vector is the zero vector (assuming non-trivial dimension reduction), they must be collinear. Thus, we can write βe = cβ0 for some non-zero scalar c. This implies that σ(βe⊤ X) = σ(β0⊤ X). Therefore, G∗Y |X = σ(βe⊤ X), completing the proof.

A.8

Proof of Theorem 5

Proof. We illustrate the proof for part (B); the proof for part (A) follows identically by dropping the additive noise h(ε) from the argument of g. Fix ε ∈ [0, 1], and let M1 , M2 ∈ M(2) pre be two models parameterized by (g1 , β1 , h1 ) and (g2 , β2 , h2 ) that induce the identical conditional distribution on the training support: 



d





g1 β1⊤ x + h1 (ε) = g2 β2⊤ x + h2 (ε) ,

∀x ∈ X .

(A.30)

By the identifiability guarantee in Theorem 6(D), it follows that β1 = β2 = β, h1 = h2 = h on [0, 1], and g1 (β ⊤ x + h(ε)) = g2 (β ⊤ x + h(ε)) almost surely for x ∈ X . Because ε ∼ U(0, 1) and h is continuous and strictly increasing (Assumption 3), the image h([0, 1]) is a closed interval [a, b] for some a < b. Fix any interior point x0 ∈ X . We can write g1 (β ⊤ x0 + t) = g2 (β ⊤ x0 + t) for all t ∈ [a, b]. This implies that g1 and g2 perfectly agree on the closed interval [β ⊤ x0 +a, β ⊤ x0 +b]. Because a closed interval contains infinitely many accumulation points, and both g1 and g2 are real analytic on the connected open set I, the identity theorem for real analytic functions (Krantz & Parks 2002) dictates that g1 = g2 everywhere on I. Next, we establish the topological boundary for extrapolation. Define the projected parameter space for this specific model as Tpre = {β ⊤ x + h(ε) | x ∈ X , ε ∈ [0, 1]}. Because X is compact, h is continuous, and the projection is linear, Tpre is a compact subset of R (see Lemma A.1.3). By assumption, Tpre ⊂ I. Define the set distance to the boundary of the analytic domain: ∆ = ds (Tpre , R \ I) =

inf

s∈Tpre ,u∈R\I

|s − u| .

Because I is open, its complement R \ I is closed. The distance between a compact set (Tpre ) and a disjoint closed set (R \ I) is strictly positive. Thus, ∆ > 0. Now, let x∗ be an extrapolation point satisfying inf x∈X ∥x − x∗ ∥2 ≤ δ for some δ > 0. Let xc ∈ X be the closest point to x∗ in the support. By the Cauchy-Schwarz inequality and the constraint ∥β∥ = 1 (Assumption 4), the distance in the projected scalar space is bounded by δ: 





β ⊤ x∗ + h(ε) − β ⊤ xc + h(ε)



= β ⊤ (x∗ − xc ) ≤ ∥β∥ ∥x∗ − xc ∥ ≤ δ.

Therefore, the projected point β ⊤ x∗ + h(ε) lies within a δ-neighborhood of Tpre . By choosing δ0 < ∆, any δ ≤ δ0 guarantees that β ⊤ x∗ + h(ε) strictly remains inside the analytic domain I for

54

all ε ∈ [0, 1]. Because we have already established that g1 = g2 everywhere on I, it follows that: g1 (β ⊤ x∗ + h(ε)) = g2 (β ⊤ x∗ + h(ε)),

∀ε ∈ [0, 1].

This establishes that D(P1 (y|x∗ ), P2 (y|x∗ )) = 0. Because this holds for any x∗ within a δneighborhood of X , taking the supremum yields UM(2) (δ) = 0 for all δ ≤ δ0 . This completes the pre proof for the fully identifiable model classes. (1)

(1) Proof for models under 2kπ-identifiability trap (M(1) pre and Mpost ). We detail the proof for Mpre ; (1) the logic applies identically to Mpost by modifying the position of the additive noise. Let M1 , M2 ∈ M(1) pre be two models that induce the identical conditional distribution on the training support. By the 2kπ-identifiability guarantee (Theorem 6(C)), the underlying parameters are not uniquely identified; instead, they satisfy h1 = h2 = h, β1 = β2 = β, and g1 (β ⊤ x + h(ε)) = g2 (β ⊤ x + h(ε)) + 2kπ for some integer k ∈ Z, almost surely for x ∈ X . Using the exact same topological arguments established previously, g1 and g2 must maintain a constant difference of 2kπ over a closed interval in the projected space. Because g1 and g2 are real analytic on the connected open set I, the function f (z) = g1 (z) − g2 (z) − 2kπ is also real analytic. Since f (z) = 0 on a closed interval, the identity theorem (Krantz & Parks 2002) requires that f (z) = 0 everywhere on I. Therefore, g1 (z) = g2 (z) + 2kπ across the entire analytic domain I. Now, consider an extrapolation point x∗ such that β ⊤ x∗ +h(ε) ∈ I. The generated circular responses for the two models at this extrapolation point are given by:

Y1∗ = g1 (β ⊤ x∗ + h(ε)) mod 2π Y2∗ = g2 (β ⊤ x∗ + h(ε)) mod 2π Substituting the analytic continuation g1 (z) = g2 (z) + 2kπ into the first equation yields: 



Y1∗ = g2 (β ⊤ x∗ + h(ε)) + 2kπ mod 2π. By the modulo 2π arithmetic, the additive 2kπ shift vanishes, and we get Y1∗ = g2 (β ⊤ x∗ + h(ε)) mod 2π = Y2∗ . Because Y1∗ = Y2∗ deterministically for all ε ∈ [0, 1], the conditional distributions match: P1 (y|x∗ ) = P2 (y|x∗ ). Consequently, D(P1 (y|x∗ ), P2 (y|x∗ )) = 0. Taking the supremum over the δ-neighborhood confirms UM(1) (δ) = 0, proving that the model is extrapolable despite being not pre fully identifiable.

Remark 1. Identifiability is not a strictly necessary condition for distributional extrapolability. For example, Shen & Meinshausen (2025) note that the class of linear pre-ANMs is not identifiable, but distributionally extrapolable.

55

B

Baseline Models and Implementation

In this section, we briefly describe the baseline models used for comparison along with their implementation strategies. As discussed in the main text, we do not consider Bayesian implementations for all models except vMQP. This is primarily due to two reasons: to avoid the excessive computational overhead of Bayesian training paradigms, and the non-availability of original open-source codes for most of the frameworks. The task-specific usage of each model for comparison in this study is mentioned in Table 1 in the main paper. 1. Circular Linear Regression (CLR): The projected multivariate linear model proposed by Presnell et al. (1998) models a circular response variable by treating it as the radial projection of an unobserved, latent bivariate normal random vector onto the unit circle. The mean of this latent Euclidean vector is formulated as a linear combination of the input covariates, and the corresponding regression weights are efficiently estimated via an Expectation-Maximization (EM) algorithm. The code implementation is taken from a public repository6 . Natively, the framework is designed to map purely linear (Euclidean) covariates to a directional response. However, to keep a fair comparison strategy, we embed the circular predictors (if any) trigonometrically into continuous sine and cosine components prior to incorporating them into the model’s design matrix. 2. Kernel CLR: Nonparametric kernel regression for a circular response (Taylor 2012, AlonsoPena & Crujeiras 2025a) extends standard local smoothing techniques to accommodate directional data. Because angular responses cannot be averaged arithmetically without artificially depending on the zero-direction, the circular response is first decomposed into its orthogonal sine and cosine components. The model then computes a locally weighted average for both components, where the assigned weights are determined by a kernel function that measures the proximity of the evaluation point to the training covariates. The final predicted mean direction is subsequently reconstructed by applying the two-argument arctangent function (atan2) to the smoothed sine and cosine estimates. We implement this nonparametric framework to natively process both linear and circular covariates simultaneously. This is achieved by utilizing a multivariate product kernel. Specifically, the joint localized weight for any given data point is constructed by multiplying individual feature kernels: a standard Euclidean kernel (such as a Gaussian kernel governed by a bandwidth parameter h) is applied to the continuous linear features, while a directional kernel (such as a von Mises kernel governed by a concentration parameter κ) is applied to the circular features. This formulation allows the model to natively map mixed-domain inputs directly to a circular response without requiring artificial embedding of the circular predictors prior to fitting. 3. Skew-Circular Beta Regressor (SCBR): The Skew-Circular Beta (SCB) regression model (Hassanzadeh 2021) provides a parametric framework for predicting multimodal and asymmetric angular responses. The model governs the circular distribution through a global location parameter µ ∈ [0, 2π) and two global shape parameters a, b > 0.5 which control the concentration and multimodality. To capture complex, observation-specific dynamics, the model dynamically predicts a directional skewness parameter λi ∈ [−1, 1] conditioned on the input covariates. This is achieved using a linear predictor ηi = xi⊤ β + bias, which is 6

https://gist.github.com/ahwillia/2941fb64e8bfe999c66291257601dfb4

56

subsequently squeezed into the valid mathematical boundaries via an arctangent link function, λi = π2 arctan(ηi ). The SCB architecture aggregates features strictly through a standard Euclidean linear predictor, and does not natively process circular covariates. We feed the sin-cos embeddings of the circular features to the model for a consistent evaluation strategy. Although Hassanzadeh (2021) presents a Bayesian scheme to estimate the model parameters, we consider two distinct non-Bayesian training paradigms: Maximum Likelihood Estimation (MLE) and Mean Circular Prediction Error (MCPE) minimization, which the author uses as a model selection criterion. Under the MLE setting, the network minimizes the negative log-likelihood (NLL) of the SCB distribution, employing strict numerical clamping within the logarithmic terms to prevent gradient explosion when angular differences approach −1. Conversely, in the MCPE setting, the model operates as a geometric point-predictor by minimizing the cosine distance between the ground truth and the theoretical mean direction derived via the distribution’s trigonometric moments. 4. LEP-mGvM: The Lifted Expectation Propagation for the multivariate Generalised von Mises (LEP-mGvM; Wu Navarro (2019)) model performs approximate Bayesian inference by mapping the N -dimensional circular target space into a relaxed 2N -dimensional Euclidean space (representing the cosine and sine components) using Gaussian site approximations. While the theoretical LEP-mGvM framework dictates that the EP moment-matching step must compute the exact trigonometric moments of the GvM distribution (a process relying on unstable infinite series of modified Bessel functions), our implementation adopts a numerically robust simplification. Specifically, we utilize a first-order von Mises approximation by discarding the second harmonic parameter (κ2 ) and employing exponentially scaled Bessel ratios, which maintains computational stability while preserving the strict annealing and damping schedules mandatory for EP convergence. Additionally, for out-of-sample predictions, the implemented model approximates the full 2N joint Gaussian posterior conditioning with a streamlined kernel projection of the learned site means. Regarding input compatibility, the theoretical Gaussian Process formulation natively ingests mixed linear and circular covariates through the construction of composite covariance functions. However, because our implementation relies exclusively on a standard Euclidean squared-exponential kernel, circular covariates must first be exogenously embedded into continuous orthogonal components (sine and cosine) prior to passing them to the model. 5. Multiple Circular Regression II (MCR2): The Multiple Circular Regression II (MCR2) model (Jha & Biswas 2017) is a geometric framework that predicts a circular response from multiple circular covariates by employing Möbius transformations. Theoretically, it linearly combines complex-valued circular predictors using convex mixing weights and rotation parameters, projects the aggregate onto the unit circle, applies a Möbius transformation, and performs a final global rotation to yield the predicted mean direction. In our PyTorch implementation, the convex weights are enforced via a softmax activation, and the first covariate’s rotation is strictly anchored to zero to guarantee theoretical identifiability. To ensure stable gradient-based optimization, the complex arithmetic of the Möbius transformation is decoupled into explicit real and imaginary components, and the network is trained by minimizing the circular distance, an objective equivalent to maximizing the von Mises log-likelihood. The MCR2 architecture is rigidly formulated for angular inputs and does not natively process continuous linear covariates. 57

6. vMQP: The von Mises Quasi-Process (vMQP) is a non-parametric Bayesian regression model specifically designed for circular data (Cohen et al. 2025). It is formulated by taking a Gaussian process over two-dimensional Euclidean random functions and conditioning it on the unit circle. Unlike alternative Gaussian process-based adaptations for directional data that rely on wrapping or radial marginalization, the vMQP yields a simple, maximumentropy density. Because the vMQP is not mathematically consistent under marginalization, it operates within a transductive learning framework where test locations are fixed during training, which motivates a fully Bayesian approach to parameter learning. To handle posterior inference efficiently, the model introduces a Stratonovich-like augmentation that linearizes the trigonometric dependencies, allowing for fast Markov Chain Monte Carlo (MCMC) sampling. We utilize the original Bayesian implementation of vMQP7 for comparison in the wind direction prediction task, where it is observed to incur the maximum training time. We thus omit it for simulation experiments and object pose estimation task.

C

Evaluation Metrics

In this section, we describe the metrics tailored for circular data used in model evaluation. Let n denote the total number of observations, θi ∈ [0, 2π) be the ground-truth circular response for the i-th observation, and θbi be the corresponding predicted mean direction. Let {θbi,m }M m=1 represent an ensemble of M predicted angles for the i-th observation generated from the trained model. The mathematical formulations for the metrics are given as follows. 1. Mean Absolute Angular Deviation (MAAD): MAAD computes the average magnitude of the shortest angular path separating the predictions from the ground truth. To handle phase wrapping at the [0, 2π) boundary, the shortest angular distance da (·, ·) is computed using the two-argument arctangent function: 













da θi , θbi = atan2 sin θi − θbi , cos θi − θbi )

.

(C.1)

The MAAD is then obtained by averaging this absolute wrapped difference over all samples (Prokudin et al. 2018): n   1X MAAD = da θi , θbi . n i=1 We report this metric in degrees for enhanced interpretability. 2. Circular Mean Directional Error (CMDE): CMDE evaluates the alignment between the predicted and true angles by leveraging the cosine of the angular difference. This metric avoids the artificial discontinuities associated with linear distance metrics evaluated on a circle. It is bounded between 0 (indicating perfect alignment) and 2 (indicating perfectly opposite directions). The formulation is given by: CMDE = 1 − 7

n   1X cos θi − θbi . n i=1

https://github.com/Yarden231/vMQP

58

3. Circular Continuous Ranked Probability Score (CRPS): To evaluate the sharpness of the predictive distributions (or ensembles), we adapt the CRPS for circular topologies. Using the shortest angular distance function (C.1), the empirical circular CRPS for a single observation i approximated by an ensemble of size M is defined as (Gneiting & Raftery 2007): CRPSi =

M M M     X X 1 X 1 da θbi,m , θi − da θbi,m , θbi,m′ . M m=1 2M (M − 1) m=1 m′ =1,m′ ̸=m

The first term quantifies the average angular distance between the ensemble members and the true observation (rewarding accuracy), while the second term quantifies the expected angular distance among the ensemble members themselves. The overall Mean Circular CRPS is simply the average across all observations: Mean Circular CRPS =

n 1X CRPSi . n i=1

4. Median Error (MedErr): The MedErr metric (Prokudin et al. 2018) provides a robust measure of central tendency for the prediction errors, mitigating the influence of extreme outliers. It is defined as the median of the shortest angular distances across all n samples. For reporting purposes, the final value is converted from radians to degrees: 



MedErr = mediani∈{1,...,N } da θi , θbi



,

where da denotes the shortest angular distance given in (C.1). 5. Accuracyθe: The ‘Threshold Accuracy’ metric quantifies the proportion of predictions that fall e It is formulated within an acceptable margin of error, defined by a predetermined threshold θ. as the mean of an indicator function over the dataset: Accθe =

N     1 X 1 da θi , θbi ≤ θe , N i=1

where 1(·) represents the indicator function, evaluating to 1 if the condition da ≤ θe is satisfied, and 0 otherwise. In our standard evaluation protocol, we set the threshold to θe = π/6 (30 degrees), following prior works (Prokudin et al. 2018, Tulsiani & Malik 2015).

D

Software Demonstration

We provide a highly modular open-source software implementation for ANGLE utilizing a PyTorch backend in Python. The source code is publicly available8 and we offer the software through our Python package anglepy9 . The primary entry point for model training is the ANGLE wrapper function, which instantiates and trains a CircularEngressor object. As mentioned in Section 6, 8 9

https://github.com/anglepy Documentation: https://anglepy.readthedocs.io

59

our implementation provides extensive flexibility regarding the network architecture, the projection of the output space, the formulation of the GCES loss, and the underlying noise injection mechanism. The input data is expected to be array-like in the form of PyTorch tensors, and training and inference can be performed in batches on both GPU and CPU.

D.1

Architectural Configurations and Sufficient Dimension Reduction

The model supports standard fully connected dense networks, which can be augmented with batch normalization (add_bn=True) and residual skip connections (resblock=True). Furthermore, the software provides native support for sufficient dimension reduction (SDR) within the neural network framework. Setting sdr=True restricts the first layer to act as a linear projection. Setting reduced_dim=1 trains a single-index model, whereas reduced_dim=b trains a multipleindex model. Upon convergence, the fitted projection matrix β ∈ Rd×b can be extracted via the beta_proj attribute. This allows the user to explicitly compute the SDR subspace β ⊤ X for downstream analysis.

D.2

Output Heads

To map unconstrained neural network outputs to the unit circle S1 ∼ = [0, 2π), the software accommodates three distinct output modes, governed by the unbounded and circular_projection arguments: • Biternion head (unbounded=False, circular_projection=’atan2’): The network’s final layer outputs a two-dimensional vector (u, v) ∈ R2 . This vector is then transformed to an angular coordinate in [0, 2π) using the two-argument arctangent function, atan2(v, u). • Scaled sigmoid (unbounded=False, circular_projection=’sigmoid’): The network applies a scaled sigmoid activation at the final layer to bound the outputs strictly within [0, 2π). • Unbounded modulo (unbounded=True): The network predicts an unconstrained scalar value y ∈ R. Before loss computation and during inference, the output is mapped to the circle via the modulo operation y (mod 2π).

D.3

Distance Metrics in GCES Loss

The core optimization objective relies on strictly proper scoring rules, specifically the GCES, which requires the specification of a distance metric on the circle. The dist_method argument allows toggling between two fundamental distances. dist_method=’chordal’ utilizes the Euclidean distance between points embedded on the unit circle. dist_method=’geodesic’ computes the shortest path along the manifold. When using the geodesic distance, the user can supply a custom (negative) kernel function via kernel_func. The kernels.py module implements the strictly positive definite kernels on the circle proposed by Gneiting (2013, Table 1).

60

D.4

Noise Injection

As a generative model, ANGLE relies on random noise to model conditional distributions. The framework allows precise control over this stochasticity. The base distribution of the noise is set via noise_dist (accepting ’gaussian’ or ’uniform’), with its scale modulated by the noise_std parameter. The noise_all_layer boolean dictates the structural causal assumption. Setting noise_all_layer=True implements a pre-ANM framework by injecting noise into the early layers or all layers. Conversely, setting it to False strictly injects noise only at the final hidden layer, yielding a post-ANM formulation.

D.5

Usage Example: Training and Inference

The following Python code snippet demonstrates how to instantiate a realistic testing scenario, train the model using a subset of the aforementioned features (e.g., SDR, biternion projection, and chordal distance), and perform both point prediction and conditional distribution estimation via sampling.

61

Listing 1: Demonstration of training and inference using the ANGLE API. 1 2

import torch from anglepy import ANGLE

3 4

5

# Assume X_train (N x d), y_train (N x 1), X_test (M x d) are torch. Tensors device = torch.device("cuda" if torch.cuda.is_available() else "cpu")

6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

# Train the ANGLE Model model = ANGLE( x=X_train, y=y_train, num_layer=3, hidden_dim=100, noise_dim=64, noise_std=1.0, noise_dist='gaussian', noise_all_layer=True, add_bn=True, resblock=False, sdr=True, reduced_dim=1, unbounded=False, circular_projection='atan2', dist_method='chordal', gamma=1.0, lr=0.05, num_epochs=500, standardize=True, device=device, verbose=False )

# Pre-ANM noise injection

# Enable SDR # Single-index model # Biternion output head # Chordal distance for GCES loss # Energy score power parameter

31 32 33

# Extract the SDR projection matrix (d x 1) beta_matrix = model.beta_proj

34 35 36

37

# Inference: Point Predictions and Probabilistic Sampling # Point estimates (e.g., conditional circular mean, median, or quantiles from 100 samples) y_pred = model.predict(X_test, target='mean', sample_size=100)

38 39 40

41

# Draw 100 conditional samples per test observation # Returns tensor of shape (M, 100) capturing the conditional distribution y_pred_ensemble = model.sample(X_test, sample_size=100)

62

E Literature Review on Estimation of Circular Conditional Functionals E.1

Conditional Mean Function

The mean regression task is the most explored in the literature for circular response. One of the earliest approaches involve modeling the circular response using a parametric link function (Fisher & Lee 1992), where the conditional distribution of the circular response is assumed to be von Mises distribution and the conditional mean function is related to the covariate through a parametric link. A relatively more popular approach was proposed by Presnell et al. (1998), which treats the circular response as a projection of unobserved response from a multivariate linear model onto S1 . Different parametrization and estimation of such projected normal distributions can be found in HernandezStumpfhauser et al. (2017), Zou & Stern (2025) among others. In a parallel approach, Lund (2002) devised the regression function through a tree-based approach, and it was utilized by Lang et al. (2020) for wind direction forecasting. Another very common method of circular regression was proposed by Kato et al. (2008), where the regression curve is expressed as a form of Möbius circle transformation. Later, the approach was further extended to various circular data applications by works like Jha & Biswas (2018), Biswas & Banerjee (2025). Recently, Fernández-Durán & Gregorio-Domínguez (2024) utilized a non-negative trigonometric sums (NNTS) model for circular data regression where the idea is to map the circular response to a vector of its first few trigonometric moments, and then exploit the hypersphere geometry of this space to convert the problem to a linear regression problem on a real line. Several other Bayesian approaches can be found in Nur Ibrahim & Xu (2024), Ye et al. (2026), and some machine learning-based modeling of the mean regression can be found in Laha & Majumdar (2022), Bruns et al. (2024). A comprehensive review of parametric models for circular variables can be found in Kim & SenGupta (2016), Mohammad et al. (2021), Jha & Biswas (2022). Nonparametric smoothing based regression approaches have also gained much popularity in estimating the conditional mean function in the mean time. In Di Marzio et al. (2013), the regression function is defined as the minimizer of the expected cosine distance between the observed response and the fitted value. Under this formulation, the regression function can be expressed as the arctangent of the ratio of two conditional trigonometric moments, which are subsequently estimated through local averaging of the sine and cosine transforms of the observed responses. Subsequently, Di Marzio et al. (2014) extended this framework using local polynomial fitting techniques for spherical data. More generally, a common nonparametric strategy in the circular data literature is to model the conditional sine and cosine components separately via smoothing procedures and then combine the resulting estimates to recover the conditional direction; see, for example, Nguyen et al. (2023), Meilan-Vila et al. (2024), Woolsey & Huang (2025). Several general additive and spline based models are also available in the literature, see McMillan et al. (2013), Hassanzadeh (2021) among others. In more recent contributions, Gottard et al. (2026) proposed modeling the circular response through the quasi-likelihood score instead of assuming the full conditional distribution, and Francisco-Fernández & Meilán-Vila (2026) utilized the Nadaraya Watson local constant estimator for a setup with circular response and a mix of continuous and categorical predictors. A comprehensive review of nonparametric regression for circular responses can be found in Alonso-Pena & Crujeiras (2025b).

63

E.2

Conditional Mode and Concentration Functions

Most proposed circular regression models are primarily concerned with modeling the conditional mean direction, with Fisher & Lee (1992) being a notable exception in additionally allowing the conditional dispersion to depend on the covariates. Likewise, comparatively little attention has been devoted to estimating the conditional modal direction, despite the fact that the mode often provides a more meaningful notion of location in multimodal or highly diffuse settings, where the mean direction may correspond to a direction that is rarely observed in practice. Applications where the estimation of the conditional mode and concentration are particularly relevant include modeling the flight orientation of nocturnal migratory songbirds as a function of wind direction and flight altitude (Sjöberg & Nilsson 2015), modeling wind direction at noon using early morning measurements (Kato & Jones 2010), modeling the onset date of acute primary angle-closure glaucoma as a function of patient age (Gao et al. 2006), and studying the timing of tropical tree growth relative to annual climatic cycles (Hogan et al. 2019). Motivated by such applications, Alonso-Pena & Crujeiras (2023), Pérez-González et al. (2026) developed nonparametric estimators for the conditional modal direction. More recently, Ameijeiras-Alonso & Gijbels (2026) proposed semiparametric modeling of the conditional concentration and modal direction through a flexible unimodal family of circular densities introduced by Ameijeiras-Alonso et al. (2022).

E.3

CDF, Quantiles, and Density Functions

Estimating cumulative distribution functions and quantiles for circular data presents several challenges that do not arise in the Euclidean setting. The fundamental difficulty stems from the periodic geometry of the circle: unlike the classical CDF on the real line, a circular CDF requires an arbitrary choice of origin or lower integration limit, and this arbitrariness propagates to downstream tasks such as quantile estimation. Consequently, circular quantiles must be defined in a manner that is invariant to the choice of origin and applicable across general circular distributions. Early contributions in this area were made by Di Marzio et al. (2012), who introduced nonparametric estimation of the circular CDF and proposed quantile estimation via inversion of the estimated CDF and through a Parzen-type smoothing approach. Later, Di Marzio, Panzera & Taylor (2016) proposed circular quantile regression using either inversion of a double-kernel CDF estimator or minimization of a smooth circular check loss analogous to the pinball loss for Euclidean quantile regression (Steinwart & Christmann 2011). In a related work, Di Marzio, Fensore, Panzera & Taylor (2016) studied nonparametric estimation of the conditional circular density. More recently, Klar et al. (2024) developed estimators of circular distribution and density functions using Fejér polynomials and additionally established validity under error-in-covariate settings. Alternative notions of circular quantiles have also been proposed through angular depth (Ley et al. 2014, Hauch & Redenbach 2024) and optimal transport formulations (Hallin et al. 2024, Bercu et al. 2024).

64

F

Tables and Figures

Table A.4: Parameter configurations for the synthetic circular data generation mechanisms. The table details the number of linear (pL ) and circular (pC ) covariates, along with their corresponding true regression coefficient vectors (β), across all 16 experimental settings. Approach

Setting

# Linear covariates

# Circular covariates

1: Arctangent link, linear relation

1.1

2

0

1.2

0

2

1.3

2

2

1.4

75

25

2.1

2

0

2.2

0

2

2.3

2

2

2.4

75

25

3.1

2

0

3.2

0

2

3.3

2

2

3.4

75

25

4.1

2

0

4.2

0

2

4.3

2

2

4.4

75

25

2: Arctangent link, nonlinear relation

3: Projected normal distribution, linear relation

4: Projected normal distribution, nonlinear relation

65

Parameters β1 = [−5.3, −7.9], β2 = [−2.8, 4.1] β1 = [13.5, −6.6, 8.1, 6.7], β2 = [−5.2, 7.4, 6.9, −5.5] β1 = [7.3, 6.6, 2.5, 1.1, 10.3, 5.9], β2 = [6.7, −8.1, −3.4, 2.7, −6.5, 5.7] β1 , β2 ∈ R125 β1,j , β2,j ∼ U(−10, 10) β1 = [−5.3, −7.9], β2 = [−2.8, 4.1] β1 = [13.5, −6.6, 8.1, 6.7], β2 = [−5.2, 7.4, 6.9, −5.5] β1 = [7.3, 6.6, 2.5, 1.1, 10.3, 5.9], β2 = [6.7, −8.1, −3.4, 2.7, −6.5, 5.7] β1 , β2 ∈ R125 β1,j , β2,j ∼ U(−10, 10) β1 = [8.3, 4.6], β2 = [6.5, −5.1] β1 = [13.5, −8.6, 12.1, 6.9], β2 = [9.2, −8.1, −5.3, 2.2] β1 = [8.3, 4.6, 2.5, 1.1, 12.3, 5.2], β2 = [6.5, −5.1, −3.4, 2.7, −6.5, 4.7] β1 , β2 ∈ R125 β1,j , β2,j ∼ U(−10, 10) β1 = [8.3, 4.6], β2 = [6.5, −5.1] β1 = [13.5, −8.6, 12.1, 6.9], β2 = [9.2, −8.1, −5.3, 2.2] β1 = [8.3, 4.6, 2.5, 1.1, 12.3, 5.2], β2 = [6.5, −5.1, −3.4, 2.7, −6.5, 4.7] β1 , β2 ∈ R125 β1,j , β2,j ∼ U(−10, 10)

Table A.5: Model performances for prediction on simulated data: mean and standard deviation of three circular metrics across all 16 data generation settings. The data generation mechanisms are summarized in Table A.4. The best and second best results are highlighted. The mean and standard deviation obtained across 50 independent training-inference runs are reported. MAAD and CRPS are reported in degrees. Model

MAAD

CMDE

CRPS

MAAD

Setting 1.1 3.370

0.007

CMDE

CRPS

MAAD

Setting 1.2 3.152

4.253

0.010

CMDE

CRPS

MAAD

Setting 1.3 3.609

3.511

0.007

CMDE

CRPS

Setting 1.4 3.172

4.345

0.012

3.187

ANGLE-Geodesic (±0.507) (±0.003) (±0.655) (±0.731) (±0.003) (±0.867) (±0.770) (±0.004) (±0.829) (±0.540) (±0.004) (±0.422) ANGLE-Chordal

3.221 2.725 3.939 0.010 2.996 2.652 0.012 0.007 3.194 0.005 4.308 3.238 (±0.318) (±0.002) (±0.399) (±0.377) (±0.001) (±0.382) (±0.320) (±0.002) (±0.315) (±0.525) (±0.004) (±0.438)

CLR

0.006 3.634 3.634 (±0.020) (±0.000) (±0.020)

2.948 0.004 3.813 5.044 0.017 5.044 2.948 0.013 3.813 (±0.032) (±0.000) (±0.032) (±0.001) (±0.000) (±0.001) (±0.003) (±0.000) (±0.003)

Kernel CLR

5.321 0.010 5.321 (±0.000) (±0.000) (±0.000)

31.444 (±0.000)

0.244 (±0.000)

31.444 (±0.000)

19.400 0.123 19.400 64.079 0.643 64.079 (±0.000) (±0.000) (±0.000) (±0.000) (±0.000) (±0.000)

SCB-NLL

33.028 0.230 33.028 (±0.069) (±0.000) (±0.069)

37.151 (±1.123)

0.277 (±0.014)

37.151 (±1.123)

46.040 0.373 46.040 44.505 0.363 44.505 (±0.687) (±0.011) (±0.687) (±0.647) (±0.009) (±0.647)

SCB-MCPE

33.352 0.238 33.352 (±0.686) (±0.008) (±0.686)

34.360 (±1.237)

0.244 (±0.010)

34.360 (±1.237)

46.441 0.378 46.441 46.441 0.378 46.441 (±0.540) (±0.005) (±0.540) (±0.540) (±0.005) (±0.540)

Lifted EP-MGvM

6.781 0.012 6.781 (±0.000) (±0.000) (±0.000)

25.661 (±0.000)

0.175 (±0.000)

25.661 (±0.000)

21.666 0.126 21.666 67.727 0.696 67.727 (±0.000) (±0.000) (±0.000) (±0.000) (±0.000) (±0.000)

MCR-II

41.277 0.352 41.277 (±0.000) (±0.000) (±0.000)

30.682 (±0.000)

0.209 (±0.000)

30.682 (±0.000)

46.074 0.419 46.074 43.826 0.389 43.826 (±1.152) (±0.013) (±1.152) (±2.180) (±0.027) (±2.180)

Setting 2.1 8.586

0.032

Setting 2.2 7.128

8.104

0.038

Setting 2.3 6.342

8.194

0.047

Setting 2.4 6.214

4.647

0.015

3.424

ANGLE-Geodesic (±3.017) (±0.025) (±1.942) (±1.184) (±0.008) (±0.876) (±1.392) (±0.014) (±0.886) (±0.587) (±0.005) (±0.464) ANGLE-Chordal

5.718 0.014 0.015 4.607 6.703 0.029 5.111 6.411 0.030 4.775 4.578 3.478 (±1.060) (±0.007) (±0.783) (±0.784) (±0.006) (±0.597) (±1.428) (±0.014) (±0.855) (±0.537) (±0.004) (±0.462)

CLR

19.281 0.144 19.281 (±0.000) (±0.000) (±0.000)

10.769 (±0.009)

0.073 (±0.000)

10.769 (±0.009)

3.865 10.422 0.072 10.422 0.015 3.865 (±0.000) (±0.000) (±0.000) (±0.002) (±0.000) (±0.002)

Kernel CLR

10.617 0.045 10.617 (±0.000) (±0.000) (±0.000)

31.988 (±0.000)

0.267 (±0.000)

31.988 (±0.000)

19.437 0.128 19.437 64.089 0.646 64.089 (±0.000) (±0.000) (±0.000) (±0.000) (±0.000) (±0.000)

SCB-NLL

28.335 0.189 28.335 (±0.179) (±0.002) (±0.179)

32.441 (±1.343)

0.219 (±0.012)

32.441 (±1.343)

33.611 0.235 33.611 37.087 0.273 37.087 (±0.477) (±0.005) (±0.477) (±1.055) (±0.012) (±1.055)

SCB-MCPE

27.815 0.186 27.815 (±3.268) (±0.041) (±3.268)

29.781 (±0.220)

0.197 (±0.003)

29.781 (±0.220)

33.259 0.229 33.259 41.401 0.325 41.401 (±3.232) (±0.040) (±3.232) (±3.365) (±0.045) (±3.365)

Lifted EP-MGvM

16.934 0.115 16.934 (±0.000) (±0.000) (±0.000)

27.560 (±0.000)

0.216 (±0.000)

27.560 (±0.000)

24.392 0.172 24.392 67.626 0.697 67.626 (±0.000) (±0.000) (±0.000) (±0.000) (±0.000) (±0.000)

MCR-II

41.435 0.346 41.435 (±0.000) (±0.000) (±0.000)

28.118 (±0.002)

0.181 (±0.000)

28.118 (±0.002)

41.959 0.357 41.959 43.495 0.383 43.495 (±0.244) (±0.003) (±0.244) (±0.828) (±0.011) (±0.828)

Setting 3.1 7.247

0.025

Setting 3.2 5.296

5.825

0.013

Setting 3.3 4.386

5.602

0.016

Setting 3.4 4.203

3.770

0.007

2.720

ANGLE-Geodesic (±0.340) (±0.002) (±0.374) (±0.319) (±0.001) (±0.427) (±0.274) (±0.001) (±0.339) (±0.366) (±0.003) (±0.236) ANGLE-Chordal

3.956 7.175 0.024 5.093 5.722 0.013 4.184 5.514 0.016 3.688 0.007 2.726 (±0.300) (±0.002) (±0.273) (±0.234) (±0.001) (±0.232) (±0.295) (±0.001) (±0.279) (±0.344) (±0.003) (±0.265)

CLR

7.306 0.026 7.306 (±0.002) (±0.000) (±0.002)

5.341 0.015 3.381 0.006 6.686 0.017 6.686 5.341 3.381 (±0.022) (±0.000) (±0.022) (±0.002) (±0.000) (±0.002) (±0.002) (±0.000) (±0.002)

Kernel CLR

0.024 7.345 7.345 (±0.000) (±0.000) (±0.000)

25.815 (±0.000)

0.172 (±0.000)

25.815 (±0.000)

20.692 0.130 20.692 68.379 0.709 68.379 (±0.000) (±0.000) (±0.000) (±0.000) (±0.000) (±0.000)

SCB-NLL

36.185 0.266 36.185 (±0.440) (±0.006) (±0.440)

30.673 (±0.143)

0.212 (±0.001)

30.673 (±0.143)

44.309 0.360 44.309 50.878 0.445 50.878 (±2.560) (±0.028) (±2.560) (±0.782) (±0.012) (±0.782)

SCB-MCPE

36.183 0.270 36.183 (±1.476) (±0.018) (±1.476)

28.712 (±0.048)

0.193 (±0.000)

28.712 (±0.048)

41.536 0.325 41.536 52.842 0.477 52.842 (±0.966) (±0.013) (±0.966) (±0.977) (±0.014) (±0.977)

Lifted EP-MGvM

8.228 0.027 8.228 (±0.000) (±0.000) (±0.000)

24.407 (±0.000)

0.130 (±0.000)

24.407 (±0.000)

21.155 0.116 21.155 85.440 0.935 85.440 (±0.000) (±0.000) (±0.000) (±0.000) (±0.000) (±0.000)

MCR-II

41.303 0.347 41.303 56.908 0.559 56.908 41.422 0.350 41.422 44.769 0.397 44.769 (±0.772) (±0.010) (±0.772) (±13.448) (±0.185) (±13.448) (±0.549) (±0.005) (±0.549) (±1.672) (±0.022) (±1.672)

Continued on next page

66

Table A.5 – continued from previous page Model

MAAD

CMDE

CRPS

MAAD

CMDE

Setting 4.1 8.017

0.034

CRPS

MAAD

Setting 4.2 6.194

0.033

7.801

CMDE

CRPS

MAAD

Setting 4.3 5.794

6.004

0.014

CMDE

CRPS

Setting 4.4 4.951

3.750

0.007

2.696

ANGLE-Geodesic (±0.478) (±0.004) (±0.513) (±0.374) (±0.002) (±0.369) (±0.715) (±0.005) (±0.529) (±0.357) (±0.003) (±0.238) ANGLE-Chordal

7.775 0.033 5.668 7.576 0.033 5.478 5.704 0.012 4.336 3.678 0.007 2.711 (±0.527) (±0.004) (±0.417) (±0.325) (±0.002) (±0.274) (±0.375) (±0.002) (±0.370) (±0.357) (±0.003) (±0.260)

CLR

8.404 0.038 8.404 (±0.003) (±0.000) (±0.003)

9.143 0.038 9.143 (±0.018) (±0.000) (±0.018)

0.005 3.163 3.163 6.472 0.017 6.472 (±0.000) (±0.000) (±0.000) (±0.002) (±0.000) (±0.002)

Kernel CLR

8.356 0.039 8.356 (±0.000) (±0.000) (±0.000)

27.293 (±0.000)

0.216 (±0.000)

27.293 (±0.000)

19.893 0.123 19.893 68.752 0.714 68.752 (±0.000) (±0.000) (±0.000) (±0.000) (±0.000) (±0.000)

SCB-NLL

28.505 0.188 28.505 (±0.293) (±0.005) (±0.293)

30.240 (±0.176)

0.213 (±0.002)

30.240 (±0.176)

34.841 0.248 34.841 48.466 0.407 48.466 (±0.422) (±0.005) (±0.422) (±0.849) (±0.013) (±0.849)

SCB-MCPE

28.628 0.189 28.628 (±0.306) (±0.004) (±0.306)

30.246 (±0.547)

0.217 (±0.005)

30.246 (±0.547)

35.607 0.259 35.607 53.406 0.485 53.406 (±4.107) (±0.048) (±4.107) (±1.074) (±0.015) (±1.074)

Lifted EP-MGvM

10.016 0.048 10.016 (±0.000) (±0.000) (±0.000)

27.938 (±0.000)

0.175 (±0.000)

27.938 (±0.000)

22.360 0.132 22.360 85.366 0.935 85.366 (±0.000) (±0.000) (±0.000) (±0.000) (±0.000) (±0.000)

MCR-II

39.088 0.318 39.088 (±0.000) (±0.000) (±0.000)

58.865 (±4.375)

0.587 (±0.061)

58.865 (±4.375)

39.400 0.315 39.400 44.962 0.401 44.962 (±0.440) (±0.007) (±0.440) (±1.724) (±0.023) (±1.724)

Table A.6: Results on the wind direction prediction task. The best and second best results are highlighted. The mean and standard deviation obtained across 50 independent training-inference runs are reported. For the proposed model using energy score loss with geodesic distance, the C 2 -Wendland kernel is used. MedErr, CRPS, and MAAD are reported in degrees.

India wind

Germany Wind

Dataset

Metric

vMQP

CLR

Kernel CLR

SCBR

Lifted EP-MGvM

MCR2

ANGLE-Geodesic

ANGLE-Chordal

Accuracy

0.85 (±0.01)

0.90 (±0.00)

0.90 (±0.00)

0.90 (±0.00)

0.90 (±0.02)

0.89 (±0.02)

0.90 (±0.03)

0.89 (±0.02)

MedErr

12.74 (±0.47)

12.66 (±0.18)

16.05 (±0.00)

13.01 (±0.02)

13.18 (±2.87)

13.28 (±0.49)

12.52 (±0.78)

12.63 (±0.87)

CRPS

12.61 (±0.11)

16.06 (±0.01)

18.12 (±0.00)

16.54 (±0.00)

15.59 (±1.53)

16.95 (±0.76)

11.47 (±0.20)

11.61 (±0.14)

MAAD

16.49 (±0.21)

16.06 (±0.01)

18.12 (±0.00)

16.54 (±0.00)

15.59 (±1.53)

16.95 (±0.76)

15.86 (±0.41)

16.29 (±0.27)

CMDE

0.06 (±0.00)

0.06 (±0.00)

0.07 (±0.00)

0.07 (±0.00)

0.07 (±0.01)

0.07 (±0.01)

0.06 (±0.00)

0.06 (±0.00)

Accuracy

0.03 (±0.00)

0.35 (±0.00)

0.68 (±0.00)

0.41 (±0.03)

0.77 (±0.00)

0.33 (±0.01)

0.72 (±0.04)

0.78 (±0.03)

MedErr

130.35 (±0.73)

43.15 (±0.01)

17.83 (±0.00)

41.41 (±1.06)

13.17 (±0.00)

41.25 (±1.25)

15.76 (±1.74)

13.04 (±1.43)

CRPS

82.63 (±2.25)

47.49 (±0.00)

28.24 (±0.00)

47.63 (±0.46)

23.24 (±0.00)

49.24 (±0.40)

19.14 (±1.95)

17.13 (±1.32)

MAAD

122.83 (±0.27)

47.49 (±0.00)

28.24 (±0.00)

47.63 (±0.46)

23.24 (±0.00)

49.24 (±0.40)

25.81 (±2.50)

23.18 (±1.60)

CMDE

1.46 (±0.00)

0.41 (±0.00)

0.20 (±0.00)

0.43 (±0.01)

0.15 (±0.00)

0.43 (±0.01)

0.18 (±0.03)

0.15 (±0.02)

67

Table A.7: Model performances on the object pose detection task. The best and second-best results are highlighted. The mean and standard deviation obtained across 50 independent training–inference runs are reported. MedErr, CRPS, and MAAD are reported in degrees.

Accuracyπ/6

Metric Model

aeroplane

bicycle

boat

Mixture-Inception

0.28

0.48

0.38

0.99

0.7

0.33

0.36

0.65

0.3

0.57

Mixture-DenseNet

0.53

0.27

0.4

0.99

0.61

0.39

0.39

0.63

0.31

0.65

Mixture-MobileNet

0.16

0.29

0.41

0.99

0.63

0.39

0.36

0.63

0.31

0.64

0.33 (±0.00)

0.27 (±0.01)

0.36 (±0.01)

0.99

Inception-v3-SCBR

0.60 (±0.01)

0.31 (±0.00)

Inception-v3-LEP-mGvM

0.60 (±0.00)

0.48 (±0.00)

0.49 (±0.00) 0.68

0.39 (±0.00)

0.68 (±0.00)

0.20 (±0.04) 0.80 (±0.00)

(±0.00)

0.58 (±0.00)

0.65 (±0.00)

Inception-v3-MCR2

0.60 (±0.01)

0.57 (±0.01)

0.44 (±0.01)

0.77 (±0.01)

0.65 (±0.01)

0.58 (±0.00)

ConvNeXt-SCBR

0.34 (±0.00)

0.28 (±0.00)

0.40 (±0.01)

0.50 (±0.00)

ConvNeXt-LEP-mGvM

0.55 (±0.00)

0.54 (±0.00)

0.48 (±0.00)

0.20 (±0.05) 0.80

ConvNeXt-MCR2

0.64 (±0.02)

0.55 (±0.04)

0.46 (±0.01)

Inception-v3-SDR-ANGLE

0.66 (±0.01)

0.59 (±0.01)

0.45 (±0.01)

0.67 (±0.01) 0.70

0.60 (±0.01) 0.61

0.47 (±0.01) 0.49

(±0.01)

(±0.01)

Inception-v3-ANGLE ConvNeXt-ANGLE

MedErr

Mean

0.66

0.69

0.53

0.65

0.64

0.54

0.66

0.67

0.51

0.62 (±0.00)

0.55 (±0.02)

0.66 (±0.00)

0.48

0.69 (±0.00)

0.73 (±0.00)

0.67 (±0.00)

0.73 (±0.00)

0.69

0.66 (±0.00)

0.68 (±0.01)

0.74 (±0.00)

0.62 (±0.01)

0.73 (±0.01)

0.67

0.39 (±0.00)

0.66 (±0.00)

0.31 (±0.00)

0.63 (±0.00)

0.66 (±0.00)

0.5

(±0.00)

0.63 (±0.00)

0.43 (±0.00)

0.66 (±0.00)

0.53 (±0.00)

0.70 (±0.00)

0.58 (±0.01) 0.71 (±0.00)

0.72 (±0.00)

0.65

0.70 (±0.02)

0.67 (±0.01)

0.58 (±0.01)

0.63 (±0.01)

0.57 (±0.07)

0.72 (±0.01)

0.59 (±0.02)

0.74 (±0.01)

0.65

(±0.00) 0.99

0.79 (±0.01) 0.80

0.65 (±0.00)

0.64 (±0.04) 0.67

0.66 (±0.04) 0.70

0.75 (±0.01)

0.68

(±0.00) 0.99

(±0.01) 0.80

(±0.00) 0.60

(±0.01) 0.67

(±0.01) 0.71

0.76 (±0.00) 0.77

0.70 (±0.01) 0.71

0.70 (±0.00)

0.66 (±0.00) 0.68

0.59 (±0.00) 0.60

(±0.01)

(±0.00)

(±0.00)

(±0.00)

(±0.00)

(±0.01)

(±0.01)

(±0.01)

(±0.00) 0.99 (±0.00) 0.99 (±0.00) 0.99 (±0.00) 0.99 (±0.00) 0.99

car

chair

table

motorbike

sofa

train

(±0.01) 0.69 (±0.01)

0.74 (±0.01) 0.76

0.70

(±0.01)

0.71

Mixture-Inception

69.89

32.55

48.24

0.96

8.9

70.64

53.62

8.51

58.83

25.38

18.52

16.31

34.36

27.63

67.03

50.35

0.62

19.88

59.3

42.16

5.9

59.84

16.42

17.63

17.45

32.02

Mixture-MobileNet

68.13

65.1

50.89

1.49

19.04

52.01

53.67

5.59

59.92

13.79

17.68

18.66

35.50

51.40 (±0.29) 34.80 (±0.00) 37.50 (±0.88) 48.71 (±0.33) 47.54 (±0.00) 34.05 (±1.29) 18.59 (±0.45) 17.39

57.07 (±0.43) 40.72 (±0.00) 41.99 (±0.53) 55.70 (±0.38) 50.85 (±0.00) 43.65 (±1.65) 21.24 (±0.64) 19.96

(±0.20)

36.05 (±0.25) 27.83 (±0.00) 29.33 (±0.10) 34.79 (±0.13) 26.07 (±0.00) 30.84 (±0.28) 18.13 (±4.24) 15.27 (±0.57)

(±0.49)

(±0.43)

21.03 (±0.50)

(±0.33)

17.66 (±0.50)

17.11 (±0.45)

11.97 (±0.36)

31.93 (±0.30) 26.92 (±0.00) 28.29 (±0.35) 31.55 (±0.11) 30.15 (±0.00) 28.67 (±0.33) 15.69 (±0.39) 14.67 (±0.28) 12.61

(±0.18)

22.80 (±0.31)

45.45 (±0.62) 43.35 (±0.00) 42.41 (±0.47) 45.15 (±0.23) 37.70 (±0.00) 44.32 (±0.83) 11.52 (±0.22) 12.05 (±0.76) 11.05

ConvNeXt-ANGLE

17.54 (±0.47)

50.77 (±0.21) 37.56 (±0.00) 39.33 (±0.70) 47.66 (±0.41) 46.52 (±0.00) 35.47 (±0.40) 15.25 (±0.21) 14.62 (±0.37) 14.32

29.22 (±0.12) 24.63 (±0.00) 24.91 (±0.66) 28.49 (±0.16) 24.85 (±0.00) 26.29 (±0.34) 13.16 (±1.06) 11.57

(±0.56)

54.85 (±1.51) 24.99 (±0.00) 28.13 (±0.46) 56.20 (±2.37) 24.56 (±0.00) 30.71 (±1.16) 6.95 (±0.22) 6.95 (±0.28) 6.12

52.10 (±0.46) 32.95 (±0.00) 32.98 (±0.42) 50.33 (±0.37) 51.45 (±0.00) 38.53 (±3.96) 18.02 (±2.81) 16.31

(±0.43)

2.91 (±0.00) 4.34 (±0.00) 3.86 (±0.10) 2.89 (±0.00) 3.39 (±0.00) 5.63 (±0.39) 2.19 (±0.31) 2.46 (±0.30)

54.97 (±0.20) 39.66 (±0.00) 38.99 (±0.12) 53.50 (±0.19) 57.75 (±0.00) 36.77 (±0.51) 23.10 (±0.30) 22.09

Inception-v3-ANGLE

65.84 (±0.52) 52.65 (±0.00) 54.47 (±0.40) 60.01 (±0.63) 54.67 (±0.00) 51.87 (±0.25) 36.29 (±1.55) 34.65 (±0.89) 30.72

(±0.62)

(±0.32)

15.47

Inception-v3-LEP-mGvM Inception-v3-MCR2 ConvNeXt-SCBR ConvNeXt-LEP-mGvM ConvNeXt-MCR2 Inception-v3-SDR-ANGLE

(±1.06)

1.64 (±0.19)

44.38 32.53 33.52 42.92 37.96 33.90 16.68 15.67

Mixture-Inception

43.38

36.25

41.75

2.82

29.49

42.96

40.77

22.32

41.99

23.35

33.54

22.47

31.76

Mixture-DenseNet

31

44.35

42.32

2.56

32.11

60.1

35.63

24.58

43.19

22.84

33.14

23.82

32.97

Mixture-MobileNet

43.8

53.33

42.22

2.95

33.03

49.58

52.75

23.86

44.06

24.91

35.64

26.35

36.04

51.40 (±0.29) 34.80 (±0.00) 37.50 (±0.88) 48.71 (±0.33) 47.54 (±0.00) 34.05 (±1.29) 23.03 (±0.16) 22.74 (±0.15) 22.55

57.07 (±0.43) 40.72 (±0.00) 41.99 (±0.53) 55.70 (±0.38) 50.85 (±0.00) 43.65 (±1.65) 27.54 (±0.22) 28.50 (±0.22) 25.36

65.84 (±0.52) 52.65 (±0.00) 54.47 (±0.40) 60.01 (±0.63) 54.67 (±0.00) 51.87 (±0.25) 36.14 (±0.42) 39.24 (±0.27) 33.70

2.91 (±0.00) 4.34 (±0.00) 3.86 (±0.10) 2.89 (±0.00) 3.39 (±0.00) 5.63 (±0.39) 2.53 (±0.07) 2.51

54.85 (±1.51) 24.99 (±0.00) 28.13 (±0.46) 56.20 (±2.37) 24.56 (±0.00) 30.71 (±1.16) 18.71 (±0.26) 20.23 (±0.21) 18.68

50.77 (±0.21) 37.56 (±0.00) 39.33 (±0.70) 47.66 (±0.41) 46.52 (±0.00) 35.47 (±0.40) 28.96 (±0.50) 29.46 (±0.28) 28.37

54.97 (±0.20) 39.66 (±0.00) 38.99 (±0.12) 53.50 (±0.19) 57.75 (±0.00) 36.77 (±0.51) 28.94 (±0.26) 30.02 (±0.18) 26.79

36.05 (±0.25) 27.83 (±0.00) 29.33 (±0.10) 34.79 (±0.13) 26.07 (±0.00) 30.84 (±0.28) 23.81 (±1.76) 21.82 (±0.12) 21.71

52.10 (±0.46) 32.95 (±0.00) 32.98 (±0.42) 50.33 (±0.37) 51.45 (±0.00) 38.53 (±3.96) 24.60 (±1.64) 22.54

45.45 (±0.62) 43.35 (±0.00) 42.41 (±0.47) 45.15 (±0.23) 37.70 (±0.00) 44.32 (±0.83) 36.82 (±0.12) 35.05 (±0.35)

31.93 (±0.30) 26.92 (±0.00) 28.29 (±0.35) 31.55 (±0.11) 30.15 (±0.00) 28.67 (±0.33) 23.45 (±0.39) 21.64 (±0.37) 20.23

(±0.21)

(±0.22)

(±0.23)

(±0.20)

(±0.24)

(±0.24)

(±0.16)

23.40 (±0.22)

29.22 (±0.12) 24.63 (±0.00) 24.91 (±0.66) 28.49 (±0.16) 24.85 (±0.00) 26.29 (±0.34) 18.67 (±1.10) 18.32 (±0.15) 16.75

Inception-v3-SCBR Inception-v3-LEP-mGvM CRPS

tvmonitor

(±0.00) 0.99

bus

Mixture-DenseNet

Inception-v3-SCBR

Inception-v3-MCR2 ConvNeXt-SCBR ConvNeXt-LEP-mGvM ConvNeXt-MCR2 Inception-v3-SDR-ANGLE Inception-v3-ANGLE ConvNeXt-ANGLE

MAAD

bottle

(±0.09) 2.70 (±0.12)

(±0.17)

(±0.33)

44.38 32.53 33.52 42.92 37.96 33.90 24.43 24.34

36.46 (±0.93)

(±0.61)

23.06

Mixture-Inception

72

50.84

62.87

3.03

38.52

80.91

64.33

27.76

70.77

31.25

41.3

31.26

47.90

Mixture-DenseNet

44.96

72.59

60.9

2.88

43.18

74.44

53.66

28.61

70.06

29.74

41.26

32.27

46.21

Mixture-MobileNet

77.42

74.02

61.66

3.24

42.8

72.09

64.42

29.3

69.18

28.6

41.32

32.69

49.73

Inception-v3-SCBR

51.40 (±0.29)

57.07 (±0.43)

65.84 (±0.52)

2.91 (±0.00)

54.85 (±1.51)

50.77 (±0.21)

54.97 (±0.20)

36.05 (±0.25)

52.10 (±0.46)

29.22 (±0.12)

45.45 (±0.62)

31.93 (±0.30)

44.38

Continued on next page

68

Table A.7 – continued from previous page Metric Model

aeroplane

bicycle

boat

bottle

bus

car

chair

table

motorbike

sofa

train

tvmonitor

34.80 (±0.00) 37.50 (±0.88) 48.71 (±0.33) 47.54 (±0.00)

40.72 (±0.00) 41.99 (±0.53) 55.70 (±0.38) 50.85 (±0.00)

52.65 (±0.00) 54.47 (±0.40) 60.01 (±0.63) 54.67 (±0.00)

4.34 (±0.00) 3.86 (±0.10) 2.89 (±0.00) 3.39 (±0.00)

24.99 (±0.00) 28.13 (±0.46) 56.20 (±2.37) 24.56 (±0.00)

37.56 (±0.00) 39.33 (±0.70) 47.66 (±0.41) 46.52 (±0.00)

39.66 (±0.00) 38.99 (±0.12) 53.50 (±0.19) 57.75 (±0.00)

27.83 (±0.00) 29.33 (±0.10) 34.79 (±0.13) 26.07

32.95 (±0.00) 32.98 (±0.42) 50.33 (±0.37) 51.45 (±0.00)

24.63 (±0.00) 24.91 (±0.66) 28.49 (±0.16) 24.85 (±0.00)

43.35 (±0.00) 42.41 (±0.47) 45.15 (±0.23) 37.70

26.92 (±0.00) 28.29 (±0.35) 31.55 (±0.11) 30.15 (±0.00)

34.05 (±1.29) 33.92 (±0.30) 32.88 (±0.29) 30.67

43.65 (±1.65) 40.14 (±0.57) 39.26 (±0.26) 36.93

51.87 (±0.25) 54.55 (±0.95) 53.92 (±0.55) 48.80

5.63 (±0.39) 3.93 (±0.23) 3.98 (±0.21)

35.47 (±0.40) 37.25 (±0.14) 37.09 (±0.26) 34.63

36.77 (±0.51) 37.57 (±0.21) 37.22 (±0.08) 35.39

38.53 (±3.96) 33.98 (±2.66) 31.80

26.29 (±0.34) 24.55 (±0.80) 23.66 (±0.25) 21.39

ConvNeXt-ANGLE

(±0.27)

(±0.41)

(±0.53)

3.41 (±0.18)

30.71 (±1.16) 25.84 (±0.26) 25.24 (±0.18) 22.40 (±0.24)

(±0.17)

Mixture-Inception

0.78

0.49

0.65

0.01

0.38

Mixture-DenseNet

0.41

0.78

0.65

0.01

Mixture-MobileNet

0.82

0.79

0.64

0.01

0.46 (±0.00)

0.53 (±0.01)

0.67 (±0.01)

0.01

Inception-v3-SCBR Inception-v3-LEP-mGvM

0.30 (±0.00)

0.37 (±0.00)

0.53 (±0.00)

Inception-v3-MCR2

0.32 (±0.01)

0.38 (±0.01)

0.53 (±0.01)

ConvNeXt-SCBR

0.42 (±0.00)

0.52 (±0.01)

0.59 (±0.01)

ConvNeXt-LEP-mGvM

0.48 (±0.00)

0.51 (±0.00)

0.57 (±0.00)

ConvNeXt-MCR2

0.28 (±0.02)

0.39 (±0.02)

0.49 (±0.00)

Inception-v3-SDR-ANGLE

0.28 (±0.00)

0.36 (±0.01)

0.54 (±0.01)

0.27 (±0.00) 0.24

0.35 (±0.00) 0.32

0.53 (±0.01) 0.46

(±0.00)

(±0.00)

(±0.01)

Inception-v3-LEP-mGvM Inception-v3-MCR2

MAAD

ConvNeXt-SCBR ConvNeXt-LEP-mGvM ConvNeXt-MCR2 Inception-v3-SDR-ANGLE

CMDE

Inception-v3-ANGLE

Inception-v3-ANGLE ConvNeXt-ANGLE

(±0.00) 30.84 (±0.28) 30.23 (±2.54) 28.07 (±0.18)

(±0.53)

(±0.18)

28.42 (±0.22)

32.30 (±0.42)

0.88

0.68

0.24

0.4

0.82

0.51

0.4

0.77

0.68

0.48 (±0.01)

0.49 (±0.00)

0.22 (±0.00)

(±0.00) 44.32 (±0.83) 39.17 (±0.10) 39.13 (±0.30)

28.67 (±0.33) 29.42 (±0.15) 27.74 (±0.12) 26.29

Mean 32.53 33.52 42.92 37.96 33.90 32.55 31.67

(±0.20)

39.76 (±0.68)

(±0.19)

30.03

0.74

0.24

0.39

0.25

0.48

0.24

0.74

0.24

0.38

0.26

0.45

0.24

0.74

0.24

0.38

0.26

0.50

0.53 (±0.00)

0.26 (±0.00)

0.46 (±0.01)

0.24 (±0.00)

0.40 (±0.00)

0.25 (±0.00)

0.40

0.36 (±0.00)

0.35 (±0.00)

0.23 (±0.00)

0.28 (±0.00)

0.19 (±0.00)

0.43 (±0.00)

0.21 (±0.00)

0.29

0.25 (±0.01)

0.37 (±0.01)

0.33 (±0.00)

0.22 (±0.00)

0.27 (±0.01)

0.17 (±0.00)

0.39 (±0.01)

0.22 (±0.00)

0.29

0.50 (±0.03)

0.46 (±0.00)

0.51 (±0.00)

0.25 (±0.00)

0.44 (±0.01)

0.23 (±0.00)

0.25 (±0.00)

0.38

0.24 (±0.00)

0.49 (±0.00)

0.60 (±0.00)

0.23 (±0.00)

0.53 (±0.00)

0.21 (±0.00)

0.40 (±0.00) 0.37 (±0.00)

0.24 (±0.00)

0.37

0.25 (±0.01)

0.32 (±0.00)

0.30 (±0.01)

0.23 (±0.00)

0.32 (±0.04)

0.17 (±0.00)

0.41 (±0.01)

0.22 (±0.00)

0.28

0.23 (±0.00)

0.34 (±0.00)

0.32 (±0.00)

(±0.00) 0.01

0.32 (±0.00) 0.29

(±0.00) 0.21

(±0.01) 0.26

0.17 (±0.00) 0.14

0.38 (±0.00) 0.37

0.23 (±0.00)

0.34 (±0.00) 0.31

0.28 (±0.03) 0.26

0.17 (±0.01)

0.22 (±0.00) 0.19

0.23 (±0.02) 0.21

(±0.00)

(±0.00)

(±0.00)

(±0.00)

(±0.00)

(±0.00)

(±0.00)

(±0.00) 0.01 (±0.00) 0.01 (±0.00) 0.01 (±0.00) 0.01 (±0.00) 0.01 (±0.00) 0.01 (±0.00) 0.01

69

(±0.00) 0.38 (±0.01)

0.21 (±0.00) 0.20 (±0.00)

0.28 0.27 0.25

Table A.8: Ablation study to evaluate the performance of different components in the proposed model on the Indian wind dataset. The negative of the kernels are taken from Gneiting (2013, Table 1). For each variant, the mean and standard deviation of the metrics are reported across 50 independent train-test runs. The best and second best results are highlighted. MedErr, CRPS, and MAAD are reported in degrees. Distance

Mode

Kernel

Output

Accuracyπ/6

MedErr

CRPS

MAAD

CMDE

Powered exponential

Modulated Scaled sigmoid atan2

0.506 ± 0.111 0.601 ± 0.025 0.583 ± 0.032

30.838 ± 12.758 23.300 ± 1.469 24.773 ± 2.355

28.314 ± 5.137 23.930 ± 0.533 24.969 ± 0.814

40.988 ± 10.284 33.319 ± 1.096 34.107 ± 1.336

0.347 ± 0.131 0.252 ± 0.011 0.258 ± 0.011

Generalized Cauchy

Modulated Scaled sigmoid atan2

0.492 ± 0.128 0.596 ± 0.028 0.556 ± 0.043

33.834 ± 17.438 23.353 ± 2.090 26.253 ± 2.940

29.117 ± 6.228 24.086 ± 0.766 25.548 ± 1.172

43.494 ± 14.470 33.573 ± 1.603 35.137 ± 1.620

0.380 ± 0.191 0.255 ± 0.016 0.266 ± 0.017

Dagum

Modulated Scaled sigmoid atan2

0.322 ± 0.076 0.463 ± 0.058 0.421 ± 0.055

52.150 ± 16.489 32.561 ± 3.952 35.461 ± 4.058

37.338 ± 4.756 27.884 ± 1.375 28.881 ± 1.575

58.019 ± 14.486 40.255 ± 2.376 41.498 ± 2.821

0.559 ± 0.201 0.321 ± 0.032 0.330 ± 0.033

Multiquadric

Modulated Scaled sigmoid atan2

0.474 ± 0.131 0.606 ± 0.020 0.590 ± 0.030

34.726 ± 16.454 22.552 ± 1.312 24.039 ± 2.186

30.373 ± 6.354 23.775 ± 0.599 24.360 ± 0.912

45.419 ± 13.377 32.939 ± 1.169 33.744 ± 1.437

0.406 ± 0.174 0.250 ± 0.014 0.256 ± 0.014

Sine power

Modulated Scaled sigmoid atan2

0.546 ± 0.093 0.600 ± 0.020 0.586 ± 0.026

27.677 ± 8.280 23.773 ± 1.836 24.890 ± 1.953

26.552 ± 4.656 23.554 ± 0.484 24.398 ± 0.783

37.190 ± 7.268 33.403 ± 1.190 34.016 ± 1.232

0.297 ± 0.090 0.249 ± 0.010 0.252 ± 0.011

Spherical

Modulated Scaled sigmoid atan2

0.347 ± 0.109 0.558 ± 0.054 0.457 ± 0.052

50.913 ± 20.542 25.493 ± 4.074 33.615 ± 4.288

35.597 ± 6.367 25.264 ± 1.407 28.206 ± 2.312

56.665 ± 17.111 35.825 ± 2.729 40.353 ± 3.806

0.544 ± 0.234 0.282 ± 0.028 0.320 ± 0.049

Askey

Modulated Scaled sigmoid atan2

0.358 ± 0.109 0.554 ± 0.057 0.462 ± 0.065

48.936 ± 19.576 25.680 ± 4.305 32.770 ± 4.848

35.737 ± 5.798 25.764 ± 1.339 27.815 ± 1.761

57.041 ± 15.566 36.398 ± 2.568 39.867 ± 2.714

0.553 ± 0.212 0.291 ± 0.025 0.315 ± 0.028

C 2 -Wendland

Modulated Scaled sigmoid atan2

0.547 ± 0.097 0.600 ± 0.018 0.588 ± 0.023

27.141 ± 7.615 23.331 ± 1.751 24.927 ± 2.029

26.047 ± 3.720 23.421 ± 0.486 24.225 ± 0.748

37.226 ± 7.359 32.985 ± 1.111 33.772 ± 1.224

0.298 ± 0.088 0.246 ± 0.010 0.249 ± 0.010

C 4 -Wendland

Modulated Scaled sigmoid atan2

0.541 ± 0.113 0.603 ± 0.022 0.595 ± 0.021

28.517 ± 11.805 22.989 ± 1.392 24.299 ± 1.739

27.149 ± 5.533 23.584 ± 0.536 24.246 ± 0.580

39.135 ± 10.713 32.906 ± 1.044 33.561 ± 1.109

0.328 ± 0.137 0.247 ± 0.011 0.250 ± 0.009

Powered exponential

Modulated Scaled sigmoid atan2

0.451 ± 0.120 0.636 ± 0.029 0.608 ± 0.026

35.858 ± 15.722 20.780 ± 1.893 22.268 ± 1.849

29.635 ± 6.077 22.108 ± 0.815 23.777 ± 0.860

43.414 ± 13.835 30.920 ± 1.167 32.897 ± 1.269

0.369 ± 0.182 0.229 ± 0.012 0.251 ± 0.014

Generalized Cauchy

Modulated Scaled sigmoid atan2

0.420 ± 0.113 0.607 ± 0.022 0.593 ± 0.027

39.099 ± 16.116 22.226 ± 1.513 23.061 ± 2.227

30.684 ± 6.301 22.970 ± 0.758 24.234 ± 0.683

45.639 ± 14.440 32.209 ± 0.931 33.579 ± 1.047

0.393 ± 0.193 0.241 ± 0.010 0.257 ± 0.011

Dagum

Modulated Scaled sigmoid atan2

0.301 ± 0.086 0.366 ± 0.031 0.372 ± 0.027

55.856 ± 19.890 39.528 ± 2.197 38.773 ± 2.181

38.068 ± 5.243 30.673 ± 0.996 29.552 ± 0.907

61.646 ± 17.354 45.573 ± 2.121 43.905 ± 1.697

0.608 ± 0.240 0.378 ± 0.029 0.353 ± 0.022

Multiquadric

Modulated Scaled sigmoid atan2

0.465 ± 0.129 0.672 ± 0.036 0.658 ± 0.029

34.633 ± 14.711 18.840 ± 2.170 19.119 ± 1.622

29.580 ± 6.322 20.698 ± 1.240 22.214 ± 0.927

42.983 ± 12.499 28.826 ± 2.034 30.367 ± 1.771

0.367 ± 0.160 0.208 ± 0.021 0.228 ± 0.020

Sine power

Modulated Scaled sigmoid atan2

0.554 ± 0.125 0.732 ± 0.043 0.665 ± 0.031

27.014 ± 10.983 16.233 ± 1.900 19.008 ± 1.857

25.174 ± 5.224 18.786 ± 1.198 21.604 ± 1.029

36.001 ± 9.907 25.509 ± 2.005 29.514 ± 1.688

0.284 ± 0.118 0.174 ± 0.019 0.217 ± 0.017

Spherical

Modulated Scaled sigmoid atan2

0.325 ± 0.088 0.484 ± 0.045 0.506 ± 0.048

50.511 ± 19.264 31.104 ± 3.022 29.524 ± 3.122

35.331 ± 5.642 26.245 ± 1.174 26.584 ± 1.333

56.105 ± 16.993 38.034 ± 2.030 37.273 ± 2.312

0.529 ± 0.235 0.292 ± 0.023 0.286 ± 0.028

Askey

Modulated Scaled sigmoid atan2

0.303 ± 0.081 0.415 ± 0.038 0.444 ± 0.036

54.639 ± 19.590 35.506 ± 2.220 33.584 ± 2.372

37.873 ± 5.146 27.932 ± 1.155 27.495 ± 0.652

59.828 ± 16.684 41.422 ± 1.923 39.308 ± 1.658

0.582 ± 0.232 0.323 ± 0.024 0.301 ± 0.017

C 2 -Wendland

Modulated Scaled sigmoid atan2

0.581 ± 0.114 0.753 ± 0.038 0.702 ± 0.036

24.504 ± 7.757 15.229 ± 1.468 17.435 ± 1.870

24.105 ± 4.472 18.180 ± 1.021 20.464 ± 1.215

33.562 ± 6.640 24.539 ± 1.711 27.494 ± 1.930

0.255 ± 0.068 0.167 ± 0.017 0.197 ± 0.020

C 4 -Wendland

Modulated Scaled sigmoid atan2

0.592 ± 0.101 0.739 ± 0.041 0.701 ± 0.033

23.940 ± 6.771 15.806 ± 1.795 17.189 ± 1.854

24.047 ± 4.219 18.460 ± 1.298 20.599 ± 1.154

34.032 ± 6.975 25.053 ± 2.038 27.715 ± 1.970

0.265 ± 0.078 0.171 ± 0.020 0.200 ± 0.020

Modulated Scaled sigmoid atan2

0.583 ± 0.059 0.602 ± 0.014 0.596 ± 0.019

25.100 ± 4.366 23.032 ± 1.255 24.012 ± 1.628

24.638 ± 2.660 22.949 ± 0.537 23.481 ± 0.472

34.387 ± 3.704 32.528 ± 0.769 33.138 ± 0.796

0.262 ± 0.043 0.242 ± 0.008 0.245 ± 0.006

Modulated Scaled sigmoid atan2

0.666 ± 0.090 0.781 ± 0.027 0.746 ± 0.029

19.418 ± 6.437 14.234 ± 1.284 14.984 ± 1.192

20.884 ± 4.080 16.948 ± 0.800 18.578 ± 0.952

28.774 ± 5.777 22.700 ± 1.259 24.813 ± 1.522

0.208 ± 0.064 0.149 ± 0.011 0.171 ± 0.016

Single-index

Geodesic

Dense

Single-index Chordal

None Dense

70

Wind Direction Distribution: India

Wind Direction Distribution: Germany

N (0/2 )

NW

N (0/2 )

NE

NW

W (3 /2)

NE

E W ( /2) (3 /2)

SW

SE

E ( /2)

SW

SE

S ()

S ()

Figure A.11: Rose plots indicating distribution of wind directions in [0, 2π) range for (A) Indian and (B) German datasets.

(A) Mean Wind Directions over India in December 2025

(B) Wind Directions over Germany 55.0° N

35.0° N 54.0° N 30.0° N

53.0° N

52.0° N

Latitude

Latitude

25.0° N

20.0° N

51.0° N

50.0° N 15.0° N 49.0° N 10.0° N

48.0° N Training Data Testing Data 70.0° E

75.0° E

80.0° E

85.0° E

Longitude

90.0° E

95.0° E

47.0° N

Training Data Testing Data 6.0° E

8.0° E

10.0° E

Longitude

12.0° E

14.0° E

Figure A.12: Training and test observations for (A) Indian and (B) German wind datasets. The train-test split is done randomly for both the datasets. The boundaries are for illustrative purposes only, and imply no political assertions.

71

India Wind Data Split for Evaluation of Extrapolability

35.0° N

30.0° N

Latitude

25.0° N

20.0° N

15.0° N

10.0° N

Training Data Testing Data 70.0° E

75.0° E

80.0° E

85.0° E

Longitude

90.0° E

95.0° E

Figure A.13: Training (70.1%) and test (29.9%) points of the Indian wind dataset for evaluation of extrapolability. The boundaries imply no political assertions and are for illustrative purposes only.

72

References Alonso-Pena, M. & Crujeiras, R. M. (2023), ‘Analyzing animal escape data with circular nonparametric multimodal regression’, The Annals of Applied Statistics 17(1), 130–152. Alonso-Pena, M. & Crujeiras, R. M. (2025a), A Review on Nonparametric Circular Regression, Springer Nature Singapore, p. 159–176. Alonso-Pena, M. & Crujeiras, R. M. (2025b), ‘A review on nonparametric circular regression’, Directional and Multivariate Statistics: A Volume in Honour of Ashis SenGupta pp. 159–176. Ameijeiras-Alonso, J. & Gijbels, I. (2026), ‘Semiparametric regression for circular response with application in ecology’, Scandinavian Journal of Statistics 53(1), 54–101. Ameijeiras-Alonso, J., Gijbels, I. & Verhasselt, A. (2022), ‘On a family of two–piece circular distributions’, Computational Statistics & Data Analysis 168, 107403. Anichini, M., De Reus, K., Hersh, T. A., Valente, D., Salazar-Casals, A., Berry, C., Keller, P. E. & Ravignani, A. (2023), ‘Measuring rhythms of vocal interactions: a proof of principle in harbour seal pups’, Philosophical Transactions of the Royal Society B: Biological Sciences 378(1875). Artemiou, A. & Li, B. (2009), ‘On principal components and regression: a statistical explanation of a natural phenomenon’, Statistica Sinica pp. 1557–1565. Bell, W. & Nadarajah, S. (2024), ‘A review of wrapped distributions for circular data’, Mathematics 12(16), 2440. Bercu, B., Bigot, J. & Thurin, G. (2024), ‘Regularized estimation of monge-kantorovich quantiles for spherical data’, arXiv preprint arXiv:2407.02085 . Bhuyan, P., Halder, S. & Jha, J. (2026), ‘Modeling Zero-Inflated Longitudinal Circular Data Using Bayesian Methods: Application to Ophthalmology’, arXiv preprint arXiv:2601.13998 . Biswas, S. & Banerjee, B. (2025), ‘Semi-parametric least-area linear-circular regression through möbius transformation’, Statistical Analysis and Data Mining: An ASA Data Science Journal 18(1), e70012. Bruns, L., Lamparter, L., Galic, M. & Jiang, X. (2024), ‘Single image estimation of cell migration direction by deep circular regression’, arXiv preprint arXiv:2406.19162 . Chen, L.-P. (2021), ‘Ultrahigh-dimensional sufficient dimension reduction with measurement error in covariates’, Statistics & Probability Letters 168, 108931. Cohen, Y., Navarro, A., Frellsen, J., Turner, R. E., Riemer, R. & Pakman, A. (2025), Bayesian Circular Regression with von Mises Quasi-Processes, in ‘28th International Conference on Artificial Intelligence and Statistics’, ML Research Press, pp. 1693–1701. Demni, H., Hubert, M., Porzio, G. C. & Rousseeuw, P. J. (2026), ‘Robust measures of dispersion for circular data with an anomaly detection rule’, arXiv preprint arXiv:2603.01237 . Deng, J., Dong, W., Socher, R., Li, L.-J., Li, K. & Fei-Fei, L. (2009), ImageNet: A large-scale hierarchical image database, in ‘2009 IEEE conference on computer vision and pattern recognition’, IEEE, pp. 248–255. 73

Di Marzio, M., Fensore, S., Panzera, A. & Taylor, C. C. (2016), ‘A note on nonparametric estimation of circular conditional densities’, Journal of Statistical Computation and Simulation 86(13), 2573– 2582. Di Marzio, M., Panzera, A. & Taylor, C. C. (2012), ‘Smooth estimation of circular cumulative distribution functions and quantiles’, Journal of Nonparametric Statistics 24(4), 935–949. Di Marzio, M., Panzera, A. & Taylor, C. C. (2013), ‘Non-parametric regression for circular responses’, Scandinavian Journal of Statistics 40(2), 238–255. Di Marzio, M., Panzera, A. & Taylor, C. C. (2014), ‘Nonparametric regression for spherical data’, Journal of the American Statistical Association 109(506), 748–763. Di Marzio, M., Panzera, A. & Taylor, C. C. (2016), ‘Nonparametric circular quantile regression’, Journal of Statistical Planning and Inference 170, 1–14. Downs, T. D. & Mardia, K. (2002), ‘Circular regression’, Biometrika 89(3), 683–698. Fan, J., Fang, C., Gu, Y. & Zhang, T. (2024), ‘Environment invariant linear least squares’, The Annals of Statistics 52(5), 2268–2292. Fernández-Durán, J. J. & Gregorio-Domínguez, M. M. (2024), ‘Regression models for circular data based on nonnegative trigonometric sums’, Journal of Statistical Planning and Inference 230, 106114. Fisher, N. I. & Lee, A. J. (1992), ‘Regression models for an angular response’, Biometrics pp. 665– 677. Fornasin, A., Rizzi, L. & Fonseca, G. (2025), ‘The effects of meteorological conditions on the circadian rhythm of births’, American Journal of Human Biology 37(8), e70120. Francisco-Fernández, M. & Meilán-Vila, A. (2026), ‘Analyzing directional errors in spatial orientation using nonparametric circular regression with mixed covariates’, arXiv preprint arXiv:2604.21498 . Gao, F., Chia, K.-S., Krantz, I., Nordin, P. & Machin, D. (2006), ‘On the application of the von mises distribution and angular regression methods to investigate the seasonality of disease onset’, Statistics in medicine 25(9), 1593–1618. Gao, H., Huang, L., Yang, Y. & Shao, X. (2026), ‘Testing equality of conditional distributions via generative models’, arXiv preprint arXiv:2606.06930 . Gill, J. & Hangartner, D. (2010), ‘Circular data in political science and how to handle it’, Political Analysis 18(3), 316–336. Gneiting, T. (2013), ‘Strictly and non-strictly positive definite functions on spheres’, Bernoulli 19(4), 1327–1349. Gneiting, T. & Raftery, A. E. (2007), ‘Strictly proper scoring rules, prediction, and estimation’, Journal of the American statistical Association 102(477), 359–378. Gorczyca, M. T. & Sefas, J. D. (2025), ‘Weighted trigonometric regression for suboptimal designs in circadian transcriptome studies’, Statistics in Medicine 44(20-22), e70201. 74

Gottard, A., Meilán-Vila, A. & Panzera, A. (2026), ‘Quasi-likelihood estimation for semiparametric circular regression models’, Biometrics 82(1). Grimit, E. P., Gneiting, T., Berrocal, V. J. & Johnson, N. A. (2006), ‘The continuous ranked probability score for circular variables and its application to mesoscale forecast ensemble verification’, Quarterly Journal of the Royal Meteorological Society: A journal of the atmospheric sciences, applied meteorology and physical oceanography 132(621C), 2925–2942. Hallin, M., Liu, H. & Verdebout, T. (2024), ‘Nonparametric measure-transportation-based methods for directional data’, Journal of the Royal Statistical Society Series B: Statistical Methodology 86(5), 1172–1196. Harvey, A., Hurn, S., Palumbo, D. & Thiele, S. (2024), ‘Modelling circular time series’, Journal of Econometrics 239(1), 105450. Hassanzadeh, F. (2021), ‘A smoothing spline model for multimodal and skewed circular responses: Applications in meteorology and oceanography’, Environmetrics 32(2), e2655. Hauch, K. & Redenbach, C. (2024), ‘Quantiles and depth for directional data from elliptically and skew-rotationally symmetric distributions’, Electronic Journal of Statistics 18(1), 2042–2063. Hernandez-Stumpfhauser, D., Breidt, F. J. & van der Woerd, M. J. (2017), ‘The general projected normal distribution of arbitrary dimension: Modeling and bayesian inference’. Hogan, J. A., McMahon, S. M., Buzzard, V., Michaletz, S. T., Enquist, B. J., Thompson, J., Swenson, N. G. & Zimmerman, J. K. (2019), ‘Drought and the interannual variability of stem growth in an aseasonal, everwet forest’, Biotropica 51(2), 139–154. Huang, H.-H., Yu, F., Li, K. & Zhang, T. (2026), ‘Fréchet sufficient dimension reduction for metric space-valued data via distance covariance’, Journal of Computational and Graphical Statistics 0(0), 1–14. Jha, J. & Biswas, A. (2017), ‘Multiple circular–circular regression’, Statistical Modelling 17(3), 142– 171. Jha, J. & Biswas, A. (2018), ‘Circular-circular regression model with a spike at zero’, Statistics in Medicine 37(1), 71–81. Jha, J. & Biswas, A. (2022), Regression models for directional variables, in ‘Directional Statistics for Innovative Applications: A Bicentennial Tribute to Florence Nightingale’, Springer, pp. 333–348. Józiak, P. (2015), ‘Conditionally strictly negative definite kernels’, Linear and Multilinear Algebra 63(12), 2406–2418. Kato, S. & Jones, M. (2010), ‘A family of distributions on the circle with links to, and applications arising from, möbius transformation’, Journal of the American Statistical Association 105(489), 249–262. Kato, S., Shimizu, K. & Shieh, G. S. (2008), ‘A circular–circular regression model’, Statistica Sinica pp. 633–645.

75

Khattak, A., Chan, P.-w., Chen, F. & Peng, H. (2023), ‘Assessing wind field characteristics along the airport runway glide slope: An explainable boosting machine-assisted wind tunnel study’, Scientific Reports 13(1), 10939. Kim, K., Song, J. & Yoo, J. K. (2026), ‘Valid asymptotic inference after sufficient dimension reduction in a single-index framework’, Statistics and Computing 36(1), 24. Kim, S. & SenGupta, A. (2016), Regressions involving circular variables: an overview, in ‘Platinum Jubilee International Conference on Applications of Statistics’, Springer, pp. 25–33. Klar, B., Milošević, B. & Obradović, M. (2024), ‘Nonparametric estimation on the circle based on fej\’er polynomials’, arXiv preprint arXiv:2411.19138 . Koning, A. J., Franses, P. H., Hibon, M. & Stekler, H. O. (2005), ‘The M3 competition: Statistical tests of the results’, International Journal of Forecasting 21(3), 397–409. Kraft, B., Stalder, S., Aeberhard, W. H., Ruiz, N. H., Meinshausen, N., Shen, X. & Gudmundsson, L. (2026), ‘Modeling uncertainty with engression: A deep generative time-series approach’, Geophysical Research Letters 53(2), e2025GL120122. Krantz, S. G. & Parks, H. R. (2002), A primer of real analytic functions, Springer Science & Business Media. Laha, A. K. & Majumdar, S. (2022), Angular-angular and linear-angular regression using ann, in ‘Directional Statistics for Innovative Applications: A Bicentennial Tribute to Florence Nightingale’, Springer, pp. 451–474. Landler, L., Ruxton, G. D. & Malkemper, E. P. (2021), ‘Advice on comparing two independent samples of circular data in biology’, Scientific reports 11(1), 20337. Lang, M. N., Schlosser, L., Hothorn, T., Mayr, G. J., Stauffer, R. & Zeileis, A. (2020), ‘Circular regression trees and forests with an application to probabilistic wind direction forecasting’, Journal of the Royal Statistical Society Series C: Applied Statistics 69(5), 1357–1374. Lee, K.-Y., Li, B. & Chiaromonte, F. (2013), ‘A general theory for nonlinear sufficient dimension reduction: Formulation and estimation’, The Annals of Statistics pp. 221–249. Ley, C., Sabbah, C. & Verdebout, T. (2014), ‘A new concept of quantiles for directional data and the angular mahalanobis depth’. Lin, W. & Kulasekera, K. (2007), ‘Identifiability of single-index models and additive-index models’, Biometrika 94(2), 496–501. Liu, J., Sun, W., Yang, H., Zeng, Z., Liu, C., Zheng, J., Liu, X., Rahmani, H., Sebe, N. & Mian, A. (2024), ‘Deep learning-based object pose estimation: A comprehensive survey’, arXiv preprint arXiv:2405.07801 . Liu, Z., Mao, H., Wu, C.-Y., Feichtenhofer, C., Darrell, T. & Xie, S. (2022), A ConvNet for the 2020s, in ‘Proceedings of the IEEE/CVF conference on computer vision and pattern recognition’, pp. 11976–11986.

76

Lund, U. J. (2002), ‘Tree-based regression for a circular response’, Communications in StatisticsTheory and Methods 31(9), 1549–1560. McMillan, G. P., Hanson, T. E., Saunders, G. & Gallun, F. J. (2013), ‘A two-component circular regression model for repeated measures auditory localization data’, Journal of the Royal Statistical Society Series C: Applied Statistics 62(4), 515–534. Meilan-Vila, A., Crujeiras, R. M. & Francisco-Fernández, M. (2024), ‘Nonparametric estimation for a functional-circular regression model’, Statistical Papers 65(2), 945–974. Mohammad, H. H., Satari, S. Z. & Yusoff, W. N. S. W. (2021), Review on circular-linear regression models, in ‘Journal of Physics: Conference Series’, Vol. 1988, IOP Publishing, p. 012108. Nguyen, T. D., Pham Ngoc, T. M. & Rivoirard, V. (2023), ‘Adaptive warped kernel estimation for nonparametric regression with circular responses’, Electronic Journal of Statistics 17(2), 4011– 4048. Ning, C., Nie, S., Li, C., Yao, Z. & Yuan, G. (2025), ‘Hybrid-driven meta-learning for wind field prediction in unmanned sailboat applications’, Ocean Engineering 341, 122771. Nonaka, S., Nobuhara, S. & Nishino, K. (2022), Dynamic 3D Gaze From Afar: Deep Gaze Estimation From Temporal Eye-Head-Body Coordination, in ‘Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition’, pp. 2192–2201. Nur Ibrahim, A. I. & Xu, G. (2024), Bayesian approach to circular-linear regression with application, in ‘AIP Conference Proceedings’, Vol. 3189, AIP Publishing LLC, p. 100010. Pathak, R. & Chakraborty, T. (2026), ‘Deep Generative Spatiotemporal Engression for Probabilistic Forecasting of Epidemics’, Transactions on Machine Learning Research . URL: https://openreview.net/forum?id=7AfAztCd5A Pérez-González, A., Cotos-Yáñez, T. R. & Crujeiras, R. M. (2026), ‘Nonparametric modal regression with missing response observations’, Computational Statistics 41(3), 58. Presnell, B., Morrison, S. P. & Littell, R. C. (1998), ‘Projected multivariate linear models for directional data’, Journal of the American Statistical Association 93(443), 1068–1077. Prokudin, S., Gehler, P. & Nowozin, S. (2018), Deep directional statistics: Pose estimation with uncertainty quantification, in ‘Proceedings of the European conference on computer vision (ECCV)’, pp. 534–551. Qiu, H., Tchetgen, E. T. & Dobriban, E. (2024), ‘Efficient and multiply robust risk estimation under general forms of dataset shift’, Annals of statistics 52(4), 1796. Quill, R., Sharples, J. J., Wagenbrenner, N. S., Sidhu, L. A. & Forthofer, J. M. (2019), ‘Modeling wind direction distributions using a diagnostic model in the context of probabilistic fire spread prediction’, Frontiers in Mechanical Engineering 5, 5. Rad, N. N., Bekker, A. & Arashi, M. (2022), ‘Enhancing wind direction prediction of south africa wind energy hotspots with bayesian mixture modeling’, Scientific Reports 12(1), 11442.

77

Ranalli, M. & Maruotti, A. (2020), ‘Model-based clustering for noisy longitudinal circular data, with application to animal movement’, Environmetrics 31(2), e2572. Rasib, M., Butt, M. A., Riaz, F., Sulaiman, A. & Akram, M. (2021), ‘Pixel level segmentation based drivable road region detection and steering angle estimation method for autonomous driving on unstructured roads’, IEEE Access 9, 167855–167867. Ren, Y., Jin, C., Xia, Y., Ke, L., Huang, L., Xue, H., Zhang, H., Guan, J. & Zhou, S. (2025), Score-based generative modeling for conditional independence testing, in ‘Proceedings of the 31st ACM SIGKDD Conference on Knowledge Discovery and Data Mining V. 2’, pp. 2410–2419. Shen, X. & Meinshausen, N. (2025), ‘Engression: extrapolation through the lens of distributional regression’, Journal of the Royal Statistical Society Series B: Statistical Methodology 87(3), 653– 677. Sjöberg, S. & Nilsson, C. (2015), ‘Nocturnal migratory songbirds adjust their travelling direction aloft: evidence from a radiotelemetry and radar study’, Biology letters 11(6). Steinwart, I. & Christmann, A. (2011), ‘Estimating conditional quantiles with the help of the pinball loss’. Steinwart, I. & Ziegel, J. F. (2021), ‘Strictly proper kernel scores and characteristic kernels on compact spaces’, Applied and Computational Harmonic Analysis 51, 510–542. Szegedy, C., Vanhoucke, V., Ioffe, S., Shlens, J. & Wojna, Z. (2016), Rethinking the inception architecture for computer vision, in ‘Proceedings of the IEEE conference on computer vision and pattern recognition’, pp. 2818–2826. Taylor, C. (2012), ‘Non-parametric Regression for Circular Responses’, Scandinavian Journal of Statistics . Tulsiani, S. & Malik, J. (2015), Viewpoints and keypoints, in ‘Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition’, pp. 1510–1519. Villarini, G. (2016), ‘On the seasonality of flooding across the continental united states’, Advances in Water Resources 87, 80–91. Waghmare, K. & Ziegel, J. (2025), ‘Proper scoring rules for estimation and forecast evaluation’, Annual Review of Statistics and Its Application 13. Wang, F. & Gelfand, A. E. (2014), ‘Modeling space and space-time directional data using projected gaussian processes’, Journal of the American Statistical Association 109(508), 1565–1580. Woolsey, N. & Huang, X. (2025), ‘Nonparametric regression for a circular response with error-incovariate’, Electronic Journal of Statistics 19(2), 4575–4636. Wu Navarro, A. K. (2019), Probabilistic Machine Learning for Circular Statistics: Models and inference using the Multivariate Generalised von Mises distribution, PhD thesis, University of Cambridge. Xiang, Y., Mottaghi, R. & Savarese, S. (2014), Beyond pascal: A benchmark for 3d object detection in the wild, in ‘IEEE winter conference on applications of computer vision’, IEEE, pp. 75–82. 78

Ye, X., Van Niekerk, J. & Rue, H. (2026), ‘A bayesian regression framework for circular models with inla’, arXiv preprint arXiv:2602.08413 . Zhang, Y., Huang, L., Yang, Y. & Shao, X. (2026), ‘Doubly robust conditional independence testing with generative neural networks’, Journal of the Royal Statistical Society Series B: Statistical Methodology 88(1), 195–220. Zheng, S., Wang, T., Lan, M. & Lin, Y. (2025), ‘A conditional distribution equality testing framework using deep generative learning’, arXiv preprint arXiv:2509.17729 . Zou, T. & Stern, H. S. (2025), ‘A dirichlet process model for directional-linear data with application to bloodstain pattern analysis’, Computational Statistics & Data Analysis 204, 108093.

79