François Coquet

Éric Lesage

Conditional inference with a complex sampling: exact computations and Monte Carlo estimations

Keywords: Auxiliary information, Conditional inference, Finite population, Inclusion probabilities, Monte Carlo methods, Sampling weights

In survey statistics, the usual technique for estimating a population total consists in summing appropriately weighted variable values for the units in the sample. Different weighting systems exit: sampling weights, GREG weights or calibration weights for example. In this article, we propose to use the inverse of conditional inclusion probabilities as weighting system. We study examples where an auxiliary information enables to perform an a posteriori stratification of the population. We show that, in these cases, exact computations of the conditional weights are possible. When the auxiliary information consists in the knowledge of a quantitative variable for all the units of the population, then we show that the conditional weights can be estimated via Monte-Carlo simulations. This method is applied to outlier and strata-Jumper adjustments.

Introduction

The purpose of this article is to give a systematic use of the auxiliary information at the estimation phase by the means of Monte Carlo methods, in a design based approach.

In survey sampling, we often face a situation where we use information about the population (auxiliary information) available only at the estimation phase. For example, this information can be provided by an administration file available only posterior to the collection stage. Another example would be the number of respondents to a survey. It is classical to deal with the non-response mechanism by a second sampling phase (often Poisson sampling conditional to the size of the sample). The size of the respondents sample is known only after the collection. This information can be compared to its counterpart estimated by the means of the sample. A significant difference typically reveals an unbalanced sample. In order to take this discrepancy into account, it is necessary to re-evaluate our estimations. In practice, two main technics exist: the model-assisted approach (ratio estimator, post-stratification estimator, regression estimator) and the calibration approach. The conditional approach we will develop in this article has been so far mainly a theoretical concept because it involves rather complex computations of the inclusion probabilities. The use of Monte-Carlo methods could be a novelty that would enable the use of conditional approach in practice. In particular, it seems to be very helpful for the treatment of outliers and strata jumpers.

Conditional inference in survey sampling means that, at the estimation phase, the sample selection is modelized by means of a conditional probability. Hence, expectation and variance of the estimators are computed according to this conditional sampling probability. Moreover, we are thus provided with conditional sampling weights with better properties than the original sampling weights, in the sense that they lead to a better balanced sample (or calibrated sample).

Conditional inference is not a new topic and several authors have studied the conditional expectation and variance of estimators, among them: [START_REF] Rao | Conditional inference in survey sampling[END_REF], [START_REF] Robinson | Conditioning ratio estimates under simple random sampling[END_REF], [START_REF] Tillé | Estimation in surveys using conditional inclusion probabilities: Simple random sampling[END_REF][START_REF] Tillé | Estimation in surveys using conditional inclusion probabilities: comlex design[END_REF] and Andersson (2004). Moreover, one can see that the problematic of conditional inference is close to inference in the context of rejective sampling design. The difference is that in rejective sampling, the conditioning event is controlled by the design, whereas, in conditional inference, the realization of the event is observed.

In section 2, the classical framework of finite population sampling and some notations are presented.

In section 3, we discuss the well-known setting of simple random sampling where we condition on the sizes of the sub-samples on strata (a posteriori stratification). This leads to an alternative estimator to the classical HT estimator. While a large part of the literature deals with the notion of correction of conditional bias, we will directly use the concept of conditional HT estimator [START_REF] Tillé | Estimation in surveys using conditional inclusion probabilities: Simple random sampling[END_REF], which seems more natural under conditional inference. A simulation study will be performed in order to compare the accuracy of the conditional strategy to the traditional one.

In section 4, the sampling design is a Poisson sampling conditional to sample size n (also called conditional Poisson sampling of size n). We use again the information about the sub-samples sizes to condition on. We show that the conditional probability corresponds exactly to a stratified conditional Poisson sampling and we give recursive formula that enables the calculation of the conditional inclusion probabilities. These results are new.

In section 5, we use a new conditioning statistic. Following [START_REF] Tillé | Estimation in surveys using conditional inclusion probabilities: Simple random sampling[END_REF][START_REF] Tillé | Estimation in surveys using conditional inclusion probabilities: comlex design[END_REF], we use the non-conditional HT estimation of the mean of the auxiliary variable to con-dition on. Whereas Tillé uses asymptotical arguments in order to approximate the conditional inclusion probabilities, we prefer to perform Monte Carlo simulations to address a non-asymptotic setting. Note that this idea of using independent replications of the sampling scheme in order to estimate inclusion probabilities when the sampling design is complex has been already proposed by [START_REF] Fattorini | Applying the Horvitz-Thompson criteion in complexe designs: A computer-intensive perspective for estimating inclusion probabilities[END_REF] and [START_REF] Thompson | Simulation-based Randomized Systematic PPS Sampling Under Substitution of Units[END_REF].

In section 6, we apply this method to practical examples: outlier and strata jumper in business survey. This new method to deal with outliers gives good results.

The context

Let U be a finite population of size N . The statistical units of the population are indexed by a label k ∈ {1, ..., N }. A random sample without replacement s is selected using a probability (sampling design) p(.). S is the set of the possible samples s. I [k∈s] is the indicator variable which is equal to one when the unit k is in the sample and 0 otherwise. The size of the sample is n(s) = |s|. Let B k = {s ∈ S, k ∈ s} = {s ∈ S, I [k∈s] = 1} be the set of samples that contain k. For a fixed individual k, let π k = p(B k ) be the inclusion probability and let d k = 1 π k be its sampling weight. For any variable z that takes the value z k on the U -unit k, the sum t z = k∈U z k is referred to as the total of z over U . t z,π = k∈s 1 π k z k is the Horvitz-Thompson estimator of the total t z .

Let x be an auxiliary variable that takes the value x k for the individual k. The x k are assumed to be known for all the units of U . Such auxiliary information is often used at the sampling stage in order to improve the sampling design. For example, if the auxiliary variable is a categorical variable then the sampling can be stratified. If the auxiliary variable is quantitative, looking for a balanced sampling on the total of x is a natural idea. These methods reduce the size of the initial set of admissible samples. In the second example, S balanced = {s ∈ S, t x,π = t x }. We wish to use auxiliary information after the sample selection, that is to take advantage of information such as the number of units sampled in each stratum or the estimation of the total t x given by the Horvitz-Thompson estimator. Let us take an example where the sample consists in 20 men and 80 women, drawn by a simple random sampling of size n = 100 among a total population of N = 200 with equal inclusion probabilities π k = 0.5. And let us assume that we are given a posteriori the additional information that the population has 100 men and 100 women. Then it is hard to maintain anymore that the inclusion probability for both men and women was actually 0.5. It seems more sensible to consider that the men sampled had indeed a inclusion probability of 0.2 and a weight of 5. Conditional inference aims at giving some theoretical support to such feelings.

We use the notation Φ(s) for the statistic that will be used in the conditioning. Φ(s) is a random vector that takes values in R q . In fact, Φ(s) will often be a discrete random vector which takes values in {1, ..., n} q . At each possible subset ϕ ⊂ Φ(S) corresponds an event A ϕ = Φ -1 (ϕ) = {s ∈ S, Φ(s) ∈ ϕ}. For example, if the auxiliary variable x k is the indicator function of a domain, say x k = 1 if the unit k is a man, then we can choose Φ(s) = k∈s I [k∈domain] = n domain the sample size in the domain (number of men in the sample). If the auxiliary variable x k is a quantitative variable, then we can choose Φ(s) = k∈s x k π k = t x,π the Horvitz-Thompson estimator of the total t x .

3 A posteriori Simple Random Sampling Stratification

Classical Inference

In this section, the sampling design is a simple random sampling without replacement(SRS) of fixed size n; S SRS = {s ∈ S, n(s) = n}; p(s) = 1/ N n and the inclusion probability of each individual k is π k = n/N . Let y be the variable of study. y takes the value y k for the individual k. The y k are observed for all the units of the sample. The Horvitz-Thompson (HT) estimator of the total t y = k∈U y k is ty,HT = k∈U y k

π k I [k∈s] . Assume now that the population U is split into H sub-populations U h called strata. Let N h = |U h |, h ∈ {1, .
.., H} be the auxiliary information to be taken into account. We split the sample s into H sub-samples s h defined by s h = s∩U h . Let n h (s) = |s h | be the size of the sub-sample s h .

Ideally, to use the auxiliary information at the sampling stage would be best. Here, a simple random stratified sampling (SRS stratified) with a proportional allocation N h n/N would be more efficient than a SRS. For such a SRS stratified, the set of admissible samples is

S SRSstratif ied = {s ∈ S, ∀h ∈ [1, H], n h (s) = N h n/N }, and the sampling design is p(s) = h∈[1,H] 1 ( N h n h ) , s ∈ S SRSstratif ied .
Once again, our point is precisely to consider setting where the auxiliary information becomes available posterior to this sampling stage

Conditional Inference

The a posteriori stratification with an initial SRS was described by [START_REF] Rao | Conditional inference in survey sampling[END_REF] and Tillé(1998). A sample s 0 of size n(s 0 ) = n is selected. We observe the sizes of the strata sub-samples:

n h (s 0 ) = k∈U h I [k∈s] , h ∈ [1, H].
We assume that ∀h, n h (s 0 ) > 0. We then consider the event:

A 0 = {s ∈ S, ∀h ∈ [1, H], n h (s) = n h (s 0 )}.
It is clear that s 0 ∈ A 0 , so A 0 is not empty.

We consider now the conditional probability: p A 0 (.) = p(./A 0 ) which will be used as far inference is concerned. The conditional inclusion probabilities are denoted

π A 0 k = p A 0 [I [k∈s] = 1] = E A 0 I [k∈s] = p [I [k∈s] = 1] ∩ A 0 /p(A 0 ).
Accordingly, we define the conditional sampling weights:

d A 0 k = 1 π A 0 k . Proposition 1.
1. The conditional probability p A 0 is the law of a stratified simple random sampling with allocation (n 1 (s 0 ), ..., n H (s 0 )),

For a unit k of the strata h: π

A 0 k = n h (s 0 ) N h and d A 0 k = N h n h (s 0 ) . Proof. |A 0 | = N 1 n 1 (s 0 ) × ... × N H n H (s 0 ) . ∀s ∈ A 0 , p A 0 (s) = 1/|A 0 |. So we have: p A 0 (s) = I [s∈A 0 ] 1 h∈[1,H] N h n h (s 0 ) = I [s∈A 0 ] * h∈[1,H] 1 N h n h (s 0 ) = h∈[1,H] I [n h (s)=n h (s 0 )] * 1 N h n h (s 0 )
and we recognize the probability law of a stratified simple random sampling with allocation (n 1 (s 0 ), ..., n H (s 0 )).

follows immediately.

Note that

