Jean-Éric Pin
email: jean-eric.pin@liafa.jussieu.fr

Pedro V Silva
email: pvsilva@fc.up.pt

A noncommutative extension of Mahler's theorem on interpolation series 1

In this paper, we prove an extension of Mahler's theorem on interpolation series, a celebrated result of p-adic analysis. Mahler's original result states that a function from N to Z is uniformly continuous for the p-adic metric d p if and only if it can be uniformly approximated by polynomial functions. We prove the same result for functions from a free monoid A * to Z, where d p is replaced by the pro-p metric, the profinite metric on A * defined by p-groups.

The aim of this paper is to give a noncommutative version of Mahler's theorem on interpolation series [START_REF] Mahler | An interpolation series for continuous functions of a padic variable[END_REF], which applies to functions from a free monoid to the set of integers. This result was first announced in [START_REF] Pin | A Mahler's theorem for functions from words to integers[END_REF] and this article is the full version of this conference paper.

The classical Stone-Weierstrass approximation theorem states that a continuous function defined on a closed interval can be uniformly approximated by a polynomial function. In particular, if a real function f is infinitely differentiable in the neighbourhood of 0, it can be approximated, under some convergence conditions, by its Taylor polynomials

k n=0 f (n) (0) n! x n
The p-adic analogue of these results is Mahler's interpolation theorem [START_REF] Mahler | An interpolation series for continuous functions of a padic variable[END_REF].

It is is usually stated for p-adic valued functions on the p-adic integers, but this full version can be easily recovered from the simpler version given in 1 Theorem 1 below. First, for each function f : N → Z, there exists a unique family a k of integers such that, for all n ∈ N,

f (n) = ∞ k=0 a k n k or, in short notation, f = ∞ k=0 a k - k .
The coefficients a k are given by the formula

a k = (∆ k f )(0)
where ∆ is the difference operator, defined by (∆f

)(n) = f (n + 1) -f (n).
In other words, the binomial functions - k form a basis of the Z-module of all functions from N to Z. The series ∞ k=0 a k k is called the Mahler's expansion of f .

Mahler's theorem states that f is uniformly continuous for the p-adic metric if and only if the partial sums of its Mahler's expansion converge uniformly to f . More precisely:

Theorem 1 (Mahler) Let ∞ k=0 a k - k
be the Mahler's expansion of a function f : N → Z. The following conditions are equivalent:

(1) f is uniformly continuous for the p-adic metric, (2) the sequence of partial sums n k=0 a k k converges uniformly to f , (3) lim k→∞ |a k | p = 0.

The most remarkable part of this theorem is the fact that any uniformly continuous function can be approximated by polynomial functions, in contrast to Stone-Weierstrass approximation theorem, which requires much stronger conditions. This nice feature makes Mahler's theorem the dream of students in mathematics: a function is equal to the sum of its Newton series if and only if it is uniformly continuous.

Our goal is to give a noncommutative extension of this result, based on the following observation: since the additive monoid N is isomorphic to the free monoid on a one-letter alphabet, the free monoid A * on a finite alphabet A can be viewed as a noncommutative generalization of the natural numbers. Accordingly, our noncommutative extension concerns functions from A * to Z. Nevertheless, generalizing Mahler's theorem to these functions requires a few preliminary questions to be solved:

(1) Extend the difference operators to word functions, (2) Extend the binomial coefficients to words, (3) Find a metric on A * which generalizes d p , (4) Find a Mahler expansion for functions from A * to Z. Fortunately, (1) can be readily solved and (2) and (3) have already be addressed by Eilenberg [START_REF] Eilenberg | Automata, Languages and Machines[END_REF]: there is a well-defined notion of binomial coefficients u v for two words u and v and the pro-p metric on A * is a natural extension of the p-adic metric. Question (4) is the topic of Section 2, which culminates with Theorem 2.2: just like in the commutative case, every function f : A * → Z has a Mahler's expansion of the form

v∈A * f, v - v
where the coefficients f, v are integers related to the difference operators.

There is also a striking similarity between Mahler's theorem and our main result, which can be stated as follows:

Theorem 2 Let v∈A * f, v - v
be the Mahler's expansion of a function f : A * → Z. The following conditions are equivalent:

(1) f is uniformly continuous for the pro-p metric, (2) the sequence of partial sums

0 |v| n f, v - v converges uniformly to f , (3) lim |v|→∞ | f, v | p = 0.
Although this result can be considered as interesting on its own right, the reader may wonder about our original motivation for proving this theorem. It was actually inspired by automata theoretic questions related to the classification of regularity-preserving functions in language theory (see [START_REF] Pin | A topological approach to transductions[END_REF][START_REF] Pin | On profinite uniform structures defined by varieties of finite monoids[END_REF] for more details). We discuss further this original motivation in Section 7.

In Section 3, we introduce the notion of Mahler polynomials and prove some of their properties. Section 4, devoted to uniform continuity, is a crucial step towards the proof of the main theorem, presented in Section 5. In Section 6, we come back to the commutative case by providing an alternative proof for Amice's Theorem, a generalization of Mahler's Theorem to functions from N k to Z.

Background and notations

In this paper, p denotes a fixed prime number.

The p-adic valuation and the p-adic metric

The p-adic valuation of a nonzero integer n is the natural number

ν p (n) = max k ∈ N | p k divides n
and its p-adic norm is the real number

|n| p = p -νp(n)
By convention, we set ν p (0) = +∞ and |0| p = 0. Note that, for all x, y ∈ Z, |x + y| p max{|x| p , |y| p } |x| p + |y| p (1.1)

|xy| p = |x| p |y| p min{|x| p , |y| p } (1.2)
The p-adic metric d p on N and Z is defined by setting, for each pair (x, y) of integers,

d p (x, y) = |x -y| p (1.3)
These definitions can be extended to N k and to Z k as follows. Given x = (x 1 , . . . , x k ) and y = (y 1 , . . . , y k ) in Z k and a positive integer r, we write x ≡ y (mod r) if, for 1 i k, x i ≡ y i (mod r). The p-adic valuation of x is

ν p (x) = min { n ∈ N | x ≡ 0 (mod p n+1 ) }
and the p-adic norm is |x| p = p -νp(x) . Finally, the p-adic metric on Z k is defined by the formula (1.3), for each pair (x, y) of elements of Z k . It is well-known that d p is actually an ultrametric, that is, satisfies the following stronger form of the triangle inequality, for all u, v, w in Z k :

d p (u, v) max{d p (u, w), d p (w, v)}

Words and languages

Let A be a finite alphabet. Recall that the free monoid A * on A is the set of all words on A, endowed with the concatenation operation. The empty word, denoted 1, is the identity of A * . The length of a word u is denoted |u|. If a is a letter of A, we denote by |u| a the number of occurrences of a in u.

A language is a subset of A * . An important class of languages is the class of recognizable (also called regular ) languages. Recognizable languages admit various characterizations, ranging from automata theory to topology or logic, but we will settle for the following algebraic definition: a language L of A * is recognizable if and only if there exists a monoid morphism ϕ from A * onto some finite monoid M such that L = ϕ -1 (ϕ(L)). It is then said that M recognizes the language L.

Binomial coefficients

Let u and v be two words of

A * . Let u = a 1 • • • a n , with a 1 , . . . , a n ∈ A. Then u is a subword of v if there exist v 0 , . . . , v n ∈ A * such that v = v 0 a 1 v 1 . . . a n v n .
For instance aaba is a subword of caacabca. Following [START_REF] Eilenberg | Automata, Languages and Machines[END_REF][START_REF] Lothaire | Combinatorics on words[END_REF], we define the binomial coefficient of u and v by setting

v u = |{(v 0 , . . . , v n ) | v = v 0 a 1 v 1 . . . a n v n }| .
Observe that if a is a letter, then v a is equal to |v| a . Further, if u = a n and v = a m , then v u = m n and hence these numbers constitute a generalization of the classical binomial coefficients. The next proposition, whose proof can be found in [START_REF] Lothaire | Combinatorics on words[END_REF]Chapter 6], summarizes the basic properties of the generalized binomial coefficients and can serve as an alternative definition.

Lemma 1.1 Let u, v ∈ A * and a, b ∈ A. Then (1) u 1 = 1, (2) u v = 0 if |u| |v| and u = v, ( 3
) ua vb = u vb if a = b u vb + u v if a = b
A third way to define the binomial coefficients is to use the Magnus automorphism of the algebra Z A of polynomials in noncommutative indeterminates in A defined by µ(a) = 1+a for all a ∈ A. One can show that, for all u ∈ A * ,

µ(u) = x∈A * u x x (1.4)
which leads to the formula

u 1 u 2 x = x 1 x 2 =x u 1 x 1 u 2 x 2 (1.5)

The pro-p metric

Let u and v be two words of A * . A p-group G separates u and v if there is a monoid morphism from A * onto G such that ϕ(u) = ϕ(v). One can show that any pair of distinct words can be separated by a p-group [START_REF] Eilenberg | Automata, Languages and Machines[END_REF]. The pro-p metric d p on A * is defined by setting

r p (u, v) = min { n | there is a p-group of order p n+1 separating u and v } d p (u, v) = p -rp(u,v)
with the usual convention min ∅ = +∞ and p -∞ = 0. It turns out that d p is an ultrametric. One can also show that the concatenation product on A * is uniformly continuous for this metric. It follows that the completion of the metric space (A * , d p ) is naturally equipped with a structure of monoid. This monoid is in fact a compact group, called the free pro-p group and denoted F p (A). When A is a one-letter alphabet, this group is of course Z p , the additive group of p-adic integers. One can define in a similar way a pro-p metric on N k or Z k . The attentive reader will notice a potential conflict of notation between the p-adic metric and the pro-p metric on Z k , which were both denoted by d p . The next lemma gets rid of this ambiguity. Proof. First, the pro-p metric on N k is the restriction to N k of the pro-p metric on Z k since every morphism from N k onto a group extends uniquely to Z k . Therefore it suffices to prove that for all u, v ∈ Z k , one has

r p (u, v) = ν p (u -v).
Since Z k is Abelian, a direct product of groups separates u and v if and only if one of its factors separates u and v. Let C r be the cyclic group of order r. In view of the structure theorem of Abelian groups, one has

r p (u, v) = min { n | C p n+1 separates u and v } ν p (u -v)
Further, for every homomorphism f :

Z k → Z, u ≡ v (mod p n+1 ) im- plies f (u) ≡ f (v) (mod p n+1 ). Therefore u ≡ v (mod p n+1 ) is a necessary condition for C p n+1 to separate u and v. Thus r p (u, v) = ν p (u -v) as claimed.
There is a nice connection [START_REF] Pin | Topologie p-adique sur les mots[END_REF] between the pro-p metric and the binomial coefficients, which comes from the characterization of the languages recognized by a p-group given by Eilenberg and Schützenberger [5, Theorem 10.1, p. 239]. Let us call a p-group language a language recognized by a p-group. Note that such a language is recognizable by definition. Proposition 1.3 A language of A * is a p-group language if and only if it is a Boolean combination of the languages

L(x, r, p) = u ∈ A * | u x ≡ r (mod p) ,
for 0 r < p and x ∈ A * .

Let us set now

r p (u, v) = min |x| x ∈ A * and u x ≡ v x (mod p) d p (u, v) = p -r p (u,v) .
It is proved in [START_REF] Pin | Topologie p-adique sur les mots[END_REF]Theorem 4.4] that d p is an ultrametric uniformly equivalent to d p . We shall use this fact in the proof of our main theorem.

The following result, which follows from [15, Theorem 2.3] provides a characterization of uniform continuity in terms of p-group languages. Let A and B be two finite alphabets. Theorem 1.4 A function f : A * → B * is uniformly continuous for d p if and only if, for every p-group language L of B * , the language f -1 (L) is a p-group language of A * .

Since p-group languages are closed under inverse of morphisms [START_REF] Eilenberg | Automata, Languages and Machines[END_REF], one gets the following corollary.

