Francesco Fedele
email: fedele@gatech.edu

Denys Dutykh
email: denys.dutykh@univ-savoie.fr

SPECIAL SOLUTIONS TO A COMPACT EQUATION FOR DEEP-WATER GRAVITY WAVES

Keywords: water waves, deep water approximation, Hamiltonian structure, travelling waves, solitons

have derived a compact form of the well known Zakharov integro-differential equation for the third order Hamiltonian dynamics of a potential flow of an incompressible, infinitely deep fluid with a free surface. Special traveling wave solutions of this compact equation are numerically constructed using the Petviashvili method. Their stability properties are also investigated. In particular, unstable traveling waves with wedge-type singularities, viz. peakons, are numerically discovered. To gain insights into the properties of these singular solutions, we also consider the academic case of a perturbed version of the compact equation, for which analytical peakons with exponential shape are derived. Finally, by means of an accurate Fouriertype spectral scheme it is found that smooth solitary waves appear to collide elastically, suggesting the integrability of the Zakharov equation.

Introduction

The Euler equations that describe the irrotational flow of an ideal incompressible fluid of infinite depth with a free surface are Hamiltonian. Their symplectic formulation was discovered by [START_REF] Zakharov | Stability of periodic waves of finite amplitude on the surface of a deep fluid[END_REF] [START_REF] Zakharov | Stability of periodic waves of finite amplitude on the surface of a deep fluid[END_REF] in terms of the free-surface elevation η(x, t) and the velocity potential ϕ(x, t) = φ(x, z = η(x, t), t) evaluated at the free surface of the fluid. Here, η(x, t) and ϕ(x, t) are conjugated canonical variables with respect to the Hamiltonian H given by the total wave energy. It is well known that the Euler quations are completely integrable in several important limiting cases. For example, in a two-dimensional (2-D) ideal fluid, unidirectional weakly nonlinear narrowband wave trains are governed by the Nonlinear Schrödinger (NLS) equation, which is integrable [START_REF] Zakharov | Exact Theory of Two-dimensional Self-focusing and Onedimensional Self-modulation of Waves in Nonlinear Media[END_REF]. Integrability also holds for certain equations that models long waves in shallow waters, in particular the Korteweg-de Vries (KdV) equation (see, for example, [START_REF] Ablowitz | The inverse scattering transform-Fourier analysis for nonlinear problems[END_REF]2,[START_REF] Drazin | Solitons: An introduction[END_REF][START_REF] Whitham | Linear and nonlinear waves[END_REF]) or the Camassa-Holm (CH) equation [4]. For these equations, the associated Lax-pairs have been discovered and the Inverse Scattering Transform unveiled the dynamics of solitons, which elastically interact under the invariance of an infinite number of time-conserving quantities [START_REF] Ablowitz | The inverse scattering transform-Fourier analysis for nonlinear problems[END_REF]2,[START_REF] Drazin | Solitons: An introduction[END_REF][START_REF] Whitham | Linear and nonlinear waves[END_REF].

