Density estimation for Compound Cox processes on hyperspheres

07/11/2017
Publication GSI2017
OAI : oai:www.see.asso.fr:17410:22549
contenu protégé  Document accessible sous conditions - vous devez vous connecter ou vous enregistrer pour accéder à ou acquérir ce document.
- Accès libre pour les ayants-droit
 

Résumé

Cox multiple scattering processes on hyperspheres are a class of doubly stochastic Poisson processes that can be used to describe scattering phenomenon in Physics (optics, micro-waves, acoustics, etc.). In this article, we present an EM (Expectation Maximization) technique to estimate the concentration parameter of a Compound Cox process with values on hyperspheres. The proposed algorithm is based on an approximation formula for multiconvolution of von Mises Fisher densities on spheres of any dimension.

Density estimation for Compound Cox processes on hyperspheres

Collection

application/pdf Density estimation for Compound Cox processes on hyperspheres Florent Chatelain, Nicolas Le Bihan, Jonathan H. Manton
Détails de l'article
contenu protégé  Document accessible sous conditions - vous devez vous connecter ou vous enregistrer pour accéder à ou acquérir ce document.
- Accès libre pour les ayants-droit

Density estimation for Compound Cox processes on hyperspheres

Média

Voir la vidéo

Métriques

0
0
287.93 Ko
 application/pdf
bitcache://53b7312af6d38a8fb071f75b1062894a9336f23d

Licence

Creative Commons Aucune (Tous droits réservés)

Sponsors

Sponsors Platine

alanturinginstitutelogo.png
logothales.jpg

Sponsors Bronze

logo_enac-bleuok.jpg
imag150x185_couleur_rvb.jpg

Sponsors scientifique

logo_smf_cmjn.gif

Sponsors

smai.png
logo_gdr-mia.png
gdr_geosto_logo.png
gdr-isis.png
logo-minesparistech.jpg
logo_x.jpeg
springer-logo.png
logo-psl.png

Organisateurs

logo_see.gif
<resource  xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
                xmlns="http://datacite.org/schema/kernel-4"
                xsi:schemaLocation="http://datacite.org/schema/kernel-4 http://schema.datacite.org/meta/kernel-4/metadata.xsd">
        <identifier identifierType="DOI">10.23723/17410/22549</identifier><creators><creator><creatorName>Nicolas Le Bihan</creatorName></creator><creator><creatorName>Jonathan H. Manton</creatorName></creator><creator><creatorName>Florent Chatelain</creatorName></creator></creators><titles>
            <title>Density estimation for Compound Cox processes on hyperspheres</title></titles>
        <publisher>SEE</publisher>
        <publicationYear>2018</publicationYear>
        <resourceType resourceTypeGeneral="Text">Text</resourceType><subjects><subject>Multiple scattering processes</subject><subject>hyperspheres</subject><subject>von Mises Fisher distribution</subject><subject>Compound Cox processes</subject><subject>characteristic function</subject><subject>parametric estimation</subject><subject>Expectation Maximization</subject></subjects><dates>
	    <date dateType="Created">Wed 7 Mar 2018</date>
	    <date dateType="Updated">Wed 7 Mar 2018</date>
            <date dateType="Submitted">Mon 15 Oct 2018</date>
	</dates>
        <alternateIdentifiers>
	    <alternateIdentifier alternateIdentifierType="bitstream">53b7312af6d38a8fb071f75b1062894a9336f23d</alternateIdentifier>
	</alternateIdentifiers>
        <formats>
	    <format>application/pdf</format>
	</formats>
	<version>37271</version>
        <descriptions>
            <description descriptionType="Abstract">Cox multiple scattering processes on hyperspheres are a class of doubly stochastic Poisson processes that can be used to describe scattering phenomenon in Physics (optics, micro-waves, acoustics, etc.). In this article, we present an EM (Expectation Maximization) technique to estimate the concentration parameter of a Compound Cox process with values on hyperspheres. The proposed algorithm is based on an approximation formula for multiconvolution of von Mises Fisher densities on spheres of any dimension.
</description>
        </descriptions>
    </resource>
.