Corollary 1.5 Every monoid morphism from A * to B * is uniformly continuous for d p .

As another consequence of Theorem 1.4, we obtain: Corollary 1.6 For every v ∈ A * , the function from A * into N which maps every word x onto the binomial coefficient x v is uniformly continuous for d p .

Proof. Let f v : A * → N be the function defined by f v (x) = x v . Since the only p-group quotients of N are the cyclic groups of order p n , a language of N recognized by a p-group is a finite union of arithmetic progressions r +p n N (r ∈ {0, . . . , p n -1}). In view of Theorem 1.4, and since the class of p-group languages is closed under union, we only need to show that f -1 v (r + p n N) is a p-group language. Now, [10, Example 4, p. 450] shows that

f -1 v (r + p n N) = x ∈ A * | x v ≡ r (mod p n )
is indeed recognized by a p-group.

Mahler's expansions

The first step to extend Mahler's theorem to functions from words to integers is to define a suitable notion of Mahler's expansion for these functions. Let f : A * → Z be a function. For each letter a, we define the difference operator ∆ a by

(∆ a f )(u) = f (ua) -f (u)
One can now define inductively an operator ∆ w for each word w ∈ A * by setting (∆ 1 f )(u) = f (u), and for each letter a ∈ A,

(∆ aw f )(u) = (∆ a (∆ w f ))(u).
For instance,

∆ aab f (u) = -f (u) + 2f (ua) + f (ub) -f (uaa) -2f (uab) + f (uaab)
The next proposition gives a direct way to compute the difference operators.

Proposition 2.1 For each word word w ∈ A * , the following relation holds:

∆ w f (u) = 0 |x| |w| (-1) |w|+|x| w x f (ux) (2.1)
Proof. We prove the result by induction on the length of w. If w is the empty word, the result is trivial. Let w be a word and let a be a letter. Then

(∆ aw f )(u) = (∆ a (∆ w f ))(u) = (∆ w f )(ua) -(∆ w f )(u)
and by the induction hypothesis applied to w,

(∆ aw f )(u) = 0 |x| |w| (-1) |w|+|x| w x f (uax) -f (ux)
On the other hand, since

aw x = w az + w z if x = az w x
otherwise we get

0 |x| |aw| (-1) |aw|+|x| aw x f (ux) = 0 |az| |w|+1 (-1) |w|+|az|+1 w az + w z f (uaz) + 0 |x| |w|+1 x / ∈aA * (-1) |w|+|x|+1 w x f (ux) = 0 |z| |w| (-1) |w|+|z| w z f (uaz) + 0 |x| |w| (-1) |w|+|x|+1 w x f (ux) = 0 |x| |w| (-1) |w|+|x| w x f (uax) -f (ux)
and thus

(∆ aw f )(u) = 0 |x| |aw| (-1) |aw|+|x| aw x f (ux)
which gives the induction step for the proof of (2.1).

By Corollary 1.6, for a fixed v ∈ A * , the function - v from A * to N which maps a word u to the binomial coefficient u v is uniformly continuous for d p . This family of functions, for v ranging over A * , is locally finite in the sense that, for each u ∈ A * , the binomial coefficient u v is null for all but finitely many words v. In particular, if (m v ) v∈A * is a family of integers, there is a well-defined function from A * to Z defined by the formula

f (u) = v∈A * m v u v
We can now state our first result, which doesn't require any assumption on f .

Theorem 2.2 Let f : A * → Z be an arbitrary function. Then there exists a unique family f, v v∈A * of integers such that, for all u ∈ A * ,

f (u) = v∈A * f, v u v (2.2)
This family is given by

f, v = (∆ v f )(1) = 0 |x| |v| (-1) |v|+|x| v x f (x) (2.3)
Proof. First observe that, according to (2.1)

(∆ v f )(1) = 0 |x| |v| (-1) |v|+|x| v x f (x) (2.4) Thus v∈A * (∆ v f )(1) u v = v∈A * |x| |v| (-1) |v|+|x| v x u v f (x) = x∈A * (-1) |u|+|x|   0 |v| |u| (-1) |v|+|u| u v v x   f (x) = f (u)
in view of the following relation from [7, Corollary 6.3.8]:

0 |v| |u| (-1) |u|+|v| u v v w = 1 if u = w 0 otherwise (2.5)
Uniqueness of the coefficients f, v follows inductively from the formula

f, u = f (u) - 0 |v| < |u| f, v u v ,
a straightforward consequence of (2.2).

The series defined by (2.2) is called the Mahler's expansion of f .

Example 2.1 Let A be the binary alphabet {0, 1}. The empty word will be exceptionally denoted by ε to avoid confusion. Let f : A * → N be the function mapping a binary word onto its value as a binary number. For instance f (010111) = f (10111) = 23. Then one sees easily that

(∆ v f ) = f + 1 if v ∈ 1{0, 1} * f otherwise (∆ v f )(ε) = 1 if v ∈ 1{0, 1} * 0 otherwise It follows that the Mahler expansion of f is f = v∈1{0,1} * - v
For instance, if u = 01001, one gets

f (u) = u 1 + u 10 + u 11 + u 100 + u 101 + u 1001 = 2+2+1+1+2+1 = 9.
The Mahler's expansion of a function f : A * → Z can be used to compute its uniform norm f p , defined by

f p = sup u∈A * |f (u)| p First observe that f p = max u∈A * |f (u)| p since |f (u)| p takes values in the set {p -n | n 0} ∪ {0}. Theorem 2.3 Let v∈A * f, v - v be the Mahler's expansion of a function f : A * → Z. Then f p = max v∈A * | f, v | p . Proof. Let us choose a word u of A * such that f p = |f (u)| p . Then f (u) = |v| |u| f, v u
v and thus we get by (1.1) and (1.2)

f p = |f (u)| p max |v| |u| | f, v | p max v∈A * | f, v | p
Conversely, let w be a word of minimal length such that

| f, w | p = max v∈A * | f, v | p The formula f, w = f (w) - 0 |v|<|w| f, v u v yields,

Mahler polynomials

A function f : A * → Z is a Mahler polynomial if its Mahler's expansion has finite support, that is, if the number of nonzero coefficients f, v is finite.

In this section, we prove in particular that Mahler polynomials are closed under addition and product. We first introduce a convenient combinatorial operation, the infiltration product. We follow the presentation of [START_REF] Lothaire | Combinatorics on words[END_REF]. Let Z A be the ring of formal power series in noncommutative indeterminates in A. Any series s is written as a formal sum s = u∈A * s, u u, a notation not to be confused with our notation for the Mahler's expansion. The infiltration product is the binary operation on Z A , denoted by ↑ and defined inductively as follows: for all u ∈ A * ,

u ↑ 1 = 1 ↑ u = u, (3.1) for all u, v ∈ A * , for all a, b ∈ A ua ↑ bv = (u ↑ vb)a + (ua ↑ v)b + (u ↑ v)a if a = b (u ↑ vb)a + (ua ↑ v)b otherwise (3.2)
for all s, t ∈ Z A ,

s ↑ t = u,v∈A * s, u t, v (u ↑ v) (3.3)
Intuitively, the coefficient u ↑ v, x is the number of pairs of subsequences of x which are respectively equal to u and v and whose union gives the whole sequence x. For instance, ab ↑ ab = ab + 2aab + 2abb + 4aabb + 2abab ab ↑ ba = aba + bab + abab + 2abba + 2baab + baba Also note that u ↑ v, u = u v . We shall need the following relation (see [7, p.131]). For all

v 1 , v 2 ∈ A * , u v 1 u v 2 = x∈A * v 1 ↑ v 2 , x u x (3.4)
Formula 3.4 leads to an explicit computation of the Mahler's expansion of the product of two functions.

Proposition 3.1 Let f and g be two functions from A * to N. Then the coefficients of the Mahler's expansion of f g are given by the formula:

f g, x = v 1 ,v 2 ∈A * f, v 1 g, v 2 v 1 ↑ v 2 , x Proof. Indeed, if f (u) = v∈A * f, v u v and g(u) = v∈A * g, v u v , then f g(u) = v 1 ,v 2 ∈A * f, v 1 g, v 2 u v 1 u v 2
and the result follows by Formula (3.4).

It is now easy to prove the result announced at the beginning of this section.

Uniform continuity

In this section, we give a simple characterization of the uniformly continuous functions from A * to Z. Uniform continuity always refers to the metric d p , but it can be replaced by d p whenever convenient in view of [START_REF] Pin | Topologie p-adique sur les mots[END_REF]Theorem 4.4].

∈ N such that, if |v| s, ν p ( f, v ) r. Setting g = |v|<s f, v - v and h = |v| s f, v - v
we get f = g + h. Further p r divides f, v for |v| s. Since g is a Mahler polynomial, it is uniformly continuous by Corollary 1.6 and so there exists t ∈ N such that d p (u, u ) p -t implies g(u) ≡ g(u ) (mod p r ) and hence f (u) ≡ f (u ) (mod p r ). Thus f is uniformly continuous. This proves the easy direction of the theorem. The key argument for the opposite direction is the following approximation result. Theorem 4.2 Let f : A * → N be a uniformly continuous function. Then there exists a Mahler polynomial P such that, for all u ∈ A * , f (u) ≡ P (u) (mod p).

Proof. We first prove the theorem for some characteristic functions related to the binomial coefficients. The precise role of these functions will appear in the course of the main proof.

Let x ∈ A * and let s be an integer such that 0 s < p. Let χ s,x : A * → N be the function defined by

χ s,x (u) = 1 if u x ≡ s (mod p) 0 otherwise Lemma 4.
3 There is a Mahler polynomial P s,x such that, for all u ∈ A * , χ s,x (u) ≡ P s,x (u) (mod p).

Proof. Let

P s,x (u) = - u x u x -1 • • • u x -(p -1) u x -s Then P s,x is a Mahler polynomial by Proposition 3.2. If u x ≡ s (mod p), then P s,x (u) ≡ 0 (mod p). If u
x ≡ s (mod p), then by Ibn Al-Haytham's theorem (also known as Wilson's theorem),

P s,x (u) ≡ -(p -1)! ≡ 1 (mod p)
It follows that P s,x (u) ≡ χ s,x (u) (mod p) in all cases.

We now prove Theorem 4.2. Since f is uniformly continuous for the metric d p , there exists a positive integer n such that if, for 0 |x| n,

u x ≡ v x (mod p) then f (u) ≡ f (v) (mod p)
It follows that the value of f (u) modulo p depends only on the residues modulo p of the family u x 0 |x| n . Let R be the set of all families r = {r x } 0 |x| n such that 0 r x < p. For 0 i < p, let R i be the set of all families r of R satisfying the following condition:

if, for 0 |x| n, u x ≡ r x (mod p), then f (u) ≡ i (mod p) (4.1)
The sets (R i ) 0 i<p are pairwise disjoint and their union is R. We claim that, for all u ∈ A * , f (u) ≡

0 i<p iP i (u) (mod p) (4.2)
where P i is the Mahler polynomial

P i = r∈R i 0 |x| n P rx,x (4.3)
First consider, for r ∈ R, the characteristic function

χ r (u) = 0 |x| n χ rx,x (u)
By construction, χ r is defined by

χ r (u) = 1 if, for 0 |x| n, u x ≡ r x (mod p) 0 otherwise
and it follows from (4.1) and from the definition of R i that

f (u) ≡ 0 i<p i r∈R i χ r (u) (mod p) (4.4) Now Lemma 4.3 gives immediately χ r (u) ≡ 0 |x| n P rx,x (u) (mod p) (4.5)
and thus (4.2) follows now from (4.3), (4.4) and (4.5). The result follows, since

P = 0 i<p iP i (u)
is a Mahler polynomial.

Theorem 4.2 can be extended as follows.

Corollary 4.4 Let f : A * → N be a uniformly continuous function. Then, for each positive integer r, there exists a Mahler polynomial P r such that, for all u ∈ A * , f (u) ≡ P r (u) (mod p r ).

