Christine Keribin

Vincent Brault

Gilles Celeux

Gérard Govaert

Estimation and Selection for the Latent Block Model on Categorical Data

Keywords: EM algorithm, Variational Approximation, Stochastic EM, Bayesian Inference, Gibbs Sampling, BIC Criterion, Integrated Completed Likelihood Algorithme EM, approximation variationnelle, Stochastique EM, inférence bayésienne, échantillonneur de Gibbs, critère BIC, critère ICL

This paper is dealing with estimation and model selection in the Latent Block Model (LBM) for categorical data. First, after providing sufficient conditions ensuring the identifiability of this model, it generalises estimation procedures and model selection criteria derived for binary data. Secondly, it develops Bayesian inference through Gibbs sampling. And, with a well calibrated non informative prior distribution, Bayesian estimation is proved to avoid the traps encountered by the LBM with the maximum likelihood methodology. Then model selection criteria are presented. In particular an exact expression of the integrated completed likelihood (ICL) criterion requiring no asymptotic approximation is derived. Finally numerical experiments on both simulated and real data sets highlight the interest of the proposed estimation and model selection procedures.

Introduction

Block clustering methods are aiming to design in a same exercise a clustering of the rows and the columns of a large array of data. These methods could be expected to be useful to summarise large data sets by dramatically smaller data sets with the same structure. Since more and more huge data sets are available, more and more block clustering methods have been proposed. Many application fields as genomic [START_REF] Jagalur | Analyzing in situ gene expression in the mouse brain with image registration, feature extraction and block clustering[END_REF] or recommendation system [START_REF] Shan | Bayesian co-clustering[END_REF] are concerned with block clustering. In particular, block clustering methods have been developed to deal with binary data present in archaeology [START_REF] Govaert | Classification croisée[END_REF] and sociology [START_REF] Wyse | Block clustering with collapsed latent block models[END_REF]. [START_REF] Inria Madeira | Biclustering algorithms for biological data analysis: A survey[END_REF] described an extensive list of block clustering methods. The checker board structure studied in this article has been considered from two point of views: determinist approaches (see for instance [START_REF] Banerjee | A generalized maximum entropy approach to bregman co-clustering and matrix approximation[END_REF][START_REF] Govaert | Algorithme de classification d'un tableau de contingence[END_REF] and modelbased approaches. The model-based view has been considered through the maximum likelihood methodology [START_REF] Govaert | Block clustering with bernoulli mixture models: Comparison of different approaches[END_REF] and through Bayesian inference [START_REF] Meeds | Nonparametric bayesian biclustering[END_REF][START_REF] Wyse | Block clustering with collapsed latent block models[END_REF]. Among them, the Latent Block Model (LBM) is attractive since it could lead to powerful representations. For instance, as illustrated in [START_REF] Govaert | Block clustering with bernoulli mixture models: Comparison of different approaches[END_REF], summarising binary tables with grey summaries derived from the conditional probabilities of belonging to a block could be quite realistic and suggestive. Moreover, this model provides a probabilistic framework to choose a relevant block clustering.

The LBM is as follows. A population of n observations described with d categorical variables of the same nature with r levels is available. Saying that the categorical variables are of the same nature means that it is possible to code them in a same (and natural) way. This assumption is needed to ensure that decomposing the data set in a block structure is making sense. Let y = (y ij , i = 1, . . . , n; j = 1, . . . , d) be the data matrix where y ij = h, 1 ≤ h ≤ r, r being the number of levels of the categorical variables.

It is assumed that there exists a partition into g row clusters z = (z ik ; i = 1, . . . , n; k = 1, . . . , g) and a partition into m column clusters w = (w j ; j = 1, . . . , d; = 1, . . . , m). The z ik s (resp. w j s ) are binary indicators of row i (resp. column j) belonging to row cluster k (resp. column cluster ), such that the random variables y ij are conditionally independent knowing z and w with parameterised density ϕ(y ij ; α k ) z ik w j . Thus, the conditional density of y knowing z and w is f (y|z, w; α) = i,j,k, ϕ(y ij ; α k ) z ik w j where α = (α k ; k = 1, . . . , g; = 1, . . . , m). Moreover, it is assumed that the row and column labels are independent: p(z, w) = p(z)p(w) with p(z) = ik π z ik k and p(w) = j ρ w j , where (π k = P (z ik = 1), k = 1, . . . , g) and (ρ = P (w j = 1), = 1, . . . , m) are the mixing proportions. Hence, the marginal density of y is a mixture density f (y; θ) = (z,w)∈Z×W p(z; π)p(w; ρ)f (y|z, w; α), Z and W denoting the sets of possible row labels z and column labels w, and θ = (π, ρ, α), with π = (π 1 , . . . , π g ) and ρ = (ρ 1 , . . . , ρ m ). The density of y can be written as

f (y; θ) = z,w i,k π z ik k j, ρ w j i,j,k, ϕ(y ij ; α k ) z ik w j . (1)
The LBM involves a double missing data structure, namely z and w, which makes statistical inference more difficult than for standard mixture model. [START_REF] Govaert | Block clustering with bernoulli mixture models: Comparison of different approaches[END_REF] show that the EM algorithm is not tractable and proposed a variational approximation of the EM algorithm (VEM) to derive the maximum likelihood (ML) estimate of θ. This VEM algorithm could give satisfactory estimates despite it is highly dependent of its initial values and has a marked tendency to provide empty clusters. Moreover, even if the parameter θ is properly estimated, and for even moderate data matrices, it is unreasonable to compute the likelihood. As a consequence, computing penalised model selection criteria such as AIC or BIC is challenging as well. Moreover, these difficulties to compute relevant model selection criteria are increased by the fact that the LBM statistical units could be defined in several different ways.

The aim of this paper is to overcome all those limitations. First, we propose algorithms aiming to avoid the above mentioned drawbacks of the VEM algorithm. Secondly, we show how it is possible to compute properly relevant model selection criteria. This is presented for the LBM for categorical data, a natural extension of the LBM for binary data, but not developed as so far. The article is organised as follows. In Section 2, the LBM is detailed for categorical data and sufficient conditions ensuring the identifiability of this model are provided. Section 3 is devoted to the presentation of a stochastic algorithm aiming to avoid the variational approximation, and it shows how Bayesian inference can be helpful to avoid empty clusters. Section 4 is concerned with model selection criteria. We first show how the integrated completed likelihood (ICL) of the LBM is closed form and derive from ICL an approximation of BIC. Section 5 is devoted to numerical experiments on both simulated and real data sets to illustrate the behaviour of the proposed estimation algorithms and model selection criteria. A discussion section proposing some perspectives for future work ends this paper.

Notation To simplify the notation, the sums and products relative to rows, columns, row clusters, column clusters and levels will be subscripted respectively by the letters i, j, k, , h without indicating the limits of variation that will be implicit as in (1). So, for instance, the sum

Model definition and identifiability

The LBM has already been defined on binary [START_REF] Govaert | Block clustering with bernoulli mixture models: Comparison of different approaches[END_REF][START_REF] Keribin | Model selection for the binary latent block model[END_REF], Gaussian [START_REF] Lomet | Sélection de modèle pour la classification croisée de données continues[END_REF]) and count data [START_REF] Govaert | Latent block model for contingency table[END_REF]). We consider here the LBM for categorical data, where the conditional distribution ϕ(y ij ; α kl ) of the outcome y ij knowing the labels z ik and w j is a multinomial distribution M(1; α k ), with parameter α k = (α h k ) h=1,...,r , where α h k ∈ (0; 1) and h α h k = 1. Using the binary indicator vector y ij = (y 1 ij , . . . , y r ij ), such that for every h = 1, . . . , r, y h ij = 1 when y ij = h and y h ij = 0 otherwise, the density per block is

ϕ(y ij ; α k ) = h (α h k ) y h ij
and the mixture density is

