A Martínez-González

C F Méndez-Barrios

S I Niculescu

J Chen

Insights in Characterizing Asymptotic Behavior for Quasipolynomials with Two Delays

Keywords: Delay Systems, Weierstrass Polynomial, Newton-Diagram, Puiseux Series

This paper proposes an analytic method to characterize the behavior of multiple critical roots of a retarded system with two delays. Expressing locally the related characteristic function, as a Weierstrass polynomial, we derive several results to analyze the stability behavior of such characteristic roots with respect to small variations on the delay parameters. The proposed results are illustrated by considering several numerical examples.

I. INTRODUCTION

It is well recognized that one of the main issues in dynamical systems are those that concern to the stability analysis (see, for instance, [START_REF] Hale | Introduction to Functional Differential Equations[END_REF]). In particular, for the case of Linear Time-Invariant (LTI) delay systems, because of the infinite dimensional state space, such an analysis is not a simple task, even in the case of fixed-delays [START_REF] Michiels | Stability, Control, and Computation for Time-Delay Systems. An Eigenvalue-Based Approach[END_REF]. Furthermore, it is well known that delay systems (and in consequence, quasipolynomials) have always infinitely many solutions (see, for further details, [START_REF] Gu | Stability of time-delay systems[END_REF] and the references therein), however in general, we will only be interested in analyzing the behavior of a critical zero of finite multiplicity. In [START_REF] Chen | An eigenvalue perturbation approach to stability analysis, part II: when will zeros of time-delay systems cross imaginary axis?[END_REF], for a general retarded LTI delay system with multiple commensurate delays, the stability behavior of the solutions has been explored in detail. More precisely, they have first, fully characterized the stability properties of such systems by proposing conditions to find the set of critical delay values, at which the system's characteristic quasi-polynomial has critical zeros on the imaginary axis. Secondly, considering the delay as a variable parameter and by adopting an operator based-approach (see, for instance, [START_REF]An eigenvalue perturbation approach to stability analysis, part I: eigenvalue series of matrix operators[END_REF]) they have expanded the solutions of the quasi-polynomial in terms of a Taylor (or Puiseux) series, allowing them to analyze the solutions behavior as the delay varies around a critical delay value. Multiple delays can also be presented in a more general form, as non-commensurate delays, that is, all delays are assuming to be independent of each other. In this case, the stability of the related quasi-polynomial becomes more complex and is less studied. For the case of solutions of multiplicity two, in [START_REF] Irofti | Some insights into the migration of double imaginary roots under small deviation of two parameters[END_REF] the authors have proposed two sectors to analyze the qualitative properties of these roots when delays are By taking into account the above arguments, the main goal of this paper is threefold. First, for multiple imaginary roots, give conditions to guarantee the existence of a convergent Puiseux (or Taylor) series. Second, extend the use of the well known Newton diagram to the case of two parameters and give a methodology to express the solutions as a Puiseux series. Finally, to characterize the solutions according to their geometric behavior.

The remaining part of the paper is organized as follows: Section II introduces some preliminary results, motivating examples and the problem formulation. Section III is devoted to the main results; specifically we present a method to compute the Weierstrass Polynomial, an algorithm to extend the Newton diagram procedure, necessary conditions to derive the generalized Puiseux series and the splitting properties around multiple roots are described. Finally, Section IV includes some numerical examples illustrating the proposed results. The contribution ends with some concluding remarks.

Notations: In the sequel, the following notations will be adopted: C is the set of complex numbers, i := √ -1. Next, R + denotes the set of positive real values. The order of a power series f (x) = i,j a i,j x i 1 x j 2 will be denoted by ord (f ) and defined as the smallest number n = i + j such that a i,j = 0. The ring of complex formal power series is denote by C[[x]], with subring C{x} of convergent power series. Finally, given two polynomials

f (z) = n j=0 a n-j z j and g(z) = m j=0 b m-j z j , the resultant of f , g is defined as R (f, g) := det               a 0 a 1 a 2 • • • a n a 0 a 1 • • • • • • a n . . . . . . . . . . . . a 0 a 1 • • • • • • a n b 0 b 1 b 2 • • • b m b 0 b 1 • • • • • • b m . . . . . . . . . . . . b 0 b 1 • • • • • • b m              
.

II. PRELIMINARIES

A. Retarded Linear Time-Invariant Systems

Consider a retarded LTI system with h a -delays τ k , described in the state-space form as:

ẋ (t) = A 0 x (t) + 2 k=1 A k x (t -τ k ) , τ k ≥ 0, (1)
or by the differential-difference equation

y (n) (t) + n-1 =0 2 k=0 a k y ( ) (t -τ k ) = 0, τ k ≥ 0, (2)
with characteristic function given by the quasi-polynomial:

f (s, τ ) = h k=0 p k (s)e -τ k s , τ k ≥ 0, (3)
where τk is a linear combination of the delays τ j such that τk = j c kj τ j , and the polynomials p k are given by

p 0 (s) = s n + n-1 =0 a 0 s , p k (s) = n-1 =0 a k s .

B. Local Representation of Analytic Functions

It it possible to reduce the analytic properties of f (x, y) to algebraic ones. To this purpose, let us consider the following result (for further details see [START_REF] Mailybaev | On the Weierstrass Preparation Theorem[END_REF]).

Theorem 1 (Weierstrass Preparation Theorem): Let f(z, x) be an analytic function vanishing at the singular point z 0 ∈ C, x 0 ∈ C n , where z = z 0 is an m-multiple root of the equation f (z, x) = 0, i.e.,

f (z 0 , x 0 ) = ∂f ∂z = • • • = ∂ m-1 f ∂z m-1 = 0, ∂ m f ∂z m = 0.
where derivatives are evaluated at (z 0 , x 0 ). Then, there exist a neighborhood U 0 ⊂ C n+1 of the point (z 0 , x 0 ) ∈ C n+1 in which the function f (z, x) can be expressed as

f (z, x) = W (z, x) b (z, x) , (4)
where W (z, x) is given by

(z -z 0 ) m + w m-1 (x) (z -z 0 ) m-1 + • • • + w 0 (x) ,
and w 0 (x),. . . ,w m-1 (x), b (z, x) are analytic functions uniquely defined by the function f (z, x), and

w i (x 0 ) = 0, b (z 0 , x 0 ) = 0. Remark 1: The holomorphic function W (z, x) = z m + w m-1 (x) z m-1 + • • • + w 0 (x) , (5)
is known as the Weierstrass polynomial (for further details on Weierstrass polynomials, see, for instance, [START_REF] Hormander | An introduction to complex analysis in several variables[END_REF]). Remark 2: It can be seen from Theorem 1, that since b (z, x) is an holomorphic non vanishing function at (0, 0), then, there must exist some neighborhood U ⊂ C n at which b(z, x) preserves the same property. Hence, based on this observation we can ensure that the roots behavior of a given quasi-polynomial f in the neighborhood U will be completely described by the roots behavior of W (x, x).

C. Newton Diagram Method

It is well known that solutions of the equation f (x, y) = 0 can be compute term by term by means of the Newton Diagram Method. Thus, in order to use such a procedure, let us introduce the following notation (for more details, see, for instance, [START_REF] Vainberg | Theory of branching of solutions of non-linear equations[END_REF]). Let f (x, y) be a pseudo-polynomial in y, i.e.,

f (x, y) = n k=0 a k (x)y k , (6)
where the corresponding coefficients are given by,

a k (x) = x ρ k ∞ r=0 a rk x r/q , (7)
a rk ∈ C, x and y are complex variables, ρ k are nonnegative rational numbers, q is an arbitrary natural number, a n (x) ≡ 0, and a 0 (x) ≡ 0.

Since by simple translation, any point on a curve can be moved to the origin, we will consider expansions of the solution of ( 6) f (x, y) = 0 around the origin, in the following form

y(x) = y 1 x 1 + y 2 x 2 + y 3 x 3 + • • • , (8)
where 1 < 2 < 3 < • • • and y 1 = 0. To determine the possible values of 1 , y 1 , 2 , y 2 , . . ., it is necessary to consider the Newton's diagram. Definition 2.1 (Newton's Diagram and Polygon): Given a pseudo-polynomial of the form [START_REF] Irofti | Some insights into the migration of double imaginary roots under small deviation of two parameters[END_REF] with coefficients given by ( 7), plot k versus ρ k for k = 0, 1, . . . , n (if a k (•) ≡ 0, the corresponding point is disregarded). Denote each of these points by π k = (k, ρ k ) and let

Π = {π k : a k (•) = 0} ,
be the set of all plotted points. Then, the set Π will be called the Newton diagram, and the Newton polygon associated with f (x, y) will be given by the lower boundary of the convex hull of the set Π.

For a given f (x, y), Fig. 1 simply illustrates Definition 2.1. Thus, the leading term of the solutions have exponents given by the slopes γ by making = γ. The coefficients are given by the non-zero solutions of the polynomial equation

P(y ) := i a 0,i y i = 0, (9)
where the sum runs over the terms satisfying ρ k + γk = ν with constant ν ∈ Q. For equations f (x, y) = 0, the Newton Diagram Method can be formalized by the following theorem ( see [START_REF] Wall | Singular points of plane curves[END_REF]). Theorem 2 (Puiseux Theorem): The equation f (x, y) = 0, with f given in formal power series such that f (0, 0) = 0, posses at least one solution in power series of the form:

x = t q , y = ∞ i=1 c i t i , q ∈ N.

D. Generalized Puiseux Series and Cones

Theorem 2 allows to find fractional power series solutions, known as Puiseux series, in the case of multiparameter algebraic equations. In this vein, when we deal with singularities of greater dimension, we must use a ring of multi-variable fractional power series. In [START_REF] Mcdonald | Fiber polytopes and fractional power series[END_REF] the author defines fractional power series ring that contains the solutions of algebraic hyper-surfaces. This can be achieved through formal power series defined in a geometric way, by taking infinite power series

∞ i=1 c ai x ai/d , where x a = x a1 1 • • • x an n ,
where the exponents a are taken from a fixed convex cones with structure related to its Newton polytopes [START_REF] Ziegler | Lectures on polytopes[END_REF].

Definition 2.2: A convex polyhedral cone is a set of the form σ = m i=1 λ i v i : λ i ∈ R + ,
where M = {v 1 , . . . , v m } is a finite set of vectors in R n . We will use fractional iterated power series of several variables as Generalized Puiseux Power Series (see [START_REF] Soto | The Newton Procedure for several variables[END_REF], [START_REF] Neumann | On Ordered Division Rings[END_REF]), denoted by K x,d . This series can be constructed by induction, taking as a bases the univariate case K x1,d and then, proceed with the field of power series in x

1/d 1 with power series coefficients in x 1/d 2 • • • x 1/d n such that K x,d = C x 1/d 1 • • • x 1/d n .

E. Motivating Examples

Even though we can reduce the analysis of a given entire function f to the study of an algebraic function P f , in this section we aim to pointing out some difficulties that arise in regarding multiparameter functions. In order to illustrate such arguments, let us consider the following motivating example. Example 2.1: Consider the following polynomial

P (z, 1 , 2 ) = z 2 + 3 1 z + 2 2 1 + 2 2 2 , (10)
where 1 and 2 are considered as perturbation parameters.

It is clear to see, that for 1 = 2 = 0, z = 0 is a root of multiplicity two. In this case, the solutions z 1,2 ( ) are not analytic at := ( 1 , 2 ) = (0, 0) = 0. Furthermore, z 1,2 ( ) does not have a unique representation as a power series which is convergent in some punctured neighborhood of the origin. In order to illustrate this assertion, let us consider the region

| 1 | < | 2 |,
in this region the solutions admit the following representation

z 1,2 ( ) = - 1 2 (3 1 ± i4 2 ) + 1 16 1 ±i 1 2 ± i 64 1 2 3 + ± i 2048
z k ( ) = -2 k-1 1 +(-1) k 4 2 2 1 +4 2 1 3 + 32 2 1 5 + O 1 2 5
.

The above arguments clearly have shown that some further considerations must be taken into account in the case of multiparameter functions. Next, as mentioned in previous sections, in the single parameter case, the Newton diagram is a powerful tool to analyze the asymptotic behavior for the solutions of pseudo -polynomials. However, in order to be able to apply such a procedure to the multiparameter case, some special situations must be taken into consideration. In order to motivate the above arguments, let us consider the following.

Example 2.2: Consider the polynomial

P (z, ) := z 5 + 1 3 2 + 2 1 2 2 z 3 + 2 1 2 2 + 3 1 2 z 2 + 4 1 2 . (11)
Clearly, z = 0 is a 5-multiple root at = (0, 0). Now, let us form the Newton diagram with respect to 1 , obtaining Π = {(0, 4), (2, 2), (3, 1), (5, 0)}, illustrated in Fig. 2-(a). The slope β 0 = 1 determines 3-solutions with respect to 1 , and coefficients that are solutions of the polynomial

P(ξ) = 2 + 2 2 ξ 2 + 3 2 ξ 3 = 0.
In this case, it is clear that the solutions cannot be easily computed. In order to compute solutions by applying the Newton procedure, we seek for a monic polynomial. Thus, with the aim of overcoming such difficulty, let us consider the change of variables (blowing-ups)

ζ; = ξ 1 = v 1 2 2 = v 2 ,
and

v 1 = y 1 v 2 v 2 = y 2 .
These change of variables enable us to avoid horizontal segments in the subsequent steps of the process. The resulting polynomial P (y 1 , y 2 ) posses the same Newton polygon, and the segment with β 0 = 1 has a monic polynomial y -5 2 P (ξ, y 2 ) = y 4 2 + y 2 ξ 2 + ξ 3 = 0. Applying Newton procedure, (see figure 2-(b)) to the above equation, we derive the fractional power series solution of P :