Proof. We prove the result by induction on r. For r = 1, the result follows from Theorem 4.2. If the result holds for r, there exists a Mahler polynomial P r such that, for all u ∈ A * , f (u) -P r (u) ≡ 0 (mod p r ). Let g = f -P r . Since g is uniformly continuous, there exists a positive integer n such that if u x ≡ v x (mod p) for |x| n, then g(u) ≡ g(v) (mod p 2r ). It follows that 1 p r g(u) ≡ 1 p r g(v) (mod p r ), and thus 1 p r g is uniformly continuous. Applying Theorem 4.2 to 1 p r g, we get a Mahler polynomial P such that, for all u ∈ A * , 1 p r g(u) ≡ P (u) (mod p)

Setting P r+1 = P r + p r P , we obtain finally f (u) ≡ P r+1 (u) (mod p r+1 ) which concludes the proof.

We now come back to the proof of Theorem 4.1. For each positive integer r, there exists a Mahler polynomial P r such that, for all u ∈ A * , f (u) ≡ P r (u) (mod p r ). Using (2.3) to compute explicitly the coefficients f -P r , v , we obtain f -P r , v ≡ 0 (mod p r ) Since P r is a polynomial, there exists an integer n r such that for all v ∈ A * such that |v| n r , P r , v = 0. It follows | f, v | p p -r and thus lim |v|→∞ | f, v | p = 0.

An extension of Mahler's theorem

Mahler's theorem is often presented as an interpolation result (see for instance [9, p. 57]). This can also be extended to functions from words to integers. Given a family of integers (c v ) v∈A * , one can ask whether there is a (uniformly) continuous function f from the free pro-p group to Z such that f (v) = c v . The answer to this question is part of our main result: (3) implies [START_REF] Amice | Interpolation p-adique[END_REF]. Formulas (1.1) and (1.2) give for all u ∈ A * and for all

n ∈ N n<|v| f, v u v p sup n<|v| |u| | f, v | p sup n<|v| | f, v | p Let ε > 0. By (3), there is an integer N such that if |v| N , then | f, v | p < ε.
It follows that for all n N and for all u ∈ A * ,

n<|v| f, v u v p sup n<|v| | f, v | p < ε.
which proves [START_REF] Amice | Interpolation p-adique[END_REF], in view of (2.1).

(2) implies (1). Theorem 4.1 shows that every Mahler polynomial is uniformly continuous. Now, if a sequence of uniformly continuous functions converges uniformly, then its limit is uniformly continuous as well. Thus (2) implies that f is uniformly continuous.

(3) implies (4). A straightforward induction shows that ∆

u (∆ v f ) = ∆ uv f holds for all u, v ∈ A * . Hence ∆ v f, u = ∆ u (∆ v f )(1) = ∆ uv f (1) = f, uv and Theorem 2.3 yields ∆ v f p = max u∈A * | ∆ v f, u | p = max u∈A * | f, uv | p (5.1)
Let ε > 0. By (3), there exists some N ∈ N such that | f, w | p < ε whenever |w| N . Thus, whenever |v| N , we get | f, uv | p < ε and so ∆ v f p < ε by (5.1). Therefore (4) holds.

(4) implies (3) follows immediately from (5.1).

(1) implies ( 5) since F p (A) (respectively Z p ) is the completion of A * (respectively Z) for d p .

(5) implies (1). Since F p (A) is compact, the continuous extension of f to F p (A) is uniformly continuous and its restriction f to A * is uniformly continuous as well.

The commutative case

A Mahler expansion can also be obtained for functions from N k into Z. First observe that the family of functions from

N k into N r 1 r 2 • • • r k | (r 1 , . . . , r k ) ∈ N k
is locally finite. Thus, given a family (m r ) r∈N k , the formula

f (n 1 , . . . , n k ) = r∈N k m r n 1 r 1 • • • n k r k defines a function f : N k → Z.
) ∈ N k , f (n 1 , . . . , n k ) = r∈N k f, r n 1 r 1 • • • n k r k
Furthermore, the coefficients are given by

f, r = r 1 i 1 =0 . . . r k i k =0 (-1) (r 1 + ... +r k )+(i 1 + ... +i k ) r 1 i 1 • • • r k i k f (i 1 , . . . , i k )
The next lemma relates commutative and noncommutative Mahler expansions. Given a finite alphabet A, we denote by c the map from A * to N |A| which maps every word u onto its commutative image c(u) = (|u| a ) a∈A .

Lemma 6.2 Let r∈N k f, r - r 1 • • • - r k be the Mahler expansion of a function f from N k into Z. Let A = {a 1 , . . . , a k } be a k-letter alphabet. Then f, r = f • c, a r 1 1 • • • a r k k for every r ∈ N k .
Proof. One has

f (n 1 , . . . , n k ) = f (c(a n 1 1 • • • a n k k )) = v∈A * f • c, v a n 1 1 • • • a n k k v = n 1 r 1 =0 . . . n k r k =0 f • c, a r 1 1 • • • a r k k n 1 r 1 • • • n k r k
The uniqueness of the Mahler expansion stated in Proposition 6.1 yields the equality f, r = f • c, a r 1 1 . . . a r k k for every r ∈ N k .

We can now deduce a new proof for the hard implication of the generalization of Mahler's Theorem to N k (see [START_REF] Amice | Interpolation p-adique[END_REF] or [18, Section 4. The opposite implication is the easy one (see [START_REF] Amice | Interpolation p-adique[END_REF]) and can be deduced from the fact that if n 1 ≡ n 2 (mod p s+νp(m!) ) then n 1 m ≡ n 2 m (mod p s ).

Details are omitted.

Back to language theory

Our original motivation was the study of regularity-preserving functions f from A * to B * , in the following sense: if L is a regular language of B * , then f -1 (L) is a regular language of A * . More generally, we are interested in functions preserving a given variety of languages V: if L is a language of V, then f -1 (L) is also a language of V. There is an important literature on the regular case [START_REF] Stearns | Regularity preserving modifications of regular expressions[END_REF][START_REF] Kosaraju | Regularity preserving functions[END_REF][START_REF] Seiferas | Regularity-preserving relations[END_REF][START_REF] Pin | Operations and transductions that preserve rationality[END_REF][START_REF] Pin | Une application de la représentation matricielle des transductions[END_REF][START_REF] Pin | A topological approach to transductions[END_REF][START_REF] Berstel | Operations preserving recognizable languages[END_REF]. A remarkable contribution to this problem can be found in [START_REF] Reutenauer | Variétés et fonctions rationnelles[END_REF], where a characterization of sequential functions preserving star-free languages (respectively group-languages) is given. A similar problem was also recently considered for formal power series [START_REF] Droste | On transformations of formal power series[END_REF].

An open problem is to characterize the functions from A * to B * preserving the p-group languages. Theorems 1.4 and 5.1 give a solution when B is a one-letter alphabet.

Theorem 7.1 Let f be a function from A * to N. Then the following conditions are equivalent:

(1) for each p-group language L of N, f -1 (L) is a p-group language of A * , (2) f is uniformly continuous for d p ,

lim |v|→∞ ∆ v f p = 0.

It would be interesting to find a suitable extension of this result for functions from A * to B * , for any alphabet B. Another challenge would be to extend other results from p-adic analysis to the word case.

Lemma 1 . 2

 12 The p-adic metric and the pro-p metric coincide on N k and on Z k .

  with the help of (1.1) and (1.2),| f, w | p max |f (w)| p , max 0 |v|<|w| | f, v | p Now, for 0 |v| < |w|, one has | f, v | p < | f, w | pby the minimality of |w| and thus max 0 |v|<|w| | f, v | p < | f, w | p . It follows that | f, w | p |f (w)| p f p and hence max v∈A * | f, v | p = f p as claimed.

Theorem 4 . 1

 41 Let v∈A * f, v - v be the Mahler's expansion of a function f : A * → Z. Then f is uniformly continuous if and only if lim |v|→∞ | f, v | p = 0. Proof. Suppose that lim |v|→∞ | f, v | p = 0 and let r ∈ N. Then there exists s