f (y; θ) = z,w i,k π z ik k j, ρ w j i,j,k, ,h (α h k ) z ik w j y h ij .
The parameter to be estimated is θ = (π, ρ, α), with g + m + (r -1)gm -2 independent components. Before considering the estimation problem, it is important to analyse the model identifiability (Frühwirth-Schnatter 2006, pp. 21-23). Obviously, LBM, as a mixture model, is not identifiable due to invariance to relabelling the blocks, but it is of no importance when concerned with Inria maximum likelihood estimation. It will be necessary to go back to this issue when concerned with Bayesian estimation. Unfortunately, it is also well-known that simple multivariate Bernoulli mixtures are not identifiable [START_REF] Gyllenberg | Non-uniqueness in probabilistic numerical identification of bacteria[END_REF], regardless of this invariance to relabelling. [START_REF] Allman | Identifiability of parameters in latent structure models with many observed variables[END_REF] set a sufficient condition to their identifiability that cannot be easily extended to LBM. We define here a set of sufficient conditions ensuring the identifiability of the Bernoulli LBM:

Theorem 1 (Binary LBM identifiability) : Consider the binary LBM. Let π and ρ be the row and column mixing proportions and A = (α kl ) the g × m matrix of Bernoulli parameters. Define the following conditions:

• C 1 : for all 1 ≤ k ≤ g, π k > 0 and the coordinates of vector τ = Aρ are distinct.

• C 2 : for all 1 ≤ ≤ m, ρ > 0 and the coordinates of vector σ = π A are distinct (where π is the transpose of π).

If conditions C 1 and C 2 hold, then the binary LBM is identifiable for n ≥ 2m -1 and d ≥ 2g -1.

Proof. The proof of this theorem is given in Appendix A.

Conditions C 1 and C 2 are not strongly restrictive since the set of vectors τ and σ violating them is of Lebesgue measure 0. Therefore, Theorem 1 asserts the generic identifiability of LBM, which is a "practical" identifiability, explaining why it works in the applications [START_REF] Carreira-Perpiñàn | Practical identifiability of finite mixtures of multivariate bernoulli distributions[END_REF]. These conditions could appear to be somewhat unnatural. However:

(i) It is not surprising that the number of row labels g (resp. column labels m) is constrained by the number of columns d (resp. rows n). In case of a simple finite mixture of g different Bernoulli products with d components, more clusters you define in the mixture, more components for the multivariate Bernoulli you need in order to ensure the identifiability: see also the following condition d > 2 log 2 g + 1 of [START_REF] Allman | Identifiability of parameters in latent structure models with many observed variables[END_REF] for simple Bernoulli mixtures, where . is the ceil function.

(ii) Conditions C 1 and C 2 are extensions of a condition set to ensure the identifiability of the Stochastic Block Model [START_REF] Celisse | Consistency of maximum-likelihood and variational estimators in the stochastic block model[END_REF], where n = d and z = w.

It follows from conditions C 1 (resp. C 2 ) that the probabilities τ k = P (y ij = 1|z ik = 1) (resp. σ = P (y ij = 1|w j = 1)) to observe an event in a cell of a row of row class k (resp. in a cell of a column of column class ) can be sorted in a strictly ascending order. Hence, contrary to what happens in the Gaussian mixture context, these conditions can be used to put a natural order on the mixture labels in a Bayesian setting. Notice that conditions C 1 and C 2 are sufficient for the identifiability, and can be proved to be not necessary for g = 2 and m = 2.

Theorem 1 is easily extended to the categorical case, where the conditions are defined on vectors τ h = A h ρ and σ h = π A h , with h = 1, . . . , r -1 and A h = (α h kl ) k=1,...,g; =1,...,m .

Model estimation

With g and m fixed, the likelihood of the model parameter L(θ) = f (y; π, ρ, α) is written

L(θ) = z,w i,k π z ik k j, ρ w j i,j,k, ,h (α h k ) z ik w j y h ij .
Even for small tables, computing this likelihood (or its logarithm) is difficult. For instance, with a data matrix 20 × 20 with g = 2 and m = 2, it requires the calculation of g n × m d ≈ 10 12 terms. In the same manner, deriving the maximum likelihood estimator with the EM algorithm is challenging. As a matter of fact, the E step requires the computation of the joint conditional distributions of the missing labels P (z ik w j = 1|y; θ (c) ) for i = 1, . . . , n, j = 1, . . . , d, k = 1, . . . , g, = 1, . . . , m with θ (c) being a current value of the parameter. Thus, the E step involves computing too many terms that cannot be factorised as for a standard mixture, due to the dependence of the row and column labels knowing the observations.

To tackle this problem, [START_REF] Govaert | Block clustering with bernoulli mixture models: Comparison of different approaches[END_REF] proposed to use a variational approximation of the EM algorithm by imposing that the conditional joint distribution of the labels knowing the observations factorizes as q

(c) zw (z, w) = q (c) z (z)q (c)
w (w). This VEM algorithm has been proved to provide relevant estimates of the LBM model in different contexts with continuous or binary data matrices [START_REF] Govaert | Block clustering with bernoulli mixture models: Comparison of different approaches[END_REF], and is directly extended to categorical data, see Appendix B. However, it presents several drawbacks illustrated in Section 5:

(i) as most variational approximation algorithms, the VEM appears to be quite sensitive to starting values, (ii) it has a marked tendency to provide solutions with empty clusters, i.e. with fewer clusters than required after a maximum a posteriori (MAP) classification rule, (iii) it can only compute a lower bound of the maximum likelihood called free energy.

A possible way to attenuate the dependence of VEM to its initial values is to use stochastic versions of EM which are not stopping at the first encountered fixed point of EM (see McLachlan and Krishnan, 2008, chap. 6). The basic idea of these stochastic EM algorithms is to incorporate a stochastic step between the E and M steps where the missing data are simulated according to their conditional distribution knowing the observed data and a current estimate of the model parameters. For the LBM, it is not possible to simulate in a single exercise the missing labels z and w and a Gibbs sampling scheme is required to simulate the couple (z, w). The SEM-Gibbs algorithm for binary data [START_REF] Keribin | Model selection for the binary latent block model[END_REF] is a simple adaptation to the LBM of the standard SEM algorithm of [START_REF] Celeux | Stochastic versions of the em algorithm[END_REF]. It is directly extended for the categorical LBM, see Appendix C.

While VEM is based on the variational approximation of the LBM, SEM-Gibbs uses no approximation, but runs a Gibbs sampler to simulate the unknown labels with their conditional distribution knowing the observations and a current estimation of the parameters. Hence, SEM-Gibbs is not increasing the loglikelihood at each iteration, but it generates an irreducible Markov chain with a unique stationary distribution which is expected to be concentrated around the ML parameter estimate. Thus a natural estimate of θ derived from SEM-Gibbs is the mean θ of (θ (c) ; c = B + 1, . . . , B + C) get after a burn-in period of length B. Moreover, as every stochastic algorithm, SEM-Gibbs is theoretically subject to label switching (see Frühwirth-Schnatter, 2006, Section 3.5.5). However this possible drawback of SEM-Gibbs does not occur in most practical situations. Numerical experiments presented in [START_REF] Keribin | Model selection for the binary latent block model[END_REF] for binary data show that SEM-Gibbs is by far less sensitive to starting values than VEM. Those results lead them to advocate initializing the VEM algorithm with the SEM-Gibbs mean parameter estimate θ to get a good approximation of the ML estimate for the latent block model.

If SEM-Gibbs can be expected to be insensitive to its initial values, there is no reason to think that it can be useful to avoid solutions with empty clusters. Bayesian inference in statistics can be regarded as a well-ground tool to regularize ML estimate in a poorly posed setting. In the LBM setting, Bayesian inference could be thought of as useful to avoid empty cluster solutions and thus to attenuate the "empty cluster" problem. In particular, for the categorical LBM, it is possible to consider proper and independent non informative prior distributions for the mixing The Bayes formula log p(θ|y) = log p(y|θ) + log p(θ) -log p(y) allows to straightforwardly define an EM algorithm for the computation of the MAP estimate:

• the E Step relies on the computation of the conditional expectation of the complete loglikelihood Q(θ, θ (c) ) as for the ML estimator:

