An algorithm for simulating Brownian increments on a sphere
AAN ALGORITHM FOR SIMULATING BROWNIAN INCREMENTS ON A SPHERE
ALEKSANDAR MIJATOVIĆ, VENO MRAMOR, AND GERÓNIMO URIBE BRAVO
Abstract.
This paper presents a novel formula for the transition density of the Brownian motion on asphere of any dimension and discusses an algorithm for the simulation of the increments of the sphericalBrownian motion based on this formula. The formula for the density is derived from an observation thata suitably transformed radial process (with respect to the geodesic distance) can be identified as a Wright-Fisher diffusion process. Such processes satisfy a duality (a kind of symmetry) with a certain coalescentprocesses and this in turn yields a spectral representation of the transition density, which can be used forexact simulation of their increments using the results of Jenkins and Spanò (2017). The symmetry thenyields the algorithm for the simulation of the increments of the Brownian motion on a sphere. We analysethe algorithm numerically and show that it remains stable when the time-step parameter is not too small. Introduction
Brownian motion is essentially a continuous time symmetric random walk and is therefore of significantimportance in science. Classically, Brownian motion is defined on a (flat) Euclidean space and this pro-cess is very well understood. For many applications, however, it is better to model the state space as acurved surface or some other manifold and then Brownian motion on the manifold is of interest. The def-inition of Brownian motion on any Riemannian manifold is possible via the Laplace-Beltrami operator asin [Hsu02] and then a typical example is Brownian motion on the sphere S ( R ) in the three dimensionalEuclidean space R . The spherical Brownian motion has been used to model fluorescent markers in cellmembranes [KDPN00], motion of bacteria [LTT08], migration of marine animals [BS98] and to study thewhole world phylogeography [Bou16], to name a few examples. Brownian motions on more general surfacesand manifold are also of interest (cf. [Far02]), but many manifolds (at least those with positive curvature)can be locally approximated by a sphere, so we henceforth only focus on a spherical Brownian motion.However, whereas the simulation of Euclidean Brownian increments is easy, as it reduces to the simulationof Gaussian random variables, the situation on a sphere is more involved. There is no closed form for thetransition density of a spherical Brownian motion, so one is usually forced to use approximations. The firstapproach is to locally approximate the sphere with a (flat) tangent space and then suitably project theGaussian increment back to the sphere. This method was for example been applied in [KDPN00], but itsdrawback is that it ignores the curvature of the sphere, so that it is a good approximation only for very smalltime-steps or very large radii. Later [NEE03] provided a better approximation for the Brownian motion onthe sphere in R and later the algorithm has been adapted to work also on the sphere in R [CEE10].In [GSS12] a different approximation is derived directly for the Brownian motion on S ( R ) . These improvedmethods not only give good results for small time-steps, but continue to provide a good approximation alsofor medium-sized time-steps. Nevertheless, all these algorithms are just approximate and as the time-stepincreases the results become less and less accurate.The main result in this article is a new representation of a transition density of the spherical Brownianmotion given in Equation (2.4) and a numerical analysis of an algorithm for the simulation of the increments a r X i v : . [ c ond - m a t . s t a t - m ec h ] D ec N ALGORITHM FOR SIMULATING BROWNIAN INCREMENTS ON A SPHERE 2 of spherical Brownian motion arising from it. The representation is a consequence of the skew-productdecomposition of the spherical Brownian motion obtained in [MMU20] and the representation of transitiondensities of Wright-Fisher diffusion processes obtained in [GL83, JS17]. An algorithm for the simulationof the increments of spherical Brownian motion is derived from the representation and is a generalizationof [MMU20, Algorithm 1] allowing general diffusion coefficients and radii of the sphere.In contrast to the existing literature in which algorithms produce samples from an approximate distribu-tion and are hence forced to use smaller time-steps for accurate results, our algorithm is exact, so that itproduces samples directly from the required distribution of the increments, and can handle arbitrary largetime-steps. In fact, its performance improves with increasing time-steps. The limiting factor of our algorithmis actually imperfect floating point arithmetic on computers. Since several quantities in the algorithm haveto be computed by complicated expression which become numerically unstable for smaller time-steps, we areactually limited with how small a time-step we are allowed to take. As such, our algorithm complements theother available methods. Moreover, in the numerical examples of Section 3 below it appears to outperformthem when the parameters are in the correct domain so we can use our algorithm.2.
Theory
We are interested in a Brownian motion on the sphere S d − ( R ) := (cid:110) (cid:126)x ∈ R d ; | (cid:126)x | = R (cid:111) , where | (cid:126)x | := (cid:113) x + · · · + x d is an absolute value of a vector (cid:126)x = ( x , . . . ,x d ) (cid:62) ∈ R d so that R > representsthe radius of the sphere and d ≥ is the dimension of the ambient Euclidean space. Hence, if we denote by ρ ( D,R ) (cid:126)y ( (cid:126)x,t ) a transition density of the spherical Brownian motion started at (cid:126)y ∈ S d − ( R ) , we are interestedin the diffusion equation:(2.1) ∂ρ ( D,R ) (cid:126)y ∂t ( (cid:126)x,t ) = D ∇ S d − ( R ) ρ ( D,R ) (cid:126)y ( (cid:126)x,t ) , ρ ( D,R ) (cid:126)y ( (cid:126)x,
0) = δ ( (cid:126)x,(cid:126)y ) , where ∇ S d − ( R ) is the Laplace-Beltrami operator on the sphere S d − ( R ) , D > is a diffusion coefficientand (cid:126)y ∈ S d − ( R ) is the initial point. Symmetry of the spherical Brownian motion allows us to deduce thefollowing fact: if A ∈ R d ⊗ R d is any orthogonal matrix, then ρ ( D,R ) A(cid:126)y ( A(cid:126)x,t ) = ρ ( D,R ) (cid:126)y ( (cid:126)x,t ) . This in particularshows that it is enough to consider a single initial point – we are going to use the north pole R · (cid:126) e d ∈ S d − ( R ) , where (cid:126) e d := (0 , . . . , , (cid:62) represents the north pole of the unit sphere. Furthermore, symmetry allows us todeduce that ρ ( D,R ) (cid:126)y ( (cid:126)x,t ) does not depend on the whole vector (cid:126)x but only depends on the geodesic distancebetween the vectors (cid:126)x and (cid:126)y. We are using the standard round metric on the sphere which is induced by the standard Euclidean scalarproduct (cid:104) (cid:126)a,(cid:126)b (cid:105) = a b + · · · + a d b d . The geodesic distance measures the shortest distance between two points(i.e. along great circles which are intersections of the sphere with two dimensional planes through the origin)and is given as a function dis : S d − ( R ) × S d − ( R ) → [0 ,Rπ ] , dis( (cid:126)x,(cid:126)y ) = R arccos( (cid:104) (cid:126)x,(cid:126)y (cid:105) /R ) =: Rθ i.e. thedistance is equal to R times the angle between the two vectors (here denoted by θ ). Usually, one looks at thestandard spherical Brownian motion which diffuses on a sphere S d − := S d − (1) of radius and correspondsto D = 1 / with initial point (cid:126) e d . Since in unit sphere coordinates relation ∇ S d − ( R ) = ∇ S d − /R holds, wecan define a new parameter τ := 2 Dt/R , which we use instead of time and then ρ ( D,R ) R(cid:126) e d ( (cid:126)x,t ) = ρ (1 / , (cid:126) e d ( (cid:126)x/R,τ ) holds. This justifies the focus on the standard spherical Brownian motion since we can extend the results to N ALGORITHM FOR SIMULATING BROWNIAN INCREMENTS ON A SPHERE 3 arbitrary radius and diffusion coefficient by rescaling and linear time-change. We will henceforth omit theparameters D = 1 / and R = 1 from the notation and simply write ρ (cid:126) e d ( (cid:126)z,τ ) for the transition density attime τ of the standard spherical Brownian motion started at (cid:126) e d . Any point (cid:126)z ∈ S d − can be represented as (cid:126)z = sin θ (cid:126)w + cos θ(cid:126) e d , where (cid:126)w ∈ S d − ⊆ R d − × { } and θ ∈ [0 ,π ] is the angle between the vectors (cid:126)z and (cid:126) e d . Due to symmetry we see that ρ (cid:126) e d ( (cid:126)z,τ ) = ρ (cid:126) e d (sin θ (cid:126)w + cos θ(cid:126) e d ,τ ) does not depend on the vector (cid:126)w but only depends on the angle θ . Hence we can and will consider theone-dimensional density ρ ( θ,τ ) := Vol( S d − (sin θ )) ρ (cid:126) e d (sin θ (cid:126)w + cos θ(cid:126) e d ,τ ) = A d − sin d − θ · ρ (cid:126) e d (sin θ (cid:126)w + cos θ(cid:126) e d ,τ ) , where A n := 2 π ( n +1) / / Γ( n +12 ) is the volume of an n -dimensional sphere S n . The additional multiplicativefactor A d − sin d − θ is included due to only considering a one-dimensional process (and its transition density)instead of the whole process on the sphere. For τ = 0 the initial distribution ρ ( θ, of this density is equalto the delta function with support at . Using (ultra)spherical harmonics and their addition formula wecan get (see [KT81, p. 339], where their Jacobi polynomials are normalised differently to our Gegenbauerpolynomials) a formal solution of the diffusion equation at time τ :(2.2) ρ ( θ,τ ) = A d − sin d − θ ∞ (cid:88) l =0 h ( l,d ) A d − C ( d/ − l (cos θ ) C ( d/ − l (1) e λ l τ/ , where C ( α ) n are Gegenbauer polynomials given by their generating series (1 − xt + t ) − α = (cid:80) ∞ n =0 C ( α ) n ( x ) t n and h ( l,d ) = (cid:0) d + l − d − (cid:1) − (cid:0) d + l − d − (cid:1) is the dimension of the space of spherical harmonics of degree l , correspond-ing to the eigenvalue λ l := − l ( l + d − . Equation (2.2) without the additional factor A d − sin d − θ (i.e.the density for the process on the sphere) has also been derived in [Cai04] where a further simplification h ( l,d ) /C ( d/ − l (1) = (2 l + d − / ( d − is used.In the special case d = 3 , i.e. for the standard unit sphere S , equation (2.2) specializes to the well-knownformula(2.3) ρ ( θ,τ ) = sin θ ∞ (cid:88) l =0 (2 l + 1) P l (cos θ ) e − l ( l +1) τ/ , where P l is the l -th Legendre polynomial. Unfortunately, formal solutions given by formulas (2.2) or (2.3)oscillate a lot even when large amount of terms are taken into account and are therefore not suitable forsimulation. One is then usually forced to use some approximative density and sampling from it.The simplest option is to locally approximate a sphere with its tangent space and use a standard Brow-nian increment which is then translated to an increment on a sphere using exponential map i.e. movingan appropriate distance along the geodesics (great circles) of a sphere. Such an approximation correspondsto density Q T ang ( θ,τ ) = N ( τ ) θ d − exp( − θ / (2 τ )) , where N ( τ ) is a normalization constant such that (cid:82) π Q T ang ( θ,τ )d θ = 1 . This approximation largely ignores the curvature of the sphere so it is a good approxi-mation only for very small values of parameter τ . Therefore, [GSS12] proposed an improved approximation.It is given by Q Approx ( θ,τ ) = N ( τ )( θ sin θ ) ( d − / exp( − θ / (2 τ )) , where N ( τ ) is a normalization constant.This new approximation continues to give good results even for an intermediate values of τ. In the case of d = 4 the same approximative density was derived in [NEE03] by using a different method.The main result in this paper is the following alternative representation of a solution of the diffusionequation, which is much more suitable for simulation. It also allows us to do exact simulation in which we N ALGORITHM FOR SIMULATING BROWNIAN INCREMENTS ON A SPHERE 4 are sampling directly from the transition density and not just from an approximative density as in the otherpreviously mentioned methods. The alternative representation of the density (2.2) is given by(2.4) ρ ( θ,τ ) = sin θ ∞ (cid:88) m =0 q d − m ( τ ) g ( d − , d − + m ) (cid:18) − cos θ (cid:19) , where g ( α,β ) ( x ) := 1B( α,β ) x α − (1 − x ) β − , x ∈ [0 , , α,β > in which B represents the Beta function so that g ( α,β ) is a density of a Beta distribution and q d − m ( τ ) := ∞ (cid:88) k = m ( − k − m b ( τ,d − k ( m ) , (2.5) b ( τ,d − k ( m ) := ( d + 2 k − d − m ) ( k − m !( k − m )! e − k ( k + d − τ/ , where the Pochhammer symbol is given by a ( x ) = Γ( a + x )Γ( a ) for a > , x ≥ − . As a first impression the representation (2.4) looks just as unwieldy, if not more, than Equation (2.2), butit does have a nice probabilistic interpretation, which makes it suitable for simulation. First, the coefficients q d − m ( τ ) are actually non-negative and they sum up to 1 (see Subsection 2.1 below). Additionally, each term inthe sum corresponds to the density of a Beta( d − , d − + m ) distribution after we have used a transformationof the form h ( x ) = arccos(1 − x ) (the inverse transformation is h − ( θ ) = − cos θ and its derivative is exactlythe factor sin θ in front of the sum). Therefore, we see that ρ ( θ,τ ) is a mixture of densities of transformedBeta distributions.This immediately yields an algorithm for the simulation from the density ρ ( θ,τ ) . If we denote by A d − ∞ ( τ ) an N -valued random variable with P [ A d − ∞ ( τ ) = m ] = q d − m ( τ ) , then we can use the standard inver-sion sampling to simulate A d − ∞ ( τ ) i.e. for a random variable U ∼ Uniform(0 , , the random variable inf (cid:110) M ∈ N ; (cid:80) Mm =0 q d − m ( t ) > U (cid:111) is distributed as A d − ∞ ( t ) . Given this value, we then simulate an appropri-ate Beta distributed random variable and after a final transformation we have found our sample from thedensity ρ ( θ,τ ) .Our goal is actually to simulate an S d − ( R ) valued random variable correctly distributed as the incrementof the Brownian motion on the sphere for an arbitrary initial point, radius and diffusion coefficient. Therefore,we have to recall symmetries of spherical Brownian motion i.e. ρ ( D,R ) A(cid:126)y ( A(cid:126)x,t ) = ρ ( D,R ) (cid:126)y ( (cid:126)x,t ) for any orthogonalmatrix A . In particular, by taking an orthogonal matrix which fixes the initial point, it implies that there isa component to the increment, which is uniformly distributed on S d − . Combining this fact with the equality ρ ( D,R ) R(cid:126) e d ( (cid:126)x,t ) = ρ (cid:126) e d ( (cid:126)x/R,τ ) , where τ = 2 Dt/R , and the aforementioned simulation from density ρ ( θ,τ ) we getthe following algorithm for the exact simulation of the increment of the spherical Brownian motion. Algorithm 1
Simulating from the transition density ρ ( D,R ) (cid:126)y ( · ; t ) of spherical Brownian motion Set R = | (cid:126)y | ,τ = 2 Dt/R Simulate M ∼ A d − ∞ ( τ ) Simulate X ∼ Beta( d − , d − + M ) Simulate Y uniformly distributed on S d − Set O ( (cid:126)y ) = I − uu (cid:62) , where u = ( (cid:126) e d − (cid:126)y/R ) / | (cid:126) e d − (cid:126)y/R | . return RO ( (cid:126)y )(2 (cid:112) X (1 − X ) Y (cid:62) , − X ) (cid:62) N ALGORITHM FOR SIMULATING BROWNIAN INCREMENTS ON A SPHERE 5
Step 4 in Algorithm 1 consists of simulating a vector N in R d − with independent standard normalcomponents and setting Y = N/ | N | . The key property of the orthogonal matrix O ( (cid:126)y ) ∈ R d ⊗ R d inAlgorithm 1 is O ( (cid:126)y ) (cid:126) e d = y/ | y | . In fact, any orthogonal matrix in R d ⊗ R d with this property would leadto an exact sample. The formula for O ( (cid:126)y ) in Algorithm 1 is chosen due to its simplicity. Additionally,if only the simulation from the density ρ ( θ,τ ) is required, we only have to do first three steps and thenrandom variable arccos(1 − X ) has the adequate law. Step 2 constitutes the main technical difficulty ofthe Algorithm 1 and all potential numerical issues stem from it. By taking R = 1 and D = 1 / we see thatthis algorithm is a more general version of [MMU20, Algorithm 1].2.1. Derivation of the alternative representation of the density in (2.4) and discussion.
We willprove that the density ρ ( θ,t ) is indeed represented as in (2.4) and hence justify the correctness of Algorithm 1.To do that we have to notice that density ρ ( θ,t ) is actually the transition density of the one-dimensional [0 ,π ] -valued process Y t = dis( Z t ,(cid:126) e d ) , where Z t denotes the standard Brownian motion on the unit sphere S d − started at (cid:126) e d . Since cos( Y t ) is nothing but the last component of the spherical Brownian motion Z t ,it is quite natural to further transform the process Y t and instead look at X t = − cos( Y t )2 . This process is [0 , -valued and surprisingly turns out to be a Wright-Fisher diffusion process, a member of a family ofprocesses which are well studied and used extensively in genetics. Computations in [MMU20] show that inour case the process X t is a Wright-Fisher diffusion process with parameters θ = d − = θ and an initialpoint y = 0 (see also [JS17, Eqs. (1),(3)]). Its transition density f ( x ; t ) is a solution of a Fokker-Planckequation ∂f ( x ; t ) ∂t = ∂ ∂x (cid:18) x (1 − x )2 f ( x ; t ) (cid:19) − ∂∂x (cid:18)(cid:18) d −
14 (1 − x ) − d − x (cid:19) f ( x ; t ) (cid:19) , f ( x ; 0) = δ ( x, . The transition densities of the Wright-Fisher diffusion processes have a spectral decomposition givenby the associated Jacobi polynomials [GS10], but such a representation is not suitable for simulation (itis essentially equivalent to representation (2.2)). Fortunately, there exists another representation of theirtransition densities, which arises from the moment duality between Wright-Fisher diffusion processes andcertain coalescent processes. This representation and its use for simulation was given in [GL83] and it isgiven explicitly in [JS17, Eq. (4)]. The transformation formula for densities immediately yields ρ ( θ,τ ) = sin θ f (cid:0) − cos θ ; τ (cid:1) and combining it with [JS17, Eq. (4)] we immediately get the representation (2.4).Returning to the Algorithm 1 we see that the steps 2 and 3 represent just a simulation from the density f ( · ; τ ) of a Wright-Fisher diffusion process and are just a particular case of [JS17, Algorithm 1]. Further-more, there is a more probabilistic description of the coefficients q d − m ( τ ) . Let n ≥ and { A d − n ( τ ); τ ≥ } be a pure death process on the non-negative integers, started at A d − n (0) = n , where the only transitions areof the form m (cid:55)→ m − and occur at a rate m ( m + d − / for each m ∈ { , . . . , n } . Then the coefficients q d − m ( τ ) can be expressed as the limit q d − m ( τ ) = P [ A d − ∞ ( τ ) = m ] = lim n →∞ P [ A d − n ( τ ) = m ] . This inparticular shows that coefficients q d − m ( τ ) are indeed non-negative and sum up to .The coefficients q d − m ( t ) are given only as an infinite series as in (2.5) and not as a closed expression, sothe exact simulation of A d − ∞ ( τ ) is not trivial, but it is possible and is achieved in [JS17]. For our purposes,sufficient accuracy is achieved simply by pre-computing the coefficients up to a certain precision and thenusing standard inversion sampling. N ALGORITHM FOR SIMULATING BROWNIAN INCREMENTS ON A SPHERE 6 Computer simulations
First, we check numerically that our alternative representation (2.4) is actually equivalent to a more well-known representation (2.2). To see that we are going to let τ = 0 . and denote by ρ d the density computedby the first terms in (2.2) and by ρ d density computed by the first terms in (2.4) where each coefficient q d − m ( τ ) is calculated using the first terms in (2.5). All the densities have been calculated in Mathematicafor d = 3 , and on Figure 1 we plot the absolute and relative error. Both errors are extremely small acrossthe interval [0 ,π ] , which confirms that both densities are equal. Additionally, increasing the dimension d seems to improve results. | ρ - ρ || ρ - ρ | × - × - × - × - × - ρ / ρ - ρ / ρ - - × - × - × - Figure 1.
Absolute errors (cid:12)(cid:12) ρ d − ρ d (cid:12)(cid:12) on the left plot and relative error ρ d /ρ d − on the rightplot for d = 3 , .Typically, when using algorithms for the simulation of the spherical Brownian motion, we need to repeat-edly get samples for a fixed dimension d and time τ. Therefore, first step is to pre-compute coefficients q d − m ( τ ) to a prescribed precision, so that we can afterwards use them in an inversion sampling of A d − ∞ ( τ ) . Since thealgorithm depends deeply on the representation (2.2), it suffers from a numerical instability as τ → . Intheory the algorithm works correctly for all values of τ , but in practice we are limited by running time and byimperfect floating number precision, and whence the potential incorrect computations of the relevant coeffi-cients q d − m ( τ ) . When τ is small we have to compute products of the form ( d +2 k − d − m ) ( k − m !( k − m )! e − k ( k + d − τ/ for large values of k and m , which means multiplying very large and very small numbers. This quicklyaccumulates numerical errors which cause the algorithm to fail. For example, using direct calculations donein Mathematica via partial sums and the equation (2.5) we calculate q (0 . to be equal to − . andeven more extremely q (0 . to be equal to . . Both of these results are clearly wrong and willcause the algorithm to fail.In order to understand better how computation of coefficients q d − m ( τ ) can go wrong, we have to look atthe terms b ( τ,d − k ( m ) from equation (2.5). We need to be able to quantify the error in our calculations,meaning we need to know how well the partial sums (cid:80) m + Kk = m ( − k − m b ( τ,d − k ( m ) approximate the coefficients q d − m ( τ ) . The partial sums form an alternating series and one is tempted to say that the partial sums witheven K are always an upper bound for q d − m ( τ ) , whereas the partial sums with odd K are always a lowerbound. Unfortunately, for this to hold, we would need b ( τ,d − k ( m ) ↓ as k → ∞ , which is not the case,although it is true once k is large enough. To see this we notice that b ( τ,d − k +1 ( m ) b ( τ,d − k ( m ) = m + k + d − k − m + 1 2 k + d k + d − e − (2 k + d − τ/ , N ALGORITHM FOR SIMULATING BROWNIAN INCREMENTS ON A SPHERE 7 which is a product of three decreasing functions in k and the limit as k → ∞ is equal to (recall that onlyterms with k ≥ m are relevant). This means that as soon as the quotient is for the first smaller than thenthe terms b ( τ,d − k ( m ) start decaying indefinitely. Therefore, for a large enough number of terms taken, theexact value q d − m ( τ ) always lies between two consecutive partial sums and they differ by just b ( τ,d − k ( m ) ,which is hence also the upper bound for the error. This means that in order to compute coefficients q d − m ( τ ) to a fixed precision ε > we simple start computing each coefficient via partial sums in an infinite seriesrepresentation (2.5) and we stop once the terms b ( τ,d − k ( m ) start decaying and additionally b ( τ,d − k ( m ) < ε holds.As we have seen above, the partial sums (cid:80) m + Kk = m ( − k − m b ( τ,d − k ( m ) eventually become arbitrarily closeto the precise value q d − m ( τ ) ∈ [0 , . However, before reaching this value they can sometimes be very far off,in some cases the difference can be of several orders of magnitude. To see this we plot in Figure 2 some ofthe terms b ( τ,d − k ( m ) and some partial sums, which are used in the computation of q d − m ( τ ) . Σ k = mm + i ( - ) k - m b k ( τ ,d - ) ( m ) b m + i ( τ ,d - ) ( m ) - Σ k = mm + i ( - ) k - m b k ( τ ,d - ) ( m ) b m + i ( τ ,d - ) ( m ) - Figure 2.
Plots of values b ( τ,d − m + i ( m ) and of partial sums (cid:80) m + ik = m ( − k − m b ( τ,d − k ( m ) for i = 0 , . . . , and parameters for the left plot: d = 3 ,τ = 0 . ,m = 2 and for the right plot: d = 3 ,τ = 0 . ,m = 13 .One of the things immediately recognisable on Figure 2 is unimodality of terms b ( τ,d − k ( m ) , which isobvious from calculations in the previous paragraph. Another noticeable thing is the size of the terms.Whereas on the left plot, where time τ is not too small, values of terms b ( τ,d − k ( m ) and partial sums arerelatively small, we can see that on the right plot, where the time τ is smaller, values become really large,in our particular case they are of order . We have to contrast this to the fact, that for the values ofparameters in the right plot the corresponding coefficient q (0 . is approximately equal to . whichis of several orders of magnitudes smaller than the terms (and partial sums) used to calculate it. Whenthe difference between the size of terms and the correct value becomes too large, numerical inaccuraciesquickly add up and numerically computed coefficients q d − m ( τ ) are far from the correct values. This problembecomes even more apparent as the time τ goes to zero. Inspecting coefficients q d − m ( τ ) suggests that thesmallest time for which calculations of coefficients still seem to go through without major numerical errorsis around τ = 0 . and the dimension d doesn’t seem to affect this as much (increasing dimension actuallyextends possible times τ for which the algorithm works correctly).We also need to analyse coefficients q d − m ( τ ) for a fixed time τ and dimension d and see which coefficientsare relevant for the simulation (i.e. which coefficients are not negligible). To get a better understanding wetherefore plot on Figure 3 some of the coefficients for various times and dimensions. N ALGORITHM FOR SIMULATING BROWNIAN INCREMENTS ON A SPHERE 8 τ = τ = τ = τ =
10 20 30 400.050.100.150.200.250.300.35 d = = = =
10 20 30 400.050.100.15
Figure 3.
Plots of values q d − m ( τ ) for m = 0 , . . . , and parameters for the left plot: d = 3 and times τ = 0 . , . , . , . and for the right plot: τ = 0 . and dimensions d = 3 , , , . Again we notice unimodality, which in this case is not directly apparent from the representation (2.5).Additionally, when time τ increases, the distribution moves to the left, which is intuitively clear from adescription via death processes. Similarly, when the dimension d increases, distribution moves to the left,but the effect is less profound.Since we are limited with how small time τ we can take i.e. τ ≥ . , numerical calculations quickly showthat the only relevant coefficients q d − m ( τ ) are those with m ≤ . All other coefficients will be essentiallyequal to zero. Of course for larger values of τ and a larger dimension d there might be even more essentiallyzero coefficients and there is no harm in setting them all equal to 0. Therefore, for a fixed time τ we willonly need to pre-calculate each coefficient q d − m ( τ ) for m ≤ . We do this calculations by taking sufficientlymany terms in representation (2.5) to get the desired accuracy. Additionally, we can stop calculations ifwe reach a coefficient q d − m ( τ ) which is already small enough and late enough so that the coefficients havealready started decreasing. Again, we can set all succeeding coefficients to . Once all relevant coefficientsare calculated, we can then use the standard inversion sampling to simulate a random variable A d − ∞ ( τ ) anduse it in Algorithm 1.As we have seen, our algorithm does not work well for small values of τ . The bottleneck of Algorithm 1is a simulation of A d − ∞ ( τ ) . Hence, if we still want to use a variant of the algorithm, we need to resort tosome kind of an approximation and one option is using the following normal approximation for A d − ∞ ( τ ) ,which was originally given in [Gri84]. Let β = ( d − t/ and η = βe β − . Then A d − ∞ ( τ ) is approximatelydistributed as a normal N (cid:16) µ ( γ,t ) , (cid:0) σ ( γ,t ) (cid:1) (cid:17) random variable, where µ ( γ,t ) = ηt and (cid:0) σ ( γ,t ) (cid:1) = ηt ( η + β ) (cid:16) ηη + β − η (cid:17) β − . Alternatively, one can resort to using the approximation Θ ∼ Q Approx ( θ,τ ) andsetting X = − cos(Θ)2 instead of steps 2 and 3. Such an approximation actually seems to give better resultsthan the normal approximation. Acknowledgements
AM and VM are supported by The Alan Turing Institute under the EPSRC grant EP/N510129/1; AMsupported by EPSRC grant EP/P003818/1 and the Turing Fellowship funded by the Programme on Data-Centric Engineering of Lloyd’s Register Foundation; VM supported by the PhD scholarship of Department ofStatistics, University of Warwick; GUB supported by CoNaCyT grant FC-2016-1946 and UNAM-DGAPA-PAPIIT grant IN114720.
N ALGORITHM FOR SIMULATING BROWNIAN INCREMENTS ON A SPHERE 9
References [Bou16] R. Bouckaert,
Phylogeography by diffusion on a sphere: whole world phylogeography , PeerJ (2016), e2406.[BS98] D.R. Brillinger and B.S. Stewart, Elephant-Seal Movements: Modelling Migration , The Canadian Journal of Statis-tics / La Revue Canadienne de Statistique (1998), no. 3, 431–443.[Cai04] J.-M. Caillol, Random walks on hyperspheres of arbitrary dimensions , Journal of Physics A: Mathematical andGeneral (2004), no. 9, 3077.[CEE10] T. Carlsson, T. Ekholm, and C. Elvingson, Algorithm for generating a Brownian motion on a sphere , Journal ofPhysics A: Mathematical and Theoretical (2010), no. 50, 505001.[Far02] J. Faraudo, Diffusion equation on curved surfaces. i. Theory and application to biological membranes , The Journalof Chemical Physics (2002), no. 13, 5831–5841.[GL83] R. C. Griffiths and W.-H. Li,
Simulating allele frequencies in a population and the genetic differentiation of popu-lations under mutation pressure , Theoretical Population Biology (1983), no. 1, 19 – 33.[Gri84] R. C. Griffiths, Asymptotic line-of-descent distributions , Journal of Mathematical Biology (1984), no. 1, 67–75.[GS10] R. C. Griffiths and D. Spanò, Diffusion processes and coalescent trees , ArXiv e-prints (2010).[GSS12] A. Ghosh, J. Samuel, and S. Sinha,
A "Gaussian" for diffusion on the sphere , EPL (Europhysics Letters) (2012),no. 3, 30003.[Hsu02] E.P. Hsu, Stochastic analysis on manifolds , Graduate Studies in Mathematics 38, vol. 38, American MathematicalSociety, 2002.[JS17] P.A. Jenkins and D. Spanò,
Exact simulation of the Wright–Fisher diffusion , Ann. Appl. Probab. (2017), no. 3,1478–1509.[KDPN00] M. M. G. Krishna, Ranjan Das, N. Periasamy, and Rajaram Nityananda, Translational diffusion of fluorescent probeson a sphere: Monte Carlo simulations, theory, and fluorescence anisotropy experiment , The Journal of ChemicalPhysics (2000), no. 19, 8502–8514.[KT81] S. Karlin and H. M. Taylor,
A second course in stochastic processes , Academic Press, 1981.[LTT08] G. Li, L.-K. Tam, and J. X. Tang,
Amplified effect of Brownian motion in bacterial near-surface swimming , Pro-ceedings of the National Academy of Sciences (2008), no. 47, 18355–18359.[MMU20] A. Mijatović, V. Mramor, and G. Uribe Bravo,
A note on the exact simulation of spherical Brownian motion ,Statistics & Probability Letters (2020), 108836.[NEE03] J. Nissfolk, T. Ekholm, and C. Elvingson,
Brownian dynamics simulations on a hypersphere in 4-space , The Journalof Chemical Physics (2003), no. 13, 6423–6432.
Department of Statistics, University of Warwick, & The Alan Turing Institute, UK
Email address : [email protected] Department of Statistics, University of Warwick, & The Alan Turing Institute, UK
Email address : [email protected] Instituto de Matematicas, Universidad Nacional Autónoma de México, México
Email address ::