Theorem 5 . 1

 51 Let v∈A * f, v - v be the Mahler's expansion of a function f : A * → Z. Then the following conditions are equivalent:(1) f is uniformly continuous for d p , (2) the sequence of partial sums0 |v| n f, v - v converges uniformly to f , (3) lim |v|→∞ | f, v | p = 0, (4) lim |v|→∞ ∆ v f p = 0,(5) f admits a continuous extension from F p (A) to Z p . Proof. (1) implies (3) follows from Theorem 4.1.

  2.4]). Given r ∈ N k , let us write |r| for |r 1 | + . . . + |r k |. Corollary 6.3 (Amice [2]) Let r∈N k f, r r 1 • • • r k be the Mahler expansion of a function f from N k into Z. Then f is uniformly continuous for d p if and only if lim |r|→∞ | f, r | p = 0. Proof. Assume that f is uniformly continuous. Let A = {a 1 , . . . , a k }. Since c : A * → N k is a morphism, it is uniformly continuous by Corollary 1.5 and so is the function f • c. By Theorem 4.1, we obtain lim |v|→∞ | f • c, v | p = 0. In particular, lim |r|→∞ | f • c, a r 1 1 • • • a r k k | p = 0 and so lim |r|→∞ | f, r | p = 0 by Lemma 6.2.

  Proposition 3.2 Mahler polynomials form a subring of the ring of all functions from A * to Z for addition and multiplication. Proof. It is clear that the difference of two Mahler polynomials is a Mahler polynomial. Further Proposition 3.1 shows that Mahler polynomials are closed under product.

  Conversely, each function from N k to Z admits a unique Mahler expansion, a result proved in a more general setting in[START_REF] Amice | Interpolation p-adique[END_REF][START_REF] Ahlswede | Approximation of continuous functions in p-adic analysis[END_REF]. Let f : N k → Z be an arbitrary function. Then there exists a unique family f, r r∈N k of integers such that, for all (n 1 , . . . , n k

	Proposition 6.1

Acknowledgement

The authors would like to thank also Daniel Barsky and Gilles Christol for pointing out several references.

 [START_REF] Ahlswede | Approximation of continuous functions in p-adic analysis[END_REF]Work supported by the project ANR 2010 BLAN 0202 02 FREC, the European Regional Development Fund through the programme COMPETE and the Portuguese Government through FCT (Fundação para a Ciência e a Tecnologia) under the project PEst-C/MAT/UI0144/2011.


Thierry Denoeux
email: thierry.denoeux@hds.utc.fr

Maximum likelihood estimation from fuzzy data using the EM algorithm

Keywords: Statistics, fuzzy data analysis, estimation, maximum likelihood principle, regression, mixture models

A method is proposed for estimating the parameters in a parametric statistical model when the observations are fuzzy and are assumed to be related to underlying crisp realizations of a random sample. This method is based on maximizing the observeddata likelihood defined as the probability of the fuzzy data. It is shown that the EM algorithm may be used for that purpose, which makes it possible to solve a wide range of statistical problems involving fuzzy data. This approach, called the Fuzzy EM (FEM) method, is illustrated using three classical problems: normal mean and variance estimation from a fuzzy sample, multiple linear regression with crisp inputs and fuzzy outputs, and univariate finite normal mixture estimation from fuzzy data.

Introduction

Recent years have seen a surge of interest in extending statistical inferential procedures to fuzzy data [START_REF] Coppi | The fuzzy approach to statistical analysis[END_REF]. Basically, two views of fuzzy data have been proposed [START_REF] Gebhardt | Fuzzy set-theoretic methods in statistics[END_REF].

The first approach assumes the data to be intrinsically fuzzy and uses the mathematical formalism of fuzzy random variables, defined as mappings from a probability space to the set of fuzzy subsets of R or R n [START_REF] Gil | Overview on the development of fuzzy random variables[END_REF][START_REF] Kruse | Statistics with vague data[END_REF][START_REF] Kwakernaak | Fuzzy random variables. Part I: definitions and theorems[END_REF][START_REF] Puri | Fuzzy random variables[END_REF], with certain measurability properties. This model fits well with a physical interpretation of fuzzy data, in which a fuzzy datum is regarded as an existing object, not necessarily connected to any underlying precise variable [START_REF] Gebhardt | Fuzzy set-theoretic methods in statistics[END_REF]. Recent work on applying this approach to estimation and hypothesis testing problems may be found, for example, in [START_REF] Colubi | Statistical inference about the means of fuzzy random variables: Applications to the analysis of fuzzy-and real-valued data[END_REF][START_REF] González-Rodríguez | Multi-sample test-based clustering for fuzzy random variables[END_REF][START_REF] González-Rodríguez | Bootstrap techniques and fuzzy random variables: Synergy in hypothesis testing with fuzzy data[END_REF].

The second approach is based on an epistemic interpretation of fuzzy data, which are assumed to "imperfectly specify a value that is existing and precise, but not measurable with exactitude under the given observation conditions" [13, page 316]. In this model, a fuzzy datum is thus seen as a possibility distribution associated to a precise realization of a random variable that has been only partially observed. This viewpoint will be adopted in this paper.

Under the latter interpretation of fuzzy data, a conceptually simple approach is to "fuzzify" standard statistical computations using Zadeh's extension principle [START_REF] Zadeh | Fuzzy sets[END_REF]. This approach makes it possible to compute fuzzy versions of classical statistics, such as the mean, variance or correlation coefficients [START_REF] Denoeux | Nonparametric rank-based statistics and significance tests for fuzzy data[END_REF][START_REF] Viertl | Univariate statistical analysis with fuzzy data[END_REF], and even to extend significance tests using fuzzy p-values [START_REF] Denoeux | Nonparametric rank-based statistics and significance tests for fuzzy data[END_REF][START_REF] Viertl | Testing hypotheses with fuzzy data: The fuzzy p-value[END_REF]. However, it usually implies solving complex nonlinear optimization problems, except in very simple cases. Furthermore, fuzzy estimates are often too complex objects to be easily handled and interpreted by end-users.

When a parametric statistic model is postulated, an alternative approach to parameter estimation based on fuzzy data relies on the extension of the likelihood function using Zadeh's definition of the probability of a fuzzy event [START_REF] Zadeh | Probability measures of fuzzy events[END_REF]. A maximum likelihood estimate (MLE) of the parameter of interest may then be defined as a crisp value maximizing the probability of observing the fuzzy data. This model was initially suggested in [START_REF] Okuda | A formulation of fuzzy decision problems with fuzzy information using probability mesures of fuzzy events[END_REF]. In the 1980's, it was studied by Gil and her colleagues in conjunction with fuzzy information systems [START_REF] Gil | An operative extension of the likelihood ratio test from fuzzy data[END_REF][START_REF] Gil | The likelihood ratio test for goodness of fit with fuzzy experimental observations[END_REF][START_REF] Gil | The minimum inaccuracy estimates in χ 2 tests for goodness of fit with fuzzy observations[END_REF]. These authors showed that the main properties of the MLE (equivariance under one-to-one transformations, asymptotic normality and asymptotic efficiency) are preserved in the fuzzy case under very general conditions [13, pages 327-333].

Although conceptually appealing and well-founded theoretically, the maximum likelihood (ML) approach to statistical inference from fuzzy data has not been widely used, essentially because of the difficulty of putting it to work in realistic statistical problems. As noted by Gebhardt et al. [13, page 329]: "In spite of these valuable properties and although the application of the maximum likelihood principle to some examples is very simple, to determine solutions in most situations becomes too complex".

In this paper, we propose to solve the above problem using the Expectation-Maximization (EM) algorithm. This algorithm, formalized by Dempster et al. [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF], provides a very general mechanism for computing MLEs from incomplete data. This method has been successfully applied to a wide range of problems involving partially observed data such as inference from censored or grouped data, finite mixture models, factor analysis, etc [START_REF] Mclachlan | The EM Algorithm and Extensions[END_REF]. As will be shown, the EM algorithm can be adapted to handle estimation problems involving fuzzy data, viewed as imperfectly specifying underlying crisp values. As will be shown, this approach makes the ML principle practically applicable to solve a wide range of statistical problems involving fuzzy data.

The rest of this paper is organized as follows. Section 2 first presents in greater detail the problem addressed in this paper, trying to clarify the considered semantics of fuzzy data. Section 3 recalls the EM algorithm and demonstrates its application in the case of fuzzy data. Section 4 then describes three applications: normal mean and variance estimation of a normal distribution, multiple regression with crisp input and fuzzy output data, and univariate normal mixture estimation. Finally, Section 5 concludes the paper.

Problem description

The problem addressed in this paper may be described as follows. We assume the existence of a random vector X, referred to as the complete data vector, taking values in a sample space X and describing the result of a random experiment. The probability density function (p.d.f.) of X is denoted by g(x; Ψ), where Ψ = (Ψ 1 , . . . , Ψ d ) ′ is a vector of unknown parameters with parameter space Ω. Although X will be generally assumed to be a continuous random vector (unless otherwise specified), g ( x; Ψ) can still be viewed as a p.d.f. in the case where X is discrete by the adoption of the counting measure.

Let x be a realization of X. If x was known exactly, we could compute the maximum likelihood estimate (MLE) of Ψ as any value maximizing the completedata likelihood function

L(Ψ; x) = g(x; Ψ). (1)
In this paper, we consider the more difficult problem where x is not observed precisely, and only partial information about x is available in the form of a fuzzy subset x of X , with Borel measurable membership function µ x : X → [0, 1]. Before tackling the problem of parameter estimation from such data, we have to examine carefully the nature of the partial information about x, in relation to the underlying random experiment. This issue will be addressed in Subsection 2.1. A generalized likelihood function will then be introduced and discussed in Subsection 2.2. Finally, another model, which leads formally to the same likelihood function with a completely different interpretation, will be discussed in Subsection 2.3.

Interpretation of fuzzy data

In the model considered here, the fuzzy observation x will be understood as encoding the observer's partial knowledge about the realization x of random vector X. In this setting, the membership function µ x is seen as a possibility distribution [START_REF] Zadeh | Fuzzy sets as a basis for a theory of possibility[END_REF][START_REF] Dubois | Possibility Theory: An approach to computerized processing of uncertainty[END_REF] interpreted as a soft constraint on the unknown quantity x. The fuzzy set x can be considered to be generated by a two-step process:

1. A realization x is drawn from X;

2. The observer encodes his/her partial knowledge of x in the form of a possibility distribution µ x .

It must be stressed that, in this model, only step 1 is considered to be a random experiment.

Step 2 implies gathering information about x and modeling this information as a possibility distribution. It is not assumed that this process may be repeated indefinitely in the same experimental conditions, i.e., it is not considered as a random experiment. Consequently, x is not considered as a realization of a fuzzy random variable in this paper. It must also be recognized that two kinds of uncertainty are involved in this setting:

• Step 1 is associated with aleatory uncertainty, which is due to the random nature of the data generation process; this uncertainty cannot be reduced before the experiment, and disappears after the experiment has been performed.

•

Step 2 is associated by epistemic uncertainty, which is related to the observer's state of knowledge. This uncertainty can sometimes be reduced by gathering additional information about x.

Example 1 Assume that n skull fragments have been found in some archaeological site, and we are interested in the volumes of these skulls. The unknown volume x i of skull i may be regarded as a realization of a random variable X i induced by random sampling from a total population of skulls (corresponding to a certain region and period of interest). The complete data x = (x 1 , . . . , x n ) is thus a realization from a random vector X = (X 1 , . . . , X n ). As only fragments are available, archaeologists have to use various indirect methods and hypotheses to assess the volume of each skull. Assume that, following a procedure reported, e.g., in [START_REF] Denoeux | Nonparametric rank-based statistics and significance tests for fuzzy data[END_REF] and [START_REF] González-Rodríguez | Fuzzy data treated as functional data: A one-way anova test approach[END_REF], two intervals are determined for each skull i:

• an interval [a i , d i ] certainly containing x i ;
• an interval [b i , c i ] containing highly plausible values for x i .

This information may be encoded as a trapezoidal fuzzy number x i = (a i , b i , c i , d i ) with support [a i , d i ] and core [b i , c i ], interpreted as a possibility distribution constraining the unknown value x i . Information about x may be represented by the joint possibility distribution

µ x (x) = µ x 1 (x 1 ) × • • • × µ x 1 (x n ). (2)
The epistemic uncertainty modeled by this possibility distribution might be reduced by carrying out further analysis of the skull fragments. This uncertainty is clearly of a different nature than the aleatory uncertainty involved in the sampling process.

We may note that, in the special case where the observed data x is crisp, it corresponds to some form of censored data as frequently encountered in statistics. Usually, censored data are handled by postulating some stochastic censoring mechanism, and maximizing the observed data likelihood using, e.g., the EM algorithm. As we will see, our approach is computationally very similar, although no random censoring mechanism is assumed.

The relationship between the unknown crisp realization x and the fuzzy observation x can be described by considering the random set view of fuzzy sets (see, e.g., [10, page 46]). In this view, x is seen as being induced by a uniform probability measure on [0, 1] and the set-valued mapping α → α x, where α x denotes the α-cut of x. If the observer knew that α = α 0 , then he/she would assert that x ∈ α 0 x. However, the observer's knowledge of α is uncertain and is described by the uniform probability distribution. In this model, µ x (u) is the observer's subjective probability that the available evidence can be interpreted as specifying an interval that contains u. We should stress here that, in the model adopted in this paper, the known subjective probability distribution on α and the unknown objective probability distribution on X are assumed to be of different natures; as such they will be treated in completely different ways in the approach described below.

Generalized likelihood function

Once x is given, and assuming its membership function to be Borel measurable, we can compute its probability according to Zadeh's definition of the probability of a fuzzy event [START_REF] Zadeh | Probability measures of fuzzy events[END_REF] (see Appendix A.1). By analogy with (1), the observed-data likelihood can then be defined as:

L(Ψ; x) = P ( x; Ψ) = X µ x (x)g(x; Ψ)dx. (3)
To understand the meaning of P ( x; Ψ) and, consequently, of L(Ψ; x), we can remark1 , as done by Höhle [START_REF] Höhle | Masse auf unscharfen Mengen[END_REF], that P ( x; Ψ) can also be written as:

P ( x; Ψ) = 1 0 P ( α x; Ψ)dα,
where α x denotes the α-cut of x (see also [10, page 52]). Using the random set view of fuzzy sets outlined in the previous section, P ( x; Ψ) can thus be seen as the average value of P ( α x; Ψ) over the random set underlying x.

The above remark can be reformulated in terms of likelihoods. If we knew that x ∈ α x, then the likelihood function would be:

L(Ψ; α x) = P ( α x; Ψ).
Assuming α to be uncertain, and our belief on α to be described by a uniform probability distribution on [0, 1], the likelihood averaged over all values of α is

1 0 L(Ψ; α x)dα = L(Ψ; x).
In the special case where the complete data x = (x 1 , . . . , x n ) is a realization of an independent identically distributed (i.i.d.) random vector X = (X 1 , . . . , X n ), and assuming the joint membership function µ x to be decomposable as in [START_REF] Colubi | Statistical inference about the means of fuzzy random variables: Applications to the analysis of fuzzy-and real-valued data[END_REF], the likelihood function (3) can be written as a product of n terms:

L(Ψ; x) = n i=1 µ x i (x)g(x; Ψ)dx, (4)
and the observed-data log likelihood is

log L(Ψ; x) = n i=1 log µ x i (x)g(x; Ψ)dx. (5)

Fuzzy Information System model

Although the model outlined in Subsection 2.1 will be adopted throughout this paper, it is interesting to compare it with another model that has been extensively used in the fuzzy statistics literature in the 1980's (see, e.g., [START_REF] Gil | The minimum inaccuracy estimates in χ 2 tests for goodness of fit with fuzzy observations[END_REF][START_REF] Gil | An operative extension of the likelihood ratio test from fuzzy data[END_REF][START_REF] Gil | The likelihood ratio test for goodness of fit with fuzzy experimental observations[END_REF]). In this alternative model, we assume the existence of fuzzy information system (FIS) , defined as a fuzzy partition F = { ξ 1 , . . . , ξ K } of X , i.e., a set of K fuzzy subsets of X verifying the orthogonality condition:

K k=1 µ ξ k (x) = 1,
where µ ξ k denotes the membership function of ξ k . When observing a realization x of X, the experimenter selects the element of F that best describes his/her perception of x. If we make the (strong) assumption that, conditionally on x, each ξ k is selected with probability µ ξ k (x), then the overall probability of selecting ξ k is

P ( ξ k ) = X µ ξ k (x)g(x; Ψ)dx = E Ψ µ ξ k (X) , (6)
assuming µ ξ k to be Borel measurable. We observe that ( 6) coincides with Zadeh's definition of the probability of a fuzzy event [START_REF] Zadeh | Probability measures of fuzzy events[END_REF]. However, in the particular model considered here, P ( ξ k ) is the objective probability (i.e., the limit frequency) of selecting ξ k as the best fuzzy description of x, based on the experimenter's perception. The usual case of categorized data is recovered as a special case when the elements of F are crisp subsets of X .

Having observed x, the likelihood function is

L(Ψ) = P ( X = x; Ψ) = X µ x (x)g(x; Ψ)dx, (7)
which is formally identical to (3), with a completely different interpretation. Consequently, the problem of finding the MLE in this model can be addressed in exactly the same way as in the model considered in Subsection 2.1, which will be adopted hereafter for clarity, unless otherwise specified. A solution to this problem is provided in the next section.

The Fuzzy EM method

The problem tackled in this section consists in maximizing the observed-data log likelihood, as defined in Section 2.2, in situations where the observed data is fuzzy and can be seen as an incomplete specification of a complete data vector x. The EM algorithm is a broadly applicable mechanism for computing MLEs from incomplete data, in situations where ML estimation would be straightforward if complete data were available [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]. In its basic form, this algorithm assumes that the observed data corresponds to a crisp set of possible values for the complete data. This algorithm thus needs to be adapted to handle fuzzy data. The EM algorithm will first be recalled in Subsection 3.1, and its application to fuzzy data, referred to as the Fuzzy EM (FEM) method, will then be presented in Subsection 3.2.

The EM algorithm

With the same notations as in Section 2, let us assume that we have a random vector X with p.d.f. g(x; Ψ). A realization x has been drawn from X, but it is incompletely observed. The observed data consists in a subset X of X such that x ∈ X. The observed-data likelihood is

L(Ψ; X) = X g(x; Ψ)dx. (8)
The EM algorithm approaches the problem of maximizing the observed-data log likelihood log L(Ψ; X) by proceeding iteratively with the complete-data log likelihood log L(Ψ; x) = log g(x; Ψ). Each iteration of the algorithm involves two steps called the expectation step (E-step) and the maximization step (M-step).

The E-step requires the calculation of

Q(Ψ, Ψ (q) ) = E Ψ (q) [log L(Ψ; X) | X] ,
where Ψ (q) denotes the current fit of Ψ at iteration q, and E Ψ (q) denotes expectation using the parameter vector Ψ (q) . The M-step then consists in maximizing Q(Ψ, Ψ (q) ) with respect to Ψ over the parameter space Ω, i.e., finding Ψ (q+1) such that

Q(Ψ (q+1) , Ψ (q) ) ≥ Q(Ψ, Ψ (q) )
for all Ψ ∈ Ω. The E-and M-steps are iterated until the difference L(Ψ (q+1) ; X) -L(Ψ (q) ; X) becomes smaller than some arbitrarily small amount.

It is proved in [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] that the observed-data likelihood L(Ψ; X) is not decreased after an EM iteration, that is, L(Ψ (q+1) ; X) ≥ L(Ψ (q) ; X)

for q = 0, 1, 2, . . .. Hence, convergence to some value L * is ensured as long as the sequence L(Ψ (q) ; X) for q = 0, 1, 2, . . . is bounded from above. As noted in [26, page 85], L * is, in most practical applications and except in pathological cases, a local maximum of the incomplete data log likelihood L(Ψ; X).

Remark 1

In [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF], Dempster et al. postulate the existence of a sample space Y, and a many-to-one mapping ϕ from X to Y. To each observed y ∈ Y thus corresponds a subset X = ϕ -1 (y) of X . If y is a realization of a random variable Y, then the set X also becomes random. To be consistent with the model described in Subsection 2.1, the set X is not considered to be random in this section. This is why Y and ϕ are not needed here. The EM algorithm remains unchanged under this interpretation.

Application to fuzzy data

Let us now assume, as we did in Section 2, that the observed data consist in a fuzzy subset x of X , with Borel measurable membership function µ x . The observed data likelihood is now given by (3), which is a direct generalization of (8). To maximize this function, we propose to adapt the EM algorithm as follows. Let the E-step now consist in the calculation of

Q(Ψ, Ψ (q) ) = E Ψ (q) [log L(Ψ; X) | x] (9) = µ x (x) log[L(Ψ; x)]g(x; Ψ (q) )dx L(Ψ (q) ; x) , (10)
where the expectation of log L(Ψ; X) is now taken with respect to the conditional p.d.f. of x given x, using parameter vector Ψ (q) : g(x| x; Ψ (q) ) = µ x (x)g(x; Ψ (q) ) µ x (u)g(u; Ψ (q) )du .

(See [START_REF] González-Rodríguez | Fuzzy data treated as functional data: A one-way anova test approach[END_REF] in Appendix A.1 for the general expression of the p.d.f. of a random variable conditionally on a fuzzy event).

The M-step is unchanged and requires the maximization of Q(Ψ, Ψ (q) ) with respect to Ψ. The proposed algorithm alternately repeats the E-and M-steps above until the increase of observed-data likelihood becomes smaller than some threshold.

To show that the above algorithm actually maximizes the observed-data likelihood (3), we may remark that it is formally equivalent to the standard EM algorithm recalled in Subsection 3.1, applied to a different statistical model with crisp observations, in which the fuzzy data play the role of known parameters2 . More precisely, let Y denote a Bernoulli random variable such that P (Y = 1|X = x) = µ x (x), and let Z = (X, Y ) be the "complete data". Its density when Y = 1 is

g(x, 1; Ψ) = µ x (x)g(x; Ψ).
Assume that we have observed Y = 1. Then the likelihood function based on this observation is

L(Ψ; Y = 1) = P (Y = 1; Ψ) = g(x, 1; Ψ)dx = µ x (x)g(x; Ψ)dx,
which is identical to (3). This likelihood may be maximized using the classical EM algorithm. The E-step of this algorithm computes

Q ′ (Ψ, Ψ (q) ) = E Ψ (q) [log L(Ψ; Z) | Y = 1] (11) = µ x (x) log[µ x (x)g(x; Ψ)]g(x; Ψ (q) )dx µ x (u)g(u; Ψ (q) )du (12) = µ x (x) log[µ x (x)]g(x; Ψ (q) )dx µ x (u)g(u; Ψ (q) )du + Q(Ψ, Ψ (q) ). ( 13
)
As the first term in the right-hand side of (13) does not depend of Ψ, it need not be calculated. Maximizing Q ′ (Ψ, Ψ (q) ) is equivalent to maximizing Q(Ψ, Ψ (q) ), which shows the formal equivalence between the computations under the two models.

The above line of reasoning shows that the Fuzzy EM method introduced in this section can be considered, from a purely formal point of view, as the application of the classical EM algorithm to a different statistical model. An immediate consequence is that the proposed procedure generates a nondecreasing sequence of observed-data likelihood values, which converges to some limit if it is bounded.

Applications

In the previous section, we have shown that the EM algorithm can be adapted to compute MLEs from fuzzy data. In this section, we will now demonstrate the application of this method to three classical estimation problems: mean and variance of univariate normal data (Subsection 4.1), multiple regression with crisp inputs and fuzzy outputs (Subsection 4.2) and univariate normal mixture estimation from fuzzy data (Subsection 4.3).

Mean and variance of univariate normal data

Let us assume that the complete data x = (x 1 , . . . , x n ) ′ is a realization of an i.i.d. random sample from a normal distribution with mean m and standard deviation σ. The observed data is supposed to take the form of a fuzzy subset x of R n with membership function

µ x (x) = n i=1 µ x i (x i ) (14)
where x = (x 1 , . . . , x n ) ′ and the x i , i = 1, . . . , n are fuzzy numbers. The complete-data p.d.f is

g(x; Ψ) = n i=1 g(x i ; Ψ),
where Ψ = (m, σ) ′ and

g(x i ; Ψ) = 1 σ √ 2π exp - (x i -m) 2 2σ 2 .
The complete-data log likelihood is thus

log L(Ψ; x) = n i=1 log g(x i ; Ψ) = - n 2 log(2π) -n log σ - 1 2σ 2 n i=1 (x i -m) 2 = - n 2 log(2π) -n log σ - 1 2σ 2 n i=1 x 2 i -2m n i=1 x i + nm 2 .
The observed-data log likelihood is given by [START_REF] Coppi | The fuzzy approach to statistical analysis[END_REF]. At iteration q + 1, the E-step of the EM algorithm requires the calculation of

Q(Ψ, Ψ (q) ) = E Ψ (q) [log L(Ψ; X)| x] = - n 2 log(2π) -n log σ - 1 2σ 2 n i=1 E Ψ (q) (X 2 i | x i ) -2m n i=1 E Ψ (q) (X i | x i ) + nm 2 .
When the x i are trapezoidal fuzzy numbers, the conditional expectations α

(q) i = E Ψ (q) (X 2 i | x i ) and β (q) i = E Ψ (q) (X i | x i
) can be computed using Equation ( 23)-( 28) in Appendix A.2.

The M-step of the EM algorithm involves maximizing Q(Ψ, Ψ (q) ) with respect to Ψ. This is easily achieved by solving the likelihood equations. From

∂Q(Ψ, Ψ (q) ) ∂m = -2 n i=1 β (q) i + 2nm,
we get

m (q+1) = 1 n n i=1 β (q) i . Now, ∂Q(Ψ, Ψ (q) ) ∂σ = - n σ + 1 σ 3 n i=1 α (q) i -2m n i=1 β (q) i + nm 2 .
Substituting m (q+1) for m in the equation

∂Q(Ψ, Ψ (q) ) ∂σ = 0
and solving for σ, we get

σ (q+1) = 1 n n i=1 α (q) i -(m (q+1) ) 2 .
The EM algorithm for this problem is summarized in Algorithm 1.

Example 2 To illustrate the above algorithm, we consider real data3 collected during the experiment reported in [START_REF] González-Rodríguez | Fuzzy data treated as functional data: A one-way anova test approach[END_REF]. In this experiment, subjects were asked about their perception of the relative length of different line segments with respect to a fixed longer segment that was used as a standard for comparison. The subjects had to describe their perception by providing the support and the core of trapezoidal fuzzy numbers. In the version of the database used in this example, there were 17 subjects, and segments with various relative lengths were presented in random order, with three repetitions for each length. We considered the n = 17 × 3 = 51 fuzzy numbers corresponding to the true relative length 61.47. Figure 1 shows the first 10 of these trapezoidal fuzzy numbers. The model described in Subsection 2.1 can be applied to this data if we assume that each reported fuzzy number x i acts a soft constraint on the observer's perceived relative length x i , considered to be a realization of a random variable X i = m + ε, where m = 61.47 and ε ∼ N (0, σ 2 ). Note that this model is different from that used in [START_REF] González-Rodríguez | Fuzzy data treated as functional data: A one-way anova test approach[END_REF], in which x i is treated as a realization of a fuzzy random variable. A detailed comparison between these two models would go beyond the scope of this paper.

Algorithm 1 EM algorithm for the estimation of the mean and standard deviation of a normal population using a fuzzy random sample. Input: x 1 , . . . , x n , ǫ > 0. Output: m, σ q ← 0 Initialization: pick m (0) and σ (0) L(Ψ (0) ; x) ← n i=1 log E Ψ (0) (µ x i (X i )) % using ( 21)-( 22) repeat % E-step 21)-( 28) As initial estimates of m and σ, we used the sample mean and standard deviations computed over the centers of cores of each trapezoidal fuzzy numbers, equal to m (0) = 60.5455 and σ (0) = 6.1423. The EM converged after 10 iterations. The stopping criterion was based on the relative change of the log likelihood, with a tolerance value ǫ = 10 -8 . The final MLEs were m = 61.0155 and σ = 5.2494.