Q(θ, θ (c) ) = IE log p(y, z, w|θ)|y, θ (c) ,
• the M Step differs in that the objective function for the maximisation process is equal to the Q(θ, θ (c) ) function augmented by the logarithm of the prior density

θ (c+1) = arg max θ Q(θ, θ (c) ) + log p(θ) .
This M Step forces an increase in the log posterior function p(θ|y) (McLachlan and Krishnan, 2008, chap. 6, p. 231). For LBM, the Bayesian approach combined with the variational approximation leads to the V-Bayes algorithm, with the following M step:

M step (V-Bayes algorithm: estimation of the posterior mode)

π (c+1) k = a -1 + s (c+1) .k n + g(a -1) , ρ (c+1) = a -1 + t (c+1) . d + m(a -1) α h k (c+1) = b -1 + i,j s (c+1) ik t (c+1) j y h ij r(b -1) + s (c+1) .k t (c+1) .
, where s (c+1) ik and t

(c+1) j are the current conditional probabilities of the labels and s

(c+1)

.k and t (c+1) .

, their sum over the lines and columns respectively, as defined in the VEM algorithm, see Appendix B. The hyperparameters a and b are acting as regularisation parameters. In this perspective, the choice for a and b is important. It appears clearly from the updating equations of this M step that V-Bayes collapses to the VEM algorithm with the uniform prior a = b = 1 and involves no regularisation. And, for Jeffreys prior (a = b = 1/2), V-Bayes has a tendency to provide more empty cluster solutions than VEM. [START_REF] Frühwirth-Schnatter | Mixtures : estimation and applications, chapter Dealing with label switching under model uncertainty[END_REF] shows the great influence of a for Bayesian inference in the Gaussian mixture context. Based on an asymptotic analysis by [START_REF] Rousseau | Asymptotic behaviour of the posterior distribution in overfitted models[END_REF] and on a thorough finite sample analysis, she advocated to take a = 4 for moderate dimension (g < 8) and a = 16 for larger dimensions to avoid empty clusters.

In Section 5, we present numerical experiments highlighting the ability of the V-Bayes algorithm to avoid the "empty cluster" cases encountered with the VEM algorithm.

However, as VEM, a V-Bayes algorithm could be expected to be highly dependent of its initial values. Thus it could be of interest to initialise V-Bayes with the solution derived from the Gibbs sampler. Obviously, since Dirichlet prior distributions are conjugate priors for the multinomial distribution, full conditional posterior distributions of the LBM parameters are closed form and Gibbs sampling is easy to implement, in such way:

Gibbs sampling for the categorical LBM 1. simulation of z (c+1) according to p(z|y, w (c) ; θ (c) ) as in SEM-Gibbs 2. simulation of w (c+1) according to p(w|y, z (c+1) ; θ (c) ) as in SEM-Gibbs 3. simulation of π (c+1) according to

D(a + z (c+1) .1 , . . . , a + z (c+1)
.g ) where z .k = i z ik denotes the number of lines in line cluster k 4. simulation of ρ (c+1) according to

D(a + w (c+1) .1 , . . . , a + w (c+1) .m
) where w . = j w j denotes the number of columns in column cluster 5. simulation of α

(c+1) k according to D(b + N 1 k (c+1) , . . . , b + N r k (c+1) )
for k = 1, . . . , g; = 1, . . . , m and with

N h k (c+1) = i,j z (c+1) ik w (c+1) j y h ij . (2)
As Gibbs sampling explores the whole distribution, it is subject to label switching. This problem can be sensibly solved for the categorical LBM by using the identifiability conditions of Theorem 1: these conditions define a natural order of row and column labels except on a set of parameters of measure zero. Hence, row and column labels are post processed separately after the last iteration, by reordering them according to the ascending values of τ h = α h ρ and σ h = π α h coordinates respectively, and for a given h.

Inria 4 Model selection

Choosing relevant numbers of clusters in a latent block model is obviously of crucial importance. This model selection problem is difficult for several reasons. First there is a couple (g, m) of number of clusters to be selected instead of a single number. Secondly penalised likelihood criteria such as AIC or BIC are not directly available since computing the maximised likelihood is not feasible. Third determining the number of statistical units of a LBM could be questionable (number of rows, number of columns, number of cells, . . .).

Fortunately, it is possible to compute the exact integrated completed loglikelihood (ICL) of the categorical LBM.

The Integrated Completed Loglikelihood (ICL)

ICL is the logarithm of the integrated completed likelihood

p(y, z, w | g, m) = Θg,m p(y, z, w | θ g,m )p(θ g,m )dθ g,m
where the missing data (z, w) have to be replaced by some values closely related to the model at hand. By taking into account the missing data, ICL is focusing on the clustering view of the model. For this very reason, ICL could be expected to select a stable model allowing to partitioning the data with the greatest evidence [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF].

As stated in Section 3, Dirichlet proper non informative priors are available for the categorical LBM and ICL can be computed without requiring asymptotic approximations. Using the conditional independence of the zs and the ws conditionally to α, π, ρ, the integrated completed likelihood can be written ICL = log p(z) + log p(w) + log p(y|z, w).

( -log Γ(n + ga) -log Γ(d + ma)

+ k log Γ(z .k + a) + log Γ(w . + a) + k, h log Γ(N h k + b) -log Γ(z .k w . + rb) .
where z .k , w . and N h k are defined as for the Gibbs sampler (see Equation 2). In practice, the missing labels z, w have to be chosen. Following [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF], they are replaced by (ẑ, ŵ) = arg max (z,w) p(z, w|y; θ), θ being the estimate of the LBM parameter computed from the V-Bayes algorithm initialised by the Gibbs sampler as described in Section 3. The maximising partition (ẑ, ŵ) is obtained with a MAP rule after the last V-Bayes step.

Penalised information criterion

BIC is an information criterion defined as an asymptotic approximation of the logarithm of the integrated likelihood [START_REF] Schwarz | Estimating the dimension of a model[END_REF]. The standard case leads to write BIC as a penalised maximum likelihood:

BIC = max θ log(p(y; θ)) - D 2 log(n),
where n is the number of statistical units and D the number of free parameters. Unfortunately, this approximation cannot be used for LBM, due to the dependency structure of the observations y.

However, a heuristic can be stated to define BIC. BIC-like approximations of each term of ICL (Equation 3) lead to the following approximation as n and d tend to infinity:

ICL(g, m) max θ log p(y, z, w; θ) - g -1 2 log n - m -1 2 log d - gm(r -1) 2 log(nd). (4)
Now, following [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF], the parameter maximizing the completed likelihood log p(y, z, w; θ) is replaced by the maximum likelihood estimator θ:

as a sum of an entropy term and BIC is classical [START_REF] Mclachlan | Finite Mixture Models[END_REF], 6.10.3) and, by analogy, we conjecture from (4) and ( 5) the following form for BIC after a straightforward factorisation of the penalty:

BIC(g, m) = log p(y; θ) - gm(r -1) + g -1 2 log n - gm(r -1) + m -1 2 log d. (6)
Notice that the maximum likelihood is not available and is replaced by a lower bound computed by variational approximation. This expression of BIC(g, m) highlights two penalty terms, one for the columns and one for the rows: the statistical units can be seen to be the rows and columns respectively. The number of row (resp. column) free parameters takes into account the proportions π (resp. ρ) and the parameter α of the categorical conditional distribution. BIC is known to be a consistent estimator of the order of mixture models with bounded log likelihoods [START_REF] Keribin | Consistent estimation of the order of mixture models[END_REF], whereas ICL is not [START_REF] Baudry | Sélection de modèle pour la classification non supervisée[END_REF]. However ICL can also be expected to be consistent for estimating the LBM order due to a specific behaviour of the conditional distribution of the labels p(z, w|y; θ). Indeed, [START_REF] Mariadassou | Convergence of the groups posterior distribution in latent or stochastic block models[END_REF] proved Inria that for a categorical LBM, this distribution concentrates on the actual configuration with large probability as soon as θ is such that α is not too far from the true value. As a consequence, relying on the MAP rule and with an estimator θ such that α is converging to the true parameter value, log p( z, w|y; θ) vanishes for true g and m. Moreover, this non positive term may add to BIC a supplementary penalty for incorrect g or m and hence leads ICL to be also consistent to estimate the order. Notice that this conjecture needs to assume the consistency of the maximum likelihood and the variational estimators for LBM. These results have not been proved yet for LBM as far as we know, and still remain an interesting challenge. But this conjecture is supported by the numerical experiments of the next section where ICL outperforms BIC for low n and d, and has the same behaviour as BIC for greater values.