ζ 1 (y 1 , y 2 ) = -y 2 y 1 + o(y 1 y 2 ), ζ 2,3 (y 1 , y 2 ) = ±iy 2/3 2 y 1 + o y 1 y 1/3 2 .

F. Problem Formulation

The present work is focused on computing the first approximation of the solution of quasi-polynomials around multiple imaginary roots. In this vein, we will focus in the following problems:

(i) compute an approximation of the associated Weierstrass polynomial; (ii) extend the Newton diagram procedure to Weierstrass polynomial of several variables; (iii) obtain conditions that allow obtaining Puiseux series solutions

s(τ 1 , τ 2 ) = c(τ 1/d 2 )τ β 1 + o(τ 1/d 1 τ 1/d 2 )
, where β = α/d and α ∈ N;

(iv) give conditions on f (s, τ ) which describes the splitting properties of its solutions s(τ 1 , τ 2 ): Regular Splitting, Completely Regular Splitting and Non Regular Splitting.

III. MAIN RESULTS

A. Computation of Weierstrass Polynomial

In [START_REF] Mailybaev | On the Weierstrass Preparation Theorem[END_REF], the authors propose a method to compute the Weierstrass polynomial for an holomorphic function. This method is based on its partial derivatives and combinatorial factors related in a recursive way. For the case of holomorphic function, f (z, x) of complex variables with x = (x 1 , x 2 ) and z = 0 a m-multiple root at (x 1 , x 2 ) = (0, 0) the computation is given as follows. The coefficients w i (5) are analytic, w i (0, 0) = 0 and can expressed as convergent power series:

w i (x 1 , x 2 ) = ∞ h1+h2=1 1 h 1 !h 2 ! w i,h x h1 1 x h2 2 ,
where h = (h 1 , h 2 ). It is not difficult to see that the coefficients w i,h can be compute by means of the following partial derivatives

w i,h = ∂ h1+h2 w i ∂x h1 2 ∂x h2 1 (0,0)
.

Since the analytic function locally satisfy f = W b, thus its partial derivatives satisfy the following recursive relations

w i,h = i j=0 α ij F j,h , (12)
F j,h = f j,h - j k=0 h +h =h c j, k; h , h w k,h b j-k,h ,
with h = 0, h = 0 and constant coefficients:

α jj = m! j!f m,0 , α ij = - m! f m,0 i-1 k=j f m+i-k,0 α kj (m + i -k)! , c (j, k; h 1 , h 2 ) = j! (j -k)! 2 s=1 (h s + h s )! h s !h s ! , and for h = 0, k = k + m, b k,h is given by k! (m+k)!   f k ,h - m-1 j=0 h +h =h c k , j; h , h w j,h b k -j,h   .
Since we are only interested in the leading terms of w i , namely a first approximation of the Weierstrass polynomial, we adopt the following notation. Definition 3.1: Let the natural numbers n

(j) i , for i ∈ {0, 1, • • • , m -1} and j = 1, 2, denote the first non zero partial derivative in (z, x 1 , x 2 ) of f , such that the following conditions hold f (0, 0, 0) = ∂ i f ∂z i = • • • = ∂ i+n (j) i -1 f ∂z i ∂τ n (j) i -1 j = 0, ∂ i+n (j) i f ∂z i ∂τ n (j) i j = 0,
with derivatives evaluated at (0, 0, 0). For n

(j) i = ∞ we have derivatives f i,(0,0) = • • • = f i,(0,n i -1) = 0, f i,(0,n i ) = 0,
evaluated at (z, x) = (0, 0, 1), with n i ∈ Z ≥0 . Leading terms of coefficients w i can be easy found up to the n (j) i and n i derivatives, as a first observation we give the following result. Proposition 1: Suppose that the Weierstrass polynomial has first non-zero partial derivative, such that

n (j) i > n (j) i+1 , 0 ≤ i < m and j = 1, 2.
Then, the leading terms of w i (x) are given by

w i (x 1 , x 2 ) = α i,i f i,(n (1) i ,0) x n (1) i 1 + α i,i f i,(0,n (2) i ) x n (2) i 2 + • • • . If n (j) i = ∞ we get w i (x 1 , x 2 ) = α i,i f i,(n i ,η) x n i 1 x η 2 + • • • . Remark 3:
There may be a case in which f i,(h1,h2) (0,0,0) = 0 ∀ h 1 , h 2 ∈ N, since w i are analytic functions, this is equivalent to w i (x) ≡ 0 for i ≤ i ≤ κ -1. Thus, according to Theorem 1 f has the following local structure:

z κ z m-κ + w m-κ (x)z m-κ-1 + • • • + w κ (x) b(z, x).
Thus, there are κ-invariant solutions z = 0 for all x. If such number κ does not exist (i.e., if such situation does not occur), then κ will be simply defined as κ := 0.

B. The Newton Diagram Method for Two Parameters

Consider the monic pseudo-polynomial f (z, x) of the form

z m + a m-1 (x 1 , x 2 )z m-1 + • • • + a 0 (x 1 , x 2 ), (13)
with a i (x) ∈ C [[x]], such that f (0, 0) = 0. The equation f = 0 can be solved by applying the Newton diagram method, this is done taking into account just one variable, say x 1 , and proceeding iteratively. We take the point π k as the order of a k in x 1 , taking x 2 as an element of C ((x 2 )).

For such a purpose, the following definition will be useful

ρ k := ord x1 (a k (x 1 , x 2 )) = ord (a k (x 1 , 1)) . (14)
Then, the Newton Polygon of f (z, x), with respect to x 1 , is defined by the lower boundary of the convex hull of the points (k, ρ k ) ∈ Π (see, Definition 2.1). In order to apply the the Newton diagram procedure, according to Section II-D, the solution z will take the following structure

z(x 1 , x 2 ) = i c i (x 2 )x i/d 1 ,
where the coefficient c i (x 2 ), is in general, given by an univariate Puiseux series in x 2 .

1) First

Step into the Newton Procedure: Let us suppose that we have determined the Newton diagram of the Weierstrass polynomial (13) of f . Since we are dealing with a monic polynomials, the Newton polygon has a finite number of segments, each one with a corresponding set of points Π ( ) and rational numbers β ≥ 0 satisfying

β 0 > β 1 > • • • > β r .
Therefore, the segments are presented in two possible ways. The first one, a Newton polygon with an horizontal segment determined by β r = 0 and the second one, where β r > 0 . In this vein, for 0 ≤ κ < m, the Newton Diagram Π is given as the set Π = Π Π {(0, ρ 0 ), . . . , (κ, 0)} {(κ, 0) , . . . , ( , ρ ), . . . , (m, 0)} .

Lets take at the first step of the process a horizontal segment with slope β r = 0. Proposition 2: Let f (z, x) be a pseudo-polynomial with the same structure as [START_REF] Ziegler | Lectures on polytopes[END_REF]. Suppose that at least one coefficient a i (x) posses order ρ i = 0. Then, the equation P(x 1 , ξ) = 0 (9) of the corresponding horizontal segment has solutions c k (x

1/d
2 ) in the form of Puiseux series. Now, at the first step of the process, the case with negative slope is considered. Proposition 3: Assume that f has the same structure than [START_REF] Ziegler | Lectures on polytopes[END_REF] and assume that the first Newton diagram posses a segment with negative slope. Then, there exist a change of variables (z, x 1 , x 2 ) → (ζ, y 1 , y 2 ) such that the polynomial P(y 2 , ξ) has Puiseux series solutions c k (y

1/d 2 ).
Hence, applying to f the change of variables z = ζ, x 1 = y a1 1 and x 2 = y a2 2 we get f (ζ, y 1 , y 2 ), which can be solved. Therefore, solutions z(x 1 , x 2 ) of f = 0 are obtained by applying the inverse change of variables to solutions ζ.

The following theorem allow us to use the iterative Newton procedure described above. Theorem 3 (See [START_REF] Soto | The Newton Procedure for several variables[END_REF].): The iteration of the classical Newton Procedure for one variable gives rise to representation of all the roots of the equation ( 13) by generalized Puiseux series with terms x a/d , d ∈ Z + , such that a belong to ndimensional, lex-positive strictly convex polyhedral cone.

C. Newton Polygon Algorithm

Let us consider the points π = ( , ρ ) ∈ Π to get the Newton polygon, obtaining a finite number of segments with slopes -β r . Now, taking as a basis the Newton procedure introduced in Section II-C, we propose the Algorithm 1. In

Algorithm 1 Auxiliary Puiseux Series Expansion Let f (s, τ ) have a critical pair such that s * = iω * is a m- multiple root at τ = (τ * 1 , τ * 1 )
. Consider the initial values as r := 0, i -1 := κ, j := κ and k := ρ κ . 1) Set E r := -k j-i : (i, ) ∈ Π, and i > j ; 2) Let β r := max E j and Π (r) := (i, ) ∈ Π : r) ; 4) Set j := i r , m r := i r -i r-1 and r = r + 1; 5) If j < m go to step 1. Otherwise the algorithm ends. the algorithm 1, κ is defined according to Remark 3. In order to used the iterated Newton diagram procedure we need to find solutions of the equation P, which can be solved by the usual Newton diagram method. For this purpose the input Π of algorithm, has points π determined by equation P = 0.