for i = 1 : n do α (q) i ← E Ψ (q) (X 2 i | x i ) % using (
β (q) i ← E Ψ (q) (X i | x i ) % using (21)-(25) end for % M-step m (q+1) ← 1 n n i=1 β (q) i σ (q+1) ← 1 n n i=1 α (q) i -(m (q+1) ) 2 L(Ψ (q+1) ; x) ← n i=1 log E Ψ (q+1) (µ x i (X i )) % using (21)-(22) D ← (L(Ψ (q+1) ; x) -L(Ψ (q) ; x))/|L(Ψ (q) ; x)| q ← q + 1 until D < ǫ m ← m (q) ; σ ← σ (q)
Figure 2 shows a contour plot of the observed data log likelihood function as well as the trajectory in parameter space from the initial parameter values (m (0) , σ (0) ) to the final MLE ( m, σ). We can check that the MLE corresponds in this case to a global maximum of the observed data log likelihood. In more complex problem, the algorithm may be trapped in local maxima. It is then necessary to start it several times with different random initial conditions.

Example 3 To illustrate experimentally the asymptotic behavior of the MLEs of m and σ under the model described in Subsection 2.3 (which lends itself easily to random simulation), we performed the following experiment. We considered i.i.d random samples of size n from the standard normal distribution. Each realization of x was fuzzified using the FIS shown in Figure 3, and the MLEs m and σ for the fuzzy sample were computed using the FEM method. For each value of n (ranging from 10 to 1000), the whole procedure was repeated 100 times and the expectation and standard error of m and σ were estimated. The results are shown in Figures 4 and5, which illustrate the convergence of both estimators towards the true parameter values. We note that similar results are obtained if each fuzzy number ξ k in the FIS is replaced by its cut at level 0.5, defining a crisp partition of the real line.

Multiple regression

Let us now consider the case of multiple regression with crisp inputs and fuzzy outputs. This is an old problem in fuzzy data analysis (see, e.g., [START_REF] Coppi | Management of uncertainty in statistical reasoning: The case of regression[END_REF][START_REF] Näther | Regression with fuzzy random data[END_REF] and references therein). As we will show in this section, the FEM method provides a new and well motivated solution to this problem, with very simple implementation.

The fuzzy data x will be assumed to have the same form as in the previous section. However, we now assume that each component x i of the complete data vector x is a realization of a normal random variable X i with mean u ′ i b and standard deviation σ, where u i = (1, u i1 , . . . , u i,p-1 ) ′ is a constant p-dimensional input vector and b is vector of p coefficients4 . The complete parameter vector is thus Ψ = (b, σ) ′ .

Denoting by U the matrix of n rows and p columns, with row i equal to u i , and assuming independence between the X i , the observed data vector x is multivariate Gaussian with mean Ub and variance σ 2 I p , where I p denotes the p-dimensional identity matrix. The complete-data log likelihood is thus:

log L(Ψ; x) = - n 2 log(2π) -n log σ - 1 2σ 2 (x -Ub) ′ (x -Ub) = - n 2 log(2π) -n log σ - 1 2σ 2 x ′ x -2b ′ U ′ x + b ′ U ′ Ub .
Taking the expectation of log L(Ψ; X) conditionally on the observed fuzzy data and using the fit Ψ (q) of Ψ to perform the E-step, we get

Q(Ψ, Ψ (q) ) = E Ψ (q) [log L(Ψ; X)| x] = - n 2 log(2π) -n log σ - 1 2σ 2 n i=1 E Ψ (q) (X 2 i | x i ) -2b ′ U ′ E Ψ (q) (X| x) + b ′ U ′ Ub . ( 15
)
As before, let us denote α

(q) i = E Ψ (q) (X 2 i | x i ), β (q) i = E Ψ (q) (X i | x i ) and
β (q) = E Ψ (q) (X| x) = (E Ψ (q) (X 1 | x 1 ), . . . , E Ψ (q) (X n | x n )) ′ .
With these notations, (15) becomes:

Q(Ψ, Ψ (q) ) = - n 2 log(2π) -n log σ - 1 2σ 2 n i=1 α (q) i -2b ′ U ′ β (q) + b ′ U ′ Ub .
The M-step requires maximizing Q(Ψ, Ψ (q) ) with respect to Ψ. This can be achieved by differentiating Q(Ψ, Ψ (q) ) with respect to b and σ, which results in

∂Q(Ψ, Ψ (q) ) ∂b = - 1 σ 2 -U ′ β (q) + U ′ Ub and ∂Q(Ψ, Ψ (q) ) ∂σ = - n σ + 1 σ 3 n i=1 α (q) i -2b ′ U ′ β (q) + b ′ U ′ Ub .
Equating these derivatives to zero and solving for b and σ, we get the following unique solution:

b (q+1) = (U ′ U) -1 U ′ β (q)
and

σ (q+1) = 1 n n i=1 α (q) i -2 b (q+1) ′ U ′ β (q) + b (q+1) ′ U ′ Ub (q+1) = 1 n n i=1 α (q) i -2 x (q+1) ′ β (q) + x (q+1) ′ x (q+1) , with x (q+1) = Ub (q+1) .

Example 4

The above algorithm was applied to the perception data reported in [START_REF] González-Rodríguez | Fuzzy data treated as functional data: A one-way anova test approach[END_REF] and already described in Example 2. We considered the responses of subject 12 for the 9 different relative lengths. There were 3 repetitions for each relative length, yielding n = 9 × 3 = 27 observations. Here the true relative length and the perceived length are taken as the independent and dependent variables, respectively. The data are shown in Figure 6. The EM algorithm was initialized as in Example 2, with the same stopping criterion. The initial values were b (0) = (-5.0942, 0.9919) ′ and σ (0) = 6.6158. The algorithm converged to the MLE b = (-3.9832, 0.9916) ′ and σ = 4.9218 in 13 iterations. Example 5 To demonstrate the interest of expressing partial information about illknown data in the form of possibility distributions, as explained in Subsection 2.1, we performed the following experiment. We generated n = 100 values u i from the uniform distribution in [0, 2], and we generated corresponding values x i using the linear regression model with b = (2, 1) ′ and σ = 0.2. To model the situation where only partial knowledge of values x 1 , . . . , x n is available, triangular fuzzy numbers x 1 , . . . , x n were generated as follows:

• For each i, a "guess" x ′ i was randomly generated from a normal distribution with mean x i and standard deviation σ i , were σ i was drawn randomly from a uniform distribution in [0, 0.5];

• x i was defined as the triangular fuzzy number with core

x ′ i and support [x ′ i - 2σ i , x ′ i + 2σ i ].
This procedure simulates the situation where the observer has only approximate knowledge of the data, and can only provide a guess x ′ i and an interval of plausible values

[x ′ i -2σ i , x ′ i + 2σ i ].
It should be emphasized that, although the x i are randomly generated for the purpose of this simulation experiment, the relationship between x i and x i is not considered to be random under the model considered here (as explained in Subsection 2.1).

Three strategies were compared for estimating the parameter vector Ψ = (b, σ) ′ :

1. Using the fuzzy data x 1 , . . . , x n (method 1). 2. Using only 0.5-cuts of the fuzzy data, i.e., the interval data 0.5 x 1 , . . . , 0.5 x n (method 2); 3. Using only the crisp guesses x ′ 1 , . . . , x ′ n (method 3).

For each of these three methods, the L 2 distance Ψ -Ψ between the true parameter vector and its MLE was computed. The whole experiment was repeated 1000 times.

Results are shown as box plots in Figure 8. As we can see, method 1 using fuzzy assessments of ill-known data and the FEM method indeed yields better estimates than those obtained using interval or crisp data.

Univariate normal mixture with common unknown variance

As a third and last illustration of the FEM method, we will consider a more complex situation where a fuzzy sample is seen as an imprecise specification of a crisp sample from a univariate normal mixture with unknown means, variance and mixing proportions. A similar problem was addressed in [START_REF] Hamdan | Mixture model clustering of uncertain data[END_REF] using the standard EM algorithm, with interval data in place of fuzzy data. As before, the observed data vector will be assumed to consist of n fuzzy numbers, denoted here w 1 , . . . , w n . Each fuzzy number w i is interpreted as a possibility distribution constraining an unknown value w i . The joint possibility distribution of the n values w 1 , . . . , w n is identified with the fuzzy subset w of R n with membership function

µ w (w) = n i=1 µ w i (w i ),
for all w = (w 1 , . . . , w n ). The n values w 1 , . . . , w n are assumed to be a realization of an i.i.d. random sample W 1 , . . . , W n from a finite normal mixture with p.d.f.

g(w; Ψ) = g i=1 π k g k (w; θ k ),
where g k (w; θ k ) is a normal p.d.f. with parameters θ k = (m k , σ k ) ′ , π k is the mixing proportion of the k-th component, g is the number of components, and

Ψ = (θ ′ 1 , . . . , θ ′ g , π 1 , . . . , π g ) ′
is the vector of parameters. Using a classical device when handling finite mixture problems using the EM algorithm, let us introduce a vector z = (z ′ 1 , . . . , z ′ n ) ′ , where z i is a vector or zeroone indicator variables such that z ik = 1 if w i arose from the k-th component of the mixture, and z ik = 0 otherwise. The complete data vector is thus x = (w ′ , z ′ ) ′ , and the complete-data p.d.f. is

g(x; Ψ) = n i=1 g(z i ; Ψ)g(w|z i ; Ψ) = n i=1 g k=1 π k z ik g k=1 g k (w; θ k ) z ik ,
from which we can deduce the expression of the complete-data log likelihood:

log L(Ψ; x) = n i=1 g k=1 z ik log π k + n i=1 g k=1 z ik log g k (w i ; θ k ) = g k=1 log π k n i=1 z ik - n 2 log(2π) -n log σ - 1 2σ 2 n i=1 g k=1 z ik (w i -m k ) 2 .
To perform the E-step of the EM algorithm, we need to compute the conditional expectation of the complete-data log likelihood conditionally on the observed data w, using the current fit Ψ (q) of the parameter vector:

Q(Ψ, Ψ (q) ) = E Ψ (q) [log L(Ψ; X) | w] = g k=1 log π k n i=1 E Ψ (q) (Z ik | w i ) - n 2 log(2π) -n log σ - 1 2σ 2 n i=1 g k=1 E Ψ (q) (Z ik W 2 i | w i ) -2m k E Ψ (q) (Z ik W i | w i ) + m 2 k E Ψ (q) (Z ik | w i ) . ( 16
)
We thus have to compute the three conditional expectations

E Ψ (q) (Z ik | w i ), E Ψ (q) (Z ik W 2 i | w i ) and E Ψ (q) (Z ik W i | w i ). Let t (q) ik = E Ψ (q) (Z ik | w i ). From Bayes' theorem, we have t (q) ik = P Ψ (q) (Z ik = 1| w i ) = P Ψ (q) ( w i |Z ik = 1)P Ψ (q) (Z ik = 1) P Ψ (q) ( w i ) = γ (q) ik π (q) k p (q) i , with γ (q) ik = µ w i (w)g k (w; θ (q) k )dw and p (q) i = g k=1 π (q) k µ w i (w)g k (w; θ (q) k )dx = g k=1 π (q) k γ (q)
ik .

Now,

E Ψ (q) (Z ik W 2 i | w i ) = E Ψ (q) (W 2 i | w i , Z ik = 1)P Ψ (q) (Z ik = 1| w i ) = ξ (q) ik t (q) ik , with ξ (q) ik = E θ (q) k (W 2 i |z ik = 1, w i ) = 1 γ (q) ik w 2 µ w i (w)g k (w; θ (q) k )dw.
Similarly,

E Ψ (q) (Z ik W i | w i ) = η (q) ik t (q) ik with η (q) ik = E θ (q) k (W i |Z ik = 1, w i ) = 1 γ (q) ik wµ w i (w)g k (w; θ (q) k )dw.
The expected complete-data log-likelihood can thus be written as:

Q(Ψ, Ψ (q) ) = g k=1 log π k n i=1 t (q) ik - n 2 log(2π) -n log σ - 1 2σ 2 n i=1 g k=1 t (q) ik ξ (q) ik -2m k η (q) ik + m 2 k . (17)
Let us now consider the M-step. Maximizing the first term on the right-hand side of [START_REF] Gil | The likelihood ratio test for goodness of fit with fuzzy experimental observations[END_REF] with respect to the mixing proportions π k , under the constraint k π k = 1 yields the classical solution

π (q+1) k = n i=1 t (q) ik g ℓ=1 n i=1 t (q) iℓ = 1 n n i=1 t (q) ik . Now, ∂Q(Ψ, Ψ (q) ) ∂m k = - 1 σ 2 - n i=1 t (q) ik η (q) ik + m k n i=1 t (q) ik .
Equating this derivative to zero and solving for m k yields:

m (q+1) k = n i=1 t (q) ik η (q) ik n i=1 t (q) ik .
Finally, the derivative of Q(Ψ, Ψ (q) ) with respect to σ is

∂Q(Ψ, Ψ (q) ) ∂σ = - n σ + 1 σ 3 n i=1 g k=1 t (q) ik ξ (q) ik -2m k η (q) ik + m 2 k .
Equating this derivative to zero and substituting m (q+1) k

for m k , we get the following fit for σ at iteration q + 1:

σ (q+1) = 1 n n i=1 g k=1 t (q) ik ξ (q) ik -2m (q+1) k η (q) ik + (m (q+1) k ) 2 .
Example 6 To illustrate the behavior of the above algorithm, we used again the perception dataset reported in [START_REF] González-Rodríguez | Fuzzy data treated as functional data: A one-way anova test approach[END_REF] and already described in Examples 2 and 4. We considered the assessments of 8 subjects for the relative lengths around 61.47, 50 and 38.3 (the relative lengths may vary slightly because of screen resolution, as explained in [START_REF] González-Rodríguez | Fuzzy data treated as functional data: A one-way anova test approach[END_REF]). The considered dataset was thus composed of n = 72 trapezoidal fuzzy sets, partitioned in three classes of equal size. A subset of the data is shown in Figure 9.

The above algorithm was applied to this dataset with g = 3. The means m k were initialized randomly using a uniform distribution in [0,100], while the initial standard deviation and proportions were set to fixed values as shown in Table 1. The EM algorithm was run 10 times, and the results corresponding to the best value of the objective criterion were retained. The MLE estimates are shown in Table 1, and the corresponding mixture density estimated is shown in Figure 9. We can see that the three relative length estimates are very close to their true values. To conclude this section, we can remark that determining the number of components in a mixture model is an important problem that has been addressed by many authors (see, e.g., [START_REF] Schwartz | Estimating the number of components in a finite mixture model[END_REF][START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF]). Some approaches are based on a Bayesian approach, which results in the addition of a regularization term to the log-likelihood criterion. Whether these or other approaches can be extended to the context of fuzzy data is an interesting question that goes beyond the scope of this paper and is left for further study.

Conclusions

We have shown that the EM algorithm can be adapted to handle estimation problems in which the observed data are fuzzy and are assumed to be related to underlying crisp data, making it possible to implement the maximum likelihood principle in this context. The proposed FEM method is very general and can be applied to a wide range of statistical problems. In this paper, it has been applied to classical parameter estimation tasks including multiple regression analysis with crisp input data and fuzzy output data, and univariate finite normal mixture estimation from fuzzy data. More complex problems such as principal component analysis [START_REF] Giordani | Two-and three-way component models for LR fuzzy data in a possibilistic framework[END_REF][START_REF] Denoeux | Principal component analysis of fuzzy data using autoassociative neural networks[END_REF], discriminant analysis, and clustering [START_REF] Côme | Learning from partially supervised data using mixture models and belief functions[END_REF][START_REF] Urso | A weighted fuzzy c-means clustering model for fuzzy data[END_REF][START_REF] Quost | Clustering fuzzy data using the fuzzy EM algorithm[END_REF] of multidimensional fuzzy data can be handled by this method as well. Another direction of research concerns the application of a similar methodology to more general types of imprecise and uncertain data such as described, e.g., by belief functions. Preliminary results in this direction have been reported in [START_REF] Denoeux | Maximum likelihood from evidential data: an extension of the EM algorithm[END_REF].

Two fuzzy events A and B in the probability space (R n , A, P ) are said to be independent if P ( A B) = P ( A)P ( B),

where A B is the fuzzy subset of R n with membership function

µ A B (x) = µ A (x) • µ B (x),
for all x in R n . The conditional probability of A given B is defined by

P ( A| B) = P ( A B) P ( B) , (19)
provided P ( B) > 0). In particular, assume that P is the probability distribution of a continuous random variable X with p.d.f. g(x). For a crisp subset A and a fuzzy subset B, Equation ( 19) becomes P (A| B) = µ A (x)µ B (x)g(x)dx µ B (u)g(u)du = A µ B (x)g(x) µ B (u)g(u)du dx.