Density estimation for Compound Cox processes on hyperspheres Florent Chatelain1 , Nicolas Le Bihan1 , Jonathan H. Manton2 1. Gipsa-lab (CNRS - UMR 5216), 2. The University of Melbourne, Dept. of Electrical and Electronic Engineering Abstract. Cox multiple scattering processes on hyperspheres are a class of doubly stochastic Poisson processes that can be used to describe scat- tering phenomenon in Physics (optics, micro-waves, acoustics, etc.). In this article, we present an EM (Expectation Maximization) technique to estimate the concentration parameter of a Compound Cox process with values on hyperspheres. The proposed algorithm is based on an approx- imation formula for multiconvolution of von Mises Fisher densities on spheres of any dimension. Keywords: Multiple scattering processes, hyperspheres, von Mises Fisher dis- tribution, Compound Cox processes, characteristic function, parametric estima- tion, Expectation Maximization. 1 Introduction We consider isotropic multiple scattering processes taking values on hyperspheres Sm−1 . The elements of Sm−1 are vectors of unit length in Rm , i.e. Sm−1 =  x ∈ Rm ; ( Pm i=1 x2 i ) 1/2 = 1 with xi, i = 1, . . . , m the components of x. Mul- tiple scattering phenomena occur in many areas of Physics, and the studied model here consists of elastic multiple scattering in any dimension. This means that the description of the wave/particle encountering multiple scattering can be performed by modeling of the direction of propagation only (no energy loss). In such situation, one can consider that the trajectory of one particle is simply a random walk on an hypersphere. In this work, we introduce a Compound Cox process model for multiple scattering in possibly dynamically varying media, and propose an estimation algorithm for the concentration parameter of the scatterers density. This density models the way a heterogeneity scatters the particle/wave, and is based on von Mises Fisher distribution (vMF). The concentration parameter of the vMF is thus a parameter of interest to describe a random medium and we propose an Expectation Maximisation procedure to infer this parameter, given a set of observations at a time t of the ouput vector in Sm−1 of the multiple scattering process. The model used in this work was introduced originally in [1]. The estimation technique here has computational advantages compared to the one proposed in [2] and is an alternate to the one studied in [3]. 2 2 von Mises Fisher random walk on Sm−1 A multiple scattering process is made up of an infinite number of contributions from random walks of finite different lengths. In this section, we give the pdf and characteristic function expressions for a random walk on hypersphere Sm−1 . These results are detailed in [1], and summarized here for use in Section 4. An isotropic random walk on Sm−1 consists of a sequence of vectors starting at x0, with element after k steps given by: xk = Rkxk−1 = RkRk−1 . . . R1x0 (1) where Ri, i = 1, . . . , k is the rotation associated with the ith step. The sequence of unit vectors x0, . . . , xk is isotropic and obeys the Markov property condition. This means that the conditional pdf f(xk|xk−1, . . . , x0) is given by: f(xk|xk−1, . . . , x0) = gk,k−1(xT k−1xk) (2) which is only a function of the cosine of the angle between xk−1 and xk. In directional statistics, the von Mises Fisher (vMF) distribution is amongst the most popular [4], thanks to its similarity with the normal distribution on the real line. Here, we consider the vMF distribution on Sm−1 , denoted Mm(µ, κ), with pdf given by: f(x; µ, κ) = κm/2−1 (2π)m/2 Im/2−1(κ) eκµT x , (3) where Iν (·) is the modified Bessel function [5, p. 374], µ ∈ Sm−1 is the mean direction and κ ≥ 0 is the concentration parameter: the larger the value of κ, the more concentrated is the distribution about the mean direction µ. The vMF random walk with concentration parameter κ is thus defined like: xk|xk−1 ∼ Mm(xk−1, κ), Indeed, this random walk is unimodal with mode µ ≡ x0 as demonstrated in [1]. In order to provide an expression of the pdf of xn, the position of the random walker after n steps on the hypersphere, we first introduce an approximation formula for the multiple convolution of vMF pdf over Sm−1 . 2.1 Convolution of vMF pdfs Consider a vMF isotropic n-step random walk of over Sm−1 , then: fn(xn; µ) = (gn,n−1 ? · · · ? g1,0) (xn) (4) where f(xk|xk−1) = gk,k−1(xT k xk−1) for k = 1, . . . , n, x0 ≡ µ is the initial di- rection and where ? represents the convolution over the double coset SO(m − 1)\SO(m)/SO(m − 1) [1]. The pdf fn(xn; µ) is unimodal (see [1] for proof and details) with mode µ and rotationally invariant with respect to µ. As it is well 3 known in directional statistics [4], the vMF distribution is not stable by convolu- tion. Here, we make use of the approximation introduced in [1, Theorem 4.1] for the multiconvolution of vMF pdf s with identical mean direction µ ∈ Sm−1 and concentration parameter κ. In the high concentration asymptotic case of large κ and small n, i.e. n/κ → 0, a n-step vMF random walk xn is distributed as Mm(µ, κ̃n) with: κ̃n = κ − 1/2 n + 1/2 (5) the equivalent concentration parameter. As a consequence, it is possible to ap- proximate the Fourier series (i.e. the characteristic function) of fn(xn; µ) based on the multiconvolution n-step random walk vMF. First, recall that a vMF pdf f ∈ L1 (SO(m − 1)\SO(m)/SO(m − 1), R) can be written as: f(x; µ, κ) = X `≥0 βm,` ˆ f`(κ)P`(µT x) (6) where the normalisation constant βm,` is given by: βm,` = 1 ωm−1 (2` + m − 2)Γ(` + m − 2) `!Γ(m − 1) for all ` ≥ 0 with ωm−1 = 2 πm/2 Γ (m/2) the area of the (m − 1)-dimensional sphere Sm−1 . The basis elements P`(µT x) are the Legendre polynomials of order ` in dimension m, taken at x and with respect to µ, the symmetry axis of f. The coefficients ˆ f`(κ) are given by [6]: ˆ f`(κ) = E [P`(µT x)] = I`+ν (κ) Iν (κ) (7) with ν = m/2 − 1 for κ > 0 and ` ≥ 0. In the case of a isotropic n-step multiconvolution, the coefficients of the Fourier series of the pdf given in (4), denoted b f⊗n ` can be approximated by: b f⊗n ` = ˜ fn ` + O n κ 3 (8) as n/κ → 0, where ˜ fn ` are the Legendre coefficients of asymptotic distribution Mm(µ, κ̃n) given in (5). Thus the Fourier series of a n-step random walk can be expressed using the asymptotic approximation introduced in this Section. In the sequel, we will make use of the high concentration asymptotic approximation to derive an expression for the pdf of a compound Cox process on Sm−1 . 3 vMF multiple scattering process on Sm−1 The study of multiple scattering processes on manifolds has been originally mo- tivated by their usage in describing the behaviour of waves/particles propagat- ing through a random medium [7, 1]. Here, we consider the case of a multiple 4 scattering process occurring in Rm , where the direction of propagation of the wave/particle is a unit vector x ∈ Sm−1 . Now assuming that during the propa- gation, the wave/particle encounters a random number of scatterers, then the di- rection of propagation after a time t, denoted xt, consists of a mixture of weighted n-steps random walks with n = 0, . . . , ∞ and with weights being the probabil- ity of having n scattering events during the elapsed time t, i.e. P(N(t) = n). The process N(t) is called the counting process. When the time between two scattering events follows an exponential law, xt is a Compound Poisson process with weights equal to e−λt (λt)n /n! and λ the Poisson intensity parameter. The parameter λ can be related to the mean free path of the random medium [7]. In the sequel, we consider the more general case where the Poisson process N(t) is no more homogeneous, i.e. its intensity parameter is a random process Λ(t). N(t) is then called a Cox process [8]. It is well known [9] that it that case, the probability P(N(t) = n) takes the form: P(N(t) = n) = Pn [fΛ(t)] = Z ∞ 0 e−λt (λt)n n! fΛ(t)(λt)dλt (9) where Pn [fΛ(t)] is the Poisson transform of the pdf fΛ(t). It is then possible to give the expression of the pdf of xt in the case where the n-steps random walks are vMFs. Using the high concentration approximation introduced in (8), for a given intensity distribution of Λ(t), in the limit of large κ, one gets: f(xt; µ, κ) ' P0 [fΛ(t)] δµ(xt) + X n≥1 Pn [fΛ(t)] f(xt; µ, κ̃n) (10) where f(· ; µ, κ̃n) is the vMF pdf of the n-step random walk, i.e. Mm(µ, κ̃n)1 . Note that this expression is valid if there exists a > 0 such that E[Λ(t)a ] < +∞. In practical applications, it is required that only a small number of scattering events that weakly pertubate the direction of propagation occured, for the random medium to be identified. The high concentration approximation used here fits exactly this framework. The high concentration approximation in (10) will be used in Section 4 for estimation purpose. 3.1 Compound Cox process with Gamma intensity distribution Thus for a given t > 0, the intensity of the Poisson process is now assumed to be Gamma distributed, i.e. Λ(t) ∼ G(rt, p) where rt > 0 is a fixed and known shape parameter and p > 0 is an unknown scale parameter. This means that the density of the mixing process reads fΛ(t)(x) = xrt−1 Γ(rt)prt e− x p , (11) 1 Here we make use of a notation abuse by expressing f(xt; µ, κ) as a sum of dirac measure and a pdf. It has to be understood the following way: when x = µ, it equals the dirac mass, and for other cases it equals the density function. 5 for all x ≥ 0. In this case the counting process N(t) for the scattering events obeys a negative binomial NB(rt, q) distribution with stopping-time rt and suc- cess probability q = p/(p + 1): Pr (N(t) = n) = Pn [fΛt ] = Γ(n + rt) n!Γ(rt) pk (p + 1)n+rt , (12) for all n ∈ N. In this case the high concentration distribution given in (10) is still valid and the Poisson transform masses are given in (12). 4 Estimation The distribution of the multiple scattering process yields a likelihood function which is a product of Fourier series. This function is quite expensive to calcu- late. In such situations, approximate Bayesian computation (ABC) is a popular likelihood-free method that can be used to perform Bayesian inference, as shown in [2], at a price of a high computational cost. However the asymptotic expression (10) allows us to model the multiple scattering process as a (infinite) mixture model of vMF distributions. As a con- sequence, it becomes quite straightforward to apply classical procedures to esti- mate the parameters of this multiple scattering process. Such procedures include expectation-maximization (EM) algorithm, or Gibbs sampling in a Bayesian framework. One benefit of these standard estimation methods is the possibility to reduce significantly the computational burden w.r.t. likelihood-free methods such as ABC. 4.1 EM algorithm In this work, we derive an EM algorithm in order to estimate the parameters of vMF compound Cox process when the intensity of the Poisson process is Gamma distributed. Observations are assumed to be made at a given time t. Let x = (x1, . . . xN ) be a sample of N independent observations in Sm−1 from a vMF multiple scattering process with known initial direction µ and whose in- tensity Λ(t) is Gamma distributed as defined in (11). The vector of the unknown parameter of the multiple scattering process to be estimated is θ = (p, κ), where p is related to the intensity of the Poisson process, and κ is the concentration parameter of the vMF distribution given in (3). In order to derive the EM equa- tions, we can introduce the following binary latent variables for all 1 ≤ i ≤ N, n ≥ 0, zi,n = ( 1 if N(t) = n for the sample xi, 0 otherwise. 6 Then the log-likelihood of the complete data (x, z), where z is the set of all the zin, can be (approximately) computed based on the high concentration distribu- tion (10) as `(θ; x, z) ' N X i=1 X n≥0 zin log [Pn [fΛt ] f(xi; µ, κ̃n)], with by convention f(xi; µ, κ̃0) = δµ(xi). E step. Given x and a current estimate of the parameter vector θ(α) =  p (α) , κ(α)  , the conditional expected value of the log-likelihood reads Q θ|θ(α)  = N X i=1 X n≥0 ti,n log [Pn [fΛt ] f(xi; µ, κ̃n)], (13) where ti,n = E  zi,n|x, θ(α)  = Pr zi,n = 1|x, θ(α)  for all 1 ≤ i ≤ N, n ≥ 0. Using Bayes rule, it comes that ti,n = Pn [fΛt ] f(xi; µ, κ̃(α) n ) P l≥0 Pl [fΛt ] f(xi; µ, κ̃(α) l ) , (14) where κ̃(α) n = κ(α)−1/2 n +1/2 as defined in (5) and Pn [fΛt ] is obtained by replacing the scale parameter p by its current value p(α) in (12). M step. In order to maximize Q θ|θ(α)  , it is useful to note that this objec- tive function is separable w.r.t. the two parameters p and κ and can thus be maximized independently. The scale parameter p is obtained by maximizing p 7→ PN i=1 P n≥0 ti,n log Pn [fΛt ], where the negative binomial masses Pn [fΛt ] are given in (12). This has the same form as a weighted maximum likelihood estimator for a negative binomial distribution. Moreover the negative binomial distribution NB(rt, p) with fixed stopping-time parameter rt describes a natural exponential family with mean parameter equal to rtp. Thus it comes directly that p(α+1) = 1 rt 1 N N X i=1 X n≥0 nti,n. (15) The concentration parameter κ is obtained by maximizing κ 7→ Q θ|θ(α)  = N X i=1 X n≥1 ti,n (log cm (κ̃n) + κ̃nµT xi) + constant, where cm(κ) = κm/2−1 (2π)m/2Im/2−1(κ) is the vMF normalizing constant, and κ̃n depends on κ as defined in (5). By differentiating this function w.r.t. κ, we obtain the 7 following score function v(κ) = N X i=1 X n≥1 ti,n −Am(κ̃n) + µT xi n , where Am(κ) = −c0 m(κ) cm(κ) = Im/2(κ) Im/2−1(κ) . For a fixed direction µ the distribution of the scalar µT y, where y is vMF distributed as Mm(µ, κ), describes a natural exponential family on the set κ > 0. Thus the associated log-likelihood is strictly concave w.r.t. κ and v(κ) is monotonically decreasing as a positively weighted sum of monotonically decreasing function. This shows that the value that max- imizes κ 7→ Q θ|θ(α)  is the unique zero of v(κ). Unfortunately there is no tractable closed-form expression of this zero. However, similarly to maximum likelihood procedures for vMF distributions [10], it is possible to perform one Newton’s iteration to improve the objective criterion Q θ|θ(α)  . Indeed such im- provement ensures, as a standard property of EM algorithm, that the marginal likelihood of the observations x is improved. This yields the following update rule κ(α+1) = κ(α) − v(κ(α) ) v0 (κ(α) ) , (16) where v0 (κ) = − P n≥1 A0 m(κ̃n) n2 PN i=1 ti,n is the derivative of the score function v(κ). Moreover A0 m(κ) can be efficiently computed using A0 m(κ) = 1 − Am(κ)2 − m−1 κ Am(κ) as shown in [4]. Note finally that the EM algorithm can be initialized by using method of moment estimates as defined in [3]. 4.2 Simulation results Several simulations have been conducted on synthetic data to evaluate the sta- tistical performances of the estimates given by the EM procedure. Note that in practice, to compute the EM equations and thus to implement the EM algorithm, we can truncate the infinite series in (14), (15) and (16) to compute the probabilities ti,n and the updated estimates. This can be done by introducing a tolerance parameter  > 0 and the associated index upper bound L = arg minL∈N PL l=0 Pl [fΛt ] > 1 − . Thus it comes that ti,n ≈ Pn [fΛt ] f(xi; µ, κ̃(α) n ) PL l=0 Pl [fΛt ] f(xi; µ, κ̃(α) l ) , (17) for 0 ≤ n ≤ L, and ti,n ≈ 0 for n > L. The biases and the standard deviations of EM estimates are reported in Table 1. These numerical results underline the accuracy of the EM estimators based on the high concentration approximation even when the ratio κ/p is not very large, e.g. κ/p ≥ 2. 8 Table 1. Performances of EM estimates (p̂, κ̂) as a function of the true parameters p and κ, from 1000 Monte-Carlo runs, for a vMF multiple scattering process with G(rt, p) mixing distribution (rt = 1). Sample size is N = 1000, and samples belong to S2 (m = 3). p 1 2 3 5 10 b p b κ b p b κ b p b κ b p b κ b p b κ κ = 20 bias 0.011 0.275 0.025 0.447 0.043 0.695 0.141 1.386 0.330 2.996 stdev 0.059 1.288 0.121 1.281 0.195 1.344 0.327 1.313 0.752 1.496 κ = 50 bias 0.007 0.425 0.022 0.624 0.042 0.975 0.151 1.914 0.437 3.433 stdev 0.059 3.011 0.125 3.063 0.193 3.241 0.337 3.564 0.752 3.966 κ = 100 bias 0.007 0.502 0.020 0.898 0.041 1.168 0.129 2.743 0.446 5.097 stdev 0.060 6.297 0.124 6.198 0.190 6.660 0.336 6.971 0.756 8.132 κ = 200 bias 0.007 1.477 0.021 2.178 0.042 2.792 0.140 5.519 0.484 9.486 stdev 0.058 12.28 0.123 12.55 0.199 13.54 0.319 13.80 0.749 15.58 References 1. Le Bihan, N., Chatelain, F., Manton, J.: Isotropic multiple scattering processes on hyperspheres. IEEE Transactions on Information Theory 62 (2016) 5740–5752 2. Chatelain, F., Le Bihan, N., Manton, J.: Parameter estimation for multiple scat- tering process on the sphere. In: IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). (2015) 3. Chatelain, F., Le Bihan, N.: von-mises fisher approximation of multiple scatter- ing process on the hypershpere. In: IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). (2013) 4. Mardia, K., Jupp, P.: Directional statistics. John Wiley & Sons Ltd (2000) 5. Abramowitz, M., Stegun, I.A.: Handbook of Mathematical Functions with Formu- las, Graphs, and Mathematical Tables. Dover Publications (1972) 6. Kent, J.: Limiting behaviour of the von Mises-Fisher distribution. Math. Proc. Camb. Phil. Soc. 84 (1978) 531–536 7. Le Bihan, N., Margerin, L.: Nonparametric estimation of the heterogeneity of a random medium using compound poisson process modeling of wave multiple scattering. Physical Review E 80 (2009) 016601 8. Lefebvre, M.: Applied stochastic processes. Springer (2006) 9. Saleh, B.: Photoelectron statistics. Springer (1978) 10. Sra, S.: A short note on parameter approximation for von mises-fisher distributions. Comput Stat 27 (2012) 177–190