β j ≡ k- j-i ; 3) Set (i r , r ) ∈ Π (r) such that i r ≥ i, ∀(i, ) ∈ Π (

D. Puiseux Series for Quasi-Polynomials with two delays

Since any critical solution (s * , τ * 1 , τ * 2 ) can always be translated to the origin by appropriate shifts s → s -

s * , τ 1 → τ 1 -τ * 1 , τ 2 → τ 2 -τ * 2
, hereinafter we will assume that (s * , τ * 1 , τ * 2 ) = (0, 0, 0). Proposition 4: Consider the following quasi-polynomial

f (s, τ 1 , τ 2 ) = p 0 (s) + p 1 (s)e -sτ1 + p 1 (s)e -sτ2 , (15)
with s = 0 a m-multiple root at τ = (0, 0) and local representation f (s, τ

) = W (τ )b(s, τ ). If n (j) i = 0 for i = 0, 1, . . . , k then, the k + 1 coefficients of the Weierstrass polynomial W satisfy w m-i (τ ) ≡ 0, i ∈ {0, 1, . . . , k}. Proposition 5: Let quasi-polynomial f (s, τ ) have a m- multiple roots s = 0 at τ = (0, 0), with associated Weierstrass polynomial W . Assume that R W, ∂W ∂s = τ a1 1 τ a2 2 U(τ 1 , τ 2 ) such that U(0, 0) = 0, where (a 1 , a 2 ) ∈ Z 2 ≥0 \ {0}, U ∈ C{τ 1 , τ 2 }.
Then, f = 0 posses m solutions given by a generalized Puiseux series. The following result gives conditions to have a regular Newton diagram. Proposition 6: Let W (s, τ 1 , τ 2 ) be the Weierstrass polynomial of a given quasi-polynomial f (s, τ 1 , τ 2 ). Assume that for a given -segment of the Newton diagram, be -β < 0 its slope with corresponding points Π ( ) = {(k 1 , ρ k1 ), (k 2 , ρ k2 ), . . . , (k s , ρ ks )}. Then, the equation P can be solved without any change of variables if the leading terms of the coefficients w ki satisfy w ki,(ρ k i ,η k i ) = 0 whenever η ki > η ks , i < s.

. Finally, using iterated Newton diagram procedure together with the Weierstrass polynomial of the quasi-polynomial f (s, τ ), we can find the leading terms of the power series solutions.

Proposition 7: Let s * = iω * be a m-multiple root of f (s, τ ) at τ * = (τ * 1 , τ * 2 )
. Assume that κ = 0 and r, β j , (i j , j ), m j and Π (j) , for j = 0, 1, . . . , r -1 are given by Algorithm 1. Then, at τ = τ * the m-zeros of f (s, τ ) can be expanded as

s jσ (τ ) = iω * +c jσ (τ 2 ) (τ 1 -τ * 1 ) βj +o |τ 1 -τ * 1 | βj |τ 2 -τ * 2 | β j ,
for j = 0, 1, . . . , r -1, σ = 0, . . . , m j and m = m 0 + • • • + m r-1 . For β j > 0 the coefficients c jσ (x 2 ) are roots of the polynomial:

P j (ξ, x 2 ) = ij k=ij-1 w k,(n (1) k ,n k ) τ n k 2 ξ k-ij-1 , k, n (1) k
∈ Π (j) , when β r-1 = 0, the coefficients are given by the solution of

P j (ξ, x 2 ) = ij k=ij-1 w k,(0,n (2) k ) τ n (2) k 2 ξ k-ij-1 , (k, 0) ∈ Π (r-1)
where n

(1)

k , n (2)
k , n k are given by the first non-zero partial derivatives of Definition 3.1; the constant terms w k,(n,η) ∈ C are computed using [START_REF] Mcdonald | Fiber polytopes and fractional power series[END_REF].

E. Splitting Properties

The main goal of this subsection is to explore some qualitative properties of the solutions s(τ 1 , τ 2 ) of the quasipolynomial f around the m-multiple critical pair (0, 0). The solution curve C ∈ C 3 , defined by the equation f = 0, is composed of m generalized Puiseux series solutions. These solutions can be arranged in r-branches, given in the following form

s jσ (τ ) = c jσ τ βj 1 τ β j 2 + o |τ 1 | βj |τ 2 | β j , (16)
with σ = 0, . . . , m j -1 and j = 0, . . . , r -1. The leading term has exponents defined by β j = p/m j , β j = q/m j , with p, q ∈ Z ≥0 . Each branch has multiplicity

m j = m 0 + m 1 + • • • + m , for some ∈ Z ≥0 such that m = m 0 + m 1 + • • • + m r-1 .
In this vein, we have the following definition. Definition 3.2: We say that there is a Completely Regular Splitting (CRS) property of the solution s * = 0 at τ * = 0 if c jσ = 0, and p • q ≤ 1, ∀j. For the Regular Splitting (RS) property, some of the coefficients c jσ for which m j = 1 may be equal to zero. In the remaining cases of the coefficient c jσ we say that Non Regular Splitting (NRS) property is present.

In this view, we have the following proposition. Proposition 8: Let f (s; τ ) be a quasi-polynomial with s * = iω a m-multiple root at (τ 1 , τ 2 ) = (τ * 1 , τ * 2 ) such that satisfy Proposition 5. Suppose that β j , m j , for j = 0, 1, • • • , r -1 are given by Algorithm 1. Then, the following relations holds:

1) if (m j •β j ) • m j •β j ≤ 1 ∀j ∈ {0, 1, . . . ,
f (s, τ ) = s 2 -2s+1 -2e -sτ1 +2πse -sτ2 +e -2sτ2 , ( 17
)
with τ = (1, π), we have a triple root at s = 0. In order to apply the proposed results, let us consider f (s, τ ) := f (s, τ 1 + 1, τ 2 + π). Now, according to the Weierstrass Preparation Theorem, the local behavior around the solution 0 of f is captured by the solutions of W (s, τ ) = s 3 + w 2 (τ ) s 2 + w 1 (τ ) s + w 0 (τ ) .

Considering Definition 3.1, we take the first partial derivatives of f . We have that n determined by

f 1,hj = (-1) j+1 2 ⇒ n (j) 1 = 1, j = 1, 2.
Similarly, for n (j) 2 we have that

f 2,h1 = -4, f 2,h2 = 4π ⇒ n (j) 2 = 1.
Applying Proposition 4 we have that w 0 ≡ 0. Proposition 1 allows to find the leading terms given by

w 1 (τ ) = 6 1 -π 3 τ 1 - 6 1 -π 3 τ 2 + • • • w 2 (τ ) = -3 3 + 4π 3 (π -1) 2 (π 3 -1) 2 τ 1 + 3 (4π -1) 2 (π 3 -1) τ 2 + • • • .
Next, according to [START_REF] Soto | The Newton Procedure for several variables[END_REF], the order with respect to τ 1 of w i (τ ) define the set of points Π = {(1, 0), (2, 0), (3, 0)}. Hence, applying Algorithm 1, we derive the results summarized in Table I. From Table I we can observe that the Newton polygon has an horizontal slope with β 0 = 0. Thus, by considering Proposition 7, we need solve the polynomial equation P 0 (ξ, τ 2 ) := ξ 2 + w 2,(0,1) τ 2 ξ + w 1,(0,1) τ 2 = 0,

= 3, κ = 1. ρ 1 = 0 r = 1, m 0 = 2, β 0 = 0 Π = {(1, 0), (2, 0), (3, 0)} Π (0) = {(1, 0), (2, 0), (3, 0)} TABLE II RESULTS SUMMARY FOR THE QUASI-POLYNOMIAL 17. Initial Data Algorithm Output Z := {z ∈ C : Pj (ξ) = 0} m = 2, κ = 0, ρ0 = 1 r = 1, m 0 = 2 β 0 = 1/2 P0(ξ) := ξ 2 + 6 1-π 3 Π = {(0, 1), (1, 1), (2, 0)} Π (0) = {(0, 1), (2, 0)} c0,σ = ± 6 1-π 3
and by Proposition 6 we know that there is no need of change of variables. Taking as input Π = {(0, 1), (1, 1), (2, 0)}, we obtain Table II from the Newton polygon algorithm. Since κ = 1, there is an invariant root at s = 0. The leading terms of solutions are given by

s 0,σ (τ ) = ± 6 1 -π 3 (τ 2 -π) 1/2 +o |τ 1 -1| 1 2 |τ 2 -π| 1 2
.

where σ ∈ {0, 1}. We can see that the solutions (0, 1, π) splits in one invariant roots s 0 and two solutions such that β 0 •m 1 = 1, thus by Proposition 8, we know that the solution possess the CRS property. Example 4.2: Lets consider the quasi-polynomial f (s, τ ) = p 0 (s) + p 1 (s)e -sτ1 + p 2 (s)e -sτ2 where

p 0 (s) = s 5 +s 4 + 4 + π 2 s 3 +2s 2 + 2 + π 2 s+2, (18a)
p 1 (s) = 1, p 2 (s) = 2s 4 +4s 2 + 2. (18b)
For (τ 1 , τ 2 ) = (π, 1), f has a double root at s = i. With the aim of applying the proposed methodology, let us shift from (i, π, 1) to the origin, obtaining f . We find that the first non-zero partial derivatives of the quasi-polynomial at (0, 0, 0) are given by

f 0,(1,0) = i ⇒ n (1) 0 = 1, f 0,(0,n) = 0 ⇒ n (2) 0 = ∞, f 1,(1,0) = 1 -iπ ⇒ n (1) 1 = 1, f 1,(0,n) = 0 ⇒ n (2) 1 = ∞.
Hence, by Proposition 1, we have that ρ j = n

(1) j for j = 1, 2. The Newton diagram is given by Π = {(0, 1), (1, 1), (2, 0)}. Table III summarizes the results deriving from Algorithm 1. Since coefficient w 1 is not over the Newton polygon, following Proposition 1, we compute the leading term of w i (τ ), as w 0 (τ ) = -2i (8 + π 2 ) + i(8 -3π) + 16e -i τ 1 + • • • .

From the algorithm output, we get a segment with negative slope β 0 = 1/2. According to Proposition 7 we need to solve: P(ξ) = ξ 2 -2i (8 + π 2 ) + i(8 -3π) + 16e -i = 0, for σ ∈ {0, 1}, the solutions are given by s 0,σ (τ ) = i ± √ 2i 3/2 (8+π 2 )+i(8-3π)+16e -i (τ 1 -π) .

Since β 0 • m 0 = 1, Proposition 8 implies that the solution (i, π, 1) has the CRS property. This behavior is illustrated in Figure 3. Example 4.3: Lets consider the special case of a model of population dynamics, see [START_REF] Irofti | Some insights into the migration of double imaginary roots under small deviation of two parameters[END_REF], with two distributed delays such that the characteristic functions is given by f (s, τ 1 , τ 2 ) = s -a 1 -e -sτ1 s -c 1 -e -sτ2 s ,

with positive parameters a = -0.214104, b = -0.996801. At (τ * 1 , τ * 2 ) ≈ (3.84003026849, 10.44866732901) the quasipolynomial (19) posses a double root at s * = i. We make a change of variable in order to shift the critical point to origin, obtaining f . The associated Weierstrass polynomial W (s, τ ) = s 2 + w 1 (τ ) + w 0 (τ ) can be computed, as a first approximation, by means of the first non zero partial derivatives using definition 3.1, given by f 0,(1,0) = 0 ⇒ n

Fig. 1 .

 1 Fig. 1. The Newton Diagram for f (x, y).

.

  Now, if instead of the previous region, we consider the region | 2 | < | 1 |, then for k ∈ {1, 2} the solutions admit the following representation

Fig. 2 .

 2 Fig. 2. Newton polygons for P (z, ) in Example 11.

=

  ∞, for h 1 = (1, 0), h 2 = (0, 1) we get the natural numbers n (j) 1

  Initial DataAlgorithm OutputZ := {z ∈ C : Pj (ξ) = 0} m = 2, κ = 0. ρ0 = 1 r = 1, m0 = 2, β0 = 1/2 P0(ξ) := ξ 2 + w 0,(1,0) Π = {(0, 1), (1, 1), (2, 0)} Π (0) = {(0, 1), (2, 0)} c0,σ = ± √ w 0,[START_REF] Hale | Introduction to Functional Differential Equations[END_REF]0)

1 Fig. 3 .

 13 Fig. 3. Root locus of quasi-polynomial f (s, τ ) (18) around (i, π, 1).

1

 1 Now we are able to compute the first terms w i (τ ) following Proposition 1.w 0 (τ ) = ατ 1 + βτ 2 + • • • , w 1 (τ ) = γτ 1 + δτ 2 + • • • ,with α = (-0.0231801 + 0.0439373i), β = (0.00397688-0.0245988i), γ = (0.0105549+0.0292551i), δ = (-0.0359858 + 0.150344i). Applying Algorithm 1, to the Newton diagram Π = Π Π = {(0, 0), (0, 1), (1, 0), (1, 1), (2, 0)}, we get TableIV. By Proposition 7 the solutions of the polynomial

  r -1}, then the solution (iω, τ * ) has the CRS property; 2)if the pairs (m k , β k ) that do not fulfill 1), satisfy β k ≥ m k ≡ 1, then the solution (iω, τ * ) has the RS property; 3)for the remaining cases of β j , the solution (iω, τ * ) has the NRS property. Remark 4: It should be mentioned that the Definition 3.2 takes into account the behavior of m -κ solutions, without considering κ invariant solutions.

	IV. NUMERICAL EXAMPLES
	In order to illustrate the proposed approach, let us consider
	the following numerical examples.
	Example 4.1: Consider the following quasi-polynomial

TABLE I RESULTS

 I SUMMARY FOR THE QUASI-POLYNOMIAL 17.

	Initial Data	Algorithm Output
	m

TABLE III RESULTS

 III SUMMARY FOR THE QUASI-POLYNOMIAL 18.

Universidad Autónoma de San Luis Potosí (UASLP), Facultad de Ingeniería, Dr. Manuel Nava No.8, San Luis Potosí, S.L.P., México jano.mtz.glz@gmail.com,cerfranfer@gmail.com

Laboratoire des Signaux et Systèmes (L2S, UMR CNRS 8506), CNRS-Supélec,

3, rue Joliot Curie, 91192, Gif-sur-Yvette, France Silviu.Niculescu@lss.supelec.fr[START_REF] Gu | Stability of time-delay systems[END_REF] Department of Electronic Engineering, City University of Hong Kong, Hong Kong, China jichen@cityu.edu.hk subject to small deviations and restricted to such sectors. In this vein, it is worth to mention that unlike the case of a single parameter, in the multiparameter case there exist some singular and unexpected behaviors (see, the motivating examples section) which have to be taken into account (see, for instance[START_REF] Monforte | Formal laurent series in several variables[END_REF]), in order that the problem is well-posed. In other words, the Puiseux type arguments cannot be extended straightforwardly from one parameter to multiparameter case.

ACKNOWLEDGMENT

The research of J. Chen was supported in part by the Hong Kong RGC under Projects F-HK006/11T and CityU 11260016, the work of Martínez-González is financially supported by CONACyT, Mexico, and the work of Méndez-Barrios has been supported by Centrale Supélec, France and partially by COPOCyT, Mexico.

m = 2, κ = 0. ρ0 = 0 r = 1, m0 = 2, β0 = 0 Π = Π Π Π (0) = {(0, 0), (1, 0), (2, 0)} P(ξ) = ξ 2 + w 1,(0,1) ξ + w 0,(1,0)

determine the leading terms of the solutions. Taking in consideration the Newton diagram Π

, we obtain Table V. Then, for σ = {0, 1}, the solutions are given by

, by Proposition 8 the solutions posses the CRS splitting property.

V. CONCLUDING REMARKS

In this paper, we have considered some issues concerning the asymptotic behavior of multiple critical roots for quasipolynomials with two delays. By means of the Weierstrass polynomial, the proposed approach is based on an iterative Newton diagram method which can be effectively applied to find the leading terms of power series solutions expressed as a generalized Puiseux series. Finally, the splitting properties of a given solution have been described by means of CRS, RS and NRS properties in order to get some insights about its geometric behavior.


Caroline Lawless

Julyan Arbel

A simple proof of Pitman-Yor's Chinese restaurant process from its stick-breaking representation

For a long time, the Dirichlet process has been the gold standard discrete random measure in Bayesian nonparametrics. The Pitman-Yor process provides a simple and mathematically tractable generalization, allowing for a very flexible control of the clustering behaviour. Two commonly used representations of the Pitman-Yor process are the stick-breaking process and the Chinese restaurant process. The former is a constructive representation of the process which turns out very handy for practical implementation, while the latter describes the partition distribution induced. Obtaining one from the other is usually done indirectly with use of measure theory. In contrast, we provide here an elementary proof of Pitman-Yor's Chinese Restaurant process from its stick-breaking representation.

Introduction