The conditional density of X given B can thus be defined as:

g(x| B) = µ B (x)g(x) µ B (u)g(u)du . (20)

A.2 Trapezoidal fuzzy events and univariate normal distribution

To illustrate the above definitions, let us assume that P is the distribution of a univariate normal random variable X with mean m, standard deviation σ and p.d.f. g(x). Let x = (a, b, c, d) be a trapezoidal fuzzy number, with membership function

µ x (x) =            x-a b-a if a ≤ x ≤ b, 1 if b ≤ x ≤ c, d-x d-c if c ≤ x ≤ d, 0 otherwise.
Denoting by g(x) the p.d.f. of X and using [START_REF] Giordani | Two-and three-way component models for LR fuzzy data in a possibilistic framework[END_REF], the probability of x can be calculated as

P ( x) = E[µ x (X)] =
which makes it possible to complete the calculation of P ( x) using [START_REF] González-Rodríguez | Bootstrap techniques and fuzzy random variables: Synergy in hypothesis testing with fuzzy data[END_REF].

Let us now compute the expectation of X given x, using the expression (20) for the conditional density of X. We have E(X| x) = µ x (x) x g(x)dx P ( x) ,

where the denominator is given by [START_REF] González-Rodríguez | Bootstrap techniques and fuzzy random variables: Synergy in hypothesis testing with fuzzy data[END_REF]. The numerator is