Numerical experiments

Relevant numerical experiments on both simulated and real data sets supporting the claims of this paper are now presented. First we analyse the ability of the Bayesian inference through Gibbs sampling to avoid the tendency of the maximum likelihood methodology through VEM or SEM-Gibbs algorithms to provide empty clusters. After a simple illustration, Monte Carlo experiments on simulated data are performed (i) to compare the behaviour of SEM-Gibbs+VEM and Gibbs+V-Bayes algorithms (ii) to analyse the influence of the hyperparameters a and b (iii) to enlighten the behaviour of model selection criteria. Then, the LBM is experimented on a real categorical data set already used by [START_REF] Wyse | Block clustering with collapsed latent block models[END_REF].

Escaping from an empty cluster solution with V-Bayes

This illustration is done with an easily separated data matrix produced by [START_REF] Lomet | Un protocole de simulation de données pour la classification croisée[END_REF] and [START_REF] Lomet | Sélection de modèle pour la classification croisée de données continues[END_REF] for the binary LBM model 1 , with 50 rows, 50 columns, (g, m) = (3, 3) clusters and equal mixing proportions. The left part of Figure 2 displays the three components of π along the iterations of the VEM algorithm while the right part of Figure 2 displays them for the V-Bayes algorithm with (a, b) = (4, 1). Both algorithms started from the same position. VEM is rapidly trapped into an empty cluster solution. But V-Bayes, for which the estimated proportions are not smaller than a-1 g(a-1) , escapes from this empty cluster solution and finally provides a satisfactory solution with the real number of row clusters g = 3.

Experiments with simulated data

Monte Carlo experiments were performed to assess the ability of Bayesian inference to avoid empty cluster solutions. A binary LBM with (g, m) = (5, 4) clusters was considered with following Bernoulli parameters

α =       ε ε ε ε 1 -ε ε ε ε 1 -ε 1 -ε ε ε 1 -ε 1 -ε 1 -ε ε 1 -ε 1 -ε 1 -ε 1 -ε      
where ε allows to produce easily (+), moderately (++) or hardly (+++) separated mixtures, see [START_REF] Govaert | Block clustering with bernoulli mixture models: Comparison of different approaches[END_REF], and with unequal ρ = (0.1 0.2 0.3 0.4) and π = (0.1 0.15 0.2 0.25 0.3) 1 The data set is available from https://www.hds.utc.fr/coclustering/doku.php. Comparing SEM-Gibbs+VEM and Gibbs+V-Bayes to avoid empty cluster solutions. For each case of separation, and each sample size (n, d) = (100, 200), (150, 150), (200,100), 500 data matrices were simulated. Both algorithms were initialized with two different numbers of clusters:

RR n°8264

(1) the true number of clusters: (g, m) = (5, 4),

(2) more than the true number of clusters: (g, m) = (8, 8).

The hyperparameters a and b for Bayesian inference were chosen in {1, 4, 16} 2 \{(1, 1)}. The results are summarised in Tables 1-2 for unequal mixing proportions. Case (a, b) = (1, 1) corresponds to the SEM-Gibbs+VEM algorithm and can be compared with Gibbs+V-Bayes algorithm on higher a and b values. It clearly appears that Bayesian inference with a = 4 or 16 and b = 1 is doing a good job to avoid empty cluster solutions. As expected, taking a > 1 is quite relevant in that purpose. On the contrary, it appears clearly that taking b > 1 is harmful. As noticed by a reviewer, one possible reason for this could be that taking b > 1 means that the prior weight in the Dirichlet is more focused on the centre of the simplex i.e. all category weights are equal. In this case, b > 1 puts most prior weight on equal proportions for success and failure in the Bernoulli distribution. Such heavy prior weight could penalize against separation of the block clusters i.e. this puts more prior weight on saying there is no difference in the block probabilities (i.e. all α kl close to 0.5), and so there could be more of a tendency for empty clusters and hence overall homogeneity in the solution. Table 2: Percentage of empty cluster solutions, for data simulated with unequal proportion models, obtained by SEM-Gibbs+VEM (cell a = 1, b = 1) and Gibbs+V-Bayes algorithms, called with (g, m) = (8, 8), for 500 matrices in the easily separated +, moderately separated ++ and hardly separated + + + cases, for three sample sizes.

Analysing the empty cluster solutions of SEM-Gibbs+ VEM. The fact that the SEM-Gibbs+VEM algorithm produces empty clusters can be thought of as beneficial: it could mean that the number of clusters has been oversized and the algorithm provides a more relevant number of clusters. To answer to this question, SEM-Gibbs+VEM was initialised with (g, m) = (8, 8) and run on 500 data sets simulated under the same conditions than previously described. Table 3 reports the number of cluster finally obtained in the hardly separated cases with n = 150, d = 150 and unequal mixing proportions. Results for other cases are available in [START_REF] Keribin | Estimation and Selection for the Latent Block Model on Categorical Data[END_REF]. In all situations, the right number of clusters is chosen less than two percent of times. Moreover, the SEM-Gibbs+VEM algorithm produces empty clusters even in case of initial 4: Number of provided clusters by SEM-Gibbs+VEM for 500 random initialisations in the separated case for g = 8 and m = 8 whereas the right repartition is g = 10 and m = 10.

Analysing the behaviour of the model selection criteria. For each g and m varying from two to eight, the Gibbs+V-Bayes algorithm was run on 50 simulated data sets with number of clusters (g, m) = (5, 4) and BIC and ICL criteria were computed. Tables 5 and6 display how many times each (g, m) was selected in the equal mixture proportion case. The ICL criterion seems to provide a better estimation of the number of clusters than BIC, especially for the smaller sizes of the data matrix. Otherwise, both criteria tend to choose less clusters than there actually are in the cases where the true number of clusters is not correctly identified. This tendency appears also on numerical experiments not reported here but available in [START_REF] Keribin | Estimation and Selection for the Latent Block Model on Categorical Data[END_REF] for mixtures with the same α parameters, but unequal mixing proportions. It deserves the following comments. First, it is worth noting that the tendency of ICL to choose fewer components than the true number of components of a mixture is well-known. The purpose of ICL is to assess the number of mixture components that leads to the best clustering and it can happen that the number of clusters is smaller than the number of mixture components. This is because a cluster may be better represented by a mixture of components than by a single component (see [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF] or [START_REF] Baudry | Combining Mixture Components for Clustering[END_REF]). For the LBM as mentioned in Section 4.2, when computing BIC, the maximum likelihood is approximated by the maximum free energy. The tendency of BIC to underestimate the true number of components could indicate that the difference between the maximum likelihood and the maximum free energy increases with the number of mixture components. As expected, for large sizes, the performances of both criteria improve and the difference between ICL and BIC decreases. This is in agreement with the conjecture presented Table 5: The frequency of the models selected by the ICL criterion in the easily separated +, moderately separated ++ and ill separated + + + cases with equal proportion for different matrix sizes.

Congressional voting in US senate

The combination Gibbs+V -Bayes was tested on UCI Congressional Voting Records data set2 recording the votes ('yes', 'no', 'abstained or absent') of 435 members of the 98th congress on 16 different key issues. This data set involves homogeneous 3-level categorical data. Table 6: The frequency of the models selected by the BIC criterion in the easily separated +, moderately separated ++ and ill separated + + + cases with equal proportion for different matrix sizes.