The Pitman-Yor process defines a rich and flexible class of random probability measures which was developed by [START_REF] Perman | Size-biased sampling of Poisson point processes and excursions[END_REF] and further investigated by [START_REF] Pitman | Exchangeable and partially exchangeable random partitions[END_REF], [START_REF] Pitman | The two-parameter Poisson-Dirichlet distribution derived from a stable subordinator[END_REF]. It is a simple generalization of the Dirichlet process [START_REF] Ferguson | A Bayesian analysis of some nonparametric problems[END_REF], whose mathematical tractability contributed to its popularity in machine learning theory [START_REF] Caron | Generalized Pólya Urn for Time-Varying Pitman-Yor Processes[END_REF], probabilistic models for linguistic applications [START_REF] Teh | A hierarchical Bayesian language model based on Pitman-Yor processes[END_REF][START_REF] Wood | The sequence memoizer[END_REF], excursion theory [START_REF] Perman | Size-biased sampling of Poisson point processes and excursions[END_REF][START_REF] Pitman | The two-parameter Poisson-Dirichlet distribution derived from a stable subordinator[END_REF], measurevalued diffusions in population genetics (Petrov, 2009, Feng and[START_REF] Feng | Some diffusion processes associated with two parameter Poisson-Dirichlet distribution and Dirichlet process[END_REF], combinatorics [START_REF] Vershik | On the Markov-Krein identity and quasiinvariance of the gamma process[END_REF][START_REF] Kerov | Coherent random allocations, and the Ewens-Pitman formula[END_REF] and statistical physics [START_REF] Derrida | Random-energy model: An exactly solvable model of disordered systems[END_REF].

Its most prominent role is perhaps in Bayesian nonparametric statistics where it is used as a prior distribution, following the work of [START_REF] Ishwaran | Gibbs sampling methods for stick-breaking priors[END_REF]. Applications in this setting embrace a variety of inferential problems, including species sampling [START_REF] Favaro | Bayesian non-parametric inference for species variety with a two-parameter Poisson-Dirichlet process prior[END_REF][START_REF] Navarrete | Some issues in nonparametric Bayesian modeling using species sampling models[END_REF][START_REF] Arbel | Bayesian nonparametric inference for discovery probabilities: credible intervals and large sample asymptotics[END_REF], survival analysis and graphical models in genetics [START_REF] Jara | Bayesian semiparametric inference for multivariate doubly-interval-censored data[END_REF][START_REF] Ni | Heterogeneous reciprocal graphical models[END_REF], image segmentation [START_REF] Sudderth | Shared segmentation of natural scenes using dependent Pitman-Yor processes[END_REF], curve estimation [START_REF] Canale | On the Pitman-Yor process with spike and slab base measure[END_REF], exchangeable feature allocations [START_REF] Battiston | A characterization of product-form exchangeable feature probability functions[END_REF] and time-series and econometrics [START_REF] Caron | Generalized Pólya Urn for Time-Varying Pitman-Yor Processes[END_REF][START_REF] Bassetti | Beta-product dependent Pitman-Yor processes for Bayesian inference[END_REF].

Last but not least, the Pitman-Yor process is also employed in the context of nonparametric mixture modeling, thus generalizing the celebrated Dirichlet process mixture model of [START_REF] Lo | On a class of Bayesian nonparametric estimates: I. Density estimates[END_REF]. Nonparametric mixture models based on the Pitman-Yor process are characterized by a more flexible parameterization than the Dirichlet process mixture model, thus allowing for a better control of the clustering behaviour [START_REF] De Blasi | Are Gibbs-type priors the most natural generalization of the Dirichlet process? Pattern Analysis and Machine Intelligence[END_REF]. In addition, see [START_REF] Ishwaran | Gibbs sampling methods for stick-breaking priors[END_REF], [START_REF] Favaro | Slice sampling σ-stable Poisson-Kingman mixture models[END_REF], [START_REF] Arbel | Stochastic approximations to the Pitman-Yor process[END_REF] for posterior sampling algorithms, Scricciolo et al. (2014), Miller and[START_REF] Miller | Inconsistency of Pitman-Yor process mixtures for the number of components[END_REF] for asymptotic properties, and [START_REF] Scarpa | Bayesian hierarchical functional data analysis via contaminated informative priors[END_REF], [START_REF] Canale | On the Pitman-Yor process with spike and slab base measure[END_REF] for spike-and-slab extensions.

The Pitman-Yor process has the following stick-breaking representation:

if v i ind ∼ Beta(1- d, α + id) for i = 1, 2, . . . with d ∈ (0, 1) and α > -d, if π j = v j j-1 i=1 (1 -v i ) for j = 1,
P(C = C) = d |C| (α) (n) α d (|C|) c∈C (1 -d) (|c|-1) , (2)
where the multiplicative factor before the product in ( 2) is also commonly (and equivalently) written as ( |C|-1 i=1 α + id)/(α) (n-1) in the literature. When the discount parameter d is set to zero, the Pitman-Yor process reduces to the Dirichlet process and the partition distribution (2) boils down to the celebrated Chinese Restaurant process (CRP, see [START_REF] Antoniak | Mixtures of Dirichlet processes with applications to Bayesian nonparametric problems[END_REF]. By abuse of language, we call the partition distribution (2) the Pitman-Yor's CRP. Under the latter partition distribution, the number of parts in a partition C of n elements, k n = |C|, grows to infinity as a power-law of the sample size, n d (see Pitman, 2003, for details). This Pitman-Yor power-law growth is more in tune with most of empirical data [START_REF] Clauset | Power-law distributions in empirical data[END_REF] than the logarithmic growth induced by the Dirichlet process CRP, α log n.

The purpose of this note is to provide a simple proof of Pitman-Yor's CRP (2) from its stick-breaking representation (1) (Theorem 2.1). This generalizes the derivation by [START_REF] Miller | An elementary derivation of the Chinese restaurant process from Sethuraman's stick-breaking process[END_REF] who obtained the Dirichlet process CRP [START_REF] Antoniak | Mixtures of Dirichlet processes with applications to Bayesian nonparametric problems[END_REF]) from Sethuraman's stickbreaking representation [START_REF] Sethuraman | A constructive definition of Dirichlet priors[END_REF]. In doing so, we also provide the marginal distribution of the allocation variables vector (3) in Proposition 2.2.

Partition distribution from stick-breaking

Suppose we make n observations, z 1 , . . . , z n . We denote the set {1, . . . , n} by [n]. Our observations induce a partition of [n], denoted C = {c 1 , . . . , c kn } where c 1 , . . . , c kn are disjoint sets and kn i=1 c i = [n], in such a way that z i and z j belong to the same partition if and only if z i = z j . We denote the number of parts in the partition C by k n = |C| and we denote the number of elements in partition j by |c j |. We use bold font to represent random variables.

We write (x) (n) = n-1 j=0 (x + j) to denote the rising factorial.

Theorem 2.1. Suppose

v i ind ∼ Beta(1 -d, α + id) for i = 1, 2, . . . , π j = v j j-1 i=1 (1 -v i ) for j = 1, 2, . . .

Let allocation variables be defined by

z 1 , . . . , z n |π = π iid ∼ π, meaning, P(z i = j | π) = π j , (3)
and C denote the random partition of [n] induced by z 1 , . . . , z n . Then

P(C = C) = d |C| (α) (n) α d (|C|) c∈C (1 -d) (|c|-1) .
The proof of Theorem 2.1 follows the lines of Miller ( 2018)'s derivation. We need the next two technical results, which we will prove in Section 3. Let C z denote the partition [n] induced by z for any z ∈ N n . Let k n be the number of parts in the partition. We define m(z) = max {z 1 , . . . , z n }, and g j (z) = #{i : z i ≥ j}.

Proposition 2.2. For any z ∈ N n , the marginal distribution of the allocation variables vector z = (z 1 , . . . , z n ) is given by

P(z = z) = 1 (α) (n) c∈Cz Γ(|c| + 1 -d) Γ(1 -d) m(z) j=1 α + (j -1)d g j (z) + α + (j -1)d . Lemma 2.3. For any partition C of [n], z∈N n 1(C z = C) m(z) j=1 α + (j -1)d g j (z) + α + (j -1)d = d |C| c∈C (|c| -d) α d (|C|) .
Proof of Theorem 2.1.

P(C = C) = z∈N n P(C = C|z = z)P(z = z) (a) = z∈N n 1(C z = C) 1 (α) (n) c∈Cz Γ(|c| + 1 -d) Γ(1 -d) m(z) j=1 α + (j -1)d g j (z) + α + (j -1)d = 1 (α) (n) c∈C Γ(|c| + 1 -d) Γ(1 -d) z∈N n 1(C z = C) m(z) j=1 α + (j -1)d g j (z) + α + (j -1)d (b) = 1 (α) (n) c∈C Γ(|c| + 1 -d) Γ(1 -d) d |C| c∈C (|c| -d) α d (|C|) (c) = 1 (α) (n) c∈C (1 -d) (|c|-1) c∈C (|c| -d) d |C| c∈C (|c| -d) α d (|C|) = d |C| (α) (n) α d (|C|) c∈C (1 -d) (|c|-1) ,
where (a) is by Proposition 2.2, (b) is by Lemma 2.3, and (c) is since Γ

(|c| + 1 -d) = (|c| -d)Γ(|c -d|).
3 Proofs of the technical results

Additional lemmas

We require the following additional lemmas.

Lemma 3.1. For a + c > 0, and b

+ d > 0, if y ∼ Beta(a, b), then E[y c (1 -y) d ] = B(a+c,b+d) B(a,b)
where B denotes the beta function.

Proof.

E[y c (1 -y) d ] = 1 0 y c (1 -y) d 1 B(a, b) y a-1 (1 -y) b-1 dy = 1 B(a, b) 1 0 y a+c-1 (1 -y) b+d-1 dy = B(a + c, b + d) B(a, b) .
Let S kn denote the set of k n ! permutations of [k n ]. The following lemma is key for proving Lemma 2.3. Lemma 3.2. For any n 1 , . . . , n kn ∈ N,

σ∈S kn kn i=1 1 a i (σ) -(k n -i + 1)d = 1 kn i=1 (n i -d)
where

a i (σ) = n σ i + n σ i+1 + • • • + n σ kn .
Proof. Consider the process of sampling without replacement k n times from an urn containing k n balls. The balls have sizes n 1 -d, . . . , n kn -d, and the probability of drawing ball i is proportional to its size n i -d. Thus for any permutation σ ∈ S kn we have that

p(σ 1 ) = n σ 1 -d n -td = n σ 1 -d a 1 (σ) -td , p(σ 2 |σ 1 ) = n σ 2 -d n -n σ 1 -(k n -1)d = n σ 2 -d a 2 (σ) -(k n -1)d , p(σ i |σ 1 , . . . , σ i-1 ) = n σ i -d n -n σ 1 -• • • -n σ i-1 -(k n -i + 1)d = n σ i -d a i (σ) -(k n -i + 1)d . Therefore, p(σ) = p(σ 1 )p(σ 2 |σ 1 ) • • • p(σ kn |σ 1 , . . . , σ kn-1 ) = kn i=1 n σ i -d a i (σ) -(k n -i + 1)d . (4)
This way, we construct a distribution on S kn . We know that σ∈S kn p(σ) = 1. Applying this to Equation ( 4) and dividing both sides by (n

σ 1 -d) • • • (n σ kn -d) = (n 1 -d) • • • (n kn -d)
gives the result.

Lemma 3.3. Let b i ∈ N for i ∈ {1, . . . , k n } and let b 0 = 0. We define bi = b 0 + b 1 + • • • + b i . Then kn i=1 b i ∈N ( α d + bi-1 ) (b i ) ( a i +α d + bi-1 ) (b i ) = ( α d ) (kn) kn i=1 ( a i d -(k n + 1 -i))
.

Proof. Let A j denote the intermediate sum

A j = kn i=j b i ∈N ( α d + bi-1 ) (b i ) ( a i +α d + bi-1 ) (b i )
. We show by induction decreasing from j = k n to j = 0 that

A j = ( α d + bj-1 ) (kn-j+1) kn i=j ( a i d -(k n + 1 -i))
.

(5)

When j = k n we have

A kn = b kn ∈N ( α d + bkn-1 ) (b kn ) ( a kn +α d + bkn-1 ) (b kn ) = b kn ∈N E[X b kn ] where X ∼ Beta( α d + bkn-1 , a kn d ). We have that b kn ∈N E[X b kn ] = E b kn ∈N X b kn = E X 1 -X = α + d bkn-1 a kn -d ,
due to Lemma 3.1, which proves the initialization for (5).

We now consider the case of an arbitrary j, greater than 0 and less than k n . By the induction hypothesis, we have that Equation (5) holds for j + 1, that is

A j+1 = ( α d + bj ) (kn-j) kn i=j+1 ( a i d -(k n + 1 -i))
.

Therefore,

A j = b j ∈N ( α d + bj-1 ) (b j ) ( a j +α d + bj-1 ) (b j ) kn i=j+1 b i ∈N ( α d + bi-1 ) (b i ) ( a i +α d + bi-1 ) (b i ) = b j ∈N ( α d + bj-1 ) (b j ) ( a j d + α d + bj-1 ) (b j ) ( α d + bj ) (kn-j) kn i=j+1 ( a i d -(k n + 1 -i))
Rearranging the rising factorials in the numerator, we can write

α d + bj-1 (b j ) α d + bj (kn-j) = α d + bj-1 (b j ) α d + bj-1 + b j (kn-j) = α d + bj-1 (b j +kn-j) = α d + bj-1 (kn-j) α d + bj-1 + k n -j (b j )
and thus factorize the terms independent of b j in order to obtain

A j = ( α d + bj-1 ) (kn-j) kn i=j+1 ( a i d -(k n + 1 -i)) b j ∈N ( α d + bj-1 + k n -j) (b j ) ( a j d + α d + bj-1 ) (b j )
.

The sum above can be rewritten, using

X ∼ Beta( α d + bj-1 + (k n -j), a j d -(k n -j)), as b j ∈N E[X b j ] = E X 1 -X = α d + bj-1 + (k n -j) a j d -(k n + 1 -j)
.

Putting this all together,

A j = ( α d + bj-1 + (k n -j)) a j d -(k n + 1 -j) ( α d + bj-1 ) (kn-j) kn i=j+1 ( a i d -(k n + 1 -i)) = ( α d + bj-1 ) (kn-j+1) kn i=j ( a i d -(k n + 1 -i))
which proves the desired result for j. By induction, this result is true for all j ∈ {1, . . . , k n }.

Letting j = 1 gives the result stated in the lemma, since b0 = b 0 = 0.

3.2 Proof of Proposition 2.2 and Lemma 2.3

Proof of Proposition 2.2. For simplicity, we fix the allocation variable vector to a value z and denote m(z) by m and g j (z) by g j . We have

P(z = z|π 1 , . . . , π m ) = n i=1 π z i = m j=1 π e j j
where e j = #{i : z i = j}. Thus,

P(z = z|v 1 , . . . , v m ) = m j=1 v j j-1 i=1 (1 -v i ) e j = m j=1 v e j j (1 -v j ) f j
where f j = #{i : z i > j}. Therefore,

P(z = z) = P(z = z|v 1 , . . . , v m )p(v 1 , . . . , v m )dv 1 • • • dv m = m j=1 v e j j (1 -v j ) f j p 1 (v 1 ) • • • p m (v m )dv 1 • • • dv m = m j=1 v e j j (1 -v j ) f j p j (v j )dv j (a) = m j=1 B(e j + 1 -d, f j + α + jd) B(1 -d, α + jd) = m j=1 Γ(e j + 1 -d)Γ(f j + α + jd)Γ(α + (j -1)d) + 1 Γ(e j + f j + α + (j -1)d) + 1)Γ(1 -d)Γ(α + jd) (b) = m j=1 Γ(e j + 1 -d) Γ(1 -d) m j=1 Γ(g j+1 + α + jd) Γ(g j + α + (j -1)d + 1) m j=1 Γ(α + (j -1)d + 1) Γ(α + jd) (c) = m j=1 Γ(e j + 1 -d) Γ(1 -d) m j=1 α + (j -1)d g j + α + (j -1)d Γ(g m+1 + α + md)Γ(α) Γ(g 1 + α)Γ(α + md) (d) = Γ(α) Γ(n + α) c∈Cz Γ(|c| + 1 -d) Γ(1 -d) m j=1 α + (j -1)d g j + α + (j -1)d
where step (a) follows from Lemma 3.1, step (b) since f j = g j+1 and g j = e j + f j , step (c) since Γ(x + 1) = xΓ(x), and step (d) since g 1 = n and g m+1 = 0.