E A 0 k∈U y k π k I [k∈s] = k∈U y k π k π A 0 k = h k∈U h y k N n h (s 0 ) nN h ,
so that the genuine HT estimator is conditionally biased in this framework. Even if, as Tillé(1998) mentioned, it is possible to correct this bias simply by retrieving it from the HT estimator, it seems more coherent to use another linear estimator constructed like the HT estimator but, this time, using the conditional inclusion probabilities.

Remark that in practice A 0 should not be too small. The idea is that for any unit k, we should be able to find a sample s such that s ∈ A 0 and k ∈ s. Thus, all the units of U have a positive conditional inclusion probability.

Definition 1. The conditional HT estimator is defined as:

ty,CHT = k∈U y k π A 0 k I [k∈s]
The conditional Horvitz-Thompson (CHT) etimator is obviously conditionally unbiased and, therefore, unconditionally unbiased.

This estimator is in fact the classical post-stratification estimator obtained from a model-assisted approach (see Särndal et Al.(1992) for example). However, conditional inference leads to a different derivation of the variance, which appears to be more reliable as we will see in next subsection.

Simulations

In this part, we will compare the punctual estimations of a total according to two strategies: (SRS design + conditional (post-stratification) estimator) and (SRS design + HT estimator). The population size is N = 500, the variable y is a quantitative variable drawn

,CHT = 1 N h k∈U h y k N h n h (s) I [k∈s] .
On Figure 1, we can see the values of μy,HT and μy,CHT for each of the 10 3 simulations. The red dots are those for which the conditional estimation is closer to the true value µ y = 2019.01 than the unconditional estimation; red dots represents 83.5% of the simulations. Moreover, the empirical variance of the conditional estimator is clearly smaller than the empirical variance of the unconditional estimator. This is completely coherent with the results obtained for the post-stratification estimator in an model-assisted approach (see Särndal et Al.(1992) for example). However, what is new and fundamental in the conditional approach, is to understand that for one fixed sample, the conditional bias and variance are much more reliable than the unconditional bias and variance. The theoretical study of the conditional variance estimation is a subject still to be developed.

Discussion

1. The traditional sampling strategy is defined as a couple (sampling design + estimator). We propose to define here the strategy as a triplet (sampling design + conditional sampling probability + estimator).

2. We have conditioned on the event:

A 0 = {s ∈ S, ∀h ∈ [1, H] n h (s) = n h (s 0 )}.
Under a SRS, it is similar to use the HT estimators of the sizes of the strata in the conditioning, that is to use Φ(s) = ( N1 (s), ..., NH (s)) t , where Nh (s) =

k∈U h I [k∈s] π k = N n n h (s). Then, A 0 = {s ∈ S, Φ(s) = Φ(s 0 )}.
We will see in Section 5 the importance of this remark.

3. The CHT estimations of the sizes of the strata are equal to the true strata sizes N h , which means that the CHT estimations, in this setting, have the calibration property for the auxiliary information of the size of the strata. Hence, conditional inference gives a theoretical framework for the current practice of calibration on auxiliary variables.

A Posteriori Conditional Poisson Stratification

Rao(1985), [START_REF] Tillé | Estimation in surveys using conditional inclusion probabilities: comlex design[END_REF] and Andersson (2005) mentioned that a posteriori stratification in a more complex setting than an an initial SRS is not a trivial task, and that one must rely on approximate procedures. In this section, we show that it is possible to determine the conditional sampling design and to compute exactly the conditional inclusion probabilities for an a posteriori stratification with a conditional Poisson sampling of size n.

Conditional Inference

Let p(s) = k∈s p k k∈s (1p k ) be a Poisson sampling with inclusion probabilities p = (p 1 , . . . , p N ) t , where p k ∈]0, 1] and s is the complement of s in U . Under a Poisson sampling, the units are selected independently.

By means of rejective technics, a conditional Poisson sampling of size n can be implemented from the Poisson sampling. Then, the sampling design is:

p(s) = K -1 1 |s|=n k∈s p k k∈s (1 -p k ),
where K = s,|s|=n k∈s p k k∈s (1p k ).

The inclusion probabilities π k = f k (U, p, n) may be computed by means of a recursive method:

f k (U, p, n) = p k 1 -p k n l∈U p l 1-p l (1 -f l (U, p, n -1)) (1 -f k (U, p, n -1))
where f k (U, p, 0) = 0. This fact was proven by [START_REF] Chen | Weighted finite poplation sampling to maximize entropy[END_REF] and one can also see [START_REF] Deville | Note sur l'algorithme de Chen[END_REF], [START_REF] Matei | Evaluation of variance approximations and estimators in unequal probability sampling with maximum entropy[END_REF], and Bondesson(2010). An alternative proof is given in Annex 1.

It is possible that the initial π k of the conditional Poisson sampling design are known instead of the p k 's. [START_REF] Chen | Weighted finite poplation sampling to maximize entropy[END_REF] have shown that it is possible to inverse the functions f k (U, p, n) by the means of an algorithm which is an application of the Newton method. One can see also [START_REF] Deville | Note sur l'algorithme de Chen[END_REF] who gave an enhanced algorithm.

Assume that a posteriori, thanks to some auxiliary information, the population is stratified in

H strata U h , h ∈ [1, H].
The size of the strata U h is known to be equal to N h , and the size of the sub-sample

s h into U h is n h (s 0 ) > 0. We consider the event A 0 = {s ∈ S, ∀h ∈ [1, H], n h (s) = n h (s 0 )}.
Proposition 2. With an initial conditional Poisson sampling of size n:

1. The probability conditional to the sub-samples sizes of the "a posteriori strata", p A 0 (s) = p(s/A 0 ), is the probability law of a stratified sampling with (independent) conditional Poisson sampling of size n h (s o ) in each stratum,

The conditional inclusion probability π

A 0 k of an element k of the strata U h is the inclusion probability of a conditional Poisson sampling of size n h (s o ) in a population of size N h . Proof. 1. For a conditional Poisson of fixe size n, a vector (p 1 , . . . , p N ) t exists, where p k ∈]0, 1], such that:

p(s) = K -1 1 |s|=n k∈s p k k∈s (1 -p k ), where K = s,|s|=n k∈s p k k∈s (1 -p k ). We remind that A 0 = {s ∈ S, ∀h ∈ [1, H], n h (s) = n h (s 0 )} Then: p(A 0 ) = K -1 p   h∈[1,H] [n h (s) = n h (s 0 )]   = K -1 h∈[1,H] p ([n h (s) = n h (s 0 )], )
where, p(.) is the law of the original Poisson sampling. Let s ∈ A 0 , then:

p A 0 (s) = p(s) p(A 0 ) = K -1 h=1,...,H k∈s h p k k∈s h (1 -p k ) K -1 h∈[1,H] p([n h (s) = n h (s 0 )]) = h=1,...,H k∈s h p k k∈s h (1 -p k ) p([n h (s) = n h (s 0 )]) = h=1,...,H k∈s h p k k∈s h (1 -p k ) s h ,|s h |=n h (s 0 ) k∈s-h p k k∈s h (1 -p k )
, which is the sampling design of a stratified sampling with independent conditional Poisson sampling of size n h (s o ) in each stratum. 2. follows immediately. We take the same population as in subsection 3.3. The sampling design is now a conditional Poisson sampling of size n = 100. The probabilities p k of the underlying Poisson design have been generated randomly, in order that k∈U p k = n and p k ∈ [0.13; 0.27]. K = 10 3 simulations were performed. Figure 2 shows that the punctual estimation of the mean of y is globally better for conditional inference. According to 77.3% of the simulations the conditional estimator is better than the unconditional estimator (red dots). The empirical variance as well is clearly better for the conditional estimator.

Definition 2. In the context of conditional inference on the sub-sample sizes of posteriori strata, under an initial conditional

Simulations

Discussion

This method allows to compute exact conditional inclusion probabilities in an "a posteriori stratification" under conditional Poisson of size n. However, one can figure out that this method can be used for any unequal probabilities sampling design, had the sampling frame been randomly sorted.

Conditioning on the Horwitz-Thompson estimator of an auxiliary variable

In the previous sections, we used the sub-sample sizes in the strata n h (s) to condition on. The good performances of this conditional approach result from the fact that the sizes of the sub-sample are important characteristics of the sample that are often used at the sampling stage. So, it was not surprising that the use of this information at the estimation stage would enhance the conditional estimators.

Another statistic that characterizes the representativeness of a sample is its HT estimator of the mean µ x (or total t x ) of an auxiliary variable. This statistic is used at the sampling stage in balanced sampling for example. So, as the sub-sample sizes into the strata, this statistic should produce good results in a conditional approach restraining the inference to the samples for which the HT estimation of µ x are equal to the value μ0 = μx,HT (s 0 ) of the selected sample s 0 .

In fact, we want the (conditional) set of the possible samples to be large enough in order that all conditional inclusion probabilities be different from zero. It is therefore convenient to consider the set of samples that give HT estimations not necessarily strictly equal to μ0 but close to μ0 . Let ϕ = [μ 0ε, μ0 + ε], for some ε > 0.

The set A ϕ of possible samples in our conditional approach will be:

A ϕ = {s ∈ S, μx,HT (s) ∈ [ µ 0 -ε, μ0 + ε]}.
The conditional inclusion probability of a unit k is:

π Aϕ k = p ([k ∈ s]/ [μ x,HT (s) ∈ [μ 0 -ε, μ0 + ε]]) = p ({s ∈ S, k ∈ s and μx,HT (s) ∈ [μ 0 -ε, μ0 + ε]}) p(A ϕ ) .
If μ0 = µ X then we are in a good configuration, because we are in a balanced sampling situation and the π Aϕ k

will certainly stay close to the π k . If µ 0 ≫ µ X say, then the sample s 0 is unbalanced, which means that in average, its units have a too large contribution x k /π k , either because they are too big (x k large) or too heavy (d k = 1 π k too large). In this case, the samples in A ϕ are also ill-balanced, because balanced on µ 0 instead of µ X : k∈s x k π k ≈ µ 0 . But conditioning on this information will improve the estimation. Indeed, the π Aϕ k will be different from the π k . For example, a unit k with a big contribution ( x k π k large) has more chance to be in a sample of A ϕ than a unit l with a small contribution. So, we can expect that π Aϕ k

> π k and π

Aϕ l < π l . And, in consequence, the conditional weight d ϕ k will be lower than d k and d ϕ l higher than d l , which will "balance" the samples of A ϕ .

Discussion:

• we can use different ways in order to define the subset ϕ. One way is to use the distribution function of Φ(s), denoted G(u) and to define ϕ as a symmetric interval:

ϕ = G -1 (max{G(Φ(s 0 )) - α 2 , 0}), G -1 (min{G(Φ(s 0 )) + α 2 , 1)} ,
where α = 5% for example.

Hence,