The Gibbs+V-Bayes algorithm was run on the 3-level data with (a, b) = (4, 1) and the ICL and BIC criteria were computed to select the number of clusters. ICL selected (g, m) = (5, 7) clusters, while BIC selected (g, m) = (4, 6). The rows and columns reorganisation derived from the block clustering for each criterion is displayed in Figure 3 where 'yes' is coloured in black, 'no' in white and 'abstained or absent' in grey. Row cluster composition by Party is shown in Table 7. Roughly speaking, the first three clusters selected with both criteria are the same. Row cluster 1 is mainly Republican while row cluster 3 is mainly Democrat. Row cluster 2 is mixed. Cluster 4 selected with BIC is split into two clusters (4 and 5) in the model selected with ICL (the six bigger abstainers are isolated in cluster 5). The first four column clusters selected with both criteria are identical. Column clusters 1 and 2 contain the issues on which Republicans voted 'yes' and Democrats 'no'. These two column clusters only differ for the vote of congress members belonging to row cluster 2. In contrast, column clusters 3 and 4 contain the issues on which Republicans voted 'no' and Democrats 'yes'. Moreover, both criteria isolated in a specific cluster the issue characterised by a high level of abstainers. Column cluster 5 selected with BIC is split into two column clusters (5 and 6) in the model selected with ICL, due to one issue more subject to abstention.

This data set was already studied by [START_REF] Wyse | Block clustering with collapsed latent block models[END_REF] where the 'no' and 'abstained or absent' were aggregated in order to work on binary data. Assuming that g and m were random, Wyse and Friel used a collapsed sampler which marginalises over the model parameters with uniform prior (i.e. (a, b) = (1, 1) in Figure 1) aiming to compute the distribution of the most visited models and the maximum a posteriori cluster membership (z, w). This maximisation is analogous to maximising the ICL criterion because : log p(z, w, g, m|y) = ICL(z, w, g, m)

+ log p(g, m) -log p(y).

(7)

Since the prior on g and m was taken to be uniform, the two algorithms should give the same optimal quadruplet (z, w, g, m). Wyse and Friel gave the a posteriori distribution for (g, m) but did not provide the maximum a posteriori cluster membership (z, w). We computed it by running their collapsed sampler with 110 000 iterations, 10 000 of them as burn-in iterations3 . This led to (g, m) = (6, 11) clusters (which is not in their 60% confidence level) with a corresponding ICL = -3565. Replicating 30 times gave variations from 10 to 14 column clusters, with 6 row clusters each time, and corresponding ICLs were computed. The best ICL value was -3554 for (g, m) = (6, 13). Running Gibbs+V-Bayes on the same binary data set and the same priors selected (g, m) = (5, 13) clusters with corresponding ICL= -3553. Notice that ICL(g, m) is quite flat around its maximum. The results are similar, requiring six times more iterations for the collapsed sampler.

Discussion

In this paper, we show that Bayesian inference through Gibbs sampling is beneficial to produce regularised estimates of the LBM for categorical data. A sufficient solution to ensure the identifiability of the LBM for categorical data has been provided. Contrary to the VEM or V-Bayes algorithms, the SEM-Gibbs algorithm and Gibbs sampling provide solutions not highly dependent of the starting values. Restricting attention to point estimation, Bayesian inference through Gibbs sampling from non informative priors produces regularised parameter estimation.

The solution providing the posterior mode of the parameters is not suffering the label switching problem and can be used as initial values of the V-Bayes algorithm, and thus avoid the empty cluster solutions which jeopardizes maximum likelihood inference for the LBM. Moreover, taking profit of the exhibited sufficient condition ensuring the identifiability of the LBM for categorical data, it is possible to define a natural order to post process the rows and columns after Gibbs sampling in order to deal with the label switching problem in a proper way. On the other hand, we select a latent block model using the integrated completed likelihood ICL which can be computed without requiring asymptotic approximations. A good perspective for future work would be to prove the conjecture stated in Section 4.2: namely ICL criterion provides a consistent estimation of the number of clusters (g, m) for the LBM, contrary to the standard mixture model for which ICL is not consistent. Otherwise, an alternative is to use the collapsed sampler of [START_REF] Wyse | Block clustering with collapsed latent block models[END_REF]. Maximising the posterior collapsed is closely related to maximising ICL for flat (g, m) priors (see equation 7). In our experiments, ICL appears to be efficient to find the (z, w, g, m) maximising the posterior log p(z, w, g, m|y) and less computationally demanding than the collapsed sampler. But, the latter is able to explore the uncertainty in the joint support of all models and corresponding missing labels. A perspective would be to analyse the variability of ICL with a random sampling of (z, w) according to the distribution p(z, w; θ). Preliminary numerical experiments on simulated data show that the variability of ICL depends on (g, m) and is clearly more important for g or m greater than the real values. It could be an interesting subject for further research.

where D k = detM k and D g = 0 as M g is a product of invertible matrices. Hence, it is possible to determine τ , and R is now known. Consequently, π is defined in a unique manner by

A π = R -1 M g R -1 .
In the same way, ρ is defined in a unique manner by considering the probabilities σ to have 1 in a column of column class and the probabilities v j to have 1 on the j first cells of the first column of y σ = IP(y ij = 1|w j = 1) and v j = IP(y 11 = 1, . . . , y j1 = 1)

To determine α, we introduce for 1 ≤ i ≤ g and 1 ≤ j ≤ m the probabilities U ij to have 1 in each i first cells of the first row and in each j first cells of the first column

U ij = IP(y 11 = 1, . . . , y 1i = 1, y 21 = 1, . . . , y j1 = 1) = k, 1 , 2 ,..., i , k1,...,kj IP(y 11 |z 1k = 1, w 1 = 1) × IP(z 1k = 1)IP(w 1 = 1) × i λ=2 IP(y 1λ |z 1k = 1, w λ λ = 1)IP(w λ λ = 1) × j η=2 IP(y η1 |z ηkη = 1, w 1 = 1)IP(z ηkη = 1) = k, IP(y 11 |z 1k = 1, w 1 = 1)IP(z 1k = 1) × IP(w 1 = 1)IP(y 1i |z 1k = 1) i-1 IP(y j1 |w 1 = 1) j-1 = k, π k τ i-1 k α k ρ σ j-1 .
These probabilities are known and we can write, with S j = (σ ) j-1 , j = 1, . . . m, = 1, . . . , m U = RA π AA ρ S , and it leads to define A = A -1 π R -1 U S -1 A -1 ρ , and hence α, in a unique manner. The proof is straightforwardly extended to the categorical case. The identification of π and ρ are obtained by considering the probabilities of y ij = 1 where 1 is the first outcome of the multinomial distribution. Then, each A h = (α h k ) k=1,...,g; =1,...,m is successively identified by considering y ij = , i = 1, . . . , m, j = 1, . . . , g.

B VEM algorithm for categorical data [START_REF] Govaert | Block clustering with bernoulli mixture models: Comparison of different approaches[END_REF] proposed to use, for the binary case, a variational approximation of the EM algorithm by imposing that the joint distribution of the labels takes the form q

(c) zw (z, w) = q (c) z (z)q (c)
w (w). To get simpler formulas, we will denote

s (c) ik = q (c) z (z ik = 1) and t (c) j = q (c)
w (w j = 1).