Proof of Lemma 2.3. As before, we denote the parts of C by c 1 , . . . , c kn , and we let k n = |C|. We denote the distinct values taken on by z 1 , . . . , z n by j 1 < • • • < j kn . We define

j 0 = b 0 = 0, b i = j i -j i-1 , and bi = b 0 + • • • + b i for i ∈ {1, . . . , k n }. We use the notation a i (σ) = n σ i + • • • + n σ kn ,

  where σ is the permutation of [k n ] such that c σ i = {ℓ :z ℓ = j i }. Then for any z ∈ N n such that C z = C, -1 ) (b i ) ( α+a i (σ) d + bi-1 ) (b i ) , because g j (z) = a i (σ) for bi-1 < j ≤ bi . Itfollows from the definition of b = (b 1 , ..., b kn ) and σ that there is a one-to-one correspondence between {z ∈ N n : C z = C} and {(σ, b) : σ ∈ S kn , b ∈ N kn }. Therefore, σ) -(k n -i + 1)d a) follows from Lemma 3.3 and step (b) follows from Lemma 3.2.

  2, . . ., and if θ 1 , θ 2 , . . . , and C is the partition of the first n integers {1, . . . , n} induced by data, then

	iid ∼ H, then the discrete random probability measure
	∞
	P =	π j δ θ j	(1)
	j=1
	is distributed according to the Pitman-Yor process, PY(α, d, H), with concentration param-
	eter α, discount parameter d, and base distribution H.
	The Pitman-Yor process induces the following partition distribution: if P ∼ PY(α, d, H),
	for some nonatomic probability distribution H, we observe data x 1 , . . . , x n |P

iid ∼ P

Acknowledgement

The authors would like to thank Bernardo Nipoti for fruitful discussions that initiated this work.


Julyan Arbel

Pierpaolo De Blasi

Igor Prünster

Stochastic approximations

Keywords: stochastic approximation, asymptotic distribution, Bayesian Nonparametrics, Pitman-Yor process, random functionals, random probability measure, stopping rule

come

Introduction

The Pitman-Yor process defines a rich and flexible class of random probability measures used as prior distribution in Bayesian nonparametric inference. It originates from the work of [START_REF] Perman | Size-biased sampling of Poisson point processes and excursions[END_REF], further investigated in [START_REF] Pitman | Exchangeable and partially exchangeable random partitions[END_REF]; [START_REF] Pitman | The two-parameter Poisson-Dirichlet distribution derived from a stable subordinator[END_REF], and its use in nonparametric inference was initiated by [START_REF] Ishwaran | Gibbs sampling methods for stick-breaking priors[END_REF]. Thanks to its analytical tractability and flexibility, it has found applications in a variety of inferential problems which include species sampling [START_REF] Lijoi | Bayesian nonparametric estimation of the probability of discovering a new species[END_REF][START_REF] Favaro | Bayesian non-parametric inference for species variety with a two-parameter Poisson-Dirichlet process prior[END_REF][START_REF] Navarrete | Some issues in nonparametric Bayesian modeling using species sampling models[END_REF], survival analysis and gene networks [START_REF] Jara | Bayesian semiparametric inference for multivariate doubly-interval-censored data[END_REF][START_REF] Ni | Heterogeneous reciprocal graphical models[END_REF], linguistics and image segmentation [START_REF] Teh | A hierarchical Bayesian language model based on Pitman-Yor processes[END_REF][START_REF] Sudderth | Shared segmentation of natural scenes using dependent Pitman-Yor processes[END_REF], curve estimation [START_REF] Canale | On the Pitman-Yor process with spike and slab base measure[END_REF] and time-series and econometrics [START_REF] Caron | Generalized pólya urn for time-varying pitman-yor processes[END_REF][START_REF] Bassetti | Beta-product dependent pitman-yor processes for bayesian inference[END_REF]. The Pitman-Yor process is a discrete probability measure

P (dx) = i≥1 p i δ ξi (dx) (1)
where (ξ i ) i≥1 are independent and identically distributed (i.i.d.) random variables with common distribution P 0 on a Polish space X , and (p i ) i≥1 are random frequencies, i.e. p i ≥ 0 and finite support approximation to the Pitman-Yor process obtained by truncating the summation in (1). A general method consists in setting the truncation level n by replacing p n+1 with 1 -(p 1 + • • • + p n ) in (1). The key quantity is the truncation error of the infinite summation (1),

R n = i>n p i = j≤n (1 -V j ), (3)
since the resulting truncated process, say P n (•), will be close to P (•) according to |P (A) -P n (A)| ≤ R n for any measurable A ⊂ X . It is then important to study the distribution of the truncation error R n as n gets large in order to control the approximation error. [START_REF] Ishwaran | Gibbs sampling methods for stick-breaking priors[END_REF] proposes to determine the truncation level based on the moments of R n . Cf. also [START_REF] Ishwaran | Exact and approximate sum representations for the Dirichlet process[END_REF]; [START_REF] Gelfand | A computational approach for full nonparametric Bayesian inference under Dirichlet process mixture models[END_REF]. In this paper we propose and investigate a random truncation by setting n such that R n is smaller than a predetermined value ∈ (0, 1) with probability one. Specifically, we define τ ( ) = min{n ≥ 1 : R n < } (4)

as the stopping time of the multiplicative process (R n ) n≥1 and, following Ghosal and van der Vaart (2017, Section 4.3.3), we call -Pitman-Yor ( -PY) process the Pitman-Yor process truncated at n = τ ( ), namely

P (dx) = τ ( ) i=1 p i δ ξi (dx) + R τ ( ) δ ξ0 (dx), (5)
where ξ 0 has distribution P 0 , independent of the sequences (p i ) i≥1 and (ξ i ) i≥1 . By construction, P is the finite stick-breaking approximation to P with the smallest number of support points given a predetermined approximation level. In fact τ ( ) controls the error of approximation according to the total variation bound

d T V (P , P ) = sup A⊂X |P (A) -P (A)| ≤ (6)
almost surely (a.s.). As such, it also guarantees the almost sure convergence of measurable functionals of P by the corresponding functionals of P as → 0, cf. Ghosal and van der Vaart (2017, Proposition 4.20). A typical application is in Bayesian nonparametric inference on mixture models where the Pitman-Yor process is used as prior distribution on the mixing measure. The approximation P can be applied to the posterior distribution given the latent variables, cf. Section 2.2 for details. In the Dirichlet process case, P has been studied by [START_REF] Muliere | Approximating distributions of random functionals of Ferguson-Dirichlet priors[END_REF]. In this setting τ ( ) -1 is Poisson distributed with parameter θ log 1/ , which makes an exact sampling of the -approximation (5) feasible. This has been implemented in the highly popular R software DPpackage, see [START_REF] Jara | Applied Bayesian non-and semi-parametric inference using DPpackage[END_REF][START_REF] Jara | DPpackage: Bayesian semi-and nonparametric modeling in R[END_REF], to draw posterior inference on the random effect distribution of linear and generalized linear mixed effect model. Finally, in Al Labadi and Zarepour (2014) a different type of finite dimensional truncation of the Pitman-Yor process based on decreasing frequencies has been proposed, see Section 5 for a discussion.

The main theoretical contribution of this paper is the derivation of the asymptotic distribution of τ ( ) as → 0 for α > 0. As (4) suggests, the asymptotic distribution of τ ( ) is related to that of R n in (3) as n → ∞. According to Pitman (2006, Lemma 3.11), the latter involves a polynomially tilted stable random variable T α,θ , see Section 2 for a formal definition. The main idea is to work with T n = -log R n so to deal with sums of the independent random variables Y i = -log(1 -V i ). The distribution of τ ( ) can be then studied in terms of the allied renewal counting process N (t) = max{n : T n ≤ t}, according to the relation τ ( ) = N (log 1/ ) + 1. The problem boils down to the derivation of an appropriate a.s. convergence of N (t) as t → ∞, which, in turn, is obtained from the asymptotic distribution of T n by showing that N (t) → ∞ a.s. as t → ∞ together with a (non standard) application of the law of large numbers for randomly indexed sequences. This strategy proves successful in establishing the almost sure convergence of τ ( ) -1 to ( T α,θ /α) -α/(1-α) as → 0. The form of the asymptotic distribution reveals how large the truncation point τ ( ) is as gets small in terms of the model parameters α and θ. In particular, it highlights the power law behavior of τ ( ) as → 0, namely the growth at the polynomial rate 1/ α/(1-α) compared to the slower logarithmic rate θ log 1/ in the Dirichlet process case. This is further illustrated by a simulation study in which we generate from the asymptotic distribution of τ ( ) by using Zolotarev's integral representation of the positive stable distribution as in [START_REF] Devroye | Random variate generation for exponentially and polynomially tilted stable distributions[END_REF]. As far as the simulation of the -PY process is concerned, exact sampling is feasible by implementing the stopping rule in (4), that is by simulating the stick breaking frequencies p j until the error R n crosses the approximation level . As this can be computationally expensive when is small, as an alternative we propose to use the asymptotic distribution of τ ( ) by simulating the truncation point first, then run the stick breaking procedure up to that point. It results in an approximate sampler of the -PY process that we compare with the exact sampler in a simulation study involving moments and mean functionals.

The rest of the paper is organized as follows. In Section 2, we derive the asymptotic distribution of τ ( ) and explain how to use it to simulate from the -PY process. Section 3 reports a simulation study on the distribution of τ ( ) and on functionals of the -PY process. In Section 4, to help the understanding and gain additional insight on the asymptotic distribution, we highlight the connections of τ ( ) with Pitman's theory on random partition structures. We conclude with a discussion of open problems in Section 5. The details of Devroye's algorithm for generating from a polynomially tilted positive stable random variable are given in Appendix A.

Theory and algorithms

Asymptotic distribution of τ ( )

In this section we derive the asymptotic distribution of the stopping time τ ( ) and show how to simulate from it. We start by introducing the renewal process interpretation which is crucial for the asymptotic results. As explained in the previous section, in order to study the distribution of τ ( ) it is convenient to work with the log transformation of the truncation error R n in (3), that is

T n = n i=1 Y i , Y i = -log(1 -V i ), (7)
with V j ind ∼ beta(1 -α, θ + jα) as in (2). Being a sum of independent and nonnegative random variables, (T n ) n≥1 takes the interpretation of a (generalized) renewal process with independent waiting times Y i . For t ≥ 0 define

N (t) = max{n : T n ≤ t}, (8)
to be the renewal counting process associated to (T n ) n≥1 , which is related to τ ( ) via τ ( ) = N (log 1/ ) + 1. Classical renewal theory pertains to iid waiting times while here there is no identity in distribution unless α = 0, i.e. the Dirichlet process case. In the latter setting, one gets Y i i.i.d.

∼ Exp(θ) so that T n has gamma distribution with scale parameter n. We immediately get from the relation {T n ≤ t} = {N (t) ≥ n} that N (t) ∼ Pois(θt) and, in turn, that τ ( ) -1 has Pois(θ log(1/ )) distribution. As far as asymptotics is concerned, T n satisfies the central limit theorem (CLT) with (T n -n/θ)/( √ n/θ) → d Z where Z ∼ N(0, 1). The asymptotic distribution of N (t) can be obtained via Ascombe theorem, cf. Gut (2013, Theorem 7.4.1), to get (N (t) -θt)/( √ θt) → d Z, as t → ∞, in accordance with the standard normal approximation of the Poisson distribution with large rate parameter.

