Page - 399 - in Differential Geometrical Theory of Statistics

Entropy2016,9, 337 derivativeγ′ij = γ ′ i(ti) is available,most of the time throughanumerical approximation. Givena quadrature formulaon [0,1]withpoints tj, j=1,. . . ,mandweightswj, j=1,. . . ,m, thedensitymay beapproximatedatx∈R2 by: d˜(x)= 1 ∑Ni=1 li N ∑ i=1 m ∑ j=1 wjK (‖x−γij‖)‖γ′ij‖. (28) where the lengths li, i=1,. . . ,Narealsoobtainedwith thesamequadraturerule: li= m ∑ j=1 wj‖γ′ij‖. Whenγ′ij is computedinanumericalway, itmaybeexpressedas: γ′ij= m ∑ k=1 w˜jkγi,k. wheretheweights w˜jkareoftenobtainedthroughtheapplicationoftheLagrangeinterpolationformula toensureexactnessonpolynomialsuptoagivendegree. Inamorecompact form, it canbewritten in matrix formas: ⎛⎜⎜⎝ γ′i1 ... γ′iq ⎞⎟⎟⎠= W˜ ⎛⎜⎝ γi1... γiq ⎞⎟⎠ where thematrix W˜ has as entries theweights w˜jk. The cost of evaluating d˜ at a single point is in o(Nm),with thekernel evaluationbeingdominant. Whendealingwithpoints inR2 orR3 and compactly-supportedkernels, asimple trickgreatlyreduces the timeneededtocompute d˜. Firstof all, thedomainof interest isdiscretizedonanevenly-spacedgrid, so that thepointsof evaluation of the density d˜ are its vertices xij, i = 1,. . . ,nx, j = 1,. . . ,ny. The grid step δx (resp. δy) in the first (resp. second) coordinate is thedifferencebetweenany twoadjacentvertices δx = xi+1,j−xi,j (resp.δy= xi,j+1−xi,j (mostof the time,δx= δy). Since theexpression(28) is linear, thecomputation canbeperformedbyaccumulatingvaluesK(‖xkl−γij‖)‖γ′ij‖ for afixed couple (i, j), where only the points xkl close enough to γij are considered. In fact, the evaluation can bewritten as a 2D discreteconvolution: d˜(xkl)= ∑ i=1,...,N,j=1,...,m wjK(‖xkl−γij‖)‖γ′ij‖. (29) When the support ofK is small compared to the overall spatial domain,much computation is saved using this procedure. Furthermore, it can be thought of as 2D filtering, so that highly efficientalgorithmscomingfromthefieldof imageprocessingcanbeapplied: inparticular, computing thedensityonagraphicsprocessingunit (GPU) is straightforwardandallowsone todecrease the computational timebyat leasta factorof ten.Whendealingwith thescaledarclength, thederivative termisnotpresent,andafactorof li appears in fromof the integral. Thediscreteversionbecomes: d˜(xkl)= ∑ i=1,...,N,j=1,...,m liwjK(‖xkl−γij‖) (30) whereγij=γi(ηj),ηjbeing incorrespondencewith tj. Pleasenote that thequadratureweightsmust be adapted to the abscissa ηj, j= 1,. . . ,m andnot to the tj, j= 1,. . . ,m. Therefore, it is advisable to resample the curves so that thepoints ηj, j= 1,. . . ,m are, for example, evenly spacedorof the Gauss–Lobatto form.Theformerwaschosenfor theexperimentsdueto itseaseof implementation, 399