Using the variational approximation, the maximisation of the loglikelihood is replaced by the maximisation of the free energy F(q zw , θ) = IE qzqw log p(y, z, w; θ) q z (z)q w (w) alternatively in q z , q w and θ (see [START_REF] Keribin | Méthodes bayésiennes variationnelles: concepts et applications en neuroimagerie[END_REF]. The difference between the maximum loglikelihood and the maximum free energy is the Kullback divergence KL(q zw ||p(z, w|y; θ)) = IE qzw log qzw(z,w) p(z,w|y;θ) . This can be extended to the categorical LBM, and leads to the following Variational EM (VEM) algorithm: E step It consists of maximising the free energy in q z and q w , and it leads to: . Note that the formulae of VEM and SEM-Gibbs are essentially the same, except that the probabilities s ik and t j are replaced with binary indicator values z ik and w j .

D Computing ICL

In this section, the exact ICL expression is derived for categorical data. Using the conditional independence of the zs and the ws conditionally to α, π, ρ, the integrated completed likelihood can be written

We now turn to the computation of p(y|z, w). Since α and (z, w) are independent, we have p(α|y, z, w) = p(y|z, w, α)p(α) p(y|z, w, α)p(α)dα , and, using the conditional independence of y ij knowing z, w and α, p(y|z, w, α) = i,j,k, h

(α h k ) y h ij z ik w j = k, h (α h k ) i,j z ik w j y h ij = k, h (α h k ) N h k . Inria Therefore p(α|y, z, w) ∝ k, p(α k ) h (α h k ) N h k ∝ k, h (α h k ) b-1 h (α h k ) N h k ∝ k, h (α h k ) b+N h k -1 .
And, the ICL criterion, presented in Section 4.1, is straightforwardly derived from equations ( 3), ( 8), ( 9) and ( 10).

Figure 1 :

 1 Figure 1: A graphical representation of the Bayesian latent block model.

)

  Using now the conjugate properties of the prior Dirichlet distributions and the conditional independence of the y ij knowing the latent vectors z and w and the LBM parameters, we get (see Appendix D for a detailed proof) ICL(g, m) = log Γ(ga) + log Γ(ma) -(m + g) log Γ(a) +mg(log Γ(rb) -r log Γ(b))

  y, z, w; θ) log p(y, z, w; θ) log p( z, w|y; θ) + log p(y; θ).Therefore ICL is the sum of the logarithm of the conditional distribution log p( z, w|y; θ), measuring the assignment confidence, and a penalised likelihood. The decomposition of ICL ICL(g, m) = log p( z, w|y; θ) + BIC(g, m)

Figure 2 :

 2 Figure 2: Evolution of π for the VEM algorithm (on left) and the V-Bayes algorithm (on right)

Figure 3 :

 3 Figure 3: Reorganisation of the data by the ICL criterion (on left) and by the BIC criterion (on right).

  = (α k1 , . . . , α km ), andψ k (y i• ; α . . . , mwhere y •j denotes the column j of the matrix y, α • = (α 1 , . . . , α g ) and φ (y •j ; α

  of p(z|y, w(c) ; θ (c) ), then simulation of z (c+1) according to p(z|y, w (c) ; θ (c) ):p(z i = k|y i• , w (c) ; θ (c) ) =π of p(w|y, z (c+1) ; θ (c) ), then simulation of w (c+1) according to p(w|y, z (c+1) ; θ (c) ).M step Denoting z .k := i z ik and w . := j w j , it leads to π

  p(y, z, w) = p(y, z, w, α, π, ρ)p(α)p(π)p(ρ)dαdπdρ = p(y|z, w, α, π, ρ)p(z, w|α, π, ρ) × p(α)p(π)p(ρ)dαdπdρ = p(y|z, w, α)p(α)dα × p(z|π)p(π)dπ p(w|ρ)p(ρ)dρ = p(z)p(w)p(y|z, w).RR n°8264Thus ICL = log p(z) + log p(w) + log p(y|z, w).Since the prior distribution of π is the Dirichlet distribution D(a, . . . , a) .1 + a, . . . , z .g + a) with the normalising factorΓ(n + ga) k Γ(z .k + a).The expression of p(z) directly follows from the Bayes theorem p. + a) Γ(d + am) .

  , . . . , b). Each term k , is the density of a non-normalised Dirichlet distribution D(b + N 1 k , . . . , b + N r k ) with the normalising factor Γ(z .k w . + rb) .k w . + rb) .

Table 1 :

 1 Percentage of empty cluster solutions, for data simulated with unequal proportion models, obtained by SEM-Gibbs+VEM (cell a = 1, b = 1) and Gibbs+V-Bayes algorithms, called with the true (g, m), for 500 matrices in a easily separated +, moderately separated ++ and hardly separated + + + cases, for three sample sizes.

	Inria

Table 3 :

 3 Number of provided clusters, for data simulated with the unequal proportion models, obtained by SEM-Gibbs+VEM for 500 matrices in the hardly separated + + + case with n = 150 and d = 150, starting with g = 8 and m = 8. numbers of clusters less than the real numbers. This algorithm was started with (g, m) = (8, 8) on 50 simulated binary matrix of sizes n = 50 and d = 50 of hardly separated (g, m) = (10, 10) clusters. Although the numbers of requested clusters are smaller than the right values, the algorithm produces empty cluster solutions almost twenty percent of times (Table4).

	H H g	H 4	H m 3 H 0	4 2	5 0	6 0	7 0	8 0
		5	1	4	4	3	2	0
		6	0	7	13 15 6	8
		7	0	3	12 26 23 21
		8	0	3	13 35 85 214
			H g H H 7	H m 6 7 H 1 2	8 18
			8		5 64 410
	Table

Table 7 :

 7 Repartition of Democrats and Republicans for each criterion.

	Inria

RR n°8264

Congressional Voting Records data set is available from http://archive.ics.uci.edu/ml/datasets/Congressional+Voting+Records. RR n°8264

The code of Wyse and Friel is available at http://www.ucd.ie/statdept/jwyse. Inria

rue Honoré d'Estienne d'Orves Bâtiment Alan Turing Campus de l'École Polytechnique 91120 Palaiseau Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgements The authors thank the reviewers and the Associate Editor for their very helpful comments which have greatly improved this paper. C. K. has been supported by LMH (Hadamard Mathematics Labex), backed by the foundation FMJH (Fondation Mathématique Jacques Hadamard).

Appendices

A Proof of Theorem 1

The identifiability of a parametric model requires that for any two different values θ = θ in the parameter space, the corresponding probability distribution IP θ and IP θ are different. We prove that, under assumptions of Theorem 1, there exists a unique parameter θ = (π, ρ, α), up to a permutation of row and column labels, corresponding to IP(y) the probability distribution function of matrix y having at least 2m -1 rows and 2g -1 columns. The proof is adapted from [START_REF] Celisse | Consistency of maximum-likelihood and variational estimators in the stochastic block model[END_REF] who set up a similar result for the Stochastic Block Model. Notice first that τ = Aρ is the vector of the probability τ k to have 1 in a cell of a row of given row class k:

As all the coordinates of τ are distinct, R defined as the g × g matrix of coefficients R ik = (τ k ) i , for 0 ≤ i < g and 1 ≤ k ≤ g, is Vandermonde, and hence invertible. Consider now u i , the probability to have 1 on the i first cells of the first row of y:

With a given IP(y), u 1 , . . . , u 2g-1 are known, and we denote u 0 = 1. Let now M be the (g + 1) × g matrix defined by M ij = u i+j-2 , for all 1 ≤ i ≤ g + 1 and 1 ≤ j ≤ g and define M i as the square matrix obtained by removing the row i from M . The coefficients of M are

We can write, with A π = Diag(π)

Now, R, unknown at this stage, can be retrieved by noticing that the coefficients of τ are the roots of the following polynomial [START_REF] Celisse | Consistency of maximum-likelihood and variational estimators in the stochastic block model[END_REF] of degree g B(x) = g k=0 (-1) k+g D k x k , Inria


Gwenaël Guillaume
email: guillaume.gwenael@neuf.fr

Pierre Aumond

Benoit Gauvreau

Guillaume Dutilleux

Application of the Transmission Line Matrix method for outdoor sound propagation modelling -Part 1: Model presentation and evaluation

The present paper deals with an original time-domain approach applied to outdoor sound propagation under meteorological effects. The Transmission Line Matrix method, based on the Huygens' principle, had already been validated over impedant grounds and complex topography. The presented formulation proposes to take into account meteorological effects (wind speed and temperature) through the relative sound speed. The necessary wavefront direction is determined through the calculation of the averaged intensity vector direction. A good agreement is found between simulations of both the transmission line matrix and parabolic equation methods. A relevant use of the method is shown in the framework of environmental acoustics and initial applications are proposed in Part 2.

Time-domain methods, Transmission Line Matrix model, Meteorological effects, Outdoor sound propagation

Introduction