A ϕ = {s ∈ S, Φ(s) ∈ G -1 (max{G(Φ(s 0 )) - α 2 , 0}), G -1 (min{G(Φ(s 0 )) + α 2 , 1)} },
and p(A ϕ ) ≤ α.

As the cdf G(u) is unknown in general, one has to replace it by an estimated cdf of Φ(s), denoted ĜK (u), computed by means of simulations.

6 Generalization: Conditional Inference Based on Monte Carlo simulations.

In this section, we consider a general initial sample design p(s) with the inclusion probabilities π k . We condition on the event A ϕ = Φ -1 (ϕ) = {s ∈ S, Φ(s) ∈ ϕ}.

For example, we can use Φ(s) = k∈s

x k π k the unconditional HT estimator of t x and ϕ = [ϕ 1 , ϕ 2 ] an interval that contains Φ(s 0 ) = k∈s 0 x k π k , the HT estimation of t x with the selected sample s 0 . In other words, we will take into account the information that the HT estimator of the total of the auxiliary variable x lies in some region ϕ.

The mathematical expression of π

Aϕ k is straightforward:

π Aϕ k = p([k ∈ s]/A ϕ ) = s p(s)1 s∈Aϕ 1 [k∈s] p(A ϕ ) .
But effective computation of the π Aϕ k 's may be not trivial if the distribution of Φ is complex. Tillé(1998) used an asymptotical approach to solve this problem when

Φ(s) = k∈s x k π k 1 [k∈s]
; he has used normal approximations for the conditional and unconditional laws of Φ.

In the previous sections, we have given examples where we were able to compute the π Aϕ k 's (and actually the p Aϕ (s)'s) exactly. In this section, we give a general Monte Carlo method to compute the π Aϕ k .

Monte Carlo

We will use Monte Carlo simulations to estimate E(1 Aϕ 1 [k∈s] ) and E(1 Aϕ ). We repeat independently K times the sample selection with the sampling design p(s), thus obtaining a set of samples (s 1 , . . . , s K ). For each simulation i, we compute Φ(s i ) and I Aϕ (s i ). Then we compute N + 1 statistics:

M Aϕ = K i=1 1 Aϕ (s i ) ∀k ∈ U, M ϕ k = K i=1 1 Aϕ (s i )1 [k∈s i ]
We obtain a consistent estimator of π Aϕ k , as K → +∞:

πAϕ k = M ϕ k /K M Aϕ /K = M ϕ k M Aϕ (1)

Point and variance estimations in conditional inference

Definition 3. The Monte Carlo estimator of the total t y is the conditional Horvitz-Thompson estimator of t y after replacing the conditional inclusion probabilities by their Monte Carlo approximations:

ty,MC = k∈s 0 1 πAϕ k y k
The Monte Carlo estimator of the variance of ty,MC is:

V( ty,MC ) = k,l∈s 0 1 πAϕ k,l y k πAϕ k y l πAϕ l (π Aϕ k,l - πAϕ k πAϕ l ),
where

πAϕ k,l = K i=1 t y 1 Aϕ (s i )1 [k∈s i ] 1 [l∈s i ] K i=1 1 Aϕ (s i )
.

Fattorini(2006) established that ty,MC is asymptotically unbiased as M Aϕ → ∞, and that its mean squared error converges to the variance of ty,HT . [START_REF] Thompson | Simulation-based Randomized Systematic PPS Sampling Under Substitution of Units[END_REF] studied the rate of convergence of the estimators πAϕ k and of the estimator ty,MC following Chebychev's inequality. Using normal approximation instead of the Chebychev's inequality gives more precise confidence intervals. We have thus a new confidence interval for πAϕ k :

p |π Aϕ k -π Aϕ k | < F -1 ((1 -α)/2) 1 4M Aϕ ≤ α,
where F is the distribution function of the normal law N (0, 1). As for the relative bias, standard computation leads to:

p | ty,CHT -ty,CHT | ty,CHT ≤ ε ≥ 1 -4 × k∈s 1 -F ǫ 1 + ǫ M Aϕ π Aϕ k ≥ 1 -4n 1 √ 2π 1 + ǫ √ M Aϕ ǫ 2 π 0 • e -M Aϕ ǫ 2 (1+ǫ) 2 π 0 , (2)
where

π 0 = min{π Aϕ k , k ∈ U }. We used the inequality 1 -F (u) ≤ 1 √ 2π e -u 2 u
which is verified for large u.

The number K of simulations is set so that K i=1 I Aϕ (s i ) reaches a pre-established M Aϕ value. Because of our conditional framework, K is a stochastic variable which follows a negative binomial distribution and we have E(K) = M Aϕ p(A ϕ ) . For instance, if p(A ϕ ) = 0.05 = 5%, with M Aϕ = 10 6 , we expect E(K) = 2.10 7 simulations.

7 Conditional Inference Based on Monte Carlo Method in Order to Adjust for Outlier and Strata Jumper

We will apply the above ideas to two examples close to situations that can be found in establishments surveys: outlier and strata jumper.

We consider an establishments survey, performed in year "n+1 ", and addressing year "n". The auxiliary information x which is the turnover of the year "n" is not known at the sampling stage but is known at the estimation stage (this information may come from, say, the fiscal administration).

Outlier

In this section, the auxiliary variable x is simulated following a gaussian law, more precisely x k ∼ N(8 000, (2 000) 2 ) excepted for unit k = 1 for which we assume that

x 1 = 50 000. The unit k = 1 is an outlier. The variable of interest y is simulated by the linear model

y k = 1000 + 0.2 x k + u k ,
where u k ∼ N(0, (500) 2 ), u k is independent from x k . The outcomes are µ x = 8 531 and µ y = 2 695.

We assume that the sampling design of the establishments survey is a SRS of size n = 20 out of the population U of size N = 100 and that the selected sample s 0 contains the unit k = 1. For this example, we have repeated the sample selection until the unit 1 has been selected in s 0 . We obtain Φ(s 0 ) = μx,HT (s 0 ) = 9 970, which is 17% over the true value µ x = 8 531 and μy,HT (s 0 ) = 3 039 (recall that the true value of µ y is 2 695).

We set Φ and ϕ as in section 5 and we use Monte Carlo simulations in order to compute the conditional inclusion probabilities πAϕ k . Each simulation is a selection of a sample following a SRS of size n = 20from the fixed population U . Recall that the value of x k will eventually be known for any unit k ∈ U .

Actually, we use two sets of simulations. The first set is performed in order to estimate the cdf of the statistic Φ(s) = μx,HT (s) which will be used to condition on. This estimated cdf will enable us to construct the interval ϕ. A last remark concerns the distribution of the statistics Φ(s) = μx,HT (s). Figure 3 shows an unconditional distribution with 2 modes and far from gaussian. This shows that in presence of outlier, we can not use the method of Tillé (1999), which assumes a normal distribution for μx,HT (s).

Strata Jumper

In this section, the population U is divided into 2 sub-populations: the small firms and the large firms. Let us say that the size is appreciated thanks to the turnover of the firm. Official statistics have to be disseminated for this 2 different subpopulations. Hence, the survey statistician has to split the population into 2 strata

We choose the interval ϕ = [8 961, 10 342] by the means of the estimated cdf of Φ(s) = μx,2,HT (s) and so that p ([μ x,2,HT (s) ∈ [8 961, 10 342]]) = α = 5%. A ϕ is then the set of the possible samples in our conditional approach:

A ϕ = {s ∈ S, μx,HT (s) ∈ [8 961, 10 342]}.

All samples in A ϕ over-estimate the mean of x.

Among the 10 6 simulations, 49 778 simulated samples (4.98%) belongs to A ϕ . 55% of them contains the strata jumper, which gives the estimated conditional inclusion probability of the strata jumper πAϕ 1 = 0.55. It is not a surprise that the strata jumper is in one sample of A ϕ over two. Indeed, its initial sampling weight d 1 = 10 000 400 = 25 is high in comparison to the weights d k = 100 20 = 5 of the other selected units of the strata U x 2 , and its contribution 25x 1 N 2 contributes to over-estimate the mean of x .

The conditional inclusion probabilities for the other units of U x 2 are comparable to their initial π k = 0.2 (see Figure 4). Figure 4 shows that sampling weights correction is here a non-monotonic function of the variable x. We point out that the usual calibration method would not be able to perform this kind of weights correction because the calibration function used to correct the weights should be monotonic.

Similarly to the outlier setting, the unconditional distribution of the statistics Φ(s) = μx,2,HT (s) has 2 modes and is far from gaussian.

Conclusion

At the estimation stage, a new auxiliary information can reveal that the selected sample is imbalanced. We have shown that a conditional inference approach can take into account this information and leads to a more precise estimator than the unconditional Horvitz-Thompson estimator in the sense that the conditional estimator is unbiased (conditionally and unconditionally) and that the conditional variance is more rigorous in order to estimate the precision a posteriori.

In practise, we recommend to use Monte Carlo simulations in order to estimate the conditional inclusion probabilities. This technic seems particularly adapted to the treatment of outliers and stratajumpers.

A

Figure 1 :

 1 Figure 1: Punctual Estimation

.

  Poisson sampling of size n, the conditional HT estimator of the total t y is:The conditional variance can be estimated by means of one of the approximated variance formulae developed for the conditional Poisson sampling of size n. See for example[START_REF] Matei | Evaluation of variance approximations and estimators in unequal probability sampling with maximum entropy[END_REF], or Andersson(2004).

Figure 2 :

 2 Figure 2: Punctual Estimation

1 = 0 . 9997 . 1 N k∈s y k πAϕ 1 leads

 10999711 More precisely, we choose the interval ϕ = [9 793, 10 110] by the means of the estimated cdf of Φ(s) = μy,HT (s) and so that p ([μ x,HT (s) ∈ [9 793, μx,HT (s 0 )]]) = α 2 = 2.5% = p ([μ x,HT (s) ∈ [μ x,HT (s 0 ), 10 110]]). A ϕ is then the set of the possible samples in our conditional approach: A ϕ = {s ∈ S, μx,HT (s) ∈ [9 793, 10 110]}. Note that p ([μ x,HT (s) ∈ [9 793, 10 110]]) = α = 5%. A ϕ typically contains samples that over-estimate the mean of x.The second set of Monte Carlo simulations consists in K = 10 6 sample selections with a SRS of size n = 20 performed in order to estimate the conditional inclusion probabilities πAϕ k . 49 782 (4.98%) simulated samples fall in A ϕ , and among them, 49 767 samples contain the outlier, which correspond to the estimated conditional inclusion probability of the outlier:πAϕ It means that almost all the samples of A ϕ contain the outlier that is mainly responsible for the over-estimation because of its large value of the variable x! The weight of the unit 1 has changed a lot, it has decreased from d k = 1 0.2 = 5 to dAϕ k = 1.0003. The conditional sampling weights of the other units of s 0 are more comparable to their initial weights d k = 5 (see Figure7.1).The conditional MC estimator μy,MC (s) = to a much better estimation of µ y : μy,MC (s 0 ) = 2 671.

