The projected normal distribution, with isotropic variance, on the 2-sphere is considered using intrinsic statistics. It is shown that in this case, the expectation commutes with the projection, and that the covariance of the normal variable has a 1-1 correspondence with the intrinsic covariance of the projected normal distribution. This allows us to estimate, after the model identification, the parameters of the underlying normal distribution that generates the data.
Spherical statisticsprojected normalparameter estimation62H1162F10This project has received partial funding from Huawei Technologies.Introduction
Directional statistics (aka. spherical statistics) form a relevant subfield of statistics in which one studies directions, rotations, and axes. A typical situation in spherical statistics is that observations are gathered on a sphere, say S2, and consequently the methods have to be adapted to non-Euclidean geometry. More generally, one can think of observations on more general compact Riemannian manifolds. Application areas are numerous as, for example, one can think of S2 representing the Earth surface, and measurements are then observations on this surface. To name just a few application areas, see [11, 22] for navigation/control in robotics, [14, 16] for modeling wind directional changes, [6] for finding lymphoblastic leukemia cells, [3] for movement of tectonic plates, [12] for modelling of protein chains, and [9, 19, 21] for radiology applications in the context of MIMO-systems. For details on spherical statistics, we refer to the monograph [5].
One of the central problems in statistics is to estimate model parameters from observations. In the spherical context, this can mean, for example, that one assumes a parametric distribution Pθ on S2 and then uses the data to estimate the unknown θ. The most commonly applied distributions on the sphere are von Mises–Fisher distributions (also called von Mises distributions in S or Kent distributions in S2) that can be viewed as the equivalent of the normal distribution on the sphere. For the parameter estimation for von Mises–Fisher distributions, see [23]. Another widely applied distribution on Sd is the projected normal distribution. That is, the distribution of X/‖X‖ for X∼N(μ,Σ). While the density function of the projected normal is well known (see [7]), the parameter estimation is much less studied due to the complicated nature of the density. In particular, the parameter estimation is studied in the circular case (S), see [15, 16, 20, 24].
In this article, we consider the problem of parameter estimation related to the projected normal distribution onto S2. In contrast to the existing literature, we do not consider estimation of the parameters of the projected normal that can be obtained via standard methodology, as the density is completely known. Instead, our aim is to extract information on the underlying normal distribution N(μ,Σ) (with support on the ambient Euclidean space R3) based solely on data consisting of projected points on S2. Obviously, we immediately run into identifiability issues if we observe only X/‖X‖ instead of X∼N(μ,Σ). First of all, it is clear that one cannot identify arbitrary shapes of Σ (for example, the distribution can be arbitrarily spread in the direction μ and this cannot be observed as all points in this direction are equally projected). For this reason, we assume, for the sake of simplicity, isotropic variance Σ=σ2I3 and assume X∼N(μ,σ2I3). However, even in this case, we can only estimate the quantities μ/‖μ‖ (the direction of the location) and σ2/‖μ‖2, the reason being that the distribution of the projection pr(aX) is invariant on a>0. This is indeed natural, as intuitively one can only estimate the direction μ/‖μ‖ from the projections (onto S2). Similarly, N(μ,σ2I3) for larger σ located in a distant μ seems similar to the normal distribution with smaller variance but located closer, if one is only observing both distributions on the surface of a sphere. We also note that estimation of the direction μ/‖μ‖ is already well known, and we claim no originality in this respect. Instead, our main contribution is the estimation of λ=σ2/‖μ‖2. For this, we study the covariance matrix of X/‖X‖ on S2 and by linking it to certain special functions and analyzing their series expansions, we show that there is a bijective mapping between λ and the covariance matrix of X/‖X‖ on S2. As the latter can be estimated by using the methods of spherical statistics, we obtain a consistent estimator for λ via the inverse mapping that can be computed, e.g., via the bisection method.
The rest of the article is organized as follows. In Section 2 we present and discuss our main results. We begin with Section 2.1 where we introduce our setup and notation. After that, we discuss the convergence of sample estimators (for mean and covariance) in the context of a general Riemannian manifold. The case of the projected normal in S2 is then discussed in Section 2.3. All the proofs are postponed to Section 3.
Main result
In this section we present and discuss our main results. First, in Section 2.1, we shall introduce some notation and clarify some terminology used in the context of manifold-valued random variables. Some of the theory relies on differential geometric concepts which are briefly summarized in Appendix A. The convergence of sample covariances on a general compact manifold is discussed in Section 2.2, while the special case of S2 and the projected normal distribution is treated in Section 2.3.
General setting
Let M⊂Rk be a smooth n-dimensional manifold. Let S∈Rk be a dense subset and let pr:S→M be the projection onto M⊆Rk, assumed smooth everywhere on S and hence almost everywhere in Rk. Let (Ω,P) be a probability space and consider a normally distributed (multivariate) random variable
X:Ω→Rk
such that
P(X∉S)=0.
Let x1,…,xL be an independently drawn sample of X. Suppose now that we observe
pr(x1),…,pr(xL).
The question is now, can one estimate the mean and covariance of this projected sample?
In order for this question to have meaning, a notion of mean and covariance is needed for manifold-valued random variables. The intrinsic mean, a.k.a. the Frechét mean (see [17]), of an absolutely continuous random variable X:Ω→M is defined as follows.
Let M be a Riemannian manifold with corresponding distance function dist and volume form dVolM. Moreover, suppose pX:M→[0,∞) is a probability density function of some absolutely continuous random variable X:Ω→M. The expected value of X is then defined by
arginfq∈M∫Mdist(q,y)2pX(y)dVolM(y)=E[X].
The argument in the infimum in (2.1) need not exist nor does it need to be unique. For example, consider X=dN(0,σ2I)∈Rn+1, then pr(X) follows the uniform distribution on Sn and all points in Sn satisfy the infimum in (2.1). This is a very important difference from the case of Rn-valued random variables. In Rn we know that if X is absolutely continuous, integrable and square integrable, then any point μ∈Rn which minimizes the least square integral
∫Rnx−μ2pXdx
is unique.
The intrinsic definition of the covariance matrix for a random variable on a Riemannian manifold is well known as well (see, e.g., [17]).
Let M be a geodesically complete Riemannian manifold, and let logq:M→TqM be the natural log map (defined a.e.). Let X be an M-valued absolutely continuous random variable with the intrinsic mean μ=E[X]∈M. Then, the covariance matrix for X is the linear map Cov(X):TμM→TμM defined by the integral
Cov(X)=∫Mlogμ(y)logμ(y)TpX(y)dVolM(y).
Convergence of sample estimators on a compact manifold
Let yℓ be a sample of L independent measurements on a compact manifold M. In order to estimate the mean of such a sample, we shall utilize the discrete version of Definition 2.1. This definition follows that of [17].
Let {yℓ}ℓ=1L be a sample of L points on a Riemannian manifold M. Then the empirical mean is defined as
arginfq∈M∑ℓ=1Ldist2(yℓ,q).
Note that since the function dist2 has good regularity, one may hope to find such an infimum by solving the (nonlinear) equation
q:∑ℓ=1Llogq(yℓ)=0.
It was shown in [8] that if M is compact and if yℓ is sampled from a distribution with a unique empirical mean x0, then the distribution of the above follows a central limit theorem. More precisely, if μ=E[yℓ] is the true mean of the underlying distribution and if μˆ is the empirical mean, it holds that
Llogμμˆ⟶dN(0,V)
for some linear map V∈L(TμM,TμM).
Let {yℓ}ℓ=1L be a sample of L points on a Riemannian manifold M with a unique empirical mean ξˆ. Then the empirical covariance of the sample is defined by
Vˆ=1L−1∑ℓ=1Llogξˆ(yℓ)logξˆ(yℓ)T
where Vˆ:TξˆM→TξˆM is a linear map.
Since the empirical mean converges by the results of [8], it remains to verify that the empirical covariance in Equation (2.2) converges to the covariance V of the limiting distribution limL→∞Llogξξˆ.
The following result shows that, in the case of isotropic covariance, the empirical covariance converges. The proof is postponed to Section 3.1.
Let{ξℓ}ℓ=1Lbe L independent identically distributed random variables on a compact geodesically complete manifold M. Suppose further they have a unique meanE[ξℓ]=μand an isotropic covarianceCov(ξℓ)=vI. Then,1L−1∑ℓ=1Llogμˆ(ξℓ)logμˆ(ξℓ)T⟶PvIwith the same rate of convergence as the empirical meanμˆof the sample{ξℓ}ℓ=1L.
As far as we know, the rate of convergence for the empirical mean on compact manifolds is not completely solved as of now. In [8] it has been shown that the empirical mean has a rate of convergence L for a large class of manifolds, but not all compact geodesically complete manifolds. However, [2] provides rates of convergence for the empirical mean for a general class of metric spaces, from which it seems that the rate of convergence L isn’t necessarily true.
Following a similar method as in [18] the update scheme for a sample of observations is given in Algorithm 1.
Estimating intrinsic average and covariance for a sample yℓ on a manifold M
Observing the projected normal in <inline-formula id="j_vmsta279_ineq_077"><alternatives><mml:math>
<mml:msup>
<mml:mrow>
<mml:mi mathvariant="double-struck">S</mml:mi>
</mml:mrow>
<mml:mrow>
<mml:mn>2</mml:mn>
</mml:mrow>
</mml:msup></mml:math><tex-math><![CDATA[${\mathbb{S}^{2}}$]]></tex-math></alternatives></inline-formula>
In this section we consider M=S2, the unit 2-sphere in R3. In this case, the projection map is simply
pr(x)=xx
and the domain of definition for pr is S=R3∖{0}.
Throughout, we shall use spherical coordinates, i.e. for the classical Cartesian coordinates (x,y,z)T in R3 we set x=cos(θ)sin(ϕ), y=sin(θ)sin(ϕ), and z=cos(ϕ). Consider two points (θ1,ϕ1)T and (θ2,ϕ2)T. Then it is a classical result (see, e.g., [10]) that the distance function may be written as
dist((θ1,ϕ1)T,(θ2,ϕ2)T)=acoscos(ϕ1)cos(ϕ2)+sin(ϕ1)sin(ϕ2)cos(θ2−θ1).
It has been shown in [7] that the projected normal, i.e. the random variable defined by pr(X) where X=dN(μ,Σ), has the probability density function
ppr(X)(θ,ϕ)=12πA3/2Σ−12exp(C)K+K2Φ(K)φ(K)+Φ(K)φ(K)
at a point u=(cos(θ)sin(ϕ),sin(θ)sin(ϕ),cos(ϕ))T, where A=uTΣ−1u, B=uTΣ−1μ, C=−12μTΣ−1μ, and K=BA−12. Moreover, the functions φ,Φ:R→R are defined as
φ(x)=exp(−x22)
and
Φ(a)=∫−∞aφ(x)dx,
respectively.
In our context, this random variable is then observed in S2. As S2 is an isometrically embedded Riemannian manifold, Definition 2.1 applies for the projected normal. Intuitively, the expectation of a projected normal distribution ought to be pr(μ). Unfortunately, for the fully general case of the projected normal, this conjecture will be false. However, by imposing an isotropic covariance Σ onto X we shall show that this intuition is indeed correct. This is the topic of the next result whose proof is postponed in Section 3.2.
Let X be a normally distributed random variable with averageμ∈R3and covariance matrix Σ, i.e.X=dN(μ,Σ)∈R3. IfΣ=σ2I3, thenE[pr(X)]=pr(μ).
The above theorem is true whenever μ is an eigenvector of Σ. However, if μ is not an eigenvector of Σ, then the statement of Theorem 2.6 is false in general. To see this, let μ=(0,0,1)T and
Σ=100010.500.51.
In this case it follows that the probability density function is not symmetric around μ, see Figure 1 below.
A normally distributed random variable X=dN(μ,Σ), with the covariance matrix Σ such that the mean μ is not an eigenvector for Σ, for which E[pr(X)]≠pr(μ)
Geometric picture showing a sphere S² and a normal random variable X with covariance matrix Σ, and mean μ. It illustrates that the projected normal random variable pr(X) might not have mean pr(μ) for a general covariance matrix Σ.
In general, the tangent space of S2 at a point q is the set of vectors orthogonal (in the Euclidean product inherited from R3) to q. Note that the initial velocity of the geodesic connecting a point μ to q has the same direction as q×μ (i.e. the cross-product of R3). The logarithm map of q centered at a point μ is thus the vector of the normalised starting velocity of the geodesic connecting μ with q times the length of the geodesic connecting μ and q, see [13]. Without loss of generality, consider μ=(0,0,1)T. Then
logμ(q)=cos(θ)sin(ϕ)sin(θ)sin(ϕ)ϕsin(ϕ),
where the basis for the vector is in a rewritten form of the basis (of TμS2)
100,010.
Hence the intrinsic covariance of pr(X) can be written as
Cov(pr(X))=∫02π∫0πcos2(θ)sin2(ϕ)cos(θ)sin(θ)sin2(ϕ)cos(θ)sin(θ)sin2(ϕ)sin2(θ)sin2(ϕ)×ϕ2sin2(ϕ)12π3/2exp(−12σ2)×K+K2Φ(K)φ(K)+Φ(K)φ(K)sin(ϕ)dϕdθ,
where K=1σcos(ϕ). Or, equivalently, by simplifying and integrating w.r.t. θ,
Cov(pr(X))=π(2π)3/2exp(−12σ2)1001×∫0πϕ2cos(ϕ)σ+cos2(ϕ)σ2+1Φ(cos(ϕ)σ)φ(cos(ϕ)σ)sin(ϕ)dϕ=:1001f(σ).
In the limit as σ→∞, it holds that ppr(X)→14π𝟙S2. The intrinsic covariance in this limit is
π2−441001.
On the other hand, in the limit as σ→0, ppr(X) converges to the point mass distribution while the covariance goes to zero.
Note that the above estimates are only for the intrinsic expectation and covariance. If we want to relate these estimates to the extrinsic expectation and covariance, we end up with the obvious problem that pr(X) and pr(aX), a>0, have the same distribution. Hence we need to choose which normal distribution in R3 corresponds to the intrinsic covariance and average. By Lemma 2.6 the average is known (up to a factor). Moreover, it turns out there is a one-to-one correspondence between the extrinsic covariance and the intrinsic covariance up to a factor a>0, and thus one may estimate the underlying covariance parameter σ using the relation (2.4) with, e.g., the bisection method.
The degeneracy of the factor a essentially means that we can only estimate the parameter σ in the case when the projected average is the true average of the R3-valued normal random variable. More generally, if X=dN(μ,σ2I3), then X‖μ‖=dNμ‖μ‖,σ2‖μ‖2I3. Consequently, we can only estimate the quantities μ‖μ‖∈S2 and σ2‖μ‖2. This means that, as expected, we can estimate the direction μ‖μ‖ of the extrinsic distribution X=dN(μ,σ2I3) but not the distance ‖μ‖. On the other hand, for the variance we can only estimate the quantity σ2‖μ‖2. This is expected as well, as for the distribution N(μ,σ2I3) located far away (‖μ‖ large) the projections onto S2 are not wide spread even if σ>0 is large.
The required one-to-one correspondence between extrinsic and intrinsic covariances is formulated in the following proposition, whose proof is presented in Section 3.3. The statement is also illustrated in Figure 2.
Let X be a normally distributed random variable with meanμ∈R3and isotropic varianceσ2I3, i.e.Xμ=dN(μμ,σ2μ2I3). ThenCov(pr(X))=vI2for0≤v<π2−44and the relationfσ2μ2:=tr(Cov(pr(X)))2=vis a bijection.
The bijection f is crucial for observing projected normal random variables. Utilizing this f, if the scalar variance of some normal random variable X=dN(μ,σ2I3) in R3 is known, then the intrinsic scalar variance of the corresponding projected normal random variable is precisely f(σ2μ2). Conversely, if some projected isotropic normal random variable has mean μ and covariance vI2, then the normal random variable X such that pr(X) has these parameters is precisely X=dN(μ,f−1(v)I3) (up to a positive scalar factor).
As a consequence of Proposition 2.7 and Theorem 2.5 we can conclude that given measurements on the sphere, we can estimate the scalar variance of an isotropically distributed normal random vector in R3. This leads to the next result that can be viewed as the main theorem of the present paper. Its proof follows directly from Theorem 2.5 and Proposition 2.7. The rate of convergence is obtained immediately from noting that [8, Theorem 2] applies to S2, and thus the empirical mean has rate of convergence L, and by Theorem 2.5 so does the empirical covariance.
Let X be a normally distributed random variable with meanμ∈R3and isotropic varianceσ2I3, i.e.X=dN(μ,σ2I3). Given independent measurements(x1,x2,…,xL)frompr(X), we can estimateλ=σ2μ2byλˆ=f−1tr(Vˆ)2whereVˆis the empirical covariance matrix given in Equation (2.2) and where f is the bijection given in Proposition2.7and defined in Equation (2.4). Moreover, it holds thatλˆ⟶Pσ2μ2asL→∞, with rate of convergenceL.
A plot of the scalar variance of pr(X), i.e. f(σ2) from Proposition 2.7, where X=dN(μ,σ2I3) and μ∈S2 is arbitrary. The red line indicates the upper bound (π2−4)/4=limσ→∞tr(Cov(pr(X)))2
Graph showing the function f(σ²)=tr(Cov(pr(X)))/2 as a blue curve and the upper limit (π²-4)/4 as a red horizontal line, with σ² on the x-axis.
We conclude this section with some simulations. Firstly, Table 1 illustrates an empirical verification of Theorem 2.5 for the case of the manifold S2. The simulations are conducted by generating L samples from X=dN(μ,σ2I3) with μ(0,0,1)T∈S2 and σ=1 and projecting these samples onto S2. The sample covariance is computed as in Algorithm 1, and its error is then compared to the theoretical covariance using the Frobenius norm. This error is then averaged over using 100 repetitions of the Monte Carlo method.
The absolute error of 100 times repeated Monte Carlo simulations of the empirical covariance, Equation (2.2), compared to the theoretical covariance of the projected normal distribution, tr(Cov(pr(X)))/2. Each simulation run uses L data points and σ=1
L
30
50
100
1000
104
105
106
error
0.013
0.0080
0.0055
0.0033
0.0015
8.4e-05
2.9e-05
Secondly, Figure 3 is an illustration of the convergence result in Theorem 2.8. The simulations are done by generating L samples from X=dN(μ,σ2I3) with μ(0,0,1)T∈S2 and σ=1. The empirical covariance Vˆ, Equation (2.2), is computed for each set of samples. Then, for each set of samples, the estimator λˆ for σ2, see Theorem 2.8, is computed by f−1(tr(Vˆ)/2). This is done with a 1000-fold repetition.
A box plot of the inferred estimator λˆ in Theorem 2.8 using the empirical covariance, Equation (2.2), for L randomly generated projected normal random variables with σ2=1 and μ=(0,0,1)T. For each L this is done with 1000 repetitions. The true underlying scalar variance, λ=f−1(tr(Cov(pr(X)))/2)=1, is shown in red color. Note that some outliers are omitted for readability of the graph
Box plot comparing estimates of lambda for 5 sample sizes (L=5 to 100), showing convergence to true lambda (red line) as sample size increases.ProofsProof of Theorem <xref rid="j_vmsta279_stat_006">2.5</xref>
First, note that
1L−1∑ℓ=1Llogμˆ(ξℓ)logμˆ(ξℓ)T
is a linear map TμˆS2→TμˆS2, and Cov(ξ) is by definition a linear map TμS2→TμS2. In order to compare Vˆ and Cov(ξ) we shall transport parallelly the linear map Vˆ:Tμˆ→Tμˆ to a linear map TμS2→TμS2. That is, we look at Pμˆ,μVˆPμ,μˆ and see how far away it is from vI2. It follows from Equation (A.1) that Pμˆ,μlogμˆ(ξℓ)logμˆ(ξℓ)TPμ,μˆ has the first order expansion
Pμˆ,μlogμˆ(ξℓ)logμˆ(ξℓ)TPμ,μˆ=logμ(ξℓ)logμ(ξℓ)T+∇logμ(μˆ)logμ(ξℓ)logμ(ξℓ)T+logμ(ξℓ)∇logμ(μˆ)logμ(ξℓ)T+O(dist(μ,μˆ)2).
Therefore,
limn→∞EPμˆ,μ1n−1∑ℓ=1nlogμˆ(ξℓ)logμˆ(ξℓ)TPμ,μˆ−vI≤limn→∞1n−1∑ℓ=1nElogμ(ξℓ)logμ(ξℓ)T−vI+limn→∞E∇logμ(μˆ)logμ(ξℓ)logμ(ξℓ)T+logμ(ξℓ)∇logμ(μˆ)logμ(ξℓ)T+limn→∞EO(dist(μˆ,μ)2)≤vI−vI+limn→∞EO(dist(μˆ,μ))
Here, according to [8, Proposition 1], it holds that dist(μ,μˆ) converges to zero almost surely. Since M is compact, it has finite mass and thus dist(μ,μˆ) converges in expectation to 0 by the dominated convergence theorem. Finally, the fact that the empirical covariance has the same rate of convergence as the empirical mean follows from the last inequality. This completes the proof.
□
Proof of Theorem <xref rid="j_vmsta279_stat_008">2.6</xref>
By the very definition, we need to show that
arginfq∈S2∫S2dist(q,y)2ppr(X)(y)dS2(y)=pr(E[X])
for X=dN(μ,σ2I3). By rescaling X with a factor of 1μ and by a rotational symmetry, we may assume μ=(−1,0,0)T without loss of generality.
We make a simple first-derivative test to find the minimum in (3.1). Let (θ1,ϕ1) be the spherical coordinates that are integrated over the sphere and let (θ2,ϕ2) be the coordinates that minimize the integral. Hence we need to show that (θ2,ϕ2)=(π,π/2). Writing the integral in spherical coordinates and using the fact that the derivative with respect to ϕ2 is zero, it follows that (θ2,ϕ2) satisfies
∂∂ϕ2∫S2dist2((θ1,ϕ1),(θ2,ϕ2))ppr(X)(θ1,ϕ1)sin(ϕ1)dθ1dϕ1=0.
Now by the Leibniz integral rule the derivative may be moved inside the integral, and it holds
∫02π∫0π2dist((θ1,ϕ1),(θ2,ϕ2))×(cos(ϕ1)sin(ϕ2)−sin(ϕ1)cos(ϕ2)cos(θ2−θ1))1−cos(ϕ1)cos(ϕ2)+sin(ϕ1)sin(ϕ2)cos(θ2−θ1)2×ppr(X)(θ1,ϕ1)sin(ϕ1)dθ1dϕ1=0.
Plugging (θ2,ϕ2)=(π,π/2) into the integral on the left-hand side yields an integral
2∫02π∫0πacos(sin(ϕ1)cos(π−θ1))cos(ϕ1)1−sin2(ϕ1)cos2(π−θ1)×ppr(X)(θ1,ϕ1)sin(ϕ1)dθ1dϕ1.
This integral equals zero, which can be seen by observing that the integrand is odd along ϕ1 around π2, since ppr(X) is symmetric around μ by construction. By similar arguments, we obtain that
∂∂θ2∫S2dist2((θ1,ϕ1),(θ2,ϕ2))ppr(X)(θ1,ϕ1)sin(ϕ1)dϕ1dθ1
reduces at (θ2,ϕ2)=(π,π/2) to
2∫02π∫0πacos(sin(ϕ1)cos(π−θ1))sin2(ϕ1)sin(π−θ1)1−sin2(ϕ1)cos2(π−θ1)×ppr(X)(θ1,ϕ1)dϕ1dθ1
which equals zero by symmetry around θ1=π. Therefore, one can conclude that μ=pr(E[X]) is a local extremum for the integral in Equation (3.1).
Next we shall argue why it is a global minimum. Note that a very similar computation will show that −μ, i.e. in spherical coordinates (θ2,ϕ2)=(0,π/2) is another local extremum. In fact, μ and −μ are the only local extremum for the integral in Equation (3.1) since it is the only two points for which the distance function is rotationally symmetric around the line in R3 which is spanned by μ. By direct computation,
ppr(X)π,π2=1(2π)3/2exp(−12σ2)1σ+1σ2Φ(1σ)φ(1σ)+Φ(1σ)φ(1σ)
and
ppr(X)0,π2=1(2π)3/2exp(−12σ2)−1σ+1σ2Φ(−1σ)φ(−1σ)+Φ(−1σ)φ(−1σ),
hence ppr(X)(μ)>ppr(X)(−μ). More generally, a similar computation shows that ppr(X)(y)>ppr(X)(Ry), if ⟨y,μ⟩>0, where R is the reflection mapping over the y,z-plane, i.e.
Rabc=−abc.
Tautologically, the distance function increases when one goes father, and points being farther away will contribute more into the integral. Therefore (θ2,ϕ2)=(π,π/2) yields a global minimum for the integral in Equation (3.1).
□
Proof of Proposition <xref rid="j_vmsta279_stat_011">2.7</xref>
Without loss of generality, we assume μ=(0,0,1)T. In this case the functions A,B,C,K inside Equation (2.3) are given by
A=1σ2,B=1σ2cos(ϕ),C=−12σ2,
and
K=1σcos(ϕ).
By (2.4) it holds that Cov(pr(X)) is isotropic and we can write
tr(Cov(pr(X)))2=1(2π)1/2exp(−12σ2)×∫0πϕ2cos(ϕ)σ+cos2(ϕ)σ2+1Φ(cos(ϕ)σ)φ(cos(ϕ)σ)sin(ϕ)dϕ.
Denote
f(σ)=exp(−12σ2)∫0πϕ2cos(ϕ)σ+cos2(ϕ)σ2+1Φ(cos(ϕ)σ)φ(cos(ϕ)σ)sin(ϕ)dϕ.
In order to obtain the claim, it suffices to prove that f(σ) is strictly increasing, from which it follows that tr(Cov(pr(X)))2 is strictly increasing in σ as well. For notational simplicity, we set x=1σ and show that f(x) is strictly decreasing in x, where now f(x), with a slight abuse of notation, is given by
f(x)=exp(−x22)×∫0πϕ2xcos(ϕ)+x2cos2(ϕ)+1Φ(xcos(ϕ))φ(xcos(ϕ))sin(ϕ)dϕ.
By differentiating in x, we get
f′(x)=exp(−x22)(I1+I2+I3),whereI1:=−x∫0πϕ2cos(ϕ)xsin(ϕ)dϕ,I2:=∫0πϕ2−x3cos2(ϕ)−x+3xcos2(ϕ)+x3cos4(ϕ)Φ(xcos(ϕ))φ(xcos(ϕ))sin(ϕ)dϕ,I3:=∫0πϕ22cos(ϕ)+x2cos3(ϕ)sin(ϕ)dϕ.
Immediately, we have that
I1=x2π24
and
I3=−π22−5π232x2.
In order to show f′(x)<0, we need to decipher I2, which is more complicated. The main idea is to show that f′(x) is analytic for all x≥0 and then to show that there is a strictly decreasing analytic function between f′(x) and 0. We begin with Lemma 3.1 showing that Φφ is a Dawson-like function and therefore analytic. In Lemma 3.2 the series expansion of I2 is integrated term-wise. The terms of the series expression for I2 are given inductively in Lemma 3.3, which are then inserted to give explicit expressions for the terms of f′(x)exp(x2/2) in Lemma 3.4. Lemma 3.5 then gives a lower bound on the odd terms of f′(x)exp(x2/2), and Lemma 3.7 gives an upper bound. By Lemma 3.8 we show that the series expression of f′(x)exp(x2/2) is eventually decreasing as a series, and in combination with Corollary 3.6 it is concluded that f′(x)exp(x2/2) is entire. The proof is finished by comparing f′(x)exp(x2/2) to a linear combination of analytic functions computed in Lemma 3.9.
It holds thatΦφ(x)=∑k=0∞1(2k+1)!!x2k+1+2π2∑k=0∞1(2k)!!x2k=:∑k=0∞dkxkwhere the series converges everywhere.
Note that Φφ(x) is very similar to the Dawson function (originally studied in [4], see also [1] for further details) given by
D−(x)=exp(x2)∫0xexp(−t2)dt,
and it has the series expansion
D−(x)=∑k=0∞2k(2k+1)!!x2k+1.
Note that
∫0xexp(−t22)dt=2∫0x/2exp(−u2)du
by the variable substitution u=t/2. Hence,
exp(x22)∫0xexp(−t22)dt=2D−(x/2)=2∑k=0∞2k(2k+1)!!x2k+12k2=∑k=0∞1(2k+1)!!x2k+1.
Moreover,
∫−∞0exp(−t22)dt=2π2
and hence
exp(x22)∫−∞0exp(−t22)dt
has the series expansion
2π2∑k=0∞1k!x2k2k=2π2∑k=0∞1(2k)!!x2k.
It follows that
Φφ(x)=exp(x22)∫−∞0exp(−t22)dt+exp(x22)∫0xexp(−t22)dt=∑k=0∞1(2k+1)!!x2k+1+2π2∑k=0∞1(2k)!!x2k
proving the claimed series representation. Finally, the convergence everywhere follows from the fact that the Dawson function converges everywhere. □
Using Lemma 3.1 it follows that I2 can be rewritten as
I2=∑k=0∞dk∫0πϕ2−x3cos2(ϕ)−x+3xcos2(ϕ)+x3cos4(ϕ)×xkcosk(ϕ)sin(ϕ)dϕ.
Each term in Equation (3.3) involves
Jm=∫0πϕ2cosm(ϕ)sin(ϕ)dϕ
which we will study next.
The sequenceJmsatisfiesJ0=π2−4, and for evenm≠0Jm=1m+1π2−4m!!(m+1)!!1+∑j=0m2−1122j+1(2j+3)2j+1j,while form∈NoddJm=π2m+1m!!(m+1)!!−1.
For J0 we observe immediately that
J0=∫0πϕ2sin(ϕ)dϕ=π2−4.
Next, let m be odd. Then integration by parts gives
Jm=−ϕ2cosm+1(ϕ)m+10π+2m+1∫0πϕcosm+1(ϕ)dϕ=−π2m+1+2m+1[m!!(m+1)!!ϕ2+ϕ∑j=0m−12cosm−2j(ϕ)sin(ϕ)m!!(m−2j)!!×(m−2j−1)!!(m+1)!!]0π−2m+1∑j=0m−12m!!(m−2j)!!(m−2j−1)!!(m+1)!!×∫0πcosm−2j(ϕ)sin(ϕ)dϕ−2m+1m!!(m+1)!!∫0πϕdϕ=−π2m+1+2m+1m!!(m+1)!!π2−1m+1m!!(m+1)!!π2=π2m+1m!!(m+1)!!−1.
Similarly, when m>0 is even, integration by parts gives
Jm=−ϕ2cosm+1(ϕ)m+10π+2m+1∫0πϕcosm+1(ϕ)dϕ=π2m+1+2m+1[ϕsin(ϕ)m!!(m+1)!!+ϕ∑j=0m2−1cosm−2jsin(ϕ)m!!(m−2j)!!×(m−1−2j)!!(m+1)!!]0π−2m+1m!!(m+1)!!∫0πsin(ϕ)dϕ−2m+1∑j=0m2−1m!!(m−2j)!!(m−1−2j)!!(m+1)!!∫0πcosm−2j(ϕ)sin(ϕ)dϕ=π2m+1−4m+1m!!(m+1)!!−4m+1m!!(m+1)!!∑j=0m2−1(m−1−2j)!!(m−2j)!!(m+1−2j)=1m+1π2−4m!!(m+1)!!1+∑j=0m2−1122j+1(2j+3)2j+1j.
This completes the proof. □
DenoteI2=∑n=1∞anxn.Then the coefficientsaksatisfya1=42π9,a2=−7π232,andak=dk−1(3Jk+1−Jk−1)+dk−3(Jk+1−Jk−1),k≥3.
Immediately I2(0)=0, so the constant term is zero. Next, we get expressions for an in terms of dn’s and Jn’s from Equation (3.3). Then straightforward calculations give
a1=d0(−J0+3J2)=2π24−π2+33π2−4231+16=42π9,a2=d1(−J1+3J3)=π24+3π2438−1=−7π232,
and, for k≥3,
ak=dk−1(3Jk+1−Jk−1)+dk−3(Jk+1−Jk−1).
□
Considering the expressions for I1, I2, and I3 gives
c0=−π22,c1=a1=42π9,c2=a2+3π232=−π28,
and
ck=akwheneverk≥3,
where ak is given in Lemma 3.3. It follows that for k≥3 odd we have
ck=dk−1(3Jk+1−Jk−1)+dk−3(Jk+1−Jk−1)=2π2(1(k−1)!![3k+2(π2−4(k+1)!!(k+2)!!×(1+∑j=0k+12−1122j+1(2j+3)2j+1j))−1k(π2−4(k−1)!!k!!(1+∑j=0k−12−1122j+1(2j+3)2j+1j))]+1(k−3)!![1k+2(π2−4(k+1)!!(k+2)!!×(1+∑j=0k+12−1122j+1(2j+3)2j+1j))−1k(π2−4(k−1)!!k!!(1+∑j=0k−12−1122j+1(2j+3)2j+1j))]),
and after simplifying we have
ck=22π(k+2)!!1+∑j=0k−12−1122j+1(2j+3)2j+1j−k+12k(k+2)kk−12.
Similarly for k≥4 even we get
ck=dk−1(3Jk+1−Jk−1)+dk−3(Jk+1−Jk−1)=π2(k−1)!!3k+2(k+1)!!(k+2)!!−1−1k(k−1)!!k!!−1+π2(k−3)!!1k+2(k+1)!!(k+2)!!−1−1k(k−1)!!k!!−1=−π2(k+2)!!.
This completes the proof. □
For allk≥3odd, setS(k)=1+∑j=0k−12−1122j+1(2j+3)2j+1j−k+12k(k+2)kk−12.ThenS(k)is increasing and satisfies22πS(k)≥π.
First note that
22πS(3)=22π1+1610−4512331=22π1+16−310=262π15≥π.
It remains to show that S(k) is an increasing function. We have
S(k+2)=1+∑j=0k−12122j+1(2j+3)2j+1j−k+32k+2(k+4)k+2k+12=S(k)+k+12k(k+2)kk−12+12k(k+2)kk−12−k+32k+2(k+4)k+2k+12=S(k)+12kkk−12−k+32k+2(k+4)k+2k+12≥S(k)+12kkk−12−12k+2k+2k+12=S(k)+(k+1)!2k+2k+12!k+32!k+3−(k+2)=S(k)+(k+1)!2k+2k+12!k+32!
and hence S(k+2)≥S(k). This completes the proof. □
Forx>0the series∑k=0∞ckxk,where the coefficientsckare determined by (3.4), is alternating.
For k even we clearly have ck<0. On the other hand, for k odd ck is positive since
ck≥π(k+2)!!
by Lemma 3.5. □
The following result shows that the sum in Lemma 3.5 can also be bounded from above.
For allk≥3odd,S(k)defined by (3.5) satisfies22πS(k)≤π2.
By elementary manipulations as in the proof of Lemma 3.5 we obtain
S(k+2)=S(k)+22π(k+3)2k+2(k+2)(k+4)k+2k+12.
From the Stirling approximation we infer
k+2k+12≤2π2k+2(k+2),
and hence
S(k+2)≤S(k)+2π(k+2)3/2.
Moreover,
S(3)=1+16−310≤1+133/2
allows us to get the estimates
limk→∞22πS(k)≤4∑k=0∞1(2k+1)3/2≤π2.
Since S is increasing by Lemma 3.5, it follows that S(k)≤π2 for all k≥3 odd. This completes the proof. □
Let the coefficientsck,k=1,2,…, be given by (3.4) and letx>0be fixed. Then the terms in the series∑k=M∞ckxkdecreases monotonically for M large enough.
Consider first the case when k is odd. Then, using the upper bound of Lemma 3.7 yields
ck+1xk+1ckxk≤π2(k+3)!!π2(k+2)!!x=(k+2)!!(k+3)!!x
which is less than one for sufficiently large k (depending on x). Similarly for k even we can use the lower bound from Lemma 3.5 to obtain
ck+1xk+1ckxk≤π2(k+3)!!π(k+2)!!x=(k+2)!!(k+3)!!πx
which again is less than one for large enough k. This yields the claim. □
By differentiating, we get
M′(x)=∑k=0∞(2k+2)x2k+1(2k+4)!!=∑k=0∞x2k+1(2k+3)!!(k+1)22k+22k+3k+1.
Now using
2k+3k+1≥22k+2k+1
gives
M′(x)≥∑k=0∞k+1(2k+3)!!x2k+1≥N(x).
Similarly, it holds that
N′(x)=∑k=0∞(2k+1)x2k(2k+3)!!=∑k=0∞x2k(2k+2)!!(2k+1)22k+2(k+2)2k+3k+1.
In this case we can use
2k+3k+1≥22k+32k+3
leading to
N′(x)≤∑k=0∞x2k(2k+2)!!k+12(k+2)(2k+3)≤M(x).
Combining the two bounds above gives us
ddxM2(x)−N2(x)=2M(x)M′(x)−2N(x)N′(x)≥M(x)N(x)−M(x)N(x)=0.
Consequently, M2−N2 is an increasing function for x≥0, which leads to
M2(x)−N2(x)≥M2(0)−N2(0)=1.
It follows that
M(x)≥1+N2(x)>N(x),
and the proof is complete. □
We are finally in the position to prove Proposition 2.7.
By Corollary 3.6 and Lemma 3.8 the series expansion (3.4) for f′(x)exp(x2/2) is convergent by the Leibniz alternating series test. Since the radius of convergence is unbounded, it is analytic, and thus we may split the series into its positive part and negative part, given by
∑k=0∞c2kx2k=:Q(x),∑k=0∞c2k+1x2k+1=:P(x).
By Lemma 3.4 it holds that Q(x)=−π2M(x) where M(x) is as in Lemma 3.9. Now Lemma 3.7 gives
c2k≤π2(k+2)!!,
and hence
Q(x)≤π2N(x),
where N(x) is as in Lemma 3.9. Therefore
f′(x)exp(x2/2)=Q(x)+P(X)≤−π2M(x)+π2N(x).
Applying Lemma 3.9 to the above inequality now gives
f′(x)exp(x22)<0,
and thus f′(x)<0 for all x≥0, where f is given by (3.2). The claim follows from this. □
Some results from differential geometry
Here we introduce some classical concepts about manifolds and their properties used in the paper.
Let M be a smooth and compact n-dimensional manifold embedded in Rk. At each point, the tangent space, TxM, is equipt with the Euclidean inner product (Riemannian metric), inherited from the ambient space Rk making M an isometrically embedded manifold as a Riemannian manifold. The Riemannian metric, denoted ·,·, induces a notion of length along smooth curves, γ:[a,b]→M, by the formula
L(γ)=∫abγ˙(t)dt.
Given two points, or equivalently a starting point and starting velocity, the curve which minimizes this distance is called a geodesic. Centered at a point x∈M, the map expx:TxM→M is called the Riemannian exponential map centered at x∈M. Given a tangent vector vx∈TxM, let γ be the unit speed geodesic such that γ(0)=x and γ˙(0)=vxvx, then
exp(vx)=γ(vx).
We shall assume that M is geodesically complete, meaning that the exponential map is defined on the whole TxM for all x∈M. Denote the inverse of expp (restricted to all points y for which there is a unique geodesic connecting y to p) by logp. The distance function, dist:M×M→[0,∞) (which is a topological metric) is then defined as dist(x,y)=logx(y). The Riemannian metric also gives a way to measure the volume (and orientation) of a parallelotope inside the tangent space, i.e. a rescaling of the determinant from linear algebra. This generalized determinant shall be referred to as dVolM, or the Riemannian volume form, which locally looks like the Lebesgue measure. By vector field, we mean a smooth assignment of a point x∈M to a tangent vector TxM, and the space of vector fields is denoted by X(M). For a given smooth function f∈C∞(M), a vector field X acts as a derivation on f in the following sense. Take a smooth curve (−ε,ε)→M such that γ(0)=p, γ˙(0)=X(p). Then the function X(f):M→R is point-wise defined by X(f)(p)=(f∘γ)′(0).
Note that for given vector fields X,Y:M→TM, the expression X,Y is a smooth function M→R. In order to differentiate vector fields, a notion of connection is required. Here we shall use the Levi-Civita connection∇:X(M)×X(M)→X(M) which is uniquely defined by the following identities:
∇fX+gYZ=f∇XZ+g∇YZ
for all smooth functions f,g∈C∞(M) and all smooth vector fields
X,Y,Z∈X(M);
∇X(gY+hZ)=g∇XY+X(g)Y+h∇XZ+X(h)Z
for all smooth functions g,h∈C∞(M) and all smooth vector fields
X,Y,Z∈X(M);
X(Y,Z)=∇X(Y),Z+Y,∇X(Z)
for all smooth vector fields X,Y,Z∈X(M);
(∇XY−∇YX)(f)=X(Y(f))−Y(X(f))
as derivations for all smooth functions f∈C∞(M) and all smooth vector fields X,Y∈X(M).
Lastly, let Pq,ℓ:TqM→TℓM denote the parallel transport from q to ℓ. That is, Pq,ℓ(vq) is the unique point τ(dist(q,ℓ))∈TℓM that satisfies the initial value problem
∇γ˙τ(t)=0,τ(0)=vq,
where γ:[0,dist(q,ℓ)]→M is the geodesic connecting q to ℓ and τ is a vector field along γ, i.e. τ(t)=τγ(t). It is worthwhile to consider the parallel transport of the Riemannian logarithm map
q⟼Pq,ℓlogq(ν),
where ν,ℓ∈M are fixed points. This map has the Taylor expansion around q=ℓ given by
Pq,ℓlogq(ν)=logℓ(ν)+∇logℓ(q)logℓ(ν)+O(dist(q,ℓ)2).
ReferencesAbramowitz, M., Stegun, I.A.: Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables vol. 55. US Government printing office (1968) MR0415956Ahidar-Coutrix, A., Le Gouic, T., Paris, Q.: Convergence rates for empirical barycenters in metric spaces: curvature, convexity and extendable geodesics. Probab. Theory Relat. Fields177(1), 323–368 (2020) MR4095017. https://doi.org/10.1007/s00440-019-00950-0Chang, T.: Spherical regression and the statistics of tectonic plate reconstructions. Int. Stat. Rev., 299–316 (1993)Dawson, H.G.: On the numerical value of ∫0hex2dx. Proc. Lond. Math. Soc.1(1), 519–522 (1897) MR1576451. https://doi.org/10.1112/plms/s1-29.1.519Fisher, N.I., Lewis, T., Embleton, B.J.J.: Statistical Analysis of Spherical Data. Cambridge University Press (1993) MR1247695Gao, F., Chia, K., Machin, D.: On the evidence for seasonal variation in the onset of acute lymphoblastic leukemia (all). Leuk. Res.31(10), 1327–1338 (2007)Hernandez-Stumpfhauser, D., Breidt, F.J., van der Woerd, M.J.: The general projected normal distribution of arbitrary dimension: Modeling and Bayesian inference. Bayesian Anal.12(1), 113–133 (2017) MR3597569. https://doi.org/10.1214/15-BA989Hotz, T., Le, H., Wood, A.T.A.: Central limit theorem for intrinsic Frechet means in smooth compact Riemannian manifolds. Probab. Theory Relat. Fields (2024) MR4771114. https://doi.org/10.1007/s00440-024-01291-3Hussien, M.T., Seddik, K.G., Gohary, R.H., Shaqfeh, M., Alnuweiri, H., Yanikomeroglu, H.: Multi-resolution broadcasting over the Grassmann and Stiefel manifolds. In: 2014 IEEE International Symposium on Information Theory, pp. 1907–1911 (2014). IEEEKells, L.M., Kern, W.F., Bland, J.R.: Plane and Spherical Trigonometry. US Armed Forces Institute (1940) MR1524230. https://doi.org/10.2307/2302987Lang, M., Feiten, W.: Mpg-fast forward reasoning on 6 dof pose uncertainty. In: ROBOTIK 2012; 7th German Conference on Robotics, pp. 1–6 (2012). VDELevitt, M.: A simplified representation of protein conformations for rapid simulation of protein folding. J. Mol. Biol.104(1), 59–107 (1976)Li, K., Frisch, D., Radtke, S., Noack, B., Hanebeck, U.D.: Wavefront orientation estimation based on progressive Bingham filtering. In: 2018 Sensor Data Fusion: Trends, Solutions, Applications (SDF), pp. 1–6 (2018). IEEENuñez-Antonio, G., Ausín, M.C., Wiper, M.P.: Bayesian nonparametric models of circular variables based on Dirichlet process mixtures of normal distributions. J. Agric. Biol. Environ. Stat.20, 47–64 (2015) MR3334466. https://doi.org/10.1007/s13253-014-0193-yNuñez-Antonio, G., Gutiérrez-Peña, E.: A Bayesian analysis of directional data using the projected normal distribution. J. Appl. Stat.32, 995–1001 (2005) MR2221902. https://doi.org/10.1080/02664760500164886Nuñez-Antonio, G., Gutiérrez-Peña, E., Escarela, G.: A Bayesian regression model for circular data based on the projected normal distribution. Stat. Model.11, 185–201 (2011) MR2849683. https://doi.org/10.1177/1471082X1001100301Pennec, X.: Intrinsic statistics on Riemannian manifolds: Basic tools for geometric measurements. J. Math. Imaging Vis.25, 127–154 (2006) MR2254442. https://doi.org/10.1007/s10851-006-6228-4Pennec, X., Arsigny, V.: Exponential barycenters of the canonical Cartan connection and invariant means on Lie groups. In: Matrix Information Geometry, pp. 123–166. Springer (2012) MR2964451. https://doi.org/10.1007/978-3-642-30232-9_7Pitaval, R., Tirkkonen, O.: Joint Grassmann-Stiefel quantization for MIMO product codebooks. IEEE Trans. Wirel. Commun.13(1), 210–222 (2013)Presnell, B., Morrison, S.P., Littell, R.C.: Projected multivariate linear models for directional data. J. Am. Stat. Assoc.93, 1068–1077 (1998) MR1649201. https://doi.org/10.2307/2669850Seddik, K.G., Gohary, R.H., Hussien, M.T., Shaqfeh, M., Alnuweiri, H., Yanikomeroglu, H.: Multi-resolution multicasting over the Grassmann and Stiefel manifolds. IEEE Trans. Wirel. Commun.16(8), 5296–5310 (2017)Sveier, A., Egeland, O.: Pose estimation using dual quaternions and moving horizon estimation. IFAC-PapersOnLine51(13), 186–191 (2018)Tanabe, A., Fukumizu, K., Oba, S., Takenouchi, T., Ishii, S.: Parameter estimation for von mises–fisher distributions. Comput. Stat.22, 145–157 (2007) MR2299252. https://doi.org/10.1007/s00180-007-0030-7Wang, F., Gelfand, A.E.: Directional data analysis under the general projected normal distribution. Stat. Methodol.10(1), 113–127 (2013) MR2974815. https://doi.org/10.1016/j.stamet.2012.07.005