Many numerical models have been developed in order to simulate meteorological effects on outdoor sound propagation. One can cite analytical formulations such as geometrical methods (e.g. ray or gaussian beam tracing approaches) [START_REF] Ostashev | Acoustics in Moving Inhomogeneous Media, E & FN SPON[END_REF][START_REF] Salomons | Computational Atmospheric Acoustics[END_REF] or numerical models based on the resolution of the wave equation. Among these latter models, the Parabolic-Equation (PE) based method has been widely used for this purpose [START_REF] Salomons | Computational Atmospheric Acoustics[END_REF][START_REF] Lihoreau | Outdoor sound propagation modeling in realistic environments: Application of coupled parabolic and atmospheric models[END_REF][START_REF] Van Renterghem | Comparison of measurements and predictions of sound propagation in a valley-slope configuration in an inhomogeneous atmosphere[END_REF][START_REF] Cotté | Estimates of the relevant turbulent scales for acoustic propagation in an upward refracting atmosphere[END_REF]. Over the last decade, with the increasing power of computational resources, time-domain methods have also been developed and applied successfully in environmental acoustics [START_REF] Chopard | Lattice boltzmann computations and applications to physics[END_REF][START_REF] Blumrich | A linearized eulerian sound propagation model for studies of complex meteorological effects[END_REF][START_REF] Ostashev | Equations for finite-difference, time-domain simulation of sound propagation in moving inhomogeneous media and numerical implementation[END_REF][START_REF] Hornikx | The extended fourier pseudospectral time-domain method for atmospheric sound propagation[END_REF]. The most popular time-domain approach is undoubtedly the finite-difference in the time-domain (FDTD) method. Dragna has investigated for instance sound propagation over a 100 m distance within a realistic context in a frequency range between 100 Hz and 2 kHz [START_REF] Dragna | Modélisation par une approche temporelle de la propagation acoustique en milieu extérieur : traitement de frontières complexes et validation sur site ferroviaire[END_REF]. An alternative time-domain method is the Transmission Line Matrix (TLM) approach [START_REF] Kagawa | Discrete Huygens' model approach to sound wave propagation[END_REF]. This model seems very promising for describing complex outdoor sound propagation yet has not been used extensively. In order to incorporate the atmospheric effects, Hofmann has proposed an interesting formulation, which is however limited to temperature effects [START_REF] Hofmann | Simulation of outdoor sound propagation with a transmission line matrix method[END_REF]. A first attempt to define a TLM scheme with unidirectional mean flow (e.g. wind field) into the TLM grid has been provided by Kagawa and his co-workers, but its method is only successful for single wind speed direction effects [START_REF] Kagawa | Discrete Huygens' modelling simulation of sound wave propagation in velocity varying environments[END_REF]. In this paper, the approach, inspired by a method proposed by Dutilleux [START_REF] Dutilleux | Applicability of TLM to wind turbine noise prediction[END_REF], consists of mod-ifying the sound speed at each point of the discretized domain as a function of the temperature and wind speed projection on the wave front direction. This paper aims to present the integration of meteorological effects in the TLM model and then its comparison with a PE model in academic cases.

The associated paper (part 2) addresses its validation by comparison with results stemming from Lannemezan's 2005 experimental campaign [START_REF] Junker | Meteorological classification for environmental acoustics -Practical implications due to experimental accuracy and uncertainty[END_REF]. The wind and temperature fields are obtained from the meso-scale meteorological model Meso-NH [START_REF] Lafore | The Meso-nH atmospheric simulation system. Part I: Adiabatic formulation and control simulations[END_REF][START_REF] Aumond | Including drag effects of canopies: Real case large-eddy simulation studies[END_REF].

The first section presents the TLM model. In section 2 the formulation taking into account wind speed and temperature effects is described. The last section proposes evaluating the ability of the TLM to treat outdoor sound propagation problems through a comparison with PE results in two academic cases.

TLM modelling

The TLM method is based on the Huygens' principle, which states that a wavefront consists of a set of secondary sources radiating spherical wavelets whose envelopes can be broken down into a new generation of secondary sources as well. Hofmann has shown the equivalence of this approach with the resolution of the discretized wave equation [START_REF] Hofmann | Simulation of outdoor sound propagation with a transmission line matrix method[END_REF]. Other authors have also derived the two-dimensional homogeneous cartesian formulation of TLM from a Lattice Boltzmann model by removing nonlinear terms, in choosing a suitable viscosity and selecting a square grid [START_REF] Martin | Low Cost Parallelizing: A Way to be Efficient[END_REF].

The TLM statement allows describing sound propagation through both a spatial and temporal discretization of the medium as well as the propa-gation phenomena. This concept is numerically conveyed by replacing the propagation medium with a transmission lines network, through which sound propagates in the form of sound pulses. Thereby, as shown in Fig. 1, each junction, or node, links N = 4 or N = 6 transmission lines to each other in two dimensions (2D) or three dimensions (3D) respectively. Thereafter, the number of dimensions is called d such that d = 2 and d = 3 for creating a 2D and 3D model respectively. An additive branch, of index N + 1, is inserted at each node of the transmission line network, in order to consider the inhomogeneities of the propagation medium (i.e. branch 7 in Fig. 1). According to the TLM concept, sound propagates in the form of pulses. Thus, incident and scattered pulses are considered at each transmission lines junction and time increment. The propagation medium is discretized by means of a uniform cartesian meshing of mesh size (∆l) d , with ∆l being the spatial step such that:

∆l λ √ d 10 , (1)
with λ the minimal wavelength of the simulation.

The scattered pulses at time increment t and node of discrete coordinates

r such that r =    (i, j) for d = 2 (i.e. in 2D), (i, j, k) for d = 3 (i.e. in 3D), (2)
are related with the incident pulses at this node at the same time iteration by the following matrix relation:

t S r = t D r × t I r , (3)
where t I r and t S r are the vectors composed of the incident pulses t I n r and scattered pulses t S n r through each transmission line n (n = 1 to N + 1) respectively. t D r is a (N + 1) × (N + 1) scattering matrix given by:

1 2 3 4 5 6 7 t S 1 (i ,j ,k ) t S 2 (i ,j ,k ) t S 3 ( i , j , k ) t S 4 ( i , j , k ) t S 5 (i,j,k) t S 6 (i,j,k) t S 7 ( i , j , k ) 1 2 3 4 5 6 7 t I 1 (i ,j ,k ) t I 2 (i ,j ,k ) t I 3 ( i , j , k ) t I 4 ( i , j , k ) t I 5 (i,j,k) t I 6 (i,j,k) t I 7 ( i , j , k ) x y z
t D r = 2 t η r + 2d            t a r 1 . 1 t η r 1 t a r . 1 t η r . . . . . 1 1 . t a r η r 1 1 . 1 t b r            , (4)
with

t a r = - t η r + 2 (d -1) 2 (5)
and

t b r = t η r -2d 2 , ( 6
)
where t η r allows locally modifying, and if needed during the simulation, the sound speed in the propagation medium. In other words, this term is used in order to model an inhomogeneous atmosphere and is defined in section 3.1.

Finally, the nodal pressure is written as a combination of all incident pulses, i.e.:

t p r = 2 t η r + 2d N -1 n=1 t I n r + t η r t I N r . (7)
In addition, the scattered pulses from nodes adjacent to node (i, j, k) at time increment t become the incident pulses to this node at the next time iteration t + ∆t, with ∆t representing the time step defined by:

∆t = ∆l √ d c 0 , (8)
with c 0 the adiabatic sound speed. This diffusion process is governed by connection laws depicted in Fig. 2 such as:

t+∆t I n r = t S m r ± n (9)
and

t+∆t I N r = t S N r , (10)
with

   n r ± n    =                         m -1 r -    if m is even,    m + 1 r +    if m is odd, (11)
and

r ± n =          (i ± 1, j, k) , for n = 1 or 2,
(i, j ± 1, k) , for n = 3 or 4, (i, j, k ± 1) , for n = 5 or 6.