Figure 7 .

 7 Figure 7.1 gives an idea of the conditional inclusion probabilities for all the units of U . Moreover, this graph shows that the correction of the sampling weights dAϕ k d k = π k πAϕ k

Figure 3 :

 3 Figure 3: Outlier, Density of Φ(s) = μx,HT (s)

1 N k∈s y k πAϕ 1 leads to a better estimation of µ y, 2 :

 112 The conditional MC estimator μy,2,MC (s) = μy,2,MC (s 0 ) = 2 649.

Figure 4 :

 4 Figure 4: Strata Jumper, Sampling Weight Corrections

Annex 1 :III

 1 Inclusion Probability with Conditional Poisson Sampling Proof. The event l∈U,l =k I [l∈s] = n -1 is independent of the events [I [k∈s] = 0] and [I [k∈s] = 1] in the Poisson model. So we can write: p l∈U,l =kI [l∈s] = n -1 = p l∈U,l =k I [l∈s] = n -1 /[I [k∈s] [l∈s] = n -1 /[I [k∈s] = 1] =k I [l∈s] = n -1   /[I [k∈s] = 0] =k I [l∈s] = n -1   /[I [k∈s] = 0] [l∈s] = n -1 /[I [k∈s] = 0] = p l∈U I [l∈s] = n -1 p [I [k∈s] = 0]/ l∈U I [l∈s] = n -1 p [I [k∈s] = 0] = p l∈U I [l∈s] = n -1 (1f k (N, p, n -[l∈s] = n /[I [k∈s] = 1] = p l∈U I [l∈s] = n p [I [k∈s] = 1]/ l∈U I [l∈s] = n p [I [k∈s] = 1] = p l∈U I [l∈s] = n f k (N, p, n) p k .So we have:f k (U, p, n) = (1f k (U, p, n -1)) p k 1p k p l∈U I [l∈s] = n -1 p l∈U I [l∈s] = n = (1f k (U, p, n -1)) p k 1p k h(U, p, n)And we can use the property k∈U f k (U, p, n) = k∈U π k = n to compute h(U, p, n) and conclude.

corresponding to the sub-populations. This may not be an easy job because the size of firms can evolve from one year to another. Here we assume that, at the time when the sample is selected, the statistician does not know yet the auxiliary information x of the turnover of the firm for the year "n", more precisely the strata the firm belongs to for the year "n". Let us assume that he only knows this information for the previous year,"n-1 ". This information is denoted by z. In practice, small firms are very numerous and the sampling rate for this strata is chosen low. On the contrary, large firms are less numerous and their sampling rate is high. When a unit is selected among the small firms but eventually happens to be a large unit of year "n", we call it a strata jumper. At the estimation stage, when the information x becomes available, this unit will obviously be transferred to strata 2 . This will bring a problem, not due to its y-value (which may well be typical in strata 2) but to its sampling weight, computed according to strata 1 (the small firms), and which will appear to be very large in comparison to the other units in strata 2 at the estimation stage.

In our simulations, the population U is split in 2 strata, by means of the auxiliary variable z: U z 1 , of size N z 1 = 10 000, is the strata of presumed small firms and U z 2 , of size N z 2 = 100, the strata of presumed large firms. The auxiliary variable x, which is the turnover of the year "n" known after collection, is simulated under a gaussian law N(8 000, (2 000) 2 ) for the units of the strata U z 2 and for one selected unit of the strata U z 1 . Let us say that this unit, the strata jumper, is unit 1. Our simulation gives x 1 = 8 002. The variable of interest y is simulated by the linear model y k = 1000 + 0.2 x k + u k , where u k ∼ N(0, (500) 2 ), u k and x k independent. We do not simulate the value of x and y for the other units of the strata U z 1 because we will focus on the estimation of the mean of y for the sub-population of large firms of year n U x 2 : µ y,2 = 1 N 2 k∈U x 2 y k . We find µ x,2 = 8 138 and µ y,2 is 2 606.

The sampling design of the establishments survey is a stratified SRS of size n 1 = 400 in U z 1 and n 2 = 20 in U z 2 . We assume that the selected sample s 0 contains the unit k = 1. In practice, we repeat the sample selection until the unit 1 (the strata jumper) has been selected.

As previously, Φ and ϕ are defined as in Section 5. We use Monte Carlo simulations in order to compute the conditional inclusion probabilities πAϕ k . A simulation is a selection of a sample with stratified SRS of size n 1 = 400 in U z 1 and n 2 = 20 in U z 2 .

We choose the statistic Φ(s) = μx,2,HT (s) in order to condition on. K = 10 6 simulations are performed in order to estimate the cdf of Φ(s) and the conditional inclusion probabilities.

Our simulations give Φ(s 0 ) = μx,2,HT (s 0 ) = 9 510, which is far from the true value µ x,2 = 8 138 and μy,2,HT (s 0 ) = 3 357 (recall that the true value of µ y,2 is 2 606).


Guillaume Blin
email: gblin@univ-mlv.fr

Romeo Rizzi
email: rrizzi@dimi.uniud.it

Stéphane Vialette
email: vialette@univ-mlv.fr

A faster algorithm for finding minimum Tucker submatrices

Keywords:

A binary matrix has the Consecutive Ones Property (C1P) if its columns can be ordered in such a way that all 1s on each row are consecutive. Algorithmic issues of the C1P are central in computational molecular biology, in particular for physical mapping and ancestral genome reconstruction. In 1972, Tucker gave a characterization of matrices that have the C1P by a set of forbidden submatrices, and a substantial amount of research has been devoted to the problem of efficiently finding such a minimum size forbidden submatrix. This paper presents a new O(∆ 3 m 2 (m∆ + n 3 )) time algorithm for this particular task for a m×n binary matrix with at most ∆ 1-entries per row, thereby improving the O(∆ 3 m 2 (mn + n 3 )) time algorithm of [M. Dom, J. Guo and R. Niedermeier, Approximation and fixed-parameter algorithms for consecutive ones submatrix problems, Journal of Computer and System Sciences, 76(3-4): [204][205][206][207][208][209][210][211][212][213][214][215][216][217][218][219][220][221] 2010 ]. Moreover, this approach can be used -with a much heavier machinery -to address harder problems related to Minimal

Introduction

A binary matrix has the Consecutive Ones Property (C1P) if its columns can be ordered in such a way that all 1s on each row are consecutive. Both deciding if a given binary matrix has the C1P and finding the corresponding columns permutation can be done in linear time [9, 17, 18, 22-24, 27, 30]. The C1P of matrices has a long history and it plays an important role in combinatorial optimization, including application fields such as scheduling [START_REF] Bartholdi | Cyclic scheduling via integer programs with circular ones[END_REF][START_REF] Hassin | Approximation algorithms for hitting objects with straight lines[END_REF][START_REF] Hochbaum | Cyclical scheduling and multi-shift scheduling: Complexity and approximation algorithms[END_REF][START_REF] Veinott | Optimal capacity scheduling[END_REF], information retrieval [START_REF] Kou | Polynomial complete consecutive information retrieval problems[END_REF], and railway optimization [START_REF] Mecke | Station location complexity and approximation[END_REF][START_REF] Mecke | Solving geometric covering problems by data reduction[END_REF][START_REF] Ruf | Set covering with almost consecutive ones property[END_REF] (see [START_REF] Dom | Algorithmic aspects of the consecutive-ones property[END_REF] for a recent survey). Furthermore, algorithmic aspects of the C1P turn out to be of particular importance for physical mapping [START_REF] Alizadeh | Physical mapping of chromosomes using unique probes[END_REF][START_REF] Christof | A branch-and-cut approach to physical mapping of chromosome by unique end-probes[END_REF][START_REF] Lu | A test for the consecutive ones property on noisy dataapplication to physical mapping and sequence assembly[END_REF] and ancestral genome reconstruction [START_REF] Adam | Common intervals and symmetric difference in a model-free phylogenomics, with an application to streptophyte evolution[END_REF][START_REF] Chauve | A methodological framework for the reconstruction of contiguous regions of ancestral genomes and its application to mammalian genome[END_REF]. (see also [10, 3-5, 13, 31] for other applications in computational molecular biology). Actually, our main motivation for studying algorithmic aspects of the C1P comes from minimal conflicting sets in binary matrices in the context of ancestral genome reconstruction [?]. A minimal conflicting set of rows in a binary matrix is a set of rows R that does not have the C1P but such that any proper subset of R has the C1P (a similar definition applies for columns). Tucker [START_REF] Tucker | A structure theorem for the consecutive 1s property[END_REF] has characterized the binary matrices that have the C1P by a set of forbidden submatrices, and the aim of this paper is to lay the foundations for efficiently computing minimal conflicting sets by presenting a new efficient algorithm for finding such minimum size forbidden Tucker submatrices [START_REF] Blin | General framework for minimal conflicting set[END_REF].

Recently, Dom et al. [START_REF] Dom | Approximation and fixed-parameter algorithms for consecutive ones submatrix problems[END_REF] investigated some natural problems arising when a matrix M does not have the C1P property (the C1P is indeed a desirable property than often leads to efficient algorithms): find a minimum-cardinality set of columns to delete such that the resulting matrix has the C1P, find a minimum-cardinality set of rows to delete such that the resulting matrix has the C1P, and find a minimum-cardinality set of 1-entries in the matrix that shall be flipped (that is, replaced by 0-entries) such that the resulting matrix has the C1P. All these problems are NP-hard even for simple instances [START_REF] Hajiaghayi | A note on the consecutive ones submatrix problem[END_REF][START_REF] Tan | The consecutive ones submatrix problem for sparse matrices[END_REF], and hence Dom et al. have focused on approximation and parameterized complexity issues. To this end, they have provided a technical solution based on efficiently detecting forbidden Tucker submatrices [START_REF] Tucker | A structure theorem for the consecutive 1s property[END_REF].

In this paper, we presents a new O(∆ 3 m 2 (m∆ + n 3 )) time algorithm for finding a a minimum size forbidden Tucker submatrix in m × n binary matrices with at most ∆ 1-entries per row, thereby improving the O(∆ 3 m 2 (mn + n 3 )) time algorithm of Dom et al. [START_REF] Dom | Approximation and fixed-parameter algorithms for consecutive ones submatrix problems[END_REF][START_REF] Dom | Algorithmic aspects of the consecutive-ones property[END_REF]. This paper is organized as follows: we first recall basic definitions in Section 2, and we then formally introduce the considered problem. In Section 3, we briefly recall the algorithm of Dom et a. Section 4 is devoted to presenting out technical improvement and we consider in Section 5 matrices with unbounded ∆.

Preliminaries

All graphs are considered as undirected. Given a graph G = (V, E), let N (v) = {u|(u, v) ∈ E} denote the neighborhood of vertex v in G. The neighborhood described above does not include v itself, and is more specifically the open neighborhood of v; it is also possible to define a neighborhood in which v itself is included, called the closed neighborhood. For any subset V ⊆ V of vertices, let G[V ] denote the subgraph of G in-duced by the vertices of V (i.e., G = (V , {(u, v) ∈ E|u, v ∈ V })). An induced cycle is a cycle that is an induced subgraph of G; induced cycles are also called chordless cycles. An asteroidal triple is an independent set of three vertices such that each pair is joined by a path that avoids the neighborhood of the third.

A binary matrix has the Consecutive Ones Property (C1P) if its columns can be ordered in such a way that all 1s on each rows are consecutive. For a matrix M , we let r i and c j stand for the i th row and the j th column of M , respectively. Let M be a m × n binary matrix. Its corresponding bipartite graph G(M ) = (V M , E M ) is defined as follows:

V M = R ∪ C, where R = {r i : 1 ≤ i ≤ m} is the set of rows of M and C = {c i : 1 ≤ i ≤ n}
is the set of columns of M , and two vertices r i ∈ R and c j ∈ C are connected by an edge if and only if M [i, j] = 1. Equivalently, M is the reduced adjacency matrix of G(M ) (i.e. the nonredundant portion of the full adjacency matrix for the bipartite graph). Actually, it will be convenient to define G(M ) as a vertex-colored bipartite graph: each r i ∈ R (row vertex) is colored black and each c j ∈ c (column vertex) is colored white. See Figure 1 for an illustration (we use capital letters for black vertices and uncapitalized letters for white vertices).

M =     c1 c2 c3 c4 c5 c6 r1 0 1 0 1 1 0 r2 0 0 0 1 1 0 r3 0 1 1 0 1 0 r4 1 0 1 0 1 1     c2 GM R1 R3 c4 c5 c3 R2 R4 c1 c6
Fig. 1. A binary matrix and its corresponding vertex-colored bipartite graph.

The following result of Tucker links the C1P for binary matrices to asteroidal triples.

Theorem 1 ([34], Theorem 6).

A binary matrix has the C1P if and only if its corresponding vertex-colored bipartite graph does not contain a white asteroidal triple, i.e. an asteroidal triple on column vertices.

Moreover, Tucker has characterized the binary matrices that have the C1P by a set of forbidden submatrices (the number of vertices in the forbidden graph defines the size of the Tucker configuration). 3 The algorithm of Dom et al.

Theorem 2 ([34], Theorem 9). A binary matrix has the C1P if and only if it contains none of the matrices

M I k , M II k , M III k (k ≥ 1), M IV and M V depicted in Figure 2.
In [START_REF] Dom | Approximation and fixed-parameter algorithms for consecutive ones submatrix problems[END_REF], Dom et al. provided an algorithm for finding a forbidden Tucker submatrix (i.e., one of

T = {M I k , M II k , M III k , M IV , M V })
in a given binary matrix. The general algorithm is as follows. For each white asteroidal triple u, v, w of G(M ), compute the sum of the lengths of three shortest paths connecting two by two u, v and w (each path has to avoid the closed neighborhood of the third vertex). Select an asteroidal triple u, v, w of G(M ) with minimum total length of the paths connecting each

O(∆mn 2 + n 3 ) O(∆ 3 m 2 (n + ∆m)) MIII k O(∆ 3 m 3 n + ∆ 2 m 2 n 2 ) O(∆mn 2 (n + ∆m)) MIV O(∆ 3 m 2 n 3 ) Algorithm of Dom et al. used verbatim MV O(∆ 4 m 2 n) Algorithm of Dom et al. used verbatim Overall O(∆ 3 m 2 (mn + n 3 )) O(∆ 3 m 2 (∆m + n 3 ))
Table 1. Comparing our results with Dom et al. [START_REF] Dom | Approximation and fixed-parameter algorithms for consecutive ones submatrix problems[END_REF][START_REF] Dom | Algorithmic aspects of the consecutive-ones property[END_REF].

pair of vertices and return the rows and columns of M that correspond to the vertices that occur along the three shortest paths. The authors proved that the returned submatrix does contain a forbidden Tucker submatrix of T but which is not necessarily of minimum size (for M III k , M IV and M V ). Indeed, since the three shortest paths may share some vertices, the sum of the lengths of the three paths is not necessarily the number of vertices in the union of the three paths. However, 1.

The main contribution of this paper is a simple O(∆ 3 m 2 (m∆ + n 3 )) time algorithm for finding a minimum size forbidden Tucker submatrix. Our technical improvement on Dom et al. [START_REF] Dom | Approximation and fixed-parameter algorithms for consecutive ones submatrix problems[END_REF][START_REF] Dom | Algorithmic aspects of the consecutive-ones property[END_REF] is based on shortest paths and two graph pruning techniques: clean and complement clean (to be defined in the next section). Graph pruning techniques were introduced by Conforti and Rao [START_REF] Conforti | Structural properties and decomposition of linear balanced matrices[END_REF]. One has to note that graph pruning technique does not always succeed in the detection of induced configurations. Indeed, in [START_REF] Bienstock | On the complexity of testing for odd holes and induced odd paths[END_REF] Bienstock gave negative results among which one can find an NP-completeness proof for the problem of deciding whether a graph contains an odd chordless cycle containing a given vertex. This negative result, which in attacking the perfect graph conjecture was useful in posing limits in what could have been a reasonable approach, also demonstrates that not everything can be done with the detection of induced configurations.

Fast detection of minimum size forbidden Tucker submatrices

In this section, we design a general algorithm that reports in polynomialtime the smallest Tucker configuration of a given matrix M if it exists. To do so, we seek for any subgraph of G(M ) that corresponds to a submatrix of T .

Let us introduce the clean and cpl clean cleaning operations. Let M be a binary matrix and

G(M ) = (V M , E M ) be the corresponding vertex- colored bipartite graph. Let v be a vertex of G(M ), then clean G(M ) (v) is the graph obtained by removing from G(M ) all neighbors of v, i.e., G(M )[V M \ N (v)]. Let v be a vertex of G(M ), then cpl clean G(M ) (v) (read complement clean))
is the graph obtained by removing from G(M ) all vertices that do not belong to the same partition nor the neighborhood of v, i.e., G(M

)[N (v) ∪ {u : color(u) = color(v)}].
For the sake of simplicity, we shall write clean G(M ) (u 1 , u 2 , . . . , u k ) for the sequence

G 1 ← clean G(M ) (u 1 ), G 2 ← clean G 1 (u 2 ), . . . G k ← clean G k-1 (u k ), Return G k ,
and cpl clean G(M ) (u 1 , u 2 , . . . , u k ) for the sequence

G 1 ← cpl clean G(M ) (u 1 ), G 2 ← cpl clean G 1 (u 2 ), . . . , G k ← cpl clean G k-1 (u k ), Return G k .
We now focus on the bipartite graphs that represent Tucker configurations (see Figure 2) and define our guessing functions. Define the function guess I (G(M )) (reps. guess II (G(M ))) that returns all possible sets of vertices {x, y, z, A, B} ⊆ V M as defined for the Tucker configuration M I k (resp. M II k ) in Figure 2. Furthermore, Define the function guess III (G(M )) that returns all possible sets of vertices {x, y, z, A} ⊆ V M as defined for the Tucker configuration M III k in Figure 2. Of particular importance, in the presented algorithms, guessed vertices will never be affected (i.e., deleted) by the clean and cpl clean functions. Proof. We apply Algorithm 1 to G(M ). Let us first prove that if G(M I k ) occurs in G(M ), then Algorithm 1 finds it. Suppose G = G(M I k ) occurs in G(M ). Then among all the guessed 5-plets x, y, z, A, B (Line 1), there should be at least one guess such that x, y, z, A, B are part of the vertices of G . By definition, G is a chordless cycle. Therefore, clean(x, A, B) preserves G since in G , [START_REF] Adam | Common intervals and symmetric difference in a model-free phylogenomics, with an application to streptophyte evolution[END_REF] x is only connected to vertices A and B, (2) A is only connected to vertices x and y, and (3) B is only connected to x and z. Therefore, looking for a shortest path p in the pruned graph between y and z after having removed A and B ensures the minimality of the returned solution.

The guessing can be done in O(m 2 ∆ 3 ) time. Indeed, once A has been identified, one can (i) select x and y among the at most ∆ neighbors of A, and (ii) identify B and one of its at most ∆ neighbors as z such that x ∈ N (B) and z / ∈ {x, y}. For each such guessing, the cleaning of x, A, B

Algorithm 1 Finding G(M I k ) in G(M )
1: for all {x, y, z, A, B} ← guess I (G(M )) do 2:

G1 ← clean G(M ) (x,A,B) 3:
Remove the vertices A and B from G1. In the resulting graph, find a shortest path p from y to z 4:

if p exists then 5:

S ← {x, y, z, A, B} ∪ {u : u ∈ p} 6:

return the induced subgraph G(M )[S] 7:

end if 8: end for 9: return None can be done in O(∆ + m) time. Finally, a shortest path between y and z can be computed in O(n+∆m) time (the pruned graph has at most m+n vertices and ∆m edges). On the whole, Algorithm 1 is O(m 2 ∆ 3 (n + ∆m)) time.

Proof. We apply Algorithm 2 to G(M ). Let us first prove that if G(M II

k ) occurs in G(M ), then Algorithm 2 finds it. Suppose G = G(M II k ) occurs in G(M ).
Then among all the guessed 5-plets x, y, z, A, B in Line 1, there must be at least one guess such that x, y, z, A, B are part of the vertices of G . By definition, in G any unguessed white node is in the neighborhood of both A and B. Thus, cpl clean(A, B) preserves G since (1) y (the only white node not in the neighborhood of B) has been guessed, and (2) z (the only white node not in the neighborhood of A) has been guessed. Moreover, in G x should be only connected to A and B. Thus, clean(x) preserves G . Finally, looking for a shortest path p in the pruned graph between y and z after having removed A and B ensures the minimality of the returned solution.

Algorithm 2 Finding G(M II k ) in G(M ) 1: for all {x, y, z, A, B} ← guess II (G(M )) do 2: G1 ← cpl clean G(M ) (A, B) 3: G2 ← cleanG 1 (x) 4:
Remove the vertices A and B from G2. In the resulting graph, find a shortest path p from y to z 5:

if p exists then 6:

S ← {x, y, z, A, B} ∪ {u : u ∈ p} 7:

return the induced subgraph G(M )[S] 8:

end if 9: end for 10: return None

The guessing can be done in O(m 2 ∆ 3 ) time. For each guessing, the cleaning/complement-cleaning of x, A, and B can be done in O(n + m) time. Finally, a shortest path between y and z can be computed in O(n + ∆m) time (the pruned graph has at most ∆ + n vertices and ∆m edges). On the whole, Algorithm 2 is O(m 2 ∆ 3 (n + ∆m)) time.

(M III k ) in G(M ) in O(m∆n 2 (n + ∆m)) time (if such a submatrix exists). Algorithm 3 Finding G(M III k ) in G(M ) Proof. 1: for all {x, y, z, A} ← guess III (G(M )) do 2: G1 ← cpl clean G(M ) (A) 3: G2 ← cleanG 1 (x) 4:
Remove the vertex A from G2. In the resulting graph, find a shortest path p from y to z.

5:

if p exists then 6:

S ← {x, y, z, A} ∪ {u : u ∈ p} 7:

return the induced subgraph G(M )[S] 8:

end if 9: end for 10: return None

We apply Algorithm 3 to G(M ). Let us first prove that if G(M III k ) occurs in G(M ), then Algorithm 3 finds it. Suppose G = G(M III k ) occurs in G(M ). Then among all the guessed 4-plets x, y, z, A in Line 1, there must be at least one guess such that x, y, z, A are part of the vertices of G . By definition, in G any unguessed white node is in the neighborhood of A. Thus, cpl clean(A) preserves G since in G y and z (the only white nodes not in the neighborhood of A) have been guessed. Moreover, in G x is only connected to A. Therefore, clean(x) preserves G . Finally, looking for a shortest path p in the pruned graph between y and z after having removed A ensures the minimality of the returned solution.

The guessing can be done in O(m∆n 2 ) time. Indeed, once A has been identified, one can (i) select x among the at most ∆ neighbors of A and (ii) then identify y and z. For each such guessing, the cleaning/complementcleaning of x and A can be done in O(n+m) time. Finally, a shortest path between y and z can be computed in O(n + ∆m) time (the pruned graph has at most ∆ + n vertices and ∆m edges). On the whole, Algorithm 3 is O(m∆n 2 (n + ∆m)) time.

As for G(M IV ) and G(M V ), a simple brute-force search yields the time complexity detailed in Lemma 4.

Lemma 4 ([16], Proposition 5.3). Let M be a m × n binary matrix with at most ∆ 1-entries per row. The smallest

G(M IV ) (resp. G(M V )) in G(M ) can be computed in O(∆ 3 m 2 n 3 ) (resp. O(∆ 4 m 2 n) time) if it exists.
We are now ready to state the main result of this paper (Table 1 compares our results with Dom et al. [START_REF] Dom | Approximation and fixed-parameter algorithms for consecutive ones submatrix problems[END_REF].). Notice that our results do not improve on the ones of Dom et al. [START_REF] Dom | Approximation and fixed-parameter algorithms for consecutive ones submatrix problems[END_REF] for large ∆. Furthermore, our results also do not improve on the ones of Dom et al.

Matrices with unbounded ∆

As mentioned in [START_REF] Dom | Approximation and fixed-parameter algorithms for consecutive ones submatrix problems[END_REF], a natural question is concerned with the complexity of the problem when the number of 1s per row is unbounded. One can distinguish two subcases: the maximum number of 1s per column is bounded (say by C) or not. In the following, we prove that using a similar approach to the one used in the preceding section, one can achieve an O(C 2 n 3 (m + C 2 n)) (resp. O(n 4 m 4 )) time complexity for the bounded (resp. unbounded) case. For Algorithm 1, the guessing can be done in O(n 3 C 2 ) time. Indeed, once x has been identified, one can select A and B among the at most C neighbors of x and then identify y and z such that y ∈ N (A), y ∈ N (B), z ∈ N (B), z ∈ N (A) and x = y = z. For each such guessing, the cleaning of x, A, and B can be done in O(C + n) time. Finally, a shortest path between y and z can be computed in O(m + Cn) time (the pruned graph has at most m + n vertices and Cn edges). We conclude that Algorithm 1 is O(C 2 n 3 (m + Cn)) time.

For Algorithm 2, the guessing can also be done in O(n 3 C 2 ). For each such guessing, the cleaning/complement-cleaning of x, A and B can be done in O(C + n) time. Finally, a shortest path between y and z can be computed in O(m + Cn) time (the pruned graphs has at most m + n vertices and Cn edges). We conclude that Algorithm 2 is O(C 2 n 3 (m + Cn)) time.

As for Algorithm 3, the guessing can be done in O(n 3 C) time. Indeed, once x has been identified, one can select A among the at most C neighbors of x and then identify y and z such that y ∈ N (A), z ∈ N (A) and x = y = z. For each such guessing, the cleaning/complement-cleaning of x and A can be done in O(C + n) time. Finally, a shortest path between y and z can be computed in O(m + Cn) time (the pruned graph has at most m + n vertices and Cn edges). We conclude that Algorithm 3 is O(Cn 3 (m + Cn)) time.

M IV and M V forbidden submatrices deserve careful consideration. Indeed, a direct application of Proposition 5.3 [START_REF] Dom | Approximation and fixed-parameter algorithms for consecutive ones submatrix problems[END_REF] for finding M IV and M V forbidden submatrices results in a O(n 6 ) time algorithm. We use the following strategy: (i) Select the central three white nodes u, v, and w (distinct from x, y and z) among the O(n 3 ) such triples, (ii) select four black nodes that are neighbors of at least one of u, v, or w among the O(C 4 ) such 4-plets, and (iii) for each such combination check in O(n) time whether every columns of the matrix M IV appears at least once in the submatrix induced by the selected rows.

A submatrix of the type M V can be found analogously: (i) select the two central white nodes (distinct from x, y, and z) among the O(n 2 ) such triples, (ii) select four black nodes that are neighbors of at least one of u and v among the O(C 4 ) such 4-plets, and (iii) for each such combination check in O(n) time whether every columns of the matrix M V appears at least once in the submatrix induced by the selected rows.

Replacing C by m in the above yields the following result. Theorem 5. Let M be m × n binary matrix. A minimum size forbidden Tucker submatrix that occurs in M can be found in O(n 4 m 4 ) time.

Fig. 2 .

 2 Fig.2. Forbidden Tucker submatrices represented as vertex-colored bipartite graphs[START_REF] Tucker | A structure theorem for the consecutive 1s property[END_REF]. Black and white vertices correspond to rows and columns, respectively.

Fig. 3 .

 3 Fig. 3. MI k Tucker configuration.

Fig. 4 .

 4 Fig. 4. MII k Tucker configuration.

Fig. 5 .

 5 Fig. 5. MIII k Tucker configuration.

Fig. 6 .

 6 Fig. 6. MIV and MV Tucker configurations.

  [START_REF] Dom | Approximation and fixed-parameter algorithms for consecutive ones submatrix problems[END_REF] for each Tucker configuration. For example, if m = n, our algorithm has time complexity O(∆ 4 m 3 ) for M I k and M II k Tucker configurations, whereas Dom et al.'s algorithm has time complexity O(∆ m 3 ).

Theorem 4 .

 4 Let M be a m × n binary matrix with at most C 1-entries per column. A minimum size forbidden Tucker submatrix that occurs inM can be found in O(C 2 n 3 (m + C 2 n)) time.Proof. We apply Algorithms 1, Algorithms 2, and Algorithms 3 to find any forbidden submatrix of type M I k , M II k or M III k . Let us now analyze the time complexity.


Laurent Aldon
email: laurent.aldon@um2.fr

Alexis Perea

2D-correlation analysis applied to in situ and operando Mössbauer spectroscopy

Keywords: two-dimensional correlation spectroscopy, analytical method, Mössbauer spectrocopy, lithium battery, positive electrode materials, LiFePO 4, FePO 4, olivine, in situ and operando measurements, data treatment

published or not. The documents may come

2D-correlation analysis applied to in situ and operando

Mossbauer spectroscopy

Introduction

The electrochemical reactions of battery materials have been extensively studied by various electrochemical and spectroscopic techniques [1,2]. We recently determined Lamb-Mössbauer factors of both FePO4 and LiFePO 4 compounds [3] in order to get precise contributions of each formed/reacted species during electrochemical process. Mössbauer is a powerful technique to follow electrochemical reactions studied by ex situ [4,5], in situ [6] or operando conditions and also for characterization of pristine materials. In some case it is possible to enriched electrode materials with a Mössbauer isotope ( 57 Fe, 119 Sn) used as a local probe in Li 4 Ti 5 O 12 [7,8,9], Li 2 Ti 3 O 7 [10,11] for instance. In other cases spectra are broadened due to a mixture of phases or an amorphization of the particles. The help of electronic structure calculation [18,19] allows non-trivial experimental attributions when mathematical ambiguities due to overlapping peaks was observed. Interpretation of both isomer shift and quadrupole splitting variations was explained with such theoretical approaches [20,21]. In the case of glassy compounds quadrupole splitting distributions have been used [22,23] based on previous works [24,25,26].

Rechargeable lithium batteries have proven to be one of the most successful solutions to achieve good electrochemical performances in terms of cycleability, low cost, environmental benefit, easy to produce and safety in handling and operation. Phosphate materials (Li 3 Fe 2 (PO 4 ) 3 [27])have been investigated in the past for their good ionics properties as solid electrolytes. The LiFePO 4 compound has been chosen as cathode material because of the requirements [28,29]. Pure LiFePO 4 (triphylite) and FePO 4 (heterosite, berlinite [30,31], high pressure [32] and amorphous [33] polymorphs) have been earlier studied for their electrochemical performances [34,35,36] and their physico-chemical properties like antiferromagnetism [37] by X-ray [38,39] and neutron [40,41] diffraction techniques, Mössbauer spectroscopy [42,43,44], Raman [45,46] and Infra-Red [47,48,49] spectroscopies, X-ray absorption [50] and X-ray photoelectron spectroscopies [51,52]. The origin of capacity loss during the first cycle was also discussed for this good candidate as cathode material [53]. Due to low ionic and electronic conductivity for LiFePO 4 , and to its ageing in presence of moisture [54,55], optimized particle size [56], mineral coating or substitutions (LiFe 1-x Mn x PO 4 [57,58,59,60] or Li 1-x Fe 1+x P 1-x Si x O 4 [61]), doping by transition metals [62], or carbon composites [63] have been proposed to overcome these issues.

In this paper we illustrate, taking advantage of 2D-correlation spectroscopy analysis, a way to analyse Mössbauer spectra recorded in operando conditions. We have used Python [64] programming langage with Numpy [65] and Pylab [START_REF]Pylab and Matplotlib (Release 0.99.1, a Graphic Library[END_REF] libraries to get 2D-correlation spectra. This concept of 2D-correlation was initially proposed by Noda [START_REF] Noda | [END_REF]68]. In general 2D-correlation spectra are obtained using a perturbation (Temperature [69,70], pH [START_REF] Dou | Science in China[END_REF], Lithium amount [72], kinetic data [73]) applied to the pristine material investigated by a given or several experimental spectroscopies. More details on this analytical method can be found here [74]. From our knowledge 2D-correlation spectroscopy analysis (2D-COSA) used for various experimental techniques (Infra-Red spectroscopy [75,76], Mass spectral data [77], Resonance light scattering [78], X-ray photoelectron spectroscopy [79]) have not yet used for Mössbauer spectroscopy. We present here an application of 2D-COSA for 57 Fe in situ Mössbauer spectroscopy during electrochemical extraction of lithium from LiFe 0.75 Mn 0.25 PO 4 with olivine structure.

Theoretical basis

Generalized 2D Correlation Spectroscopy

A set of m experimental Mössbauer spectra are recorded as a function of a perturbation variable x (e.g., Li amount extracted/inserted during electrochemical tests) ranging from x min to x max . The Li variation between two collected spectra is x=(x max -x min )/(m-1). This continuous variation induces a broadening of the experimental spectra. Each raw spectrum is folded and normalized giving a set of m perturbationinduced dynamic spectra, is then transformed into a set of 2D correlation spectra by a form of cross-correlation analysis. We first define the reference spectrum of the system to be the averaged spectrum given by :

(1)

Then we need the dynamic spectra obtained by :

(2)

The intensity of 2D correlation spectrum X(v 1 ,v 2 ) reprensents the quantitative measure of a comparative similarity or dissimilarity of spectral intensity variations measured at two different velocities, v 1 and v 2 , along with a fixed interval. In a Mössbauer spectrum v i corresponds to the position of a Lorentzian absorption line. In order to simplify numerical manipulation, X(v 1 ,v 2 ) is treated as a complex number :

(3) comprising two orthogonal components, known respectively as the synchronous S(v 1 ,v 2 ) and the

asynchronous A(v 1 ,v 2 ) 2D correlation intensities.
Using Python programming language, S(v 1 ,v 2 ) is calculated as follows :

S=dot(y_dyn.conj().transpose(),y_dyn)/float(m-1)

The intensity of the synchronous 2D correlation spectrum represents the simultaneous changes of two separate spectral intensity variations measured at v 1 and v 2 within the interval between x min and x max .

The autocorrelation spectrum a(v i ), defined by the diagonal of S matrix, presents peaks called autopeaks, whose magnitude is always positive. It represents the overall extent of spectral intensity variation at the specific velocity v i within the interval between x min and x max . Crosspeaks located at the off-diagonal positions of a synchronous 2D spectrum represent simultaneous changes of spectral intensities observed at two differents velocities v 1 and v 2 . Such a synchronized change suggests the possible existence of a coupled origin of the spectral intensity variations. The sign of synchronous crosspeak becomes positive if the spectral intensities of the two velocities corresponding to the coordinates of the crosspeak are either increasing or decreasing together within the observation interval.

However, the negative sign of crosspeaks indicates that one of the spectral intensities in increasing while the other is decreasing.

The asynchronous 2D correlation spectrum A(v 1 ,v 2 ) is calculated in Python as follows:

A=dot(y_dyn.conj().transpose(),dot(Njk,spectres_dyn))/float(m-1)

In this previous expression N jk is the Hilbert-Noda matrix calculated as:

Experimental

Synthesis procedure LiFe 0.75 Mn 0.25 PO 4 pristine compound was obtained by solid-state reaction. The precursors:

Li 2 CO 3 (Sigma-Aldrich, 99%),, FeC 2 O 4 .2H 2 O (Labosi, 99.99%), Mn II (CH 3 COO) 2 (Sigma-Aldrich, 98%)

and NH 4 H 2 PO 4 (Acros Organics, 99.99%) were first ball milled during 1h30 and then thermally treated in a tube furnace under argon flow for 8 hours at 800°C.

X-ray diffraction pattern

Electrochemical measurements

For Mössbauer measurements in operando conditions, we used a designed cell allowing to use it in a transmission configuration. The cell has been designed in such a way that gamma rays can go through the Be windows which are the current collector, the sample under study is in contact with an aluminum foil avoiding corrosion of the Be windows at high potential, and the lithium foil used as reference electrode. Olivine has been used as actived cathode material with carbon black (80:20) versus Li used as anode with a Whatman paper as separator using 1M LiPF 6 in EC:PC:DMC. We used a rate of C/40, corresponding to 1 Li extracted in 40h with a potential window ranging from 3.2V to 4.3V.

Results and discussion

For LiFe 0.75 Mn 0.25 PO 4 sample, XRD measurements suggest pure phase within the detection limit of the apparatus as shown in figure 1. Our pristine material present XRD lines indexed in the Pnma space group with a set of cell parameters a=10.356(6) Å, b=6.025(4) Å and c=4.703(3) Å. Refined cell parameters have been obtained by pattern matching using FullProf program [80]. These values are in agreement with those published ealier [81].

In figure 2 we have represented the 57 Fe Mössbauer spectrum of the pristine material. The main contribution (91%) with an isomer shift (IS) of 1.20(1)mm/s and a quadrupole splitting (QS) of 2.96(1)mm/s correspond to Fe 2+ in the olivine structure. These values are in agreement with those found in the litterature [43]. Concerning the Fe 3+ contribution detected in our Mössbauer spectra, the hyperfine parameters are IS=0.31(4)mm/s and QS=0.85(8)mm/s. These values are not very closed to those of Li x FePO 4 (OH) x (taverite) [55] the reaction product obtained from LiFePO 4 ageing under moisture. The presence of such an impurity is not excluded since about 9% of Fe 3+ is present and only 90% of the Li can be extracted by electrochemical reaction as can be seen in Figure 3. An other explanation of the Fe 3+ Li are used to oxidize Fe 2+ to Fe 3+ and Mn 2+ to Mn 3+ , respectively. These values are rather close to the expected one : 0.75×0.91 = 0.683 to oxidize Fe 2+ and 0.25×0.91 = 0.227 to oxidize Mn 2+ , in a rather good agreement with the stoechiometric ratio Fe:Mn of 0.75:0.25.

We have represented in figure 4 the 57 Fe Mössbauer spectra obtained during electrochemical Li extraction as mentioned above. We clearly see the evolution of the absorption lines upon the perturbation, i.e. the Li content. In spectrum labeled 1, we retrieve the doublet of the olivine (LiFe 2+ 0.75 Mn 2+ 0.25 PO 4 ) constituted of the two peaks located at v 1 =-0.28mm/s and v 2 =2.68mm/s, corresponding to the hyperfine parameters IS=1.2mm/s and QS=2.96mm/s. After extraction of x=0.68

Li (spectrum labeled 13), we observe the doublet of the lithiophilite (Li 0.25 Fe 3+ 0.75 Mn 2+ 0.25 PO 4 ) constituted of the two peaks located at v 1 =-0.22mm/s and v 2 =0.98mm/s, corresponding to the hyperfine parameters IS=0.38mm/s and QS=1.20mm/s. At the end of lithium extraction (spectrum labeled 18), we observe the doublet of the heterosite (Fe 3+ 0.75 Mn 3+ 0.25 PO 4 ) constituted of the two peaks located at v 1 =-0.36mm/s and v 2 =1.18mm/s, corresponding to the hyperfine parameters IS=0.41mm/s and QS=1.54mm/s.

In figure 5 we have plotted the averaged spectrum wich give an indication of the regions affected by Li extraction. We see again three main peaks located at about [-0.2, 1.1 and 2.6mm/s]. This averaged spectrum is needed to calculate both 2D-correlation synchronous (Figure 6) and asynchronous (Figure 7) spectra.

In the 2D-correlation synchronous spectrum (figure 6 The transmission intensity is slightly correlated due to the presence of small positive cross-peak with [c,g] coordinates in the 2D-correlation synchronous spectrum. In figure 7, cross-peak located at [c,g] presents a null intensity in the 2D-correlation asynchronous spectrum. In such analysis, positive (negative) values means that variation at a given velocity occurs before (after) the perturbation. In the context of Mössbauer spectroscopy, a close-to-zero value in the 2D-correlation asynchronous spectrum suggests that variation observed at 2 given velocities are coincidental. Hence, these two positions belongs to the same doublet of a specie defined by an isomer shift, IS=(v 1 +v 2 )/2 and a quadrupole splitting, QS = v 2 -v 1 . Using this kind of analysis, we propose an attribution given in table 1. In table 2, we have reported the hyperfine parameters deduced from the 2D-correlation analysis and those obtained from a classical fitting procedure. Our 2D-correlation becomes a new tool for analyses of rather complicated Mössbauer spectra presenting some overlaps between absorption lines and small variation from a spectrum to an other. Another point to emphasize is that 2D-correlation applied to Mössbauer spectroscopy allows to extract both isomer shift and quadrupole splitting of the species using the whole set of spectra instead of using data separately.

The autocorrelation spectrum plotted in figure 8 shows the positions of these characteristic values.

These sets of coordinates are used to represent the relative transmission variation upon Li content as shown in figure 9. Hence, the presence of Mn 3+ in the neighborhood of Fe 3+ induces a higher quadrupole splitting due to the second neighbor shell. This can be related to change induced in the local symmetry of the Fe 3+ site due the Jahn-Teller effect on Mn 3+ (3d 4 ). We also observed this by in situ XRD measurement (not the purpose of this paper) and our results are rather in a good agreement with those previously published by Yamada et al [63].

Relative transmission at velocity [g] continously increases from

Concerning the Mössbauer investigations, we have shown by in situ and operando measurements that 57 Fe is directly sensitive to oxidation state of the iron (Fe 2+ /Fe 3+ ) during the first stage of the Li extraction as expected in figure 3. We have also shown that 57 Fe is indirectly sensitive the the second neighbor shell through variations of Mn oxidation state. The quadrupole splitting is ranging from 1.20mm/s when Fe 3+ are surrounded by Fe 3+ and Mn 2+ sharing the M2 site of the olivine structure, to 1.54mm/s when Fe 3+ are surrounded by Fe 3+ and Mn 3+ .

Conclusion

2D-Correlation technique can be used in other interesting application like catalysis, in situ corrosion

measurements, f-Lamb-Mössbauer factor determination for instance.

  or 0 if j=k The intensity of an asynchronous 2D correlation spectrum represents sequential, but not coincidental, changes of spectral intensities measured separately at v 1 and v 2 . Crosspeaks develop only if the intensity varies out of phase with each other. The sign of asynchronous crosspeak provide useful information on sequential order of events observed by the Mössbauer spectroscopy along the lithium content. If the signs of synchronous and asynchronous crosspeaks are the same, the intensity change at v 1 occurs before v 2 . If the signs of synchronous and asynchronous crosspeaks are different, the intensity change at v 1 occurs after v 2 .

LiFe 0. 75 40 ( 1

 75401 Mn 0.25 PO 4 pristine compound has been characterized by X-ray Power Diffraction with a Philips -2 diffractometer using CuK radiation (=1.5418 Å) and a nickel filter.Crystallographic reference data for LiFePO 4 (triphylite) has been taken from Powder diffraction File (JCPDS #81-1173, a=10.332 Å, b=6.010 Å, c=4.692 Å). These values and the orthorhombic Pnma space group have been used in the initial stages of refinement of the cell parameters for our LiFe 0.75 Mn 0.25 PO 4 compound.57 Fe Mössbauer measurementsMössbauer spectra have been recorded in the constant acceleration mode and in transmission geometry on a standard Mössbauer spectrometer composed of devices from Ortec and WissEl. A 57 Co(Rh) source with a nominal activity of 10 mCi have been used. Each Mössbauer spectrum has been collected every 128minutes with a deintercalation uptake x=0.053 Li corresponding to an electrochemical rate of C/Li extracted in 40 hours). The source has been always kept at RT. Calibration of the velocity scale was done using absorption lines of a 25µm iron foil. Isomer shifts are given relative to -Fe.

2D-

  Correlation spectroscopy applied to Mössbauer spectroscopy during electrochemical studies of electrode materials in real conditions (in situ and operando) can be considered as a new tool for the Mössbauer community. This analytical technique gives valuable information on the electrochemical mechanisms involved in extraction (insertion) of lithium from (in) electrode materials. Compared to other spectroscopies, Mössbauer spectroscopy is a complementary technique since this electrode material is based on Fe 2+ /Fe 3+ redox system, 57 Fe and gives valuable information at atomic scale in the bulk material, while XPS and Raman are more dedicated to surface analysis. With the help of 2D-correlation analysis, Mössbauer spectroscopy can be used in operando conditions, i.e. during Li insertion (extraction) in (from) electrode materials tested in real battery with the presence of current collector (Be windows), separator and Li foil. Hence Mössbauer effect measurements yield the evidence for the oxidation processes taking place in the course of lithium extraction.

Figure 1 .

 1 Figure 1. X-ray diffraction pattern of pristine LiFe0.75Mn0.25PO4. Cell parameters determined by LeBail method[80] using the Pnma space group are a=10.356(6) Å, b=6.025(4) Å, c=4.703(3) Å.

Figure 2 .

 2 Figure 2. 57Fe Mössbauer spectrum recorded at room temperature. Experimental (open circle) and calculculated spectrum (solid line) are represented. A small contribution (in red) due to Fe3+ impurity (see text for explanation) is also shown for clarity.

Figure 3 .

 3 Figure 3. Electrochemical curve recorded during lithium extraction from LiFe0.75Mn0.25PO4. The numbered labels correspond to the file number of the recorded 57Fe Mössbauer spectra in operando conditions.

Figure 4 .

 4 Figure 4. 57Fe Mössbauer spectra recorded in operando conditions during Li extraction from pristine LiFe0.75Mn0.25PO4. Labels 1 and 18 correspond to LiFe0.75Mn0.25PO4 and Fe0.75Mn0.25PO4 respectively. Colored lines [a-g] correspond to cross-correlation peaks positions observed in the synchronous spectra. The evolution of the intensity for those specific positions are discussed in the text.

Figure 5 .

 5 Figure 5. Averaged spectrum determined from the mean of the 18 Mössbauer spectra recorded during extraction of Li. The main absorption peak located at negative velocity is due to the overlap between contributions of both Fe2+ and Fe3+. The peak located near 1mm/s correspond to Fe3+ contribution and the peak at 2.6mm/s correspond to Fe2+ contribution.

Figure 6 .

 6 Figure 6. 2D-correlation synchronous spectrum obtained from experimental 57Fe spectra recorded during electrochemical lithium extraction from LiFe0.75Mn0.25PO4. Positive (light green) and negative (dark green) crosspeaks are observed. Positive crosspeaks have the following coordinates : [a,f], [d,e] and [c,g], while negative crosspeaks have mainly the following coordinates: [a,g], [e,g] and [f,g].

Figure 7 .

 7 Figure 7. 2D-correlation asynchronous spectrum obtained from experimental 57Fe spectra recorded during electrochemical lithium extraction from LiFe0.75Mn0.25PO4. Positive crosspeaks (light green) negative crosspeaks (dark green) suggest variation at a given velocity occurs before (after) for an other velocity.

Figure 8 .

 8 Figure 8. Autocorrelation spectrum determined from the 2D-correlation synchronous spectrum. Features previously observed in both synchronous and asynchronous correspond to either close to maxima in the autocorrelation spectrum [a, d, e, f] or minimum [c]. This latter peak is more broadened than those centered at about 2.7mm/s. It contains some meaningful information as can be read in the text.

Figure 9 .

 9 Figure 9. Evolution of the relative transmission for various velocities [a-g] vs. file number. A decrease in the relative transmission correspond to an increase of absorption at a given velocity.

Figure 1 . 4 Figure 3 .

 143 Figure 1. X-ray diffraction pattern of pristine LiFe 0.75 Mn 0.25 PO 4 . Cell parameters determined by LeBail

Figure 4 .

 4 Figure 4. 57 Fe Mössbauer spectra recorded in operando conditions during Li extraction from pristine

Figure 5 .

 5 Figure 5. Averaged spectrum determined from the mean of the 18 Mössbauer spectra recorded during

Figure 6 .

 6 Figure 6. 2D-correlation synchronous spectrum obtained from experimental 57 Fe spectra recorded

Figure 7 .

 7 Figure 7. 2D-correlation asynchronous spectrum obtained from experimental 57 Fe spectra recorded

Figure 8 .

 8 Figure 8. Autocorrelation spectrum determined from the 2D-correlation synchronous spectrum.

Figure 9 .

 9 Figure 9. Evolution of the relative transmission for various velocities [a-g] vs. file number. A decrease

  contrubution may due to Fe 2 O 3 . From Mössbauer data recorded at various temperature, it has been shown that anisotropy constant can be determined[82]. Hence for -Fe 2 O 3 57 Fe Mössbauer spectra present rather a doublet instead of a sextet at room temperature when particles are nanosized. This doublet presents hyperfine parameters ranging from IS=0.32mm/s QS=0.98mm/s (D=4nm) to IS=0.38mm/s, QS=0.44mm/s (D=18nm) depending on the mean particle size. So, in our case, both IS and QS could be some Fe 2 O 3 formed at high temperature. From linear correlation between QS and nanoparticle size we found a mean particle size of D=7±2nm for our Fe 2 O 3 impurity. Recently, Rho et al. have mentioned[83] the formation of both Fe 2 P (singlet+doublet) and FeP (IS=0.6(3)mm/s and QS=0.8(6)mm/s) as impurity in LiFePO 4 . The electrochemical curve shown in figure 3. presents two main regions. Starting for [x=1 Li,V=3.2V] in Li x Fe 0.75 Mn 0.25 PO 4 we can extract about x=0.68 Li (Mössbauer spectra#1 to 13) on a mean voltage of about 3.48V. This voltage is known to be Fe 2+ /Fe 3+ redox potential. Then from [x=0.32Li,V=3.88V] to [x=0.10,V=4.3V] we observe a mean voltage at about 4.02V, corresponding to Mn 2+ /Mn 3+ redox potential (Mössbauer spectra#14 to 18). Assuming that 0.91Li have been used (9% of nonelectrochemically active impurity observed in the Mössbauer spectrum, fig. 2), x=0.68 Li and x=0.23

  ), we see autopeaks with positive contribution (diagonal represented in figure8) and crosspeaks with positive and negative intensities depending on the respective evolution of the intentities of two given velocities. Positive and negative crosspeaks have

the following coordinates: (

[a,f], [d,e] and [c,g]) and ([a,g], [e,g

] and [f,g]) respectively. In Table

1

, we have reported the characteristic values of these coordinates.

  file #1 to #13 or #14 and then still constant up to the end of the reaction (file#18). This unambiguous position corresponds to the oxidation process of Fe 2+ to Fe 3+ in Li x Fe 0.75 Mn 0.25 PO 4 . From x=1 (file #1) to x=0.34 (between file #13 and #14), we can write the redox process as the following:Now, focussing our attention on relative transmission at velocities [d,e,f], we can have an idea of the mechanism involved during the next steps of the Li extraction. In fact, it corresponds in some way to what happens in the Fe 3+ local environment. It is clear that any variation detected in the Mössbauer spectra may be due to oxidation of Mn 2+ to Mn

	LiFe 2+	0.75 Mn 2+	0.25 PO 4 -yLi → Li 1-y Fe 2+	0.75-y Fe 3+	y Mn 2+ 0.25 PO 4 .	(0<y<0.75) (i)
	Li 0.25 Fe 3+	0.75 Mn 2+ 0.25 PO 4 -(y-0.75) Li → Li 1-y Fe 3+	0.75 Mn 2+	1-y Mn 3+	y-0.75 PO 4	(0.75<y<1)	(ii)

3+

. Quadrupole splitting is known to be very sensitive to second neighbors effect, and to the charge of the ions. In figure

9

we can see that the relative transmission at [e] decreases from file#1 to file#13 and still constant until the end of the Li extraction.

While the relative transmission at [d] decreases from file#1 to file#13 and then increases from file#14 to file#18. These observations suggest the following process:

TABLES . Table 1 .

 .1 Velocities of the main features [a-g] deduced from crosspeaks analysis of both 2D-correlation synchronous and asynchronous spectra and autocorrelation spectrum. Some attributions are given as an help for the discussion.

	Features	a	b	c	d	e	f	g
	Velocity	-0.52	-0.36	-0.22	-0.02	1.02	1.18	2.68
	(mm/s)
	Attribution		v 1 (Fe III 2 )	v 1 (Fe II )	v 1 (Fe III 1 )	v 2 (Fe III 1 )	v 2 (Fe III 2 )	v 2 (Fe II )

Table 2 .

 2 Comparison of isomer shift and quadrupole splitting ( calc and  calc in mm/s) deduced from 2Dcorrelation analysis and those obtained by least-square fitting procedure ( fit and  fit in mm/s).

		v 1	v 2	 calc	 calc	 fit	 fit
	Fe II	-0.22(c)	2.68(g)	1.23	2.90	1.20	2.96
	Fe III 1	-0.02(d)	1.02(e)	0.52	1.04	0.38	1.20
	Fe III 2	-0.36(a)	1.18(f)	0.41	1.54	0.41	1.54

Acknowledgements

The authors would like to thank SAFT Company, Bordeaux and CNRS (contract SAFT/CNRS N°029888) for the financial support through the Ph.D. grant of Alexis Perea and the Région Languedoc Roussillon for the Mössbauer platform.

(see text for explanation) is also shown for clarity.