Figure 1 :

 1 Figure 1: Ten trapezoidal fuzzy numbers considered in Example 2. The heights have been jittered to better separate the membership functions.

Figure 2 :

 2 Figure 2: Contour plot of the observed data log likelihood function and trajectory in parameter space (broken line) from the initial parameter values (m (0) , σ (0) ) (o) to the final MLE ( m, σ) (x), for the data of Example 2.

Figure 3 :

 3 Figure 3: Fuzzy information system used to encode the simulated fuzzy data of Example 3.

Figure 4 :

 4 Figure 4: Mean plus or minus one standard deviation (over 100 trials) of m for a fuzzy sample of size n from a standard normal distribution, as a function of n (Example 3).

Figure 5 :

 5 Figure 5: Mean plus or minus one standard deviation (over 100 trials) of σ for a fuzzy sample of size n from a standard normal distribution, as a function of n (Example 3).

Figure 6 :

 6 Figure 6: Data of Example 4, and fit of the regression line. The support and core of each fuzzy number are represented by solid and bold line segments, respectively. Jitter has been added to the u i to avoid superposition of fuzzy numbers.

Figure 7 :

 7 Figure 7: Example of randomly generated data for the experiment of Example 5, with the fit of the regression line. The points (u i , x i ) and (u i , x ′ i ) are represented by crosses and circles, respectively. The vertical segments represent the supports of the triangular fuzzy numbers x i .

Figure 8 :

 8 Figure 8: Box plots of the distributions of Ψ -Ψ over 1000 trials for the three methods and the experiment described in Example 5. Boxes whose notches do not overlap indicate that the medians of the two groups differ at the 5% significance level.

Figure 9 :

 9 Figure 9: Subset of the data of Example 6 (with 5 fuzzy numbers in each class), and normalized estimated mixture p.d.f. (bold line). The heights of fuzzy membership functions have been jittered for better visualization.

  b * ) -Φ(a * )) + Φ(c * ) -Φ(b * ) + d dc (Φ(d * ) -Φ(c * )) -1 dc d c xg(x)dx,[START_REF] González-Rodríguez | Bootstrap techniques and fuzzy random variables: Synergy in hypothesis testing with fuzzy data[END_REF] where Φ denotes the c.d.f. of the standard normal distribution, and x * denotes (xm)/σ for all x. It is easy to show that Φ(b * ) -Φ(a * )) ,