In the general Pitman-Yor case α > 0, the waiting times Y i are no more identically distributed. More importantly, generalizations of the CLT such as the Lindeberg-Feller theorem do not apply for T n , hence we cannot resort to Anscombe's theorem to derive the asymptotic distribution of N (t) and, in turn, of τ ( ). Nevertheless, the limit exists but is not normal as stated in Theorem 1 below. To this aim, let T α be a positive stable random variable with exponent α, that is E(e -sTα ) = e -s α , and denote its density by f α (t). A polynomially tilted version of T α is defined as the random variable T α,θ with density proportional to t -θ f α (t), that is

f α,θ (t) = Γ(θ + 1) Γ(θ/α + 1) t -θ f α (t), t > 0. (9)
The random variable T α,θ is of paramount importance in the theory of random partition structures associated to the frequency distribution of the Pitman-Yor process, see Section 4 for details. In particular, the convergence of R n can be expressed in terms of T α,θ . In Theorem 1 the a.s. limit of log N (t) as t → ∞ is obtained from that of T n = -log R n as n → ∞ by showing that N (t) → ∞ a.s. as t → ∞ and by an application of the law of large numbers for randomly indexed sequences.

Theorem 1. Let N (t) be defined in (7)-( 8) and T α,θ be the random variable with density in (9). Then t -(1/α -1) log N (t) + log α → a.s. log T α,θ as t → ∞.

Proof. By definition (8), the renewal process N (t) is related to the sequence of renewal epochs

T n through {T n ≤ t} = {N (t) ≥ n}. ( 10
) Since N (T n ) = n, we have T N (t) = T n when t = T n , thus 0 = t-T N (t) for t = T n . Moreover, since N (t) is increasing, when T n < t < T n+1 , N (T n ) < N (t) < N (t) + 1, hence T N (t) < t < T N (t)+1 , i.e. 0 < t -T N (t) < T N (t)+1 -T N (t) = Y N (t)+1
. Together the two relations above yield