t+ ∆ t I 1 (i ,j ,k ) t+ ∆ t I 2 (i ,j ,k ) t + ∆ t I 3 ( i , j , k ) t + ∆ t I 4 ( i , j , k ) t+∆t I 5 (i,j,k) t+∆t I 6 (i,j,k) t + ∆ t I 7 ( i , j , k ) t S 2 (i -1, j, k) t S 1 (i + 1, j, k) t S 4 ( i , j - 1 , k ) t S 3 ( i , j + 1 , k ) t S 6 (i,j,k-1) t S 5 (i,j,k+1) t S 7 ( i , j , k ) (i -1, j, k) (i + 1, j, k) (i, j -1, k) (i, j + 1, k) (i, j, k -1) (i, j, k + 1)
x y z Boundaries are implemented in the TLM model at a distance ∆l/2 from the nearest node in order to ensure the synchronism of sound pulses. They can be characterized by a reflection coefficient in pressure [START_REF] Kagawa | Discrete Huygens' model approach to sound wave propagation[END_REF] or by an impedance boundary condition [START_REF] Guillaume | Time-domain impedance formulation for transmission line matrix modelling of out-door sound propagation[END_REF]. In addition, absorbing layers can be introduced in order to model an unbounded propagation medium as depicted in Ref. [START_REF] Guillaume | A simple absorbing layer implementation for transmission line matrix modeling[END_REF].

TLM formulation for acoustic propagation in a meteorological field

Regarding outdoor sound propagation, the temperature and wind fields combine to produce local variations of sound speed. The literature review above has shown that the implementation of thermal gradients is straightforward in TLM since sound speed is defined at the local level, whereas a suitable implementation for wind gradients is still underdeveloped. In order to address the general case of outdoor sound propagation, both thermal and wind effects must be taken into account simultaneously.

In order to for wind speed gradients in the TLM model, the approach developed in this section satisfies the previous requirement since it is based on the so-called effective sound speed. This approach requires knowledge of the local direction of the wavefront. A proof of concept of this approach for the TLM method has been proposed by one of the authors in Ref. [START_REF] Dutilleux | Applicability of TLM to wind turbine noise prediction[END_REF].

Effective sound speed

The effective sound speed is expressed as follows [START_REF] Salomons | Computational Atmospheric Acoustics[END_REF]:

t c effr = γ R t T r + t u r • t n r ( 13
)
where γ = 1.41 and R = 287 J.kg -1 .K -1 , representing the specific heat ratio and the perfect gas constant for air respectively. t T r is the temperature, t u r is the wind speed vector and t n r is the wavefront direction. The sound speed at a point r of the TLM grid verifies [START_REF] Guillaume | Application of the TLM method for the sound propagation modeling in urban area (Application de la méthode TLM à la modélisation de la propagation acoustique en milieu urbain)[END_REF]:

t c TLMr = c 0 2 t η r + 2d (14)
So, if one sets t c effr = t c TLMr , solving for the inhomogeneity parameter η leads to:

t η r = 2 √ γ R t T r + t u r • t n ∆l 2 ∆t 2 -2d. ( 15
)
In order to ensure that η 0, the adiabatic sound speed c 0 must verify:

c 0 γ R T max + u max ( 16
)
where the subscript max is computed over the whole simulated domain and simulation duration.

Wavefront direction

Vector n defines the direction of the wavefront. This information is not readily available in the TLM approach, yet it can be approximated by a quantity relating the instantaneous sound intensity vector:

t I r = t p r • t v r ( 17
)
where t v r is the local instantaneous velocity defined in [START_REF] Kagawa | Discrete Huygens' model approach to sound wave propagation -Reverberation in a room. Sound source identification and tomography in time reversal[END_REF]:

t v r = 1 ρ 0 c 0 ( t S m r -t S n r ) ( 18
)
where m is even (see Eq. 11).

A purely local definition of the sound intensity is not suitable as an approximation for n because it is influenced by the reactive intensity which does not propagate. In order to eliminate the reactive intensity, a possible approach is to average intensity over space, for instance over a onewavelength-size domain. As a first approximation and since meteorological effects (e.g. wind-induced refraction) are a long-range phenomenon, this averaging approach is only valid far away from the source. Thus, the average intensity vector corresponds with the active intensity vector and is as follows:

t I avg,r = 1 (2a + 1) d rn+a rn-a t I rn , (19)
where a sets the span of the average around r n .

From Eq. ( 19), one defines:

t n r = I avg,r I avg,r . (20)
The angle θ between the wavefront and the wind directions can be defined in 2D as:

t θ r = t u r • t n r || t u r | | × || t n r | | . (21)
Fig. 3 illustrates the wave direction and the angle θ obtained without averaging (i.e. the vector t n r corresponds with the normalized instantaneous intensity vector) and by averaging the variables with a spatial span a = 10.

The reactive intensity contribution appears clearly when no averaging is applied, i.e. where the sign of the sound pressure changes. Spatial averaging significantly improves the determination of the wavefront propagation direction.

Comparison with PE results

In order to evaluate the ability of the proposed formulation to take into account meteorological effects (through wind and temperature fields) in the TLM model, 2D simulation results are compared with results obtained with the Parabolic Equation (PE) model [START_REF] Lihoreau | Outdoor sound propagation modeling in realistic environments: Application of coupled parabolic and atmospheric models[END_REF], which is widely used for this purpose.

Three simulation results are presented. The source is located alternatively at 0.15 m, 2.05 m and 5.05 m from the ground. First, the sound attenuation levels are compared over a line of receivers located at a height of 2.05 m and at a fixed frequency of 200 Hz. Then, the attenuation spectra are given at 150 m from the source. The upper frequency bound is fixed in order to limit the computation cost.

The simulation presented in Fig. 4 is performed over a rigid ground through an atmosphere characterized by a linear vertical gradient of temperature (dT /dz = 0.35 • C m -1 ) and a wind field set to zero. The sound attenuations levels are very similar. The discrepancy in the first 50 m can be explained by the PE limit at short distances from the source [START_REF] Salomons | Computational Atmospheric Acoustics[END_REF]. The depth peak differences are due to the sensitivity and numerical uncertainties of the two models. This academic case allows us to validate the effective sound speed approach without any air flow (c effr = √ γ R T r ).

Since the meteorological effets tend to increase ground effects in the case of downward refraction conditions, the second simulation proposes to consider an impedant ground defined by the Miki impedance model [START_REF] Miki | Acoustical properties of porous materials -Modifications of Delany-Bazley models[END_REF] with an air flow resistivity σ = 100 kN.s.m -4 . In addition, the effective sound speed field is calculated for a homogeneous temperature field (T = 20 • C) and for a realistic vertical wind gradient of 0.2 s -1 . Thus, it allows evaluating the In the second part of this paper, TLM simulations are performed in realistic cases (with measured meteorological fields as input data) using the approach presented in this part. Thus, numerical sound predictions are compared with experimental data obtained under various meteorological conditions.

Figure 1 :

 1 Figure 1: Representation of the incident (left) and the scattered (right) pulses at node (i, j, k) in 3D.

Figure 2 :

 2 Figure 2: Representation of the connection laws through the main transmission lines (n = 1to N ) at the node (i, j, k) in 3D.

Figure 3 :

 3 Figure 3: Computation of the wave propagation direction and of the angle θ. Representation of (a) the wavefront, and of the wave propagation direction and of the angle θ obtained from (b)-(c) the instantaneous acoustic intensity and (d)-(e) the averaged acoustic intensity with a spatial span a = 10.

Figure 4 :

 4 Figure 4: Comparison of the TLM model simulations (lines) with the PE ones (dots) for the sound propagation over a rigid ground with a linear temperature gradient such as dT /dz = 0.35 • C m -1 . The source is located at 0.15 m (blue), 2.05 m (green) and 5.05 m (red) from the ground.

Figure 5 :

 5 Figure 5: Comparison of the TLM model simulations (lines) with the PE ones (dots) for the sound propagation over a rigid ground with a linear wind gradient such as dU/dz = 0.2 s -1 . The source is located at 0.15 m (blue), 2.05 m (green) and 5.05 m (red) from the ground.