2 + (m 2 + 26 ) 2 +

 22262 (x)dx,[START_REF] Kruse | Statistics with vague data[END_REF] which can be computed using[START_REF] Hamdan | Mixture model clustering of uncertain data[END_REF] and b a x 2 g(x)dx = σ 2 √ 2π a * exp -a * 2 2 b * expσ 2 ) (Φ(b * ) -Φ(a * )) .(25)Finally, let us computeE(X 2 | x) = µ x (x) x 2 g(x)dx P ( x) . (The numerator isµ x (x)x 2 g(x)dx = m 3 (Φ(b * ) -Φ(a * )) .(28)

Table 1 :

 1 Initial parameter values and MLEs with g = 3 components for the data of

	Example 6.
	parameter	m 1	m 2	m 3	σ	π 1	π 2	π 3
	initial value 89.2054 12.7208 78.3867	10	0.3333 0.3333 0.3333
	MLE	60.9896 37.4851 49.9999 2.1915 0.2940 0.3283 0.3777

We thank an anonymous referee for bringing this point to our attention.

This argument was suggested to the author by an anonymous referee.

The dataset can be downloaded from the webpage of the SMIRE research group at http://bellman.ciencias.uniovi.es/SMIRE/Perceptionsdata.html.

To be consistent with our previous notations, we have to depart from the usual convention in regression analysis, where the dependent and independent variables are denoted by Y and x, respectively. We hope that the reader will not be confused by this change of notation.

Acknowledgement

The author thanks the members of SMIRE group and, in particular, Ana Colubi, María Ángeles Gil and Gil González-Rodríguez, for making the perception dataset publicly available. He also expresses his thanks to the two anonymous referees for their constructive and very helpful comments.

A Probability of fuzzy events

The notion of probability was extended to fuzzy events by Zadeh [START_REF] Zadeh | Probability measures of fuzzy events[END_REF]. In this appendix, we recall the main definitions and we derive some of the formula used in the paper.

A.1 Basic definitions

Let (R n , A, P ) be a probability space in which A is the σ-field of Borel sets and P is a probability measure on R n . Then a fuzzy event in R n is a fuzzy subset A of R n whose membership function µ A is Borel measurable. The probability of A is defined as the expectation of µ A with respect to P : P ( A) = µ A (x)dP. [START_REF] Giordani | Two-and three-way component models for LR fuzzy data in a possibilistic framework[END_REF]


Vincent Bombardier
email: bombardier@cran.uhp-nancy.fr

Laurent Wendling

Emmanuel Schmitt
email: schmitt@cran.uhp-nancy.fr

FEATURE SELECTION BASED ON IMPORTANCE AND INTERACTION INDEXES Hierarchical Fuzzy Rule Classifier Application

Keywords: Feature Selection, Pattern Recognition, Fuzzy Rules

This paper proposed an extension of an iterative method to select suitable features for pattern recognition context. The main improvement is to replace its iterative step with another criterion based on importance and interaction indexes, providing suitable feature reduced set. This new scheme is embedded on a hierarchical fuzzy rule classification system. At last, each node gathers a set of classes having a similar aspect. The aim of the proposed method is to automatically extract an efficient subset of suitable features for each node. A selection of features is given. The associated criterion is directly based on importance index and assessment of positive and negative interaction between features. An experimental study, made in a wood defect recognition industrial context, shows the proposed method is efficient to producing significantly fewer rules.

INTRODUCTION

In many pattern recognition applications, a feature selection scheme is fundamental to focus on most significant data while decreasing the dimensionality of the problem under consideration. The information to be extracted from the images is not always trivial, and to ensure that the maximum amount of information is obtained, the number of extracted features can strongly increase. The feature selection area of interest consists in reducing the problem dimension. It can be translated as being an optimization problem where a subset of features is searched in order to maximize the classification rate of the recognition system.

Because of specific industrial context, there are many constraints. One constraint is the necessity of working with very small training data sets (sometimes, there is only one or two samples for a specific class because of its rareness). Another difficulty is to respect the real time constraint in the industrial production system. So, low complexity must be kept for the recognition model. Such a classification problem has been relatively poorly investigated in the early years [START_REF] Abdulhady | Performance of neural classifiers for fabric faults classification[END_REF], [START_REF] Yang | Fabric defect classification using wavelet frames and minimum classification error training[END_REF], [START_REF] Murino | Statistical classification of raw textile defects[END_REF]. Thus, this work takes place on a "small scale" domain according to (Kudo 200), [START_REF] Zhang | Feature selection using Tabu Search method[END_REF] definition because of the weak number of used features. The Fuzzy Rule Iterative Feature Selection (FRIFS) method proposed in [START_REF] Schmitt | Improving Fuzzy Rule Classifier by Extracting Suitable Features from Capacities with Respect to the Choquet Integral[END_REF], is based on the analysis of a training data set in three steps. It combines an original Fuzzy Rule Classifier [START_REF] Bombardier | Fuzzy rule classifier: Capability for generalization in wood color recognition[END_REF] and feature selection associated to capacity learning has been proposed. This approach allows reducing the dimensionality of the problem while keeping a high recognition rate and improving the system interpretability by discarding weak parameters. First a reference set of features is set and the associated average recognition rate is kept to check the next training. Then an iterative global feature selection process is performed and can be roughly split into two steps:

Step 1. From the previous set of features, an interactivity process is applied to determine the less representative features.

Step 2. Generate the recognition model without the first less representative features and test it. The reached recognition rate is stored to be compared to the previous step one.

And so on while minimization is on the way. A rollback step provides the remaining set of suitable features. This model has shown its ability to efficiently detect fibre defects in industrial application [START_REF] Schmitt | Improving Fuzzy Rule Classifier by Extracting Suitable Features from Capacities with Respect to the Choquet Integral[END_REF].

Nonetheless in many industrial applications, it is interesting to group together classes or products following their specificities. A new selection scheme based on hierarchical description is proposed. New sets of features are assessed by a new criterion calculated from both importance and interaction indexes. In this paper, we propose to restrict the iterative part of FRIFS Method (Steps 1 and 2) by directly selecting a pertinent subset of features. This approach is applied to an industrial pattern recognition problem. The aim is to identify wood singularities. Comparisons with well-known selection methods like Sequential Backward Feature Selection (SBFS) or Sequential Forward Feature Selection (SBFS) methods, attest of the good behaviour of our method.

SELECTION OF SUITABLE SUBSET OF FEATURES

Hierarchical description

Feature selection based on Choquet integral [START_REF] Grabisch | Classification by fuzzy integral -performance and tests[END_REF] provided a suitable set of features to the Fuzzy Rules Classifier [START_REF] Bombardier | Fuzzy rule classifier: Capability for generalization in wood color recognition[END_REF].

As said previously, the rules are then obtained by learning. This principle is extended here to process with a hierarchical description associated to classes.

For each node, a set of suitable features is extracted from an analysis of their degree of importance combined with their level of positive and negative interactions (see next section). The selected subset directly depends on the recognition rates achieved by considering independently each node with the FRC and so on in order to process with the whole description.

Capacity indexes

Once the fuzzy measure is learned, it is possible to interpret the contribution of each decision criterion in the final decision. Several indexes can be extracted from the fuzzy measure, helping to analyze the behavior of DC [START_REF] Grabisch | The application of fuzzy integral in multicriteria decision making[END_REF]. The importance of each criterion is based on the definition proposed by Shapley in game theory [START_REF] Shapley | A value for n-person games. Contributions to the Theory of Games[END_REF]. Let a fuzzy measure μ and a criterion D i be considered:

( ) ( ) ( ) [ ] ∑ ∑ = ⊆ - = - ∪ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - = t T D N T i n t i i T D T t n n D | | \ 1 , 0 1 1 1 , μ μ μ σ (3)
The Shapley value can be interpreted as a weighted average value of the marginal contribution μ(T∪D i ) -μ(T) of criterion Di alone in all combinations.

The interaction index, also called the Murofushi and Soneda index [START_REF] Murofushi | Techniques for reading fuzzy measures(iii): interaction index[END_REF], [START_REF] Rendek | Extraction of Consistent Subsets of Descriptors using Choquet Integral[END_REF] represents the positive or negative degree of interaction between two Decision Criteria. If the fuzzy measure is non-additive then some sources interact. The marginal interaction between Di and Dj, conditioned to the presence of elements of combination T ⊆ X\D i D j is given by:

( ) ( ) ( ) ( )( ) T n t t n D D I j i j i D D N T D D j i ∑ ⊆ Δ - - - = \ ! 1 ! ! 2 , μ μ (4) With:
( ) ( ) (

) ( ) ( ) ( )

j i j i D D D T D T T D D T T j i - - - - + ∪ = Δ μ μ μ μ μ
And so on, considering any pair (D i ,D j ) with i ≠ j. Obviously the index are symmetric, i.e I(μ,D i D j )=I(μ,D j D i ). A positive interaction index for two DC D i and D j means that the second one reinforces the importance of one decision criterion.

In other words, both DC are complementary and their combined use betters the final decision. The behaviour is given by the value of the index. A negative interaction index indicates that the sources are antagonist.

Node features

The aim is to find a suitable subset of features for all the nodes. Each node is assumed to be a cluster of classes having similar characteristics. Considering the whole set of features it is obvious that the set of suitable parameters for one node could be not the same for another one. The processing time relies on the cardinality of the considered hierarchical description. A selection scheme is introduced to decrease processing complexity while quickly focusing suitable sets of features. The method relies on a combination of both indexes previously described.

A Shapley value property is Σ i=1,n σ(μ, D i )= 1. Generally values are multiplied by a number of decision criteria n=|N|. Hence, a DC with an importance index value less that 1 can be interpreted as a low impact in the final decision. Otherwise an importance index greater than 1 describes an attribute more important than the average that is

( ) { } 1 , / ≥ × ∈ = i i D n N D F μ σ
, the weaker:

{ } F N F \ =
Among them, some decision criteria may impact the recognition by having positive interactions with selected important criteria. Then, the DCs having an interaction of order 2 greater to the mean are selected:

( ) ( ) ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ Δ ≥ Δ ⊆ = ∑ ∑ ≠ ≠ ⊆ j i D v t F D D v t j i i j v t D D F D D F F D G / , 2 2 , | | 1 , | | 1 / μ μ At last final set G F N ∪ = '
is sent to the Fuzzy Rule Classifier to determine if it improves or not the recognition using the same learning samples. If the recognition rate is better than previous epoch one, the process is run again with N' instead of N. Otherwise the previous set of features is kept (rollback step). This process is run until all the nodes are studied.

EXPERIMENTAL STUDY

The choice of a fuzzy logic-based method for our application in the wood defect detection field could be justified by three main reasons. Firstly, the singularities to be recognized are intrinsically fuzzy (gradual transition between clear wood and knots for instance). The features extracted from the images are thus uncertain (but precisely calculated) and the use of fuzzy logic allows taking this into account.

Secondly, the customer expresses his needs under a nominal form; the output classes are thus subjective and often not separated (non strict boundary between the classes representing a small knot and the class representing a large knot). Finally, the customer needs and the human operator experience are subjective and mainly expressed in natural language.

Wood Defect Recognition Case

The results presented in this section are based on a real set of samples collected from a wood industrial case. The objective of this application is to develop a vision system for singularity identification on wooden boards used to estimate the quality of the final products.

During the image segmentation, a set of features is calculated on the achieved regions to provide a characteristic vector to the recognition step. This set is composed with geometrical (SURF, MIN_AXIS, MAJ_AXIS…) and topological (C1, C2, C3…) features but we cannot make them explicit because of confidentiality. We can note those calculated features are rather basic and simple due to the real time industrial constraints. So, associated values are redundant and often contradictory bringing noise to the final decision

Results and Discussion

The classification is done with the features issued from the segmentation stage performed by the industrialist. The database is composed of 877 samples divided in nine classes of wood singularities (called nuodo muerto, grieta, medula, resina …).

The learning database, used to compute the learning recognition rates, is composed of 250 samples. The generalisation rates are obtained with the generalisation database constituted by the 627 remaining samples. These databases are relatively heterogeneous, for instance, the 250 samples of the learning database are composed of 8, 56 7, 7, 18, 47, 5, 93, 9 samples of the nine classes. The fuzzy inference engine consists in a single inference where all the features are in input of the model and all the classes to recognize in output. Let us consider the hierarchical description provided in Figure 1. automatic feature selection methods used as references [START_REF] Pudil | Floating search methods in feature selection[END_REF]. The third one is an expert method where the feature sets have been manually determined.

Each selected feature set is provided for input to the FRC. This classifier is used in a hierarchical version rather than a single node version. Here, the structure is composed of three nodes as shown in figure 1. The main advantages are to reduce the total number of rules of the system, as the rule set of each node is smallest and also to enhance the recognition rate. Table 1 shows the learning rates (LR), the generalisation rates (GR), the number of features (#Param.) and the number of generated rules (#Rules) obtained for five methods.

Despite the proposed method provided better results considering separately each node, the recognition rates obtained by SFFS, SBFS and the proposed method (PM) are comparable. Scores between 97% and 98% are reached for the generalization step considering the whole database. The expert selection (EXP) gives rise to a rate around 96%. The difference is insignificant.

However, the number of parameters extracted by our method is the lowest assuming the better interpretability of the model. The numbers of rules are carried out using expert (EXP: 3775 rules) and the Proposed Method (PM: 1375 rules). Other methods lead to a rule number higher than 350 000.

CONCLUSIONS

An enhancement of a Fuzzy Rule Iterative Feature Selection method has been presented. It allows the decreasing of learning time processing while focusing on relevant samples. A suitable set of features is obtained. Then, FRIFS method is adapted to handle with a hierarchical fuzzy inference system. Industrial real-data tests show the efficiency of the proposed method. The recognition rate is similar to other methods but the number of features significantly decreases and thus the number of rules too. Actually, the extension of our model to provide selection of parameters per class is under consideration

Figure 1 :

 1 Figure 1: Example of hierarchical description.Tests aim to reduce the dimensionality by removing non-efficient features per node until reaching an interpretable model while keeping a "good" recognition rate. Three selection methods are used: SBFS, SFFS and EXP. Two of them are

Table 1 .

 1 Recognition rates and Generated rules.