Another important limiting case of the Euler equations for the free-surface of an ideal flow is that considered by [START_REF] Zakharov | Stability of periodic waves of finite amplitude on the surface of a deep fluid[END_REF][START_REF] Zakharov | Statistical theory of gravity and capillary waves on the surface of a finite-depth fluid[END_REF]. By means of a third order expansion of H in the wave steepness, he derived an integro-differential equation in terms of canonical conjugate Fourier amplitudes, which has no restrictions on the spectral bandwidth. To derive the Zakharov (Z) equation, fast non-resonant interactions are eliminated via a canonical transformation that preserves the Hamiltonian structure [17,[START_REF] Zakharov | Statistical theory of gravity and capillary waves on the surface of a finite-depth fluid[END_REF][START_REF] Gramstad | Hamiltonian form of the modified nonlinear Schrödinger equation for gravity waves on arbitrary depth[END_REF]. The integrability of the Z equation is still an open question, but the fully nonlinear Euler equations are nonintegrable [7,8]. Indeed, non-integrability can be easily proven by considering the terms of the perturbation series of the Hamiltonian in powers of the wave steepness limited on their resonant manifolds. To have non integrability, it is enough to prove that at least one of these amplitudes is nonzero. In this regard, [7,8] conjectured that the Z equation for unidirectional water waves (2-D) is integrable since the nonlinear fourth-order term of the Hamiltonian vanishes on the resonant manifold, leaving only trivial wave-wave interactions, which just cause nonlinear frequency shifts of the separate Fourier modes. They also pushed their analysis to the next order and proved that the effective fifth-order term does not vanish on the corresponding resonant manifold. Thus, the Euler equations of the free-surface hydrodynamics are not integrable in two dimensions [8].

Recently, [START_REF] Dyachenko | Compact Equation for Gravity Waves on Deep Water[END_REF] [9] realized that the trivial resonant quartetinteractions that occur in the 2-D free surface dynamics could be further removed by a canonical transformation. As a result, the Z equation drastically simplifies to the compact form

ib t = Ωb + i 8 b * (b x ) 2 x -b * x (b 2 ) x x - 1 4 bK{|b x | 2 } -b x K{|b| 2 } x , (1.1
η = ω 0 2g a + c.c., (1.2)
where a is a conjugate canonical variable that relates to b via a non-trivial canonical transformation (see [9]), and c.c. denotes complex conjugation. To leading order in the wave steepness a = b.

In this study, we explore (1.1), hereafter referred to as cDZ, for analytical studies and numerical investigations of special solutions in the form of solitary waves. The present paper is organized as follows. First, we derive the envelope equation associated to cDZ and then smooth ground states and traveling waves are numerically computed using the Petviashvili method ([26], see also [18]). Further, unstable traveling waves with wedgetype singularities, viz. peakons, are also numerically identified. To gain insights into the properties of singular traveling waves we consider a perturbed version of the local compact equation, where the non-local terms involving the operator K are neglected in (1.1). For this academic case we are able to derive analytical peakons with exponential shape. Finally, the interaction of smooth traveling waves is numerically investigated by means of an accurate Fourier-type pseudo-spectral scheme.

Envelope dynamics and special solutions

Consider the following ansatz for wave trains in deep water

b(X, T ) = ε 2g ω 0 a 0 B(X, T )e i(k 0 x-ω 0 t) , (2.1)
where B is the envelope of the carrier wave e i(k 0 x-ω 0 t) , and X = εk 0 (xc g t), T = ε 2 ω 0 t, with k 0 = ω 2 0 g and ω 0 as characteristic wavenumber and frequencies. The small parameter ε = k 0 a 0 is a characteristic wave steepness associated to an amplitude a 0 and c g is the wave group velocity in deep water. From (1.2), the leading order wave surface η is given in terms of the envelope B as η(X, T ) = εa 0 B(X, T )e i(k 0 x-ω 0 t) + c.c.

(2.2) Using (2.1) in (1.1) yields the cDZ envelope equation

iB T = Ω ε B + i 4 B * S((SB) 2 ) + iB * (SB) 2 -2S B|SB| 2 - ε 2 BK{|SB| 2 } -S SBK{|B| 2 } , (2.3)
where S = ε∂ X + i. The approximate dispersion operator Ω ε is defined as follows

Ω ε := 1 8 ∂ XX + i 16 ε∂ XXX - 5 128 ε 2 ∂ XXXX + 7i 256 ε 3 ∂ XXXXX ,
where o(ε 3 ) dispersion terms are neglected. Equation (2.3) admits three invariants, viz. the action A, momentum M and the Hamiltonian H given, respectively, by

H = R B * Ω ε B + i 4 |SB| 2 [B(SB) * -B * SB] - ε 2 |SB| 2 K(|B| 2 ) dX, (2.4) A = R B * B dX, M = R i B * SB -B(SB) * dX. (2.5)
Expanding the operator S in terms of ε, (2.3) can be written in the form of a generalized derivative NLS equation as

iB T = Ω ε B + |B| 2 B -3iε|B| 2 B X - ε 2 BK{|B| 2 } + ε 2 N 2 (B) + ε 3 N 3 (B) = 0, (2.6)
where B X denotes differentiation with respect to X,

N 2 (B) = - 3 2 B * (B X ) 2 + B|B X | 2 -|B| 2 B XX + 1 2 B 2 B * XX + i 2 BK(B * B X -BB * X ) + B X K|B| 2 + BK|B| 2 X , and
N 3 (B) = - i 2 |B X | 2 B X + i 2 B XX (B * B X -BB * X )- i 2 BB X B * XX - 1 2 BK|B X | 2 -B X K|B| 2 X .
To leading order (2.6) reduces to the NLS equation.

2.1.

A Hamiltonian Dysthe equation. Hereafter, we will study a special case of (2.6), which is a symplectic version of the temporal Dysthe equation ( [START_REF] Dysthe | Note on a modification to the nonlinear Schrödinger equation for application to deep water[END_REF], see also [31,[START_REF] Stiassnie | Note on the modified nonlinear Schrödinger equation for deep water waves[END_REF]16,[START_REF] Trulsen | Frequency downshift in three-dimensional wave trains in a deep basin[END_REF]). Keeping terms to O(ε) in (2.6) yields

iB T = 1 8 ∂ XX + iε 16 ∂ XXX B + |B| 2 B -3iε|B| 2 B X - 1 2 εBK|B| 2 , (2.7)
hereafter referred to as cDZ-Dysthe. Note that the original temporal Dysthe equation is not Hamiltonian since expressed in terms of multiscale variables, which are usually non canonical (see, for example, [START_REF] Fedele | Hamiltonian form and solitary waves of the spatial Dysthe equations[END_REF]). This is common in mathematical physics. Typically, the dynamics is governed by partial differential equations expressed in terms of physicallybased variables, which are not usually canonical. A transformation to new variables is needed in order to unveil the desired structure explicitly (see, for example, [27]). This is the case for the equations of motion for an ideal fluid: in the Eulerian description, they cannot be recast in a canonical form, whereas in a Lagrangian frame the Hamiltonian structure is revealed by Clebsch potentials (see, for example, [27,[START_REF] Morrison | Hamiltonian description of the ideal fluid[END_REF]). Moreover, multiplescale perturbations of differential equations expressed in terms of non-canonical variables typically lead to approximate equations that do not maintain the fundamental conserved quantities, as the hydrostatic primitive equations on the sphere [22], where energy and angular momentum conservation are lost under the hydrostatic approximation.

As an example, consider the finite dimensional system of an harmonic oscillator in the classical canonical variables q(t) (coordinate) and p(t) (momentum). This admits the canonical form

q = ∂H ∂p = p, ṗ = - ∂H ∂q = -q, (2.8)
where ṗ denotes time derivative, and H = (q 2 + p 2 )/2. The flow in the phase-space is 'incompressible' since the divergence vanishes:

∂ q ∂q + ∂ ṗ ∂p = 0.
The transformation z = q + ip is canonical and (2.8) transforms to

ż = -i ∂H ∂z * = -iz, ż * = i ∂H ∂z = iz * , (2.9)
where H = |z| 2 /2, and z * is the complex conjugate of z. It is straightforward to prove that the gauge transformation z = we iα|w| 2 , with α as a free parameter, is also canonical, and (2.9) remains unchanged in the new variables w and w * . On the other hand, if one considers the coordinate change

Q = q 1 + αq 2 , P = p, (2.10)
then (2.8) transforms to the noncanonical form

Ṗ = -Q √ 1 + αP 2 -2αP 2 , Q = P √ 1 + αP 2 . (2.11)
This flow does not preserve volume as (2.8) does, nonetheless equations (2.11) are those of a disguised harmonic oscillator obtained via the noncanonical change of variables (2.10) (see also [START_REF] Morrison | Hamiltonian description of the ideal fluid[END_REF]). The original Dysthe equation [START_REF] Dysthe | Note on a modification to the nonlinear Schrödinger equation for application to deep water[END_REF] shares the same roots as (2.11). They both come from a noncanonical transformation of a Hamiltonian system. For the spatial version of the Dysthe equation, canonical variables have been identified by means of a gauge transformation, and the hidden Hamiltonian structure is unveiled [START_REF] Fedele | Hamiltonian form and solitary waves of the spatial Dysthe equations[END_REF]. The more general canonical transformation proposed by [9] yields the symplectic form (2.7) of the Dysthe equation [START_REF] Dysthe | Note on a modification to the nonlinear Schrödinger equation for application to deep water[END_REF].

In the following, insights into the underlying dynamics of the cDZ equation and associated Dysthe equation are to be gained if we construct some special families of solutions in the form of ground states and traveling waves of the envelope B, often just called solitons or solitary waves, as described below.

2.2. Ground states and traveling waves. To begin, consider the cDZ equation (2.3). We construct numerically ground states and traveling waves (TWs) of the envelope B of the form B(X, T ) = F (X -cT )e -iωT , where the function F (•) is in general complex, ω and c are dimensionless frequency and velocity of the TW with respect to a reference frame moving with the group velocity c g . In physical space the true frequency and velocity are given, respectively, by ω = ε 2 ω 0 ω and c = εcω 0 /k 0 . After substituting this ansatz into the governing equation (2.3) we obtain the following nonlinear steady problem for F (in the moving frame X -cT ) LF = N (F ), where L = ωic -Ω ε and N (F ) denotes nonlinear terms of the cDZ of (2.3). In order to solve this equation we use the Petviashvili method [26, [START_REF] Pelinovsky | Convergence of Petviashvili's iteration method for numerical approximation of stationary solutions of nonlinear wave equations[END_REF]18,[START_REF] Yang | Nonlinear Waves in Integrable and Nonintegrable Systems[END_REF]. This numerical approach has been successfully applied in deriving TWs of the spatial version of the Dysthe equation [START_REF] Fedele | Hamiltonian form and solitary waves of the spatial Dysthe equations[END_REF]. Schematically, the iteration takes the following form

F n+1 = S γ L -1 • N (F n ), S = F n , L • F n F n , N (F n ) ,
where S is the so-called stabilyzing factor and the exponent γ is usually defined as a function of the degree of nonlinearity p (p = 3 for the cDZ equation). The rule of thumb prescribes the following formula γ = p p-1 . The scalar product is defined in the L 2 space.

The inverse operator L -1 can be efficiently computed in the Fourier space. To initialize the iterative process, one can use the analytical solution to the leading order NLS equation (see the associated Dysthe equation (2.7)). We point out that this method can be very efficiently implemented using the Fast Fourier Transform (FFT) (see, for example, [START_REF] Frigo | The Design and Implementation of FFTW3[END_REF]). Without loosing generality, hereafter we just consider the leading term of the dispersion operator, viz. Ω ε = 1 8 ∂ XX , since the soliton shape is only marginally sensitive to the higher order dispersion terms as shown in Figure 1. We also observed that the Petviashvili method does not converge if we retain the full dispersion operator. Figure 2 shows the action A of the ground states (c = 0, but moving with the group velocity in the physical frame) as a function of the frequency ω computed via the Petviashvili iteration and numerical continuation for ε = 0.20. The stability of a ground state is investigated numerically by means of an accurate Fourier-type spectral scheme [START_REF] Boyd | Chebyshev and Fourier Spectral Methods[END_REF]32], see also [START_REF] Fedele | Hamiltonian form and solitary waves of the spatial Dysthe equations[END_REF]. We found that smooth ground states are stable in agreement with the criterion formulated by [34], since dA dω > 0 (see also [START_REF] Zakharov | Wave turbulence in one-dimensional models[END_REF][START_REF] Yang | Nonlinear Waves in Integrable and Nonintegrable Systems[END_REF]). Further, we notice that an abrupt reduction in the action A occurs at a critical frequency ω c (ε) and solitons with wedge-type singularities, viz. peakons, bifurcate from a smooth solitary waves as shown in the right panel of Figure 2. Clearly, this bifurcation can be interpreted as a possible indication of the non-existence of smooth solitons above the critical threshold ω c . As one can see, as ω increases the soliton shape tends to become asymmetric and steeper until the critical threshold ω c = 0.85, above which the smooth solitary waves bifurcate to peakons, which are unstable in agreement with [34], since dA dω < 0. Peakons also bifurcate from smooth solitons of given frequency ω as the steepness ε increases, as clearly seen in Figure 3. Note that in both cases ground states grow asymmetrically before bifurcating to smaller peakons, and so do travelling waves (c = 0). However, such bifurcation is not observed in the 2-D NLS or Dysthe equations (see, for example, [START_REF] Fedele | Hamiltonian form and solitary waves of the spatial Dysthe equations[END_REF]), as clearly seen in Figures 567where we report the dependence of the three invariants A, M and H on the parameters ω, c and ε for the cDZ and cDZ-Dysthe equations respectively.

Note that stable and elastic peakons have been discovered in a special limit of the integrable CH equation [4]. On the contrary, from numerical investigations of the cDZ equation peakons appear unstable. The left panel of Figure 4 shows the numerically converged peakon obtained via the Petviashvili scheme using N ∼ 1.5 × 10 6 Fourier modes for ε = 0.20.

As the frequency ω increases, peakons bifurcate at smaller steepness ε as clearly seen in the right panel of Figure 4. Indeed, for ground states we observed numerically that ω c (ε) ∼ ε -2 (see left panel of Figure 8), and the same scaling holds also for travelling waves (c = 0). In the physical frame (x, t), owing to the scaling T = ε 2 ω 0 t, bifurcation occurs at a well defined frequency ω = ε 2 ω c ≈ 0.037ω 0 , which is independent of ε. Further, in the domain X both the peakon amplitude a p = |B|| X=0 -and the associated slope of the wedge singularity (corner)

s p = 1 a p ∂ X |B|| X=0 -= cot(θ/2)
scales as ε -3/2 , θ being the interior angle of the peakon's singularity (see central and right panels of Figure 7). In the physical frame x, from (2.2) the slope s p becomes

s p = εa 0 ∂ X |B|| X=0 -∂ x X = ε 3 a p s p ≈ cot(170 • /2),
which is independent of ε and corresponds to an angle θ ≈ 170 • . For a meaningful comparison of this result with the almost-highest solitary wave theory of Longuet-Higgins and Fox [20, 21, 19], which predicts an angle of θ ≈ 120 • , higher order corrections to (2.2) hidden within the full canonical transformation of [9] should be accounted for. This is currently under study and will be discussed elsewhere. However, it is also important to mention that the cDZ equation is derived assuming weak nonlinearities, so the strongly nonlinear peakon could be a non-physical artifact, whereas the almost-highest wave satisfies the full Euler equations and thus accounts for nonlinearities of all orders. The derivation of an analytical solution for peakons of the non-local cDZ equation (2.6) is a challenge. To gain more insights into the properties of these singular solutions we consider the academic case of a perturbed version of the local cDZ, viz.

iB T = 1 8 B XX + i 4 B * S((SB) 2 ) + (1 + δ)iB * (SB) 2 -2S B|SB| 2 , (2.12)
Here, the non-local terms are dropped simply because peakons with exponential shape can be found analytically as function of the free parameter δ. Indeed, the ansatz

B(X, T ) = a p e i(κX-ωpT ) e -α 2 |X-cpT | , (2.13)
satisfies (2.12) for

κ = - 8 + δ 8ε , α 2 = |δ| 8ε , ω p = - 4 + δ 32ε 2 , c p = 8 + δ 32ε . (2.14)
The amplitude a p is still unknown and it rules the existence of peakons. Indeed, they arise as a balance between the dispersion 1 8 B XX and the O(ε 2 ) term of the nonlinearity i 4 B * S((SB) 2 ) in (2.12). These two terms are interpreted in distributional sense because they give rise to Dirac delta functions that must vanish by properly chosing the amplitude a p , thus satisfying the differential equation (2.12) in the sense of distributions. As a result

a p = 1 ε - 2 δ .
(2.15)

Clearly, peakons exist for δ < 0 and for δ = 0 they reduce to periodic wave solutions of the local cDZ. Thus, peakons can bifurcate from periodic waves of the local cDZ as δ changes sign, and ground states occur for δ = -8. We note that both the NLS and Dysthe equations do not have higher nonlinear terms similar to i 4 B * S((SB) 2 ) and thus they cannot support peakons. We conclude that wedge-type singular solutions are special features of the cDZ equation. The left panel of Figure 9 shows the remarkable agreement of a numerical peakon (δ = -3/2 and ε = 0.3) computed using the Petviashvili scheme (Fourier modes N ∼ 1.5 × 10 6 ) and the corresponding analytical solution (2.13). Moreover, the amplitude of the numerical peakon agrees with the associated theoretical counterpart (2.15) as function of ε (see right panel of Figure 9). Finally, note that the analytical peakon (2.13) has a slope that scales as ε -1 , in contrast to the ε -3/2 scaling observed numerically for the full cDZ singular solution. In the physical frame x, the interior angle of the wedge singularity is independent of ε and so is the peakon frequency ω p , in agreement with the numerical cDZ peakons. This is an indirect evidence of the validity of our numerical results for the cDZ equation. Finally, we note that as ε approaches zero, (2.12) tends to the NLS equation, and all the peakon parameters approach infinite values as an indication of the non-existence of singular solutions in the NLS limit.

Traveling wave interactions

Hereafter, we investigate the collision of smooth traveling waves of the cDZ equation (2.3) by means of a highly accurate Fourier-type pseudo-spectral method. We will first briefly describe the adopted spectral approach and then discuss the numerical results. 3.1. Numerical method description. We rewrite (2.3) in the operator form

B t + L • B = N (B), (3.1)
where L = i 1 8 ∂ XX and

N (B) = 1 4 B * S((SB) 2 ) + iB * (SB) 2 -2S B|SB| 2 + εi 2 BK{|SB| 2 } -S SBK{|B| 2 } .
We solve (3.1) numerically by applying the Fourier transform in the space variable X.

The transformed functions will be denoted by B = F {B}. We recall that the symbol of the non-local term is equal to |k|, and that of L is i 1 8 k 2 , k being the Fourier transform parameter. The nonlinear terms are computed in physical space, while space derivatives ∂ X and the nonlocal operator K = H[∂ x ] are computed in Fourier space. For example, the

F {B 2 B X } = F { F -1 ( B) 2 • F -1 {ik B}},
and all nonlinear terms are treated in a similar way. We note that the usual 4/3 rule is applied for anti-aliasing since we have to deal with cubic nonlinearities [32,[START_REF] Clamond | A fast method for fully nonlinear water-wave computations[END_REF]14].

In order to improve the stability of the space discretization procedure, we can integrate exactly the linear terms. This is achieved by making a change of variables [23,14]: Finally, the resulting system of ODEs is discretized in space by the Verner's embedded adaptive 9(8) Runge-Kutta scheme [35]. The step size is chosen adaptively using the socalled H211b digital filter [START_REF] Dyachenko | Toward an integrable model of deep water[END_REF]29] to meet the prescribed error tolerance, set as of the order of machine precision.

Ŵt = e (t-t 0 )L • N e -(t-t 0 )L • Ŵ , Ŵ (t) ≡ e (t-t 0 )L • B(t), Ŵ (t 0 ) = B(t 0 ).

Numerical results.

In all the performed simulations the accuracy has been checked by following the evolution of invariants (2.4), (2.5). From a numerical point of view the cDZ equation becomes gradually stiffer as the steepness parameter ε increases, or if higher order dispersion terms are included. Consequently, the conservation of invariants H, A and M might be degraded. Nevertheless, even in the stiffest situations a decent accuracy was assured by choosing a sufficiently large number of Fourier modes and the dispersion operator Ω ε = 1 8 ∂ XX . For example, for ε = 0.10 the number N = 16384 of Fourier modes were sufficient to achieve conservation of the invariants close to ∼ 10 -13 . As an application, consider the interaction between two solitary waves traveling in opposing directions with the same speed (ε = 0.10). The plot of Figure 10 shows that the two solitons emerge out of the collision with the same shape, but a slight phase shift. The interaction appears elastic as clearly seen from the plot of Figure 11, which reports the initial and final shapes of the two solitary waves. Further, the interaction of two traveling waves with a ground state appears also elastic (see . This suggests the integrability of the cDZ equation (2.3) in agreement with the recent results of [START_REF] Dyachenko | Collision of two breathers at surface of deep water[END_REF]. However, the associated Hamiltonian version of the Dysthe equation (2.7) does not support elastic collisions as shown in Figure 14.

Conclusions

Special solutions of the compact Zakharov equation in the form of traveling waves are numerically constructed using the Petviashvili method. The stability of ground states agrees with the Vakhitov-Kolokolov criterion. Further, unstable ground states with wedgetype singularities, viz. peakons, are numerically identified bifurcating from smooth ground states. As an academic case, we considered a perturbed version of the local form of the cDZ equation, for which an analytical solution of peakons with exponential shape is derived. Finally, by means of an accurate Fourier-type pseudo-spectral scheme, it is also shown that smooth solitary waves appear to collide elastically, suggesting the integrability of the compact Zakharov equation, but not that of the associated Hamiltonian version of the Dysthe equation.

  ) where b is a canonical variable, b t , b x denote partial derivatives with respect to t and x and the symbols of the pseudo-differential operators Ω and K = H[∂ x ] (H being the Hilbert transform) are given, respectively, by g|k| and |k|, k being the Fourier transform parameter. Further, b relates to the wave surface η as

Figure 1 .

 1 Figure 1. cDZ equation: Convergence of the ground state B j (ε = 0.19, ω = 0.1, c = 0) with respect to the approximate dispersion relation that retains terms up to second, third, fourth and fifth orders (j = 2, . . . , 5). Left: soliton shapes; right: normalized differences |B j | -|B 2 | for j = 3, 4, 5.

15 D 3 E , ω = 1 Figure 2 .

 15312 Figure 2. cDZ equation: (left) Action A of ground states of the envelope B as function of the frequency ω and (right) bifurcation of a peakon from smooth ground states (ε = 0.2, c = 0).

A 28 Figure 3 .

 283 Figure 3. cDZ equation: (left) Action A of ground states of the envelope B as function of the steepness ε and (right) bifurcation of a peakon from smooth ground states (ω = 0.4, c = 0).

Figure 4 .

 4 Figure 4. cDZ peakon: (left) numerical solution via the Petviashvili scheme (ε = 0.2, c = 0, number of Fourier modes N ∼ 1.5 × 10 6 ) and (right) action of ground states of the envelope B as function of the steepness ε for different values of the frequency ω.

Figure 5 .Figure 6 .

 56 Figure 5. cDZ equation: Action A, momentum M and Hamiltonian H dependence on the frequency ω for ground states of the envelope B (ε = 0.23, c = 0).

Figure 7 .

 7 Figure 7. cDZ equation: Action A, momentum M and Hamiltonian H dependence on the velocity c for travelling waves of the envelope B(ε = 0.12, ω = 0.90).

Figure 8 .

 8 Figure8. cDZ peakon: dependence of the critical frequency ω c , peakon slope s p and amplitude a p on the steepness parameter ε. This result is obtained for ground states of the envelope B. Note that the same scalings hold true also for travelling waves.

Figure 9 .

 9 Figure9. Left: numerical peakon (solid) of the perturbed local cDZ equation (2.12) and associated analytical solution (2.13) (dash) for ε = 0.30, δ = -1.5; Right: amplitude of the numerical peakon against its theoretical counterpart (2.15) as function of ε.

Figure 10 .

 10 Figure 10. cDZ envelope equation: collision of two smooth travelling waves (ε = 0.10).

Figure 11 .

 11 Figure 11. cDZ envelope equation: Initial (1) shape and (2) after collision of the two interacting traveling waves of Figure 10.

Figure 12 .

 12 Figure 12. cDZ envelope equation: Elastic collision of two traveling waves with a ground state (ε = 0.10).

Figure 13 .

 13 Figure 13. cDZ envelope equation: Initial (1) shape and (2) after collision of the three solitons of Figure 12.

Figure 14 .

 14 Figure 14. cDZ-Dysthe envelope equation: Initial (1) shape and (2) after collision of two traveling waves and a ground state (ε = 0.10).

Acknowledgements D. Dutykh acknowledges the support from French Agence Nationale de la Recherche, project MathOcéan (Grant ANR-08-BLAN-0301-01). F. Fedele acknowledges the travel support received by the Geophysical Fluid Dynamics (GFD) Program to attend part of the summer school at the Woods Hole Oceanographic Institution in August 2011 and 2012. F. Fedele thanks Prof. Gregory Eyink for useful discussions on weak solutions of partial differential equations. The authors are also grateful to Profs. Mike Banner, Didier Clamond, Taras Lakoba, Paul Milewski, and Jianke Yang for useful discussions on the subject of nonlinear waves.


Alexandre Gramfort
email: alexandre.gramfort@inria.fr

Gaël Varoquaux

Bertrand Thirion

Christophe Pallier

Decoding Visual Percepts Induced by Word Reading with fMRI

Keywords: fMRI, supervised learning, decoding, word, reading, visual cortex, retinotopy

Word reading involves multiple cognitive processes. To infer which word is being visualized, the brain first processes the visual percept, deciphers the letters, bigrams, and activates different words based on context or prior expectation like word frequency. In this contribution, we use supervised machine learning techniques to decode the first step of this processing stream using functional Magnetic Resonance Images (fMRI). We build a decoder that predicts the visual percept formed by four letter words, allowing us to identify words that were not present in the training data. To do so, we cast the learning problem as multiple classification problems after describing words with multiple binary attributes. This work goes beyond the identification or reconstruction of single letters or simple geometrical shapes [1], [2] and addresses a challenging estimation problem, that is the prediction of multiple variables from a single observation, hence facing the problem of learning multiple predictors from correlated inputs.

I. INTRODUCTION

Supervised learning can be used to predict from functional MRI (fMRI) volumes one or multiple variables describing a stimulus [START_REF] Mitchell | Learning to decode cognitive states from brain images[END_REF]. This procedure, often called decoding [START_REF] Kamitani | Decoding the visual and subjective contents of the human brain[END_REF], is a key element in answering a question that naturally arises in fMRI studies of cognitive processes: given the observation of brain activity, what can be concluded about the mental state of the subject in a particular task [START_REF] Poldrack | Can cognitive processes be inferred from neuroimaging data?[END_REF]? Decoding methods have been extensively used to study the visual system, from texture orientation [START_REF] Kamitani | Decoding the visual and subjective contents of the human brain[END_REF], to high-level object recognition [START_REF] Haxby | Distributed and Overlapping Representations of Faces and Objects in Ventral Temporal Cortex[END_REF]. However, to inform about the full visual processing stream, a decoding experiment must go beyond classifying objects in well-separated predefined categories or predicting shape regardless of content: it should capture the brain's ability to represent rich stimuli. Written words are a typical case of such complex stimulus with high-level cognitive content [START_REF] Mitchell | Predicting Human Brain Activity Associated with the Meanings of Nouns[END_REF].

Learning predictive models from fMRI is challenging because of the high dimensionality of the data, comprising tens of thousands of voxels and the limited number of measurements, seldom more than a few hundreds of fMRI volumes. When predicting complex objects like words, it is even more challenging as the prediction space is large: we would like to predict a word among thousands of candidates, some of which have not been presented during the training sessions. This problem is also known as zero shot learning [START_REF] Palatucci | Zero-shot learning with semantic output codes[END_REF].

Reading involves various cognitive tasks, from the processing of the visual percept, assembling letters, bigrams then forming the word, and eventually disambiguating different words based on context or prior expectation like word frequency [START_REF] Dehaene | The neural code for written words: a proposal[END_REF]. A machine able to do a similar inference using fMRI data recorded while a subject reads the word should first capture the processing in the early visual areas. That is the problem addressed by this contribution using four letter words. Similar to what was done in [START_REF] Miyawaki | Visual image reconstruction from human brain activity using a combination of multiscale local image decoders[END_REF], one needs to predict letters. However, as we aim at predicting words, the stimulus should contain multiple letters rather than random patterns as in [START_REF] Miyawaki | Visual image reconstruction from human brain activity using a combination of multiscale local image decoders[END_REF]. To be able to predict unseen words, we propose a learning strategy that works with binary attributes as explained below.

Notations: We write vectors in bold, a ∈ R n , matrices with capital bold letters, A ∈ R n×n . We denote by A F = i,j a 2 ij the 2 Frobenius norm of a matrix.

II. THE DECODING OF VISUAL WORDS

We are interested in predicting words presented at fixed positions from fMRI volumes. Let us denote the target variable to be predicted as y ∈ Y and the data as x ∈ R P , where P is the number of voxels in the fMRI volume. Each x corresponds to one stimulus presentation. A training set with N volumes then forms a matrix X ∈ R N ×P .

Prediction in a space Y formed by words is challenging especially if one wants to be able to predict words not present in the training set. If one denotes the predictive function f : R P → Y, f is obtained by minimizing the empirical risk:

min f 1 N N i=1 L(y i , f (x i )) + λ Ω(f ), λ ≥ 0 (1)
where the parameter λ balances the data fit and the regularization Ω(f ). To go further we need a model for f and to define a loss function L. A four letter word as displayed on the screen during the experiment. Each letter is assembled from 12 lines (6 horizontal and 6 vertical).

To render the problem tractable we propose to describe a word by a set of binary attributes starting with visual descriptors. By displaying words as in Fig. 1, each letter is characterized by a set of lines. Discarding letters using diagonal lines (e.g. R), 12 lines suffice. The visual percept of a four letter word can then be described by 48 binary attributes. We thus cast the problem of word prediction as K = 48 binary classification problems with shared data. The space Y now equals {0, 1} K .

A natural idea is then to try to predict each binary variable independently, which boils down to minimizing a hamming loss which in practice can be achieved using a convex upper bound such as K independent hinge or logistic loss functions. If one denotes f k the function predicting for attribute k, 1 ≤ k ≤ K and if we further assume a linear model, as classically done with fMRI data, it leads to f k (x) = x T w k , which can be estimated with:

( ŵk ) k = argmin (w k ) k 1 N N i=1 K k=1 L(y i , x T i w k )+λ Ω(W), λ ≥ 0 (2) where W ∈ R P ×K .
We use a logistic loss as in [START_REF] Ryali | Sparse logistic regression for whole-brain classification of fMRI data[END_REF], [START_REF] Michel | Total variation regularization for fMRI-based prediction of behaviour[END_REF] as it gives not only a prediction of which class an observation is associated with, but also a probability for the prediction. Given an fMRI volume, the predictive framework gives the probability of each attribute composing a word, p(y k i |x i ). By further assuming independence of each attribute, which is clearly a coarse model approximation, it leads to a probabilistic output for a full word y ∈ {0, 1} K :

p(y|x i ) = K k=1 p(y k i |x i ) . (3)
This procedure provides a means to rank possible candidate words given the data. In order to go further and take into account the correlation structure between the different attributes, we also investigated a two-steps approach where a second logistic regression model was fitted on a new set of variable formed by XW ∈ R N ×K . This approach can be seen as a hierarchical decoder able to capture some correlation between attributes in order to obtain more reliable probability estimates in Eq. ( 3). Result with both models are detailed in the following section.

III. RESULTS

Before detailing results we present the experimental setup.

Paradigm: Each word was presented on the screen for 3 s at a flickering frequency of 15 Hz. A 5 s rest interval was inserted between each word presentation. The subject was asked to fixate a colored cross at the center of the screen. Each session comprised 46 words including 6 verbs. To ensure that subjects were reading, they were asked to report with a button press when a verb was presented on the screen. Repetitions corresponding to verbs were then removed from the analysis. Six sessions were recorded, leading to 240 different words used in the analysis. A partial retinotopic mapping using rotating wedges was performed to get a phase map able to delineate low level visual areas [START_REF] Sereno | Borders of multiple visual areas in human revealed by functional magnetic resonance imaging[END_REF].

MRI acquisition: MRI data were recorded on 3 subjects with a 3 Tesla Siemens Magnetom TrioTim scanner. An echo-planar imaging (EPI) scan was used to acquire functional images over the entire occipital lobe and along the temporal lobe (TR, 2400 ms; TE, 30 ms; flip angle, 81 • ; FOV, 192×192 mm; voxel size, 1.5×1.5×2 mm; 26 slices).

Data preprocessing: Data were corrected for motion artifacts. Volumes were not normalized, although SPM8 was used to segment tissues on an anatomical T1 image and to compute a normalization transform. Using this normalization, the Harvard-Oxford (HO) probabilistic atlas was warped back into subject-space and used to restrict the analysis. Analysis was done in the Lateral Occipital Cortex superior division, Lateral Occipital Cortex inferior division, Supracalcarine Cortex and Occipital Pole. The analysis was then further restricted by intersecting the HO areas with the estimated mask of the gray matter (threshold at probability > 10%). This lead to analyses on approximately P = 35 000 voxels on average across subjects. Volumes were slightly smoothed (FWHM=3mm).

Data exploration with retinotopy and univariate statistics: Retinotopic mapping was computed for each subject on the volume. For visualization purposes the phase map was projected to the cortical mesh using FreeSurfer. The flattened portion of the occipital cortex for subject 3 is presented in Fig. 2. One can observe the primary visual cortices corresponding to the largest clusters of active vertices along the calcarine fissure.

Supervised learning was performed on coefficients obtained by a standard GLM with one regressor per word presentation as well as regressors for motion and drifts with a cutoff frequency of 128 s. This resulted in one volume per word. Before running the decoding procedure, a T-test was used to detect the voxels capturing the effect of each attribute. By assigning to each voxel with a p-value below 0.001 for at least one attribute the index of the attribute leading to the smallest p-value, it is possible to map each voxel to one attribute, in the present case, one line among the 48 forming a word. Results are presented in Fig. 3. We observe activations mainly in V1 with the expected symmetry between left (resp. right) visual field and right (resp. left) hemisphere. This procedure recovers the retino- topic organization (Fig. 2) confirming a good fixation by the subject. One also observe that the lines at the periphery of the visual field yield smaller cortical activation than foveal regions. This is known as the cortical magnification factor of V1. It also suggests that a weaker fMRI signal is expected from peripheral lines, which is confirmed by further analysis using supervised learning. Note that these observations are reproducible across the 3 subjects.

Decoding results: The ultimate goal of this work is to identify the word presented to the subject in the scanner. Before addressing this hard problem, we first investigate if one can predict reliably the presence or absence of a line on the screen. To do this, we use logistic regression with a leave-one-session-out cross-validation (6 folds where each fold is a session). For computational reasons, all logistic models were fit after a univariate T-test to select the 100 most informative voxels within the training sessions. A value of 100 voxels was set as an upper bound on the expected brain volume involved in the prediction of each attribute. To control for overfitting, regularization was a squared 2 norm: Ω(W) = W 2 F . All computations were performed using the scikit-learn package [START_REF] Pedregosa | Scikit-learn: Machine Learning in Python[END_REF]. One of the challenges when working on real words is that each line does not appear on the screen half of the time. In other words, the classes are unbalanced and the number of occurrences is different for each line. This makes classification accuracy impossible to compare between lines. To alleviate this difficulty, we propose to convert the classification scores to pvalues. Assuming each prediction is drawn from a Bernoulli distribution with parameter 0 < p < 1, which corresponds to the class balance, it is possible to associate a p-value to a classification score. Prediction p-values for subject 3 are presented in Fig. 4. We can observe, that lines at the center of the screen are much more reliably predicted than lines at the periphery. This is consistent with data shown in Fig. 3. The cortical magnification factor affects the fMRI spatial resolution at the periphery of the visual field.

Prediction of words: Following Eq. ( 3), from the probabilities obtained by the logistic regression models for each line, one can compute the likelihood of a word given the data. The likelihood can then be used to rank words in a corpus. Results in terms of ranking are presented in Fig. 5. Given a rank k, this plot shows the percentage of time the true word was in the first k words. To quantify the overall performance of the method, we compute an Area under the curve (AUC) index, which would be at 0.5 should the algorithm perform at chance level. Using a Bernoulli distribution as above, it is also possible to assign a pvalue to a ranking score as illustrated in Fig. 5. From this result one can observe that our prediction algorithm performs significantly above chance. Results with the 2-level logistic Figure 6. Word prediction scores using a ranking metric for the 3 subjects using a 2-level logistic model. AUC obtained consistently increase with respect to the single logistic model. regression model, aiming at capturing correlation between visual attributes, are reported in Fig. 6. It can be observed that this hierarchical model consistently improve the results on the 3 subjects.

To further assess the statistical significance of our prediction, we used a similar testing scenario as in [START_REF] Mitchell | Predicting Human Brain Activity Associated with the Meanings of Nouns[END_REF]. The procedure consists in removing two words from the training set and testing if the algorithm trained on the remaining words can predict which fMRI volume corresponds to which of the two left out words. Results are presented in Table I.

IV. CONCLUSION

In this work we have shown that it is possible to predict the visual percept induced by a complex stimulus, such as a four letter word, using fMRI data. Experimental results confirmed well known neuroscience findings about the retinotopic organization of the primary visual areas. This turned out to be an issue when working with multi-letter words as peripheral letters were not very accurately predicted by V1. Working with multiple letters is nevertheless required when training an algorithm to extract relevant information from fMRI data during the reading process. Modeling correlations between predicted targets using a hierarchical predictive model achieves better prediction scores. More elaborate attributes about the word such as the semantic content as investigated in [START_REF] Mitchell | Predicting Human Brain Activity Associated with the Meanings of Nouns[END_REF], [START_REF] Pereira | Generating text from functional brain images[END_REF] will be used in further work to improve prediction.

Figure 1 .

 1 Figure 1.

Figure 2 .

 2 Figure 2. Retinotopic phase maps of subject 3 projected onto the flattened cortical mantel of the occipital lobe (left and right hemispheres). The cut in the mantel corresponds the calcarine fissure along the medial wall. The colormap represent the polar angle of the visual field, shifted by -π/2 and constrained to the [-π, π] interval.

Figure 3 .

 3 Figure 3. The 48 visual bars displayed with a color code (blue to red from left to right) with corresponding voxels yielding the maximum responses (threshold on T-test p < 0.001). Results were projected on the flattened cortical mantel around the calcarine fissure (left and right hemispheres).

Figure 4 .

 4 Figure 4. Prediction scores for each visual bar (S3). As lines have different class balance (some appear more often then others) the classification scores are converted to p-values assuming each prediction is drawn from a Bernoulli distribution with appropriate parameter. What is represented here is -log 10(p value ) meaning that red corresponds to p-values below 0.01.

Figure 5 .

 5 Figure 5. Word prediction scores using a ranking metric for the 3 subjects.

ACKNOWLEDGMENT

This work was supported by the ANR grant ViMAGINE ANR-08-BLAN-0250-02 and ANR grant IRMGroup, ANR-10-BLAN-0126-02.


Franc ¸ois Caron

Yee Whye Teh

Bayesian nonparametric models for ranked data

Keywords: choice models, generalized Bradley-Terry model, Plackett-Luce model, gamma process, Markov Chain Monte Carlo Modèles de choix, modèle de Bradley-Terry généralisé, modèle de Plackett-Luce, processus de gamma, méthodes de

We develop a Bayesian nonparametric extension of the popular Plackett-Luce choice model that can handle an infinite number of choice items. Our framework is based on the theory of random atomic measures, with the prior specified by a gamma process. We derive a posterior characterization and a simple and effective Gibbs sampler for posterior simulation. We develop a time-varying extension of our model, and apply it to the New York Times lists of weekly bestselling books.

Introduction

Data in the form of partial rankings, i.e. in terms of an ordered list of the top-m items, arise in many contexts. For example, in this paper we consider datasets consisting of the top 20 bestselling books as published each week by the New York Times. The Plackett-Luce model [START_REF] Luce | Individual choice behavior: A theoretical analysis[END_REF][START_REF] Plackett | The analysis of permutations[END_REF] is a popular model for modeling such partial rankings of a finite collection of M items. It has found many applications, including choice modeling [START_REF] Luce | The choice axiom after twenty years[END_REF], sport ranking [START_REF] Hunter | MM algorithms for generalized Bradley-Terry models[END_REF], and voting [START_REF] Gormley | Exploring voting blocs with the Irish electorate: a mixture modeling approach[END_REF]. [START_REF] Diaconis | Group representations in probability and statistics[END_REF]Chap. 9] provides detailed discussions on the statistical foundations of this model.

In the Plackett-Luce model, each item k ∈ [M ] = {1, . . . , M } is assigned a positive rating parameter w k , which represents the desirability or rating of a product in the case of choice modeling, or the skill of a player in sport rankings. The Plackett-Luce model assumes the following generative story for a top-m list ρ = (ρ 1 , . . . , ρ m ) of items ρ i ∈ [M ]: At each stage i = 1, . . . , m, an item is chosen to be the ith item in the list from among the items that have not yet appeared, with the probability that ρ i is selected being proportional to its desirability w ρi . The overall probability of a given partial ranking ρ is then:

P (ρ) = m i=1 w ρi M k=1 w k - i-1 j=1 w ρj . ( 1
)
with the denominator in (1) being the sum over all items not yet selected at stage i.

In many situations the collection of available items can be very large and potentially unknown. In this case, a nonparametric approach can be sensible, where the pool of items is assumed to be infinite and the model allows for the possibility of items not observed in previous top-m lists to appear in new ones. In this paper we propose such a Bayesian nonparametric Plackett-Luce model. Our approach is built upon recent work on Bayesian inference for the (finite) Plackett-Luce model and its extensions [START_REF] Gormley | A grade of membership model for rank data[END_REF][START_REF] Guiver | Bayesian inference for Plackett-Luce ranking models[END_REF][START_REF] Caron | Efficient Bayesian inference for generalized Bradley-Terry models[END_REF]. Our model assumes the existence of an infinite pool of items {X k } ∞ k=1 , each with its own rating parameter, {w k } ∞ k=1 . The probability of a top-m list of items, say (X ρ1 , . . . , X ρm ), is then a direct extension of the finite case [START_REF] Luce | Individual choice behavior: A theoretical analysis[END_REF]:

P (X ρ1 , . . . , X ρm ) = m i=1 w ρi ∞ k=1 w k - i-1 j=1 w ρj . (2)
To formalize the framework, a natural representation to encapsulate the pool of items along with their ratings is using an atomic measure:

G = ∞ k=1 w k δ X k (3)
Using this representation, note that the top item X ρ1 in our list is simply a draw from the probability measure obtained by normalizing G, while subsequent items in the top-m list are draws from probability measures obtained by first removing from G the atoms corresponding to previously picked items and normalizing. Described this way, it is clear that the Plackett-Luce model is basically a partial size-biased permutation of the atoms in G [START_REF] Patil | Diversity as a concept and its implications for random communities[END_REF], and the existing machinery of random measures and exchangeable random partitions [START_REF] Pitman | Combinatorial stochastic processes[END_REF] can be brought to bear on our problem. In particular, in Section 2 we will use a gamma process as the prior over the atomic measure G. This is a completely random measure [START_REF] Kingman | Completely random measures[END_REF] with gamma marginals, such that the corresponding normalized probability measure is a Dirichlet process. We will show that with the introduction of a suitable set of auxiliary variables, we can characterize the posterior law of G given observations of top-m lists distributed according to [START_REF] Plackett | The analysis of permutations[END_REF]. A simple Gibbs sampler can then be derived to simulate from the posterior distribution. In Section 3 we develop a time-varying extension of our model and derive a simple and effective Gibbs sampler for posterior simulation. In Section 4 we apply our time-varying Bayesian nonparametric Plackett-Luce model to the aforementioned New York Times bestsellers datasets, and conclude in Section 5.

RR n°8140

A Bayesian nonparametric model for partial ranking

We start this section by briefly describing a Bayesian approach to inference in finite Plackett-Luce models [START_REF] Caron | Efficient Bayesian inference for generalized Bradley-Terry models[END_REF], and taking the infinite limit to arrive at the nonparametric model. This will give good intuitions for how the model operates, before we rederive the same nonparametric model more formally using gamma processes. Throughout this paper we will suppose that our data consists of L partial rankings, with ρ = (ρ 1 , . . . , ρ m ) for ∈ [L]. For notational simplicity we assume that all the partial rankings are length m.

Finite Plackett-Luce model with gamma prior

Suppose we have M choice items, with item k ∈ [M ] having a positive desirability parameter w k . A partial ranking ρ = (ρ 1 , . . . , ρ m ) can be constructed generatively by picking the ith item ρ i at the ith stage for i = 1, . . . , m, with probability proportional to w ρ i as in [START_REF] Luce | Individual choice behavior: A theoretical analysis[END_REF]. An alternative Thurstonian interpretation, which will be important in the following, is as follows: For each item k let z k ∼ Exp(w k ) be exponentially distributed with rate w k . Thinking of z k as the arrival time of item k in a race, let ρ i be the index of the ith item to arrive (the ith smallest value among (z k ) M k=1 ). The resulting probability of ρ can then be shown to still be [START_REF] Luce | Individual choice behavior: A theoretical analysis[END_REF]. In this interpretation (z k ) can be understood as latent variables, and the EM algorithm can be applied to derive an algorithm to find a ML parameter setting for (w k ) M k=1 given multiple partial rankings. Unfortunately the posterior distribution of (z k ) given ρ is difficult to compute directly, so we instead consider an alternative parameterization: Let Z i = z ρ i -z ρ i-1 be the waiting time for the ith item to arrive after the i -1th item (with z ρ 0 defined to be 0). Then it can be shown that the joint probability is:

P ((ρ ) L =1 , (Z i ) L,m =1,i=1 |(w k ) M k=1 ) = L =1 m i=1 w ρ i exp -Z i M k=1 w k - i-1 j=1 w ρ j (4)
Note that the posterior of (Z i ) m i=1 is simply factorized with

Z i |ρ, w ∼ Exp( M k=1 w k - i-1
j=1 w ρ j ), and the ML parameter setting can be easily derived as well. Taking a further step, we note that a factorized gamma prior over (w k ) is conjugate to (4), say w k ∼ Gamma( α M , τ ) with hyperparameters α, τ > 0. Now Bayesian inference can be carried out either with a VB EM algorithm, or a Gibbs sampler. In this paper we shall consider only Gibbs sampling algorithms. In this case the parameter updates are of the form

w k |(ρ ), (Z i ), (w k ) k =k ∼ Gamma α M + n k , τ + L =1 m i=1 δ ik Z i (5)
where n k is the number of occurrences of item k among the observed partial rankings, and δ ik = 0 if there is a j < i with ρ j = k and 1 otherwise. These terms arise by regrouping those in the exponential in [START_REF] Hunter | MM algorithms for generalized Bradley-Terry models[END_REF].

A nonparametric Plackett-Luce model can now be easily derived by taking the limit as the number of choice items M → ∞. For those items k that have appeared among the observed partial rankings, the limiting conditional distribution (5) is well defined since n k > 0. For items that did not appear in the observations, (5) becomes degenerate at 0. Instead we can define w * = k:n k =0 w k to be the total desirability among all infinitely many previously unobserved items, and show that

w * |(ρ ), (Z i ), (w k ) k:n k >0 ∼ Gamma α, τ + L =1 m i=1 Z i (6)
The Gibbs sampler thus alternates between updating (Z i ), and updating the ratings of the observed items (w k ) k:n k >0 and of the unobserved ones w * . This nonparametric model allows us to estimate the probability of seeing new items appearing in future partial rankings in a consistent manner. While intuitive, this RR n°8140

derivation is ad hoc in the sense that it arises as the infinite limit of the Gibbs sampler for finite models, and is unsatisfying as it did not directly capture the structure of the underlying infinite dimensional object, which we will show in the next subsection to be a gamma process.

A Bayesian nonparametric Plackett-Luce model

Let X be a measurable space of choice items. A gamma process is a completely random measure over X with gamma marginals. Specifically, it is a random atomic measure of the form (3), such that for each measurable subset A, the (random) mass G(A) is gamma distributed. Assuming that G has no fixed atoms (that is, for each element x ∈ X we have G({x}) = 0 with probability one) and that the atom locations {X k } are independent of their masses {w k }, it can be shown that such a random measure can be constructed as follows: each X k is iid according to a base distribution H (which we assume is non-atomic with density h(x)), while the set of masses {w k } is distributed according to a Poisson process over R + with intensity λ(w) = αw -1 e -wτ where α > 0 is the concentration parameter and τ > 0 the inverse scale. We write this as G ∼ Γ(α, τ, H). Under this parametrization, we have that

G(A) ∼ Gamma(αH(A), τ ).
Each atom X k is a choice item, with its mass w k > 0 corresponding to the desirability parameter. The Thurstonian view described in the finite model can be easily extended to the nonparametric one, where a partial ranking (X ρ 1 . . . X ρ m ) can be generated as the first m items to arrive in a race. In particular, for each atom X k let z k ∼ Exp(w k ) be the time of arrival of X k and X ρ i the ith item to arrive. The first m items to arrive (X ρ 1 . . . X ρ m ) then constitutes our top-m list, with probability as given in [START_REF] Plackett | The analysis of permutations[END_REF]. Again reparametrizing using inter-arrival durations, let

Z i = z ρ i -z ρ i-1 for i = 1, 2, . . . (with z ρ0 = 0).
Then the joint probability is:

P ((X ρ i ) m i=1 , (Z i ) m i=1 |G) = P ((z ρ 1 . . . z ρ m ), and z k > z ρ m for all k ∈ {ρ 1 , . . . , ρ m }) (7)
= m i=1 w ρ i e -wρ i zρ i k ∈{ρ i } m i=1 e -w k zρ m = m i=1 w ρ i exp -Z i ∞ k=1 w k - i-1 j=1 w ρ j
Marginalizing out (Z i ) m i=1 gives the probability of (X ρ i ) m i=1 in (2). Further, conditional on ρ it is seen that the inter-arrival durations Z 1 . . . Z m are mutually independent and exponentially distributed:

Z i |(X ρ i ) m i=1 , G ∼ Exp ∞ k=1 w k - i-1 j=1 w ρ j (8)
The above construction is depicted on Figure 1(left). We visualize on right some top-m lists generated from the model, with τ = 1 and different values of α.

Posterior characterization

Consider a number L of partial rankings, with the th list denoted

Y = (Y 1 . . . Y m ) , for ∈ [L].
While previously our top-m list (X ρ1 . . . X ρm ) consists of an ordered list of the atoms in G. Here G is unobserved and (Y 1 . . . Y m ) is simply a list of observed choice items, which is why they were not expressed as an ordered list of atoms in G. The task here is then to characterize the posterior law of G under a gamma process prior and supposing that the observed partial rankings were drawn iid from the nonparametric Plackett-Luce model given G. Re-expressing the conditional distribution (2) of Y given G, we have:

P (Y |G) = m i=1 G({Y i }) G(X\{Y 1 . . . Y i-1 }) (9)
RR n°8140 As before, for each , we will also introduce a set of auxiliary variables Z = (Z 1 . . . Z m ) (the interarrival times) that are conditionally mutually independent given G and Y , with:

Z i |Y , G ∼ Exp(G(X\{Y 1 , . . . , Y i-1 })) (10)
The joint probability of the item lists and auxiliary variables is then (c.f. ( 7)):

P ((Y , Z ) L =1 |G) = L =1 m i=1 G({Y i }) exp(-Z i G(X\{Y 1 , . . . , Y i-1 })) (11)
Note that under the generative process described in Section 2.2, there is positive probability that an item appearing in a list Y appears in another list Y with = . Denote the unique items among all L lists by X * 1 . . . X * K , and for each k = 1, . . . , K let n k be the number of occurrences of X * k among the item lists. Finally define occurrence indicators

δ ik = 0 if ∃j < i with Y j = X * k ; 1 otherwise. (12)
i.e. δ ik is the indicator of the occurence that item X * k does not appear at a rank lower than i in the th list. Then the joint probability under the nonparametric Plackett-Luce model is:

P ((Y , Z ) L =1 |G) = K k=1 G({X * k }) n k × L =1 m i=1 exp(-Z i G(X\{Y 1 , . . . , Y i-1 })) = exp -G(X) i Z i K k=1 G({X * k }) n k exp -G({X * k }) i (δ ik -1)Z i (13)
Taking expectation of ( 13) with respect to G using the Palm formula gives:

Theorem 1 The marginal probability of the L partial rankings and auxiliary variables is:

P ((Y , Z ) L =1 ) = e -ψ( i Z i ) K k=1 h(X * k )κ n k , i δ ik Z i ( 14
)
where ψ(z) is the Laplace transform of λ,

ψ(z) = -log E e -zG(X) = R + λ(w)(1 -e -zw )dw = α log 1 + z τ (15)
and κ(n, z) is the nth moment of the exponentially tilted Lévy intensity λ(w)e -zw :

κ(n, z) = R + λ(w)w n e -zw dw = α (z + τ ) n Γ(n) (16)
Details are given in the appendix. Another application of the Palm formula now allows us to derive a posterior characterisation of G:

Theorem 2 Given the observations and associated auxiliary variables (Y , Z ) L =1 , the posterior law of G is also a gamma process, but with atoms with both fixed and random locations. Specifically,

G|(Y , Z ) L =1 = G * + K k=1 w * k δ X * k ( 17
)
where G * and w * 1 , . . . , w * K are mutually independent. The law of G * is still a gamma process,

G * |(X , Z ) L =1 ∼ Γ(α, τ * , h) τ * = τ + i Z i (18)
while the masses have distributions,

w * k |(Y , Z ) L =1 ∼ Gamma n k , τ + i δ ik Z i (19)

Gibbs sampling

Given the results of the previous section, a simple Gibbs sampler can now be derived, where all the conditionals are of known analytic form. In particular, we will integrate out all of G * except for its total mass w * * = G * (X). This leaves the latent variables to consist of the masses w * * , (w * k ) and the auxiliary variables (Z i ). The update for Z i is given by [START_REF] Patil | Diversity as a concept and its implications for random communities[END_REF], while those for the masses are given in Theorem 2:

Gibbs update for Z i :

Z i |rest ∼ Exp w * * + k δ ik w * k (20)
Gibbs update for w * k :

w * k |rest ∼ Gamma n k , τ + i δ ik Z i (21)
Gibbs update for w * * :

w * * |rest ∼ Gamma α, τ + i Z i (22)
Note that the auxiliary variables are conditionally independent given the masses and vice versa. Hyperparameters of the gamma process can be simply derived from the joint distribution in Theorem 1. Since the marginal probability of the partial rankings is invariant to rescaling of the masses, it is sufficient to keep τ fixed at 1. As for α, if a Gamma(a, b) prior is placed on it, its conditional distribution is still gamma:

Gibbs update for α: α|rest ∼ Gamma a + K, b + log 1 + i Z i τ ( 23
)
Note that this update was derived with w * * marginalized out, so after an update to α it is necessary to immediately update w * * via ( 22) before proceeding to update other variables.

RR n°8140

3 Dynamic Bayesian nonparametric ranking models

In this section we develop an extension of the Bayesian nonparametric Plackett-Luce model to model time-varying rankings, where the rating parameters of items may change smoothly over time and reflected in a changing series of rankings. Given a series of times indexed by t = 1, 2, . . ., we may model the rankings at time t using a gamma process distributed random measure G t as in Section 2.2, with Markov dependence among the sequence of measures (G t ) enabling dependence among the rankings over time.

Pitt-Walker dependence model

We will construct a dependent sequence (G t ) which marginally follow a gamma process Γ(α, τ, H) using the construction of [START_REF] Pitt | Constructing stationary time series models using auxiliary variables with applications[END_REF]. Suppose G t ∼ Γ(α, τ, H). Since G t is atomic, we can write it in the form:

G t = ∞ k=1 w tk δ X tk (24)
Define a random measure C t with conditional law:

C t |G t = ∞ k=1 c tk δ X tk c tk |G t ∼ Poisson(φ t w tk ) (25)
where φ t > 0 is a dependence parameter. Using the same method as in Section 2.3, we can show:

Proposition 3 Suppose the law of G t is Γ(α, τ, H).
The conditional law of G t given C t is then:

G t = G * t + ∞ k=1 w * tk δ X tk (26)
where G * t and (w * tk ) ∞ k=1 are all mutually independent. The law of G * t is given by a gamma process, while the masses are conditionally gamma,

G * t |C t ∼ Γ(α, τ + φ t , H) w * tk |C t ∼ Gamma(c tk , τ + φ t ) (27)
The idea of [START_REF] Pitt | Constructing stationary time series models using auxiliary variables with applications[END_REF] is to define the conditional law of G t+1 given G t and C t to coincide with the conditional law of G t given C t as in Proposition 3. In other words, define

G t+1 = G * t+1 + ∞ k=1 w t+1,k δ X tk (28)
where G * t+1 ∼ Γ(α, τ + φ t , H) and w t+1,k ∼ Gamma(c tk , τ + φ t ) are mutually independent. If the prior law of G t is Γ(α, τ, H), the marginal law of G t+1 will be Γ(α, τ, H) as well when both G t and C t are marginalized out, thus maintaining a form of stationarity. Further, although we have described the process in order of increasing t, the joint law of G t , C t , G t+1 can equivalently be described in the reverse order with the same conditional laws as above. Note that if c tk = 0, the conditional distribution of w t+1,k will be degenerate at 0. Hence G t+1 has an atom at X tk if and only if C t has an atom at X tk , that is, if c tk > 0. In addition, it also has atoms (those in G * t+1 ) where C t does not (nor does G t ). Finally, the parameter φ t can be interpreted as controlling the strength of dependence between G t+1 and G t . Indeed it can be shown that

E[G t+1 |G t ] = φ t φ t + τ G t + τ φ t + τ H. (29)
Another measure of dependence can be gleaned by examining the "lifetime" of an atom. Suppose X is an atom in G 1 with mass w > 0. The probability that X is an atom in C 2 with positive mass is 1 -exp(-φ 1 w), in which case it has positive mass in G 2 as well. Conversely, once it is not an atom, it will never be an atom in the future since the base distribution H is non-atomic. The lifetime of the atom is then the smallest t such that it is no longer an atom. We can show by induction that: (details in appendix)

Proposition 4 The probability that an atom X in G 1 with mass w > 0 is dead at time t is given by

P (G t ({X}) = 0|w) = exp(-y t|1 w)
where y t|1 can be obtained by the recurrence y t|t-1 = φ t-1 and y t|s-1 = y t|s φs-1 φs-1+τ +y t|s .

Posterior characterization and Gibbs sampling

Assume for simplicity that at each time step t = 1, . . . , T we observe one top-m list Y t = (Y t1 , . . . , Y tm ) (it trivially extends to multiple partial rankings of differing sizes). We extend the results of the previous section in characterizing the posterior and developing a Gibbs sampler for the dynamical model. Since each observed item at time t has to be an atom in its corresponding random measure G t , and atoms in G t can propagate to neighboring random measures via the Pitt-Walker dependence model, we conclude that the set of all observed items (through all times) has to include all fixed atoms in the posterior of G t . Thus let X * = (X * k ), k = 1, . . . , K be the set of unique items observed in Y 1 , . . . , Y T , let n tk ∈ {0, 1} be the number of times the item X * k appears at time t, and let ρ t be defined as Y t = (X * ρ1 , . . . , X * ρm ). We write the masses of the fixed atoms as w tk = G t ({X * k }), while the total mass of all other random atoms is denoted w t * = G t (X\X * ). Note that w tk has to be positive on a random contiguous interval of time that includes all observations of X * k -it's lifetime-but is zero outside of the interval. We also write c tk = C t ({X * k }) and c t * = C t (X\X * ). As before, we introduce, for t = 1, . . . , T and i = 1, . . . , m, latent variables

Z ti ∼ Exp w t * + K k=1 w tk - i-1 j=1 w tρj (30)
Each iteration of the Gibbs sampler then proceeds as follows (details in appendix). The latent variables (Z ti ) are updated as above. Conditioned on the latent variables (Z ti ), (c tk ) and (c t * ), we update the masses (w tk ), which are independent and gamma distributed since all likelihoods are of gamma form. Note that the total masses (G t (X)) are not likelihood identifiable, so we introduce an extra step to improve mixing by sampling them from the prior (integrating out (c tk ), (c t * )), scaling all masses along with it.

Directly after this step we update (c tk ), (c t * ). We update α along with the random masses (w t * ) and (c t * ) efficiently using a forward-backward recursion. Finally, the dependence parameters (φ t ) are updated.

Continuous time formulation using superprocesses

The dynamic model described in the previous section is formulated for discrete time data. When the time interval between ranking observations is not constant, it is desirable to work with dynamic models evolving over continuous-time instead, with the underlying random measures (G t ) defined over all t ∈ R, but with observations at a discrete set of times t 1 < t 2 < • • • . Here we propose a continuous-time model based on the Dawson-Watanabe superprocess [START_REF] Watanabe | A limit theorem of branching processes and continuous state branching processes[END_REF][START_REF] Dawson | Stochastic evolution equations and related measure processes[END_REF] (see also [START_REF] Ethier | The transition function of a measure-valued branching diffusion with immigration[END_REF][START_REF] Mena | On a construction of Markov models in continuous time[END_REF][START_REF] Feng | Poisson-Dirichlet Distribution and Related Topics[END_REF][START_REF] Cox | A theory of the term structure of interest rates[END_REF]). This is a diffusion on the space of measures with the gamma process Γ(α, τ, H) as its equilibrium distribution. It is defined by a generator with ξ parametrizing the rate of evolution. Figure 2 gives a sample path, where we see that it is continuous but non-differentiable. For efficient inference, it is desirable to be able to integrate out all G t 's except those G t1 , G t2 , . . . at observation times. An advantage to using the Dawson-Watanabe superprocess is that, the conditional distribution of G ts given G ts-1 is remarkably simple [START_REF] Ethier | The transition function of a measure-valued branching diffusion with immigration[END_REF]. In particular it is simply given by the discrete-time process of the previous section with dependence parameter φ ts|ts-1 = τ e τ ξ(ts-t s-1 ) -1

L = ξ G(dX) ∂ 2 ∂G(X) 2 + α H(dX) ∂ ∂G(X) -τ G(dX) ∂ ∂G(X) RR n°8140
. Thus the inference algorithm developed previously is directly applicable to the continuous-time model too.

Experiments

We apply the discrete-time dynamic Plackett-Luce model to the New York Times bestsellers data. These consist of the weekly top-20 best-sellers list from June 2008 to April 2012 in various categories. We consider here the categories paperback nonfiction (PN) and hardcover fiction (HF), for which respectively 249 and 916 books appear at least once in the top-20 lists over the 200 weeks. We consider that the correlation parameter φ t = φ is constant over time, and assign flat improper priors p(α) ∝ 1/α and p(φ) ∝ 1/φ. In order to take into account the publication date of a book, we do not consider books in the likelihood before their first appearance in a list. We run the Gibbs sampler with 10000 burn-in iterations followed by 10000 samples. Mean normalized weights for the more popular books in both categories are shown in Figure 3.

The model is able to estimate the weights associated to each book that appeared at least once, as well as the total weight associated to all other books, i.e. the probability that a new book enters at the first rank in the list, represented by the black curve. Moreover, the Bayesian approach enables us to have a measure of the uncertainty on the weights. The hardcover fiction category is characterized by rapid changes in successive lists, compared to the paperback nonfiction. This is quantified by the estimated value of the parameter φ, which are respectively 85 ± 20 and 140 ± 40 for PN and HF. The estimated values of the shape parameter α are 7 ± 1.5 and 2 ± 1 respectively.

Discussion

We have proposed a Bayesian nonparametric Plackett-Luce model for ranked data. Our approach is based on the theory of atomic random measures, where we showed that the Plackett-Luce generative model corresponds exactly to a size-biased permutation of the atoms in the random measure. We characterized the posterior distribution, and derived a simple MCMC sampling algorithm for posterior simulation. Our approach can be see as a multi-stage generalization of posterior inference in normalized random measures [START_REF] James | Posterior analysis for normalized random measures with independent increments[END_REF][START_REF] Griffin | Posterior simulation of normalized random measure mixtures[END_REF][START_REF] Favaro | MCMC for normalized random measure mixture models[END_REF], and can be easily extended from gamma processes to general completely random measures.

We also proposed dynamical extensions of our model for both discrete and continuous time data, and applied it to modeling the bestsellers' lists on the New York Times. Our dynamic extension may be useful for modeling time varying densities or clusterings as well. In our experiments we found that our model is insufficient to capture the empirical observation that bestsellers often start off high on the lists and tail off afterwards, since our model has continuous sample paths. We adjusted for this by simply not including books in the model prior to their publication date. It may be possible to model this better using models with discontinuous sample paths, for example, the Orstein-Uhlenbeck approach of [START_REF] Griffin | The Ornstein-Uhlenbeck Dirichlet process and other time-varying processes for Bayesian nonparametric inference[END_REF] where the process evolves via a series of discrete jump events instead of continuously.

A Proof of Theorem 1

The marginal probability ( 14) is obtained by taking the expectation of ( 13) with respect to G. Note however that ( 13) is a density, so to be totally precise here we need to work with the probability of infinitesimal neighborhoods around the observations instead, which introduces significant notational complexity. To keep the notation simple, we will work with densities, leaving it to the careful reader to verify that the calculations indeed carry over to the case of probabilities.

P ((Y , Z ) L =1 ) =E P ((Y , Z ) L =1 |G) =E e -G(X) i Z i K k=1 G({X * k }) n k e -G({X * k }) i (δ ik -1)Z i
The gamma prior on G = ∞ j=1 w j δ Xj is equivalent to a Poisson process prior on N = ∞ j=1 δ (wj ,Xj ) defined over the space R + × X with mean intensity λ(w)h(x). Then,

=E   e -wN (dw,dx) i Z i K k=1 ∞ j=1 w n k j 1(X j = X * k )e -wj i (δ ik -1)Z i  
Applying the Palm formula for Poisson processes to pull the k = 1 term out of the expectation,

= E   e -w(N +δ w * 1 ,x * 1 )(dw,dx) i Z i K k=2 ∞ j=1 w n k j 1(X j = X * k )e -wj i (δ ik -1)Z i   × (w * 1 ) n k h(X * 1 )e -w * 1 i (δ i1 -1)Z i λ(w * 1 )dw * 1 =E   e -wN (dw,dx) i Z i K k=2 ∞ j=1 w n k j 1(X j = X * k )e -wj i (δ ik -1)Z i   × h(X * 1 ) (w * 1 ) n k e -w * 1 i δ i1 Z i λ(w * 1 )dw * 1 RR n°8140
Now iteratively pull out terms k = 2, . . . , K using the same idea, and we get:

=E e -G(X) i Z i K k=1 h(X * k ) (w * k ) n k e -w * k i δ ik Z i λ(w * k )dw * k =e -ψ( i Z i ) K k=1 h(X * k )κ n k , i δ ik Z i (31)
This completes the proof of Theorem 1.

B Proof of Theorem 2

Let f : X → R be measurable with respect to H. Then the characteristic functional of the posterior G is given by:

E[e -f (x)G(dx) |(Y , Z ) L =1 ] = E[e -f (x)G(dx) P ((Y , Z ) L =1 |G)] E[P ((Y , Z ) L =1 |G)] (32)
The proof is essentially obtained by calculating the numerator and denominator of (32). The denominator is already given in Theorem 1. The numerator is obtained using the same technique with the inclusion of the term e f (x)G(dx) , which gives:

E e -f (x)G(dx) P ((Y , Z ) L =1 |G) =E e -(f (x)+ i Z i )G(dx) K k=1 h(X * k ) (w * k ) n k e -w * k (f (X * k )+ i δ ik Z i ) λ(w * k )dw * k
By the Lévy-Khintchine Theorem (using the fact that G has a Poisson process representation N ),

= exp -(1 -e -w(f (x)+ i Z i ) )λ(w)h(x)dwdx × K k=1 h(X * k ) (w * k ) n k e -w * k (f (X * k )+ i δ ik Z i ) λ(w * k )dw * k (33)
Dividing the numerator (31) by the denominator (33), the characteristic functional of the posterior G is:

E e -f (x)G(dx) |(Y , Z ) L =1 = exp -(1 -e -wf (x) )e -i Z i λ(w)h(x)dwdx × K k=1 h(X * k ) e -f (X * k ) (w * k ) n k e -w * k i δ ik Z i λ(w * k )dw * k (w * k ) n k e -w * k i δ ik Z i λ(w * k )dw * k (34)
Since the characteristic functional is the product of K + 1 terms, we see that the posterior G consists of K + 1 independent components, one corresponding to the first term above (G * ), and the others corresponding to the K terms in the product over k. Substituting the Lévy measure λ(w) for a gamma process, we note that the first term shows that G * is a gamma process with updated inverse scale τ * . The kth term in the product shows that the corresponding component is an atom located at X * k with density

(w * k ) n k e -w * k i δ ik Z i λ(w * k )
; this is the density of the gamma distribution over w * k in Theorem 2. This completes the proof.

RR n°8140

C Proof of Proposition 4

We have P

(G t (X 1k ) = 0|w t-1,,k ) = exp(-φ t-1 w t-1,k ) Assume that P (G t (X 1k ) = 0|w sk ) = exp(-y t|s w sk ) then P (G t (X 1k ) = 0|w s-1,k ) = exp(-y t|s w sk )p(w sk |w s-1,k )dw sk = c s-1,k exp(-y t|s w sk )p(w sk |c s-1,k )p(c s-1,k |w s-1,k )dw sk = c s-1,k exp -c s-1,k log 1 + y t|s φ s-1 + τ p(c s-1,k |w s-1,k ) = exp -y t|s φ s-1 φ s-1 + τ + y t|s w s-1,k

D Gibbs sampler for the dynamic nonparametric Plackett-Luce model

For ease of presentation, we assume that φ t takes the same value φ at each time step. The Gibbs sampler will iterate between the following steps The steps are now fully described.

1.a) Sample (G t (X)) given (α, φ) We have G 1 (X)|α ∼ Gamma(α, τ ) and for t = 1, . . . , T -1 G t+1 (X) ∼ Gamma(α + M t , τ + φ)
where M t ∼Poisson(φG t (X)). The weights (w t , w t * ) are then appropriately rescaled.

RR n°8140 If w t+1,k = 0, we only have two possible moves: c tk = 0 or c tk = 1, given by the following probabilities

P (c tk = 0|w t+1,k = 0, w tk ) = exp(-φw tk ) exp(-φw tk ) + φw tk exp(-φw tk )(τ + φ) = 1 1 + φw tk (τ + φ) P (c tk = 1|w t+1,k = 0, w tk ) =
φw tk exp(-φw tk )(τ + φ) exp(-φw tk ) + φw tk exp(-φw tk )(τ + φ) = φw tk (τ + φ) 1 + φw tk (τ + φ) Note that the above Markov chain is not irreducible, as the probability is zero to go from a state (c tk > 0, w t+1,k > 0) to a state (c tk = 0, w t+1,k = 0), even though the posterior probability of this event is non zero in the case item k does not appear after time t. We can resolve that by the following procedure, that uses a backward forward recursion.

Assume that item k does not appear after time step τ + k (the same procedure applies if item k does not appear before time step τ - k ). Then we can sample jointly the whole sequence (w k,t , c k,t ) t=τ k +1,...,T using the following backward forward recursion. Let

Figure 1 :

 1 Figure 1: Bayesian nonparametric Plackett-Luce model. Left: G and U = k u k δ X k where u k = -log(z k ). The top-3 ranking is (ρ 1 , ρ 2 , ρ 3 ). Right: Visualization of top-5 rankings with rows corresponding to different rankings and columns to items sorted by size biased order. A lighter shade corresponds to a higher rank. Each figure is for a different G, with α = .1, 1, 3.

Figure 2 :

 2 Figure 2: Sample path drawn from the Dawson-Watanabe superprocess. Each colour represents an atom, with height being its (varying) mass. Left shows (G t ) and right (G t /G t (X)), a Fleming-Viot process.

Figure 3 :

 3 Figure 3: Mean normalized weights for paperback nonfiction (left) and hardcover fiction (right). The black lines represent the weight associated to all the books that have not appear in the top-20 lists.

φx t+1 1 +• 4 )δ tik w tk (45) 5 )

 145 φ + x t+1We have, for k = 1, . . . , K and t = τ + k c tk |(Z, φ, w tk ) ∼ Poisson 1 + φ 1 + φ + x t φw tk (36)w t+1,k |c tk , Z ∼ Gamma (c k,t , τ + φ + x t+1 )(37)RR n°81402.a) Sample α given (Z, φ)We can sample from the full conditional which is given byα|(Z, γ, φ) ∼ Gamma (a + K, b + y 1 + log(1 + x 1 ))(38)where x 1 and y 1 are obtained with the following recursion t = T -1, . . . , 1x t = m k=1 Z tk + φx t+1 1 + φ + x t+1 y t = y t+1 -log 1 + φ 1 + φ + x t+1 . 2.b) Sample (c * , w * ) given (Z, φ, α)We can sample from the full conditional which is given byw 1 * |(Z, φ, α) ∼ Gamma (α, τ + x 1 )(41)where x 1 is defined above. Then for t = 2, . . . , T , letc t-1 * |(Z, φ, α, w t-1 * ) ∼ Poisson 1 + φ 1 + φ + x t φw t-1 * w t * |c t-1 * , Z, α ∼ Gamma (α + c t-1 * , τ + φ + x t )3) Sample (w, w * ) given (Z, α, c, c * , φ)For each time step t = 1, . . . , T• For each item k = 1, . . . , K, samplew tk |c t-1,k , c tk , Z t ∼ Gamma n tk + c t-1,k + c tk , τ + 2φ + m i=1 δ tik Z ti(42)if c tk + c t-1,k + n tk > 0, otherwise, set w tk = 0. The occurence indicator δ tik is defined asδ tik = 0 if ∃j < i with Y tj = X * k ; Samplethe total mass w t * |c t * , c t-1 * , Z t , α ∼ Gamma α + c t * + c t-1 * , τ + 2φ + Sample Z given (w, w * ) For t = 1, . . . , T and i = 1, . . . m, sample Z ti |w, w * ∼ Exp w t * + K k=1 Sample φ given w, w * , α, φ We sample φ using a MH step. Propose φ = φ exp(σε) where σ > 0 and ε ∼ N (0, 1). And accept it with probability t+1 * | φ, w t * ) p(w t+1 * |φ, w t * ) K k=1 p(w t+1,k | φ, w tk ) p(w t+1,k |φ, w tk ) (46) RR n°8140 RESEARCH CENTRE BORDEAUX -SUD-OUEST 351, Cours de la Libération Bâtiment A 29 33405 Talence Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr

  Consider first the sampling of c 1:T . We have, for t = 1, . . . , T and k = 1, . . . , K p(c tk |w tk , w t+1,k ) ∝ p(c tk |w tk )p(w t+1,k |c tk ) where p(c tk |w tk ) = Poisson(c tk ; φw tk ) Gamma(w t+1,k ; c tk , τ + φ) if w tk > 0 Hence we can have the following MH update. If w t+1,k > 0, then we necessarily have c tk > 0. We sample c * tk ∼zPoisson(φw tk ) where zPoisson(φw tk ) denotes the zero-truncated Poisson distribution and accept c * tk w.p. min 1, Gamma(w t+1,k ; c * tk , τ + φ) Gamma(w t+1,k ; c tk , τ + φ)

	and
	p(w t+1,k |c tk ) =	δ 0 (w t+1,k )	if w tk = 0

1.b) Sample (c, c * ) given (w, w * , φ, α)

RR n°8140

Acknowledgements

YWT thanks the Gatsby Charitable Foundation for generous funding.

RR n°8140