0 ≤ t -T N (t) < Y N (t)+1 . (11)
From Lemma 3.11 of [START_REF] Pitman | Combinatorial stochastic processes[END_REF] and an application of the continuous mapping theorem (see Theorem 10.1 in [START_REF] Gut | Probability : a graduate course[END_REF] the asymptotic distribution of T n is obtained as

T n -(1/α -1) log n + log α → a.s. log T α,θ as n → ∞.
Now we would like to take the limit with respect to n = N (t) as t → ∞, that is apply the law of large numbers for randomly indexed sequence (see Theorem 6.8.1 in [START_REF] Gut | Probability : a graduate course[END_REF]. To this aim, we first need to prove that N (t) → a.s. ∞ as t → ∞. Since N (t) is non decreasing, by an application of Theorem 5.3.5 in [START_REF] Gut | Probability : a graduate course[END_REF], it is sufficient to prove that N (t) → ∞ in probability as t → ∞, that is P(N (t) ≥ n) → 1 as t → ∞ for any n ∈ N. But this is an immediate consequence of the inversion formula (10). We have then established that

T N (t) -(1/α -1) log N (t) + log α → a.s. log T α,θ as t → ∞.
To conclude the proof, we need to replace T N (t) with t in the limit above. Note that, from (11

), |t -T N (t) | ≤ Y N (t)
+1 so it is sufficient to show that the upper bound goes to zero a.s.. Actually, by a second application of Theorem 6.8.1 in [START_REF] Gut | Probability : a graduate course[END_REF] it is sufficient to show that Y n → a.s. 0 as n → ∞. This last result is established as follows. Recall that

Y j = -log(1 -V j ) for V j ind ∼ beta(1 -α, θ + jα). For > 0, P(1 -V n < e -) = e - 0 Γ(θ + nα + 1 -α) Γ(θ + nα)Γ(1 -α) v θ+nα-1 (1 -v) -α dx ≤ (1 -e -) -α Γ(1 -α) Γ(θ + nα + 1 -α) Γ(θ + nα) e -(θ+nα) θ + nα = (1 -e -) -α Γ(1 -α) (θ + nα) -α e -(θ+nα) 1 + O 1 θ + nα (12)
where in equality ( 12) we have used Euler's formula

Γ(z + α)/Γ(z + β) = z α-β 1 + (α -β)(α + β -1) 2z + O(z -2 )
for z → ∞, see [START_REF] Tricomi | The asymptotic expansion of a ratio of gamma functions[END_REF]. Since P(Y n > ) = P(1 -V n < e -), ( 12) implies that P(Y n > ) is exponentially decreasing in n and, in turn, that n≥1 P(Y n > ) < ∞. An application of Borel-Cantelli Lemma yields Y n → a.s. 0 and the proof is complete.

The asymptotic distribution of τ ( ) is readily derived from Theorem 1 via the formula τ ( ) = N (log 1/ ) + 1 and an application of the continuous mapping theorem. The proof is omitted.

Theorem 2. Let τ ( ) be defined in (4) and T α,θ be the random variable with density in (9).

Then τ ( ) -1 ∼ a.s. ( T α,θ /α) -α/(1-α) as → 0.
In order to sample from the asymptotic distribution of τ ( ), the key ingredient is random generation from the polynomially tilted stable random variable T α,θ . Following [START_REF] Devroye | Random variate generation for exponentially and polynomially tilted stable distributions[END_REF], we resort to Zolotarev's integral representation, so let A(u) be the Zolotarev function

A(x) = sin(αx) α sin((1 -α)x) 1-α sin(x) 1 1-α , x ∈ [0, π]
and Z α,b , α ∈ (0, 1) and b > -1 be a Zolotarev random variable with density given by

f (x) = Γ(1 + bα)Γ(1 + b(1 -α)) πΓ(1 + b)A(x) b(1-α) , x ∈ [0, π].
According to Theorem 1 of [START_REF] Devroye | Random variate generation for exponentially and polynomially tilted stable distributions[END_REF], for G a a gamma distributed random variable with shape a > 0 and unit rate,

T α,θ d = A(Z α,θ/α ) G 1+θ(1-α)/α 1-α α
so that random variate generation simply requires one gamma random variable and one Zolotarev random variable. For the latter, rejection sampler can be used as detailed in [START_REF] Devroye | Random variate generation for exponentially and polynomially tilted stable distributions[END_REF]. See Algorithm 3 in Appendix A.

Simulation of the -PY process

Given α, θ, and a probability measure P 0 on X , an -PY process can be generated by implementing the stopping rule in the definition of τ ( ), cf. ( 4). The algorithm consists in a while loop as follows:

Algorithm 1 (Exact sampler of -PY)

1. set i = 1, R = 1 2. while R ≥ : generate V from beta(1 -α, θ + iα). set p i = V R, R = R(1 -V ), i = i + 1 3. set τ = i, R τ = R 4. generate τ + 1 random variates ξ 0 , ξ 1 , . . . , ξ τ from P 0 5. set P (dx) = τ i=1 p i δ ξi (dx) + R τ δ ξ0 (dx)
When is small, the while loop happens to be computationally expensive since conditional evaluations at each iteration slow down computation, and memory allocation for the frequency and location vectors cannot be decided beforehand. In order to avoid these pitfalls and make the algorithm faster, one should generate the stopping time τ ( ) first, and the frequencies up to that point later. We propose to exploit the asymptotic distribution of τ ( ) in Theorem 2 as follows:

Algorithm 2 (Approximate sampler of -PY)

1: generate T d = T α,θ 2: set τ ← 1 + ( T /α) -α/(1-α) 3. for i = 1, . . . , τ : generate V i from beta(1 -α, θ + iα). set p i = V i i-1 j=1 (1 -V j ) 4. set R τ = 1 - τ i=1 p i = τ i=1 (1 -V j ) 5: generate τ + 1 random variates ξ 0 , ξ 1 , . . . , ξ τ from P 0 6: set P (dx) = τ i=1 p i δ ξi (dx) + Rδ ξ0 (dx)
Algorithm 2 is an approximate sampler of the -PY process (while Algorithm 1 is an exact one) since it introduces two sources of approximations. First, through the use of the asymptotic distribution of τ ( ). Second, through Step 3 since the V i 's are not generated according to the conditional distribution given τ ( ), rather unconditionally. Finding the conditional distribution of V i , or an asymptotic approximation thereof, is not an easy task and is object of current research. In terms of the renewal process interpretation in ( 7)-( 8), the problem is to generate the waiting times Y i = -log(1 -V i ), i = 1, . . . , n, from the conditional distribution of the renewal epochs (T 1 , . . . , T n ) given N (t) = n for t = -log 1/ .

A typical use of samples from the Pitman-Yor process we have in mind is in infinite mixture models. In fact, the discrete nature of the Pitman-Yor process makes it a suitable prior on the mixing distribution. Algorithm 1 or Algorithm 2 can be then applied to approximate a functional of the posterior distribution of the mixing distribution. In such models, the process components can be seen as latent features exhibited by the data. Let P denote such a process, n denote the sample size and X 1:n = (X 1 , . . . , X n ) be an exchangeable sequence from P , that is

X 1:n |P i.i.d. ∼ P . Variables X 1:n are latent variables in a model conditionally on which observed data Y 1:n come from: Y j |X j ind ∼ f ( • |X j )
where f denotes a kernel density. Actually, independence is not necessary here and applications also encompass dependent models such as Markov chain transition density estimation. In order to deal with the infinite dimensionality of the process, a strategy is to marginalize it and to draw posterior inference with a marginal sampler. Since draws from a marginal sampler allows to make inference only on posterior expectations of the process, for more general functionals of P , in the form of ψ(P ), one typically needs to resort to an additional sampling step. Exploiting the composition rule

L(ψ(P )|Y 1:n ) = L(ψ(P )|X 1:n ) × L(X 1:n |Y 1:n ) this
additional step boils down to sampling P conditional on latent variables X 1:n . At this stage, recalling the conditional conjugacy of the Pitman-Yor process is useful. Among X 1:n , there are a number k ≤ n of unique values that we denote by X * 1:k . Let n * 1:k denote their frequencies. Then the following identity in distribution holds

P |X 1:n = k j=1 q j δ X * j + q k+1 P * ,
where, independently, (q 1 , . . . , q k , q k+1 ) ∼ Dirichlet(n * 1 -α, . . . , n * k -α, θ + αk) and P * is a Pitman-Yor process of parameter (α, θ +αk), see Corollary 20 of [START_REF] Pitman | Some developments of the Blackwell-MacQueen urn scheme[END_REF]. Thus sampling from L(P |X 1:n ), hence from L(ψ(P )|X 1:n ), requires sampling the infinite dimensional P * . Cf. Ishwaran and James (2001, Section 4.4). For the sake of comparison, the conjugacy of the Dirichlet process similarly leads to the need of sampling an infinite dimensional process, where P |X 1:n takes the form of a Dirichlet process. As already noticed, the truncation of the Dirichlet process is very well understood, both theoretically and practically. The popular R package DPpackage [START_REF] Jara | Applied Bayesian non-and semi-parametric inference using DPpackage[END_REF][START_REF] Jara | DPpackage: Bayesian semi-and nonparametric modeling in R[END_REF] makes use of the posterior truncation point τ * ( ), as defined in ( 5), but here with respect to the posterior distribution of the process. Thus, it satisfies τ * ( ) -1 ∼ Pois((θ + n) log(1/ )), where θ + n is the precision of the posterior Dirichlet process. Adopting here similar lines for the Pitman-Yor process, we replace P * by the truncated process

P * P * (dx) = τ * ( ) i=1 p * i δ ξi (dx) + R τ * ( ) δ ξ0 (dx),
cf. equation (5). Here (p * i ) i≥1 are defined according to (2) with θ + αk in place of θ, i.e. V j ind ∼ beta(1 -α, θ + α(k + j)). Hence, according to Theorem 2 we have τ * ( ) -1 ∼ a.s. ( T α,θ+αk /α) -α/(1-α) , as → 0 hence Algorithm 2 can be applied here.

3 Simulation study 3.1 Stopping time τ ( ) According to Theorem 2, the asymptotic distribution of τ ( ) changes with , α and θ. For illustration, we simulate τ ( ) from Steps 1.-2. in Algorithm 2 using Devroye's sampler, cf. Algorithm 3 in Appendix A. In Figure 1 we compare density plots obtained with 10 4 iterations with respect to different combinations of , α and θ. The plot in the left panel shows how smaller values of result in larger values of τ ( ). In fact, as → 0, τ ( ) increases proportional to 1/ α/(1-α) . Note also that ( T α,θ /α) -α/(1-α) is nonnegative for T α, < α/ , which happens with high probability when is small. As for α, the plot in the central panel shows how τ ( ) increases as α gets large. In fact, it is easy to see that ( T α,θ /α) -α/(1-α) is increasing in α when T α, < e 1-α α/ , which also happens with high probability when is small, so the larger α, the more stick-breaking frequencies are needed in order to account for a prescribed approximation error . Finally, the plot in the right panel shows that the larger θ, the larger τ ( ). In fact, by definition, the polynomial tilting makes T α,θ stochastically decreasing in θ. In order to illustrate the rate of convergence in Theorem 2, we compare next the exact distribution of τ ( ) with the asymptotic one. To do so, we repeat the following experiment several times: we simulate τ ( ) from Steps 1.-3. in Algorithm 1, then we compare the empirical distribution of ( /α) α (τ ( ) -1) 1-α with T -α α,θ , the latter corresponding to the αdiversity of the PY process, see Section 4 for a formal definition. In Table 1 we report the Kolmogorov distance together with expected value, median, first and third quartiles for both the exact and the asymptotic distribution obtained with 10 4 iterations. This is repeated for α = 0.5, θ = {0, 1, 10} and = {0.10, 0.05, 0.01}. As expected, as we decrease , the Kolmogorov distance gets smaller to somehow different rates according to the parameter choice. The derivation of convergence rates is left for future research.

Functionals of the -PY process

In the case that P is defined on X ⊆ R, the total variation bound (6) implies that |F (x) -F (x)| < almost surely for any x ∈ R, where F and F are the cumulative distribution functions of P and P . Also, measurable functionals ψ(P ) such as the mean µ = xP (dx) can be approximated in distribution by the corresponding functionals ψ(P ). For illustration, we set X = [0, 1] and P 0 the uniform distribution on [0, 1]. For given α and θ, we then compare the distribution under P with that under the -PY process P for F (1/2), F (1/3) and µ = xP (dx).

As for the distribution of the finite dimensional distributions F (1/2) and F (1/3) under the full process P , we set α = 0.5 so to exploit results in [START_REF] James | On the posterior distribution of classes of random means[END_REF]. According to their Proposition 4.7, the finite dimensional distributions of P when α = 0.5 are given by f (w 1 , . . . , w n-1 ) = (

n i=1 p i )Γ(θ + n/2) π (n-1)/2 Γ(θ + 1/2) w -3/2 1 • • • w -3/2 n-1 (1 - n-1 i=1 w i ) -3/2 A n (w 1 , . . . , w n-1 ) θ+n/2 for any partition A 1 , . . . , A n of X with p i = P 0 (A i ) and A n (w 1 , . . . , w n-1 ) = n-1 i=1 p 2 i w -1 i + p 2 n (1- n-1 i=1 w i ) -1 . Direct calculation shows that F (1/2) has beta distribution with parameters (θ + 1/2, θ + 1/2) while F (1/3) has density f (w) = 2 √ π 9 θ Γ(θ + 1) Γ(θ + 1/2) (w(1 -w)) θ-1/2 (1 + 3w) θ+1 .
As for the mean functional µ = xP (dx), the distribution under the full process P is approximated by simulations by setting a deterministic truncation point sufficiently large. As for the distribution under P , we use both Algorithm 1 and Algorithm 2.

In Figure 2 we compare the density plots of F (1/2) for = {0, 1.0.05, 0.001} and θ = {0, 10} under P with the beta density under P so to illustrate that the two distributions get close as gets small. As for F (1/3) and µ = xP (dx), in Tables 2 and3 we report the Kolmogorov The density under the Pitman-Yor process is the black dotted curve. The parameter α is fixed equal to 0.5, θ is equal to 0 on the first row and 10 on the second row, while is respectively equal to {0.10, 0.05, 0.01} in the left, center and right columns.

distance between P and P for the two sampling algorithms, together with expected value, median, first and third quartiles. For each case and each parameter configuration, we have sampled 10 4 trajectories from the -PY process and 10 4 trajectories from the Pitman-Yor process in the case of µ = xP (dx). As expected, the Kolmogorov distances are generally larger, still close, when using Algorithm 2 versus Algorithm 1 due to the approximate nature of the former.

Computation time

In this section, we provide a concrete justification of the computational advantage of using Algorithm 2 versus Algorithm 1. We simulate 10 4 -PY iterations by using Algorithm 1 and Algorithm 2 for different combinations of the α and θ parameters and of the error threshold. In Table 4 (resp. Al1 Al2 Al1 Al2 PY Al1 Al2 PY Al1 Al2 PY Al1 Al2 PY 0 0.10 16.29 16.48 0.33 0.33 0.33 0.04 0.01 0.04 0.20 0.16 0.20 0.60 0.64 0.59 0 0.05 11.53 12.52 0.33 0.33 0.33 0.05 0.01 0.04 0.20 0.17 0.20 0.58 0.63 0.59 0 0.01 5.49 5.60 0.34 0.33 0.33 0.04 0.03 0.04 0.21 0.19 0.20 0.59 0.61 0.59 1 0.10 3.08 5.65 0.33 0.33 0.33 0.14 0.12 0.14 0.29 0.28 0.28 0.49 0.50 0.49 1 0.05 1.34 3.11 0.33 0.33 0.33 0.14 0.13 0.14 0.28 0.28 0.28 0.48 0.50 0.49 1 0.01 0.56 0.89 0.33 0.34 0.33 0.14 0.14 0.14 0.28 0.29 0.28 0.49 0.49 0.49 10 0.10 3.10 3.81 0.33 0.33 0.33 0.25 0.25 0.26 0.32 0.32 0.32 0.40 0.41 0.40 10 0.05 1.41 1.38 0.33 0.33 0.33 0.26 0.26 0.26 0.32 0.32 0.32 0.40 0.40 0.40 10 0.01 0.75 0.65 0.33 0.33 0.33 0.26 0.26 0.26 0.33 0.32 0.32 0.40 0.40 0.40 Table 3: Simulation study on µ = xP (dx)

d K Mean 25% Median 75% θ Al1 Al2 Al1 Al2 PY Al1 Al2 PY Al1 Al2 PY Al1
Al2 PY 0 0.10 1.60 3.57 0.50 0.50 0.50 0.36 0.34 0.36 0.50 0.50 0.50 0.64 0.67 0.65 0 0.05 0.94 2.72 0.50 0.50 0.50 0.35 0.34 0.36 0.50 0.50 0.50 0.64 0.66 0.65 0 0.01 1.18 2.10 0.50 0.50 0.50 0.36 0.35 0.36 0.50 0.50 0.50 0.64 0.65 0.65 1 0.10 1.61 3.18 0.50 0.50 0.50 0.40 0.39 0.40 0.50 0.50 0.50 0.60 0.61 0.60 1 0.05 1.28 2.32 0.50 0.50 0.50 0.40 0.40 0.40 0.50 0.50 0.50 0.59 0.60 0.60 1 0.01 1.12 0.57 0.50 0.50 0.50 0.40 0.41 0.40 0.50 0.50 0.50 0.60 0.60 0.60 10 0.10 2.81 4.18 0.50 0.50 0.50 0.45 0.46 0.46 0.50 0.50 0.50 0.55 0.55 0.54 10 0.05 1.78 1.28 0.50 0.50 0.50 0.46 0.46 0.46 0.50 0.50 0.50 0.54 0.54 0.54 10 0.01 2.01 1.09 0.50 0.50 0.50 0.46 0.46 0.46 0.50 0.50 0.50 0.54 0.54 0.54 Summary statistics for F (1/3) (Table 2) and µ = xP (dx) (Table 3) using Algorithm 1 (Al1) and Algorithm 2 (Al2) to sample from the -PY process. The Kolmogorov distance (d K ) is between the cumulative distribution functions with respect to the Pitman-Yor (PY) process (multiplied by a factor of 100). The parameter values are α = 0.5, θ ∈ {0, 1, 10} and ∈ {0.10, 0.05, 0.01}.

realization of the -PY process including frequencies and locations, while by support point, we mean that we divide the total time by the number of support points τ ( ) + 1. In order to account for the computational task required per iteration, the expected stopping time E[τ ( )] is also reported. Both tables illustrate that our proposed approach is faster than Algorithm 1 when the -PY is composed of about 20 support points or more. The more support points, the faster Algorithm 2 is compared to Algorithm 1. This disadvantage of the former for small numbers of support points comes from the fixed cost of initially generating a random variable with the same distribution as T α,θ . Conversely, as the number of support points increases, this fixed cost is largely counterbalanced by the fast vector-sampling of a prescribed size, which is in contrast with Algorithm 1 while loop whose cost increases with the number of support points. This can be seen in Table 5 where the actual sampling time per support point is essentially increasing for Algorithm 1 and decreasing for Algorithm 2. With the parameter configurations tested, Algorithm 2 can be up to 90 times faster Algorithm 1 for α = 0.6, θ = 10 and = 0.01. The random variable T α,θ in Theorem 1 plays a key role in the Pitman-Yor process, in particular for its link with the α-diversity of the process. The α-diversity is defined as the almost sure limit of n -α K n where K n denotes the (random) number of unique values in the first n terms of an exchangeable sequence from P in (1). According to Theorem 3.8 in [START_REF] Pitman | Combinatorial stochastic processes[END_REF], n -α K n ∼ a.s. (T α,θ ) -α , in particular, for θ = 0, T -α α has a Mittag-Leffler distribution with p-th moment Γ(p + 1)/Γ(pα + 1), p > -1. According to Pitman (2006, Lemma 3.11, eqn (3.36)), the asymptotic distribution of the truncation error R n can be derived from that of K n to get R n ∼ a.s. α(T α,θ ) -1 n 1-1/α as n → ∞. The proof relies on Kingman's representation of random partitions [START_REF] Kingman | The representation of partition structures[END_REF] together with techniques set forth by [START_REF] Gnedin | Notes on the occupancy problem with infinitely many boxes: general asymptotics and power laws[END_REF]. In the proof of Theorem 1 the asymptotic distribution of T n = -log R n is a direct consequence of the above by an application of the continuous mapping theorem.

When θ = 0 it is possible to give an interpretation of the asymptotic distribution of R n in terms of the jumps of a stable subordinator. In this case the weights of P can be represented as the renormalized jumps of a stable subordinator, with T α denoting the total mass. Denote the (unormalized) jumps as (J i ) i≥1 in decreasing order and as ( Ji ) i≥1 when in size-biased order,

T α = i≥1 J i = i≥1 Ji , and T α R n = i>n Ji .
By the asymptotic distribution of R n , n 1/α-1 i>n Ji → a.s. α as n → ∞. That is, once properly scaled, the small jumps of the stable subordinator (in size-biased random order), interpreted as the " dust ", converge to the " proportion " α. This is reminiscent to the number of singletons which is asymptotically (n → ∞) a α proportion of the number of groups in a sample of size n, see Lemma 3.11, eqn (3.39), of [START_REF] Pitman | Combinatorial stochastic processes[END_REF].

Regenerative random compositions and Anscombe's theorem

We review next the connections of the counting renewal process N (t) defined in ( 7)-( 8) and the theory of regenerative random compositions. The reader is referred to the survey of Gnedin (2010) for a review. Recall that, when α = 0 (Dirichlet process case), V i i.i.d.

∼ beta(1, θ) in the stick-breaking representation (2), and in turns

Y i = -log(1 -V i ) i.i.d.
∼ Exp(θ) and T n = -log R n ∼ Gamma(n, θ). By direct calculus, N (t) ∼ Pois(θt) so that τ ( ) -1 = N (log 1/ ) has Pois(θ log 1/ ) distribution. More generally, the stick-breaking frequencies (p i ) i≥1 correspond to the gaps in [0, 1] identified by the multiplicative regenerative set R ⊂ (0, 1) consisting of the random partial sums 1 -R k = i≤k p i . The complement open set R c = (0, 1)/R can be represented as a disjoint union of countably many open intervals or gaps,

R c = ∞ k=0 (1 - R k , 1 -R k+1 ), R 0 = 1.
A random composition of the integer n into an ordered sequence κ n = (n 1 , n 2 , . . . , n k ) of positive integers with j n j = n can be generated as follows: independently of R, sample U 1 , U 2 , . . . from the uniform distribution on [0, 1] and group them in clusters by the rule: U i , U j belong to the same cluster if they hit the same interval. The random composition of κ n corresponds then to the record of positive counts in the left-to-right order of the intervals. The composition structure (κ n ) is called regenerative since for all n > m ≥ 1, conditionally given the first part of κ n is m, if the part is deleted then the remaining composition of n -m is distributed like κ n-m . The regenerative set R corresponds to the closed range of the multiplicative subordinator {1 -exp(-S t ), t ≥ 0}, where S t is the compound Poisson process with Lévy intensity ν(dy) = θe -θy dy. Since the range of S t is a homogeneous Poisson point process on R + with rate θ, R is an inhomogeneous Poisson point process N (dx) on [0, 1] with Lévy intensity ν(dx) = θ/(1 -x)dx so that, for t = log 1/ ,

N (log 1/ ) = N [0, 1 -] ∼ Pois(λ), λ = 1- 0 θ 1 -x dx = θ log 1/
as expected. Suppose now that (V i ) i≥1 are independent copies of some random variable V on [0, 1], not necessarily beta(1, θ) distributed. The corresponding random composition structure has been studied in [START_REF] Gnedin | The Bernoulli sieve[END_REF]; [START_REF] Gnedin | The Bernoulli sieve revisited[END_REF] as the outcome of a Bernoulli sieve procedure. We recall here the relevant asymptotic analysis. Let µ = E(-log(1 -V )) and σ 2 = Var(-log(1 -V )), equal respectively to 1/θ and 1/θ 2 in the DP case, respectively. If those moments are finite, by the CLT,

T n -nµ √ nσ → d Z, as n → ∞,
where Z ∼ N(0, 1), and, by means of Anscobe's Theorem, one obtains that

N (t) -t/µ σ 2 t/µ 3 → d Z, as t → ∞.
It turns out that the normal limit of N (log n) corresponds to the normal limit of K n ,

K n -log n/µ σ 2 log n/µ 3 → d Z, as n → ∞
provided that E(-log V ) < ∞. To see why, consider iid random variables X 1 , X 2 , . . . with values in N such that

{X i = k} = {U i ∈ (1 -R k-1 , 1 -R k )}. Hence P(X 1 = k|R) = p k .
We then have that K n = #{k : X i = k for at least one i among 1, . . . , n}. Define M n = max{X 1 , . . . , X n }.

For U 1,n ≤ U 2,n ≤ . . . ≤ U n,n denoting the order statistics corresponding to the uniform variates U 1 , . . . , U n , we have

M n = min{j : 1 -R j ≥ U n,n } = min{j : T j ≥ E n,n } upon transformation x → -log(1 -x)
, where E n,n is the maximum of an iid sample of size n from the standard exponential distribution. Since N (t) = max{n : [START_REF] Gnedin | The Bernoulli sieve revisited[END_REF] proves the equivalence

T n ≤ t} = min{n : T n ≥ t} -1 we have M n -1 = N (E n,n ).
M n -b n a n → d X ⇐⇒ N (log n) -b n a n → d X
where X is a random variable with a proper and non degenerate distribution with a n > 0, a n → ∞ and b n ∈ R. A key fact exploited in the proof is that, from extreme-value theory, E n,n -log n has an asymptotic distribution of Gumbel type. That M n can be replaced by K n in the equivalence relation above follows from the fact that M n -K n , the number of integers k < M n not appearing in the sample X 1 , . . . , X n , is bounded in probability when E(-log V ) < ∞, see Proposition 5.1 in [START_REF] Gnedin | The Bernoulli sieve revisited[END_REF].

Back to the Pitman-Yor process case, by Theorem 1 we have n -α/(1-α) N (log n) → d (T α,θ /α) -α/(1-α) while n -α K n → a.s. (T α,θ ) -α . So we see that N (log n) and K n do not have the same asymptotic behavior as in the α = 0 case. By using the fact that

P(X 1 > n|(p i )) = R n , R n ∼ a.s. αn -(1-α)/α T -1 α,θ ,
and the fact that, conditional on (p i ) i≥1 , M n belongs to the domain of attraction of Fréchet distribution, Pitman and Yakubovich (2017, Theorem 6.1) establishes that

P(M n ≤ xn α/(1-α) ) → E exp -αT -1 α,θ x -(1-α)/α
so we see that N (log n) and M n do not have the same asymptotic behavior as in the α = 0 case, although they share the same growth rate n α/(1-α) . Finally, the non correspondence of the asymptotic distribution of M n and K n suggests that the behavior of M n -K n is radically different with respect to the α = 0 case.

Discussion

In this paper we have studied stochastic approximations of the Pitman-Yor process consisting in the truncation of the sequence of stick-breaking frequencies at a random stopping time τ ( ) that controls the accuracy of the approximation in the total variation distance by . We name this finite dimensional approximation the -Pitman-Yor process. We have derived the asymptotic distribution of τ ( ) as goes to zero and we have advanced its use to devise a sampling scheme that generates the stopping time first, and then the frequencies up to that point. The simulations in Section 3 show that the proposed sampler proves computationally very efficient in the moderate to large stopping time regime (for approximately τ ( ) ≥ 20). The asymptotic distribution illustrates how large the stopping time is as the approximation error gets small in terms of the prior parameters θ and α. In particular, it shows that the distribution of τ ( ) in the Dirichlet process case is not recovered in the limit α → 0 in Theorem 2. In fact, in the Dirichlet process case τ ( ) grows at a logarithmic rate in 1/ while in Pitman-Yor case it grows at the polynomial rate α/(1-α) and the first regime is not recovered by letting α approach 0 in the second regime. We have also drawn important connections with the theory of random partition structures developed by Jim Pitman and coauthors which highlight the relationship of the the stopping time τ ( ) with the number K n of unique values in a sample of size n from the Pitman-Yor process.

We have left as open problem for future research the study of the conditional distribution of the stick-breaking frequencies given the stopping time. In the Dirichlet process case one can exploit the renewal process interpretation to generate exactly from this conditional distribution. In fact, when α = 0, the sequence (-log R n ) n≥1 corresponds to the jump times of a Poisson process and the conditional distribution of the jumps given the number of jumps at time t can be described in terms of the ordered statistics of i.i.d. uniform random variates on (0, t). The case α > 0 does not seem to be easily tractable, as it would be if the counting process associated to τ ( ) were a mixed sample process or, equivalently, a Cox process, cf. Grandell (1997, Section 6.3).

It would be also interesting to compare the accuracy of our finite dimensional approximation of the Pitman-Yor process to the one proposed in Al Labadi and Zarepour (2014). The latter is based on a representation of the frequencies in decreasing order, cf. Pitman and Yor (1997, Proposition 22). Al Labadi and Zarepour (2014) compare the accuracy of their approximation scheme to a stick-breaking truncation at a number n of stick-breaking frequencies that matches the number of frequencies used in their scheme. Not surprisingly, their approximation is superior since it generates weights in decreasing order, specially when α is large. In contrast, Theorem 2 describes precisely how large the truncation threshold n should be as α gets large for a given approximation level , cf. the center panel of Figure 1. It also underlines that the approximation deteriorates for fixed n and increasing α, which is coherent with the findings in Al [START_REF] Al Labadi | On simulations from the two-parameter poisson-dirichlet process and the normalized inverse-gaussian process[END_REF]. A fair comparison with their approach can only be done for a given nominal approximation error, but unfortunately the authors did not provide a precise assessment of it. The number of stick-breaking frequencies needed to match the approximation accuracy of Al Labadi and Zarepour (2014) would be de facto larger due to the non monotonicity. However, since the stopping rule (4) adapts to the size of α, we do not expect the accuracy of our approximation scheme to deteriorate for α large. As for computation time, the techniques used by Al Labadi and Zarepour (2014) in order to obtain decreasing frequencies are computational heavy. Their average computing time for α = 0.5 is about 2.30 seconds/iteration with 10 4 locations. This amounts to 0.23 milli seconds/support point, which is 1000 times slower than the computing time for our Algorithm 2 in the parameter configuration α = 0.5, θ = 10 and = 0.01, equal to 0.23 microsecond/support point. It would be interesting to investigate what are the consequences in terms of computation time per iteration for a given approximation error.

A Random generation of T α,θ

Let Y be a standard Exp(1) random variable. Note that P((Y /T α ) α > x) = ∞ 0 P(Y > x 1/α t)f α (t)dt = ∞ 0 exp[-x 1/α t]f α (t)dt = E e -x 1/α Tα = e -x = P(Y > x) so we have Y = d (Y /T α ) α . For r < α, E(Y -r/α ) < ∞, so we find that E(Y -r/α ) = E(T r α )E(Y -r ) and

E(T r α ) = E(Y -r/α ) E(Y -r ) = Γ(1 -r/α) Γ(1 -r) . (13)
The normalizing constant in f α,θ (t) is ∞ 0 t -θ f α (t)dt = E(T -θ α ), so set r = -θ and note that -θ < α. Let G a be a gamma random variable with shape a > 0 and unit rate. Simple moment comparisons using (13) yield the distributional equality G 1+θ/α d = (G 1+θ /T α,θ ) α , which, however, does not provide a way to generate from T α,θ . For this we resort to [START_REF] Devroye | Random variate generation for exponentially and polynomially tilted stable distributions[END_REF]. This Gaussian upper bound suggests a simple rejection sampler for sampling Zolotarev random variates. Following [START_REF] Devroye | Random variate generation for exponentially and polynomially tilted stable distributions[END_REF], it is most efficient to adapt the sampler to the value of σ. If σ ≥ √ 2π, rejection from a uniform random variate is best. Otherwise, use a normal dominating curve as suggested in the bound above. The details are given below.

Figure 1 :

 1 Figure 1: Density plot for the asymptotic approximation of τ ( ) based on 10 4 values under the following parameter configurations. Left: ∈ {0.10, 0.05, 0.01}, α = 0.4, θ = 1. Center: α ∈ {0.4, 0.5, 0.6}, θ = 1, = 0.1. Right: θ ∈ {0, 1, 10}, α = 0.25, = 0.05.

Figure 2 :

 2 Figure2: Density plots for the random probability F (1/2) using the Algorithm 2 (in red solid curve) and Algorithm 1 (in blue dashed curve) to sample from the -PY process. The density under the Pitman-Yor process is the black dotted curve. The parameter α is fixed equal to 0.5, θ is equal to 0 on the first row and 10 on the second row, while is respectively equal to {0.10, 0.05, 0.01} in the left, center and right columns.

  ) for Algorithm 1 (Al1) and Algorithm 2 (Al2) based on 10 4 iterations, and expected stopping time n = E[τ ( )]. The parameter values are α ∈ {0.4, 0.5, 0.6}, θ ∈ {0, 1, 10} and ∈ {0.10, 0.05, 0.01}. 4 Connections with random partition structures 4.1 α-diversity and asymptotic distribution of R n

  First we recall how to generate a Zolotarev random variable Z α,b for α ∈ (0, 1) and b = θ/α > -1. LetC = Γ(1 + bα)Γ(1 + b(1 -α)) πΓ(1 + b)andB(u) = A(u) -(1-α) = sin(u) sin(αu) α sin((1 -α)u) 1-α .A simple asymptotic argument yields the value B(0) = α -α (1 -α) -(1-α) . Then f (x) = C B(x) b , 0 ≤ x ≤ π. The following bound holds f (x) ≤ CB(0) b e -x 22σ 2 , with σ 2 = 1 bα(1 -α) .

Algorithm 3 (

 3 Sampler of T α,θ )1. set b = θ/α and σ = bα(1 -α) 2. if σ ≥ √ 2π: then repeat: generate U ∼ Unif(0, π) and V ∼ Unif(0, 1). set X ← U , W ← B(X). until V ≤ (W/B(0)) b else repeat: generate N ∼ N(0, 1) and V ∼ Unif(, 1). set X ← σ|N |, W ← B(X). until X ≤ π and V e -N 2 /2 ≤ (W/B(0)) b 3. generate G d = G 1+b(1-α)/α 4. set T ← 1/(W G 1-α ) 1/α5. return T

Table 1 :

 1 Summary statistics for the asymptotic distribution (As) and exact distribution (Ex) of τ ( ) at the scale of the α-diversity based on 10 4 values. The Kolmogorov distance (d

		d K	Mean	25%	Median	75%
	θ		As Ex	As Ex	As Ex	As Ex
	0 0.10	3.42 1.06 1.05 0.45 0.45 0.89 0.89 1.61 1.55
	0 0.05	2.17 1.10 1.08 0.45 0.45 0.95 0.95 1.64 1.58
	0 0.01	1.73 1.14 1.11 0.45 0.45 0.97 0.95 1.64 1.60
	1 0.10	4.79 2.24 2.14 1.55 1.48 2.14 2.10 2.86 2.76
	1 0.05	2.38 2.25 2.20 1.55 1.52 2.17 2.14 2.86 2.79
	1 0.01	1.40 2.26 2.25 1.57 1.54 2.19 2.19 2.87 2.85
	10 0.10 11.93 6.39 6.07 5.69 5.40 6.34 6.06 7.04 6.72
	10 0.05	6.12 6.39 6.24 5.70 5.56 6.34 6.22 7.05 6.88
	10 0.01	1.93 6.40 6.37 5.71 5.70 6.34 6.34 7.05 7.00

K ) is between the empirical cumulative distribution function of the sample from the exact distribution and the asymptotic one (multiplied by a factor of 100). The parameter values are α = 0.5, θ ∈ {0, 1, 10} and ∈ {0.10, 0.05, 0.01}.

Table 5

 5

), we report the average computing time 1 per iteration (resp. per support point) for Algorithm 1 and Algorithm 2. By iteration, we mean a full

Table 2 :

 2 Simulation study on F (1/3)

	d K	Mean	25%	Median	75%
	θ

Table 4 :

 4 Computing time (ms) per iteration

		α = 0.4	α = 0.5	α = 0.6
	θ	Al1 Al2 n Al1 Al2	n Al1 Al2	n
	0 0.10 0.01 0.20 5 0.02 0.04 11 0.11 0.05	38
	0 0.05 0.01 0.04 8 0.04 0.04 21 0.20 0.06 105
	0 0.01 0.04 0.04 20 0.36 0.05 101 15.10 0.26 1163
	1 0.10 0.03 0.19 17 0.07 0.05 31 0.23 0.07	92
	1 0.05 0.06 0.06 26 0.13 0.06 61 0.80 0.12 258
	1 0.01 0.18 0.09 73 0.80 0.12 301 27.75 0.57 2877
	10 0.10 0.22 0.15 121 0.61 0.10 211 2.11 0.18 567
	10 0.05 0.45 0.10 191 1.52 0.15 421 9.22 0.37 1603
	10 0.01 1.93 0.20 558 13.24 0.48 2101 760.68 4.01 17911

Table 5 :

 5 Computing time (µs) per support pointAverage computing time per iteration (in milli second in Table4) and per support point (in microsecond in Table5

		α = 0.4	α = 0.5	α = 0.6
	θ	Al1 Al2 n Al1 Al2	n Al1 Al2	n
	0 0.10 1.92 38.46 5 1.82 3.64 11 2.91 1.32	38
	0 0.05 1.30 5.22 8 1.90 1.90 21 1.91 0.57 105
	0 0.01 1.95 1.95 20 3.56 0.50 101 12.98 0.22 1163
	1 0.10 1.81 11.45 17 2.26 1.61 31 2.50 0.76	92
	1 0.05 2.33 2.33 26 2.13 0.98 61 3.10 0.46 258
	1 0.01 2.45 1.23 73 2.66 0.40 301 9.65 0.20 2877
	10 0.10 1.82 1.24 121 2.89 0.47 211 3.72 0.32 567
	10 0.05 2.35 0.52 191 3.61 0.36 421 5.75 0.23 1603
	10 0.01 3.46 0.36 558 6.30 0.23 2101 42.47 0.22 17911

The experiments were conducted on an Intel Core i5 processor (3.1 GHz) computer.