For a moment, nothing happened. Then, after a second or so, nothing continued to happen." -Douglas Adams, The Hitchhiker's Guide to the Galaxy Voilà qui résume de manière assez fidèle les premiers instants du processus ayant abouti à ces remerciements. Ajoutez même quelques heures de néant à tout cela pour être au plus proche de la réalité. Rien de nouveau sous le soleil, mais l'exercice des remerciements est difficile. Difficile parce que personnel. Difficile parce que stylistique. Difficile parce que ce sont souvent les premières lignes sur lesquelles vont se poser les yeux du lecteur. Parfois même les seules. Ne mentez pas, ne rougissez pas, c'est quelque chose que tout le monde fait. À défaut d'être originaux, j'espère qu'ils transmettront le moins partiellement et erronément possible toute la reconnaissance que j'aurais souhaité qu'ils expriment.

Stéphane, merci. Merci de m'avoir ouvert les portes de cette aventure que rien ne prédestinait à être graphique et arborescente, mais qui l'a pourtant été. Elle fut également passionnante, souriante et profondément formatrice. Il m'est littéralement impossible de concevoir un environnement doctoral plus épanouissant et plus agréable que celui que tu as pu m'offrir. Un savant mélange de liberté et de conseils avisés, dispensé avec générosité et disponibilité. Tes qualités sont si nombreuses qu'il m'est difficile de leur faire justice ici. La perspective embarrassante d'une tentative pédestre et maladroite m'invite à en rester là, mais que ce manque d'audace ne soit qu'une preuve supplémentaire de mon admiration.

Quelques mots en anglais, puisque c'est outre-Manche que l'histoire a commencé. Thank you Michael, for taking me in as a Master student and for having started all this, in a way. The Imperial College has been a special place for me. Je remercie Marina Meilă et Étienne Birmelé d'avoir sacrifié une partie de leur été pour rapporter cette thèse. J'espère que la tâche ne fut pas trop pénible. Je remercie également Christophe Giraud et Steffen Lauritzen d'avoir accepté de faire partie de mon jury. Plus que d'un laboratoire, c'est d'une communauté dont j'ai fait partie durant mes années de thèse. Une communauté éminemment scientifique et profondément humaine. Merci à chacun d'entre vous. Pour ces pauses café, que l'éclectisme des conversations rend si vivantes et imprévisibles. Pour ces moments passés, au Vieux Chêne et ailleurs, à travailler sur une définition claire du bobo, du hipster et du parisien, pour finalement toujours tomber dans v les pièges de l'induction. Pour ces quiz, où la bonne humeur vole la vedette à la mauvaise foi et à la culture. Que la tradition perdure. Ces instants cristallisent pour moi un esprit que je ne peux qu'espérer retrouver ailleurs, aussi loin que les vents me portent. Vous êtes vraiment de belles personnes.

Un grand merci au Bureau des Pauses, lieu de vie (et de travail) où l'entropie est étrangement hétérogène et où les idées les plus déjantées ont une certaine tendance à ne pas rester abstraites. Anna et Marie, je suis assez fier de ce que l'on a accompli, et je ne suis pas malheureux d'avoir partagé ce bureau avec vous. J'eus été bien démuni, sans votre soutien, devant ces pages web où il faut confirmer et confirmer encore. Petit conseil pour la fin : commencez dès maintenant à écrire vos remerciements. Merci à Pierre G. de nous avoir apporté cette saine émulation et de nous rappeler quotidiennement que l'humilité est la plus haute des qualités. Je remercie également nos voisins et chaperons pour leur bienveillance permanente.

À toutes les personnes ayant contribué à cette thèse de par la qualité des moments où je n'y travaillais pas, merci. Et elles sont nombreuses. Elles font avancer la recherche. Elles rient, elles dansent, elles chantent. Elles foulent les planches. Elles habitent loin. Ou pas. Elles rugissent comme tout le monde. Elles sont d'une candeur entière. Elles lézardent au soleil. Elles engueulent la pauvre Suzanne. Elles adorent (perdre) les paris. Elles mangent du pâté de campagne. Elles vont toujours bien. Elles poussent des cris (très) aigus. Elles sont au moins aussi folles que moi. Elles ont de nombreux petits-enfants. Elle n'aiment pas être bousculées. Elles ont des perches à selfie. Elles ne savent pas faire un duck face. Elles naissent

A network is a system of interconnected entities. These entities can be of various natures. They might for instance share their feelings on the Internet along with pictures of cats, in which case the resulting network is often qualified of social. But they might as well exchange information through copper wires and optical fibres, or via chemical compounds. The bottom line is that networks can be found at all scales, spanning the entire planet or the confined spaces of a cell. Once these networks have been observed, studying their properties is of crucial interest, as it can unravel deep knowledge about underlying mechanisms. Some might say it could even lead to the answer to life, the universe, and everything. But that is not our point.

In the problems we are interested in, the network is not observed. We just have access to snapshots giving the individual states of the involved entities. The idea is to retrieve the hidden network organizing the system, that we assume to exist, from these snapshots. Let us consider the following thought experiment as an example. You put ten imaginary people in a room. As the mastermind behind this Gedankenexperiment, you give these people instructions before leaving the imaginary room. They have to choose an order between them. Then, the first person tosses a coin. If it is a heads, this person sits otherwise they stand. The next person tosses a coin. If it is a heads, they behave like the person before. Otherwise, an other coin toss decides if they stand or sit, and so on. When you get back in the room, you can only see the posture of each person. Nonetheless, if this procedure is repeated many times with the same order, you might finally be able to make a good guess on the order your imaginary friends had chosen, based on the different configurations that you have seen. This is called network inference. Now, imagine that the procedure is repeated a hundred times, and that participants are allowed to change their order three times during the whole experiment. The changes can only be made between different runs of the procedure. Your goal is now to determine when the three changes occurred. This is called change-point 1 detection. These are basically the two problems that provided an occupation for yours truly during the passed few years.

Graphical Models & Bayesian Inference

From a more formal point of view, networks find a natural representation in mathematical objects called graphs. A graph is made of a set of vertices, that are meant to represent the entities of interest, and a collection of edges linking these vertices. These edges can either be directed or not. In the example of our thought experiment, each vertex symbolises a participant, and directed edges are drawn between consecutive participants in the chosen order. Network inference supposes that, while being impossible to observe directly, the underlying graph is visible through the distribution of the quantities that we actually observe. Thus, the probabilistic models deployed in such circumstances have to take that into account. This leads to the development of the so-called graphical models. The rationale behind these models is to depict the conditional independence properties of a multivariate distribution by means of a graph. The graph acts as a more easy to handle proxy encoding abstract probabilistic properties in a visual manner. The model used in the example above belongs to one of the simplest classes of graphical models called Markov chains. Knowing the posture of a person makes the knowledge about all the previous participants in the order irrelevant to predict the posture of the following participants. The general framework of graphical models is meant to extend this reasoning to arbitrary dependence structures.

Performing network inference in graphical models requires to consider the graph itself as a parameter in a larger model, that includes the graph, the distribution of the observations and the observations themselves. In a Pirandellian twist, this model is itself a graphical model, whose graphical representation is given below.

Graph Distribution Observations

Our choice was to consider this model from a Bayesian point of view. This decision was motivated by the fact that it would be inconvenient to deal with uncertainty on a graph in a frequentist setting because of their discrete nature. The posterior distributions provided by Bayesian inference are more practical than point estimates in this case. It was therefore necessary to specify prior distributions on graphs and distributions. Network inference then boils down to computing the posterior distribution on the space of graphs, conditionally to the observations. Said like this, it does not seem like much, but we are actually dealing with a distribution on a frighteningly large discrete set. For instance, considering all possible undirected graphs on p = 10 vertices, we are faced with a set of size 2 p(p-1)/2 = 2 45 ≈ 3.5 • 10 13 , and usual problems might involve many more vertices. Performing exact inference under these circumstances would not be possible in a time amounting to less than a few geological eras.

Exact Inference through Algebra

Whenever exact inference is out of reach, sampling methods are one road to go down. As an example, let us imagine that we are given a probability distribution on the set of undirected graphs and that we want to evaluate the probability associated to the set of graphs with one connected component. We can try to compute this probability directly, but this might turn out to be a quite tedious calculation if we consider more than a handful of vertices. Ruling out exact computation, an estimate of this probability can be obtained by sampling a lot of graphs from the distribution and looking at the proportion of these graphs that were in our set of interest, namely the graphs with one connected component. If we are sufficiently thorough in our exploration of graphs, we get a sketch of the distribution landscape that is accurate enough to produce a good estimate of the desired integral. This is called Monte Carlo integration. In cases where direct sampling is intractable, the roundabout way is to build a Markov chain whose stationary distribution is the targeted distribution. After a burnin period, samples of this chain can be used as surrogates for independent samples of the initial distribution. This is a brief description of the general principle behind Markov Chain Monte Carlo (MCMC) methods. Taking random walks in vast, jagged spaces, you might wander infinitely before getting anywhere, or find yourself caught in some local topographic feature. In spite of these potential pitfalls, MCMC algorithms have been extensively and successfully used in a variety of contexts, including graphical models.

We made the decision to follow a different path, shaded by trees, and taking short-cuts through algebraic loopholes. Our roadmap for graph-related inference was to focus on situations where exact calculations were tractable, with spanning trees as a common thread. We define such graphs as connected undirected graphs without any cycle. They can alternatively be seen as the sparsest connected graphs, or the most connected graphs without cycle. The basic principle of our approach on network inference was to replace a partial exploration of the full set of graphs by an exhaustive exploration of the subset made of spanning trees. As can be seen in Figure 1, the latter only represent a small fraction of undirected graphs. But in many practical cases, it is actually not necessary to allow all possible structures when trying to extract a network from the data. The Occam's razor principle states that assumptions introduced to explain a phenomenon must not be multiplied beyond necessity, and sparsity is a feature that is often sought-after in network inference. Restricting our attention to the bare spanning trees might therefore not be as ludicrous as it could seem at first sight.

As lovely as woodland paths can be, the main appeal of spanning trees, as far as we are concerned, lies in their specific algebraic properties. In the 19th century, Arthur Cayley and Gustav Kirchhoff both looked into the combinatorial analysis of spanning trees and showed that the number of trees embedded in a given graph could be computed in polynomial time, as the determinant of the Laplacian matrix associated to this graph. This result is either referred to as Cayley's formula, Kirchhoff's theorem or, in a less eponymous manner, as the Matrix-Tree theorem. In particular, it states that there are p p-2 spanning trees on p vertices. An extension of this theorem to weighted trees can be derived, under the assumption that the weight of a tree can be expressed as the product of the weights of its edges. Summation over spanning trees can therefore be performed in a reasonable amount of time and that mainly motivated our decision to use them for network inference.

Other situations involving summation over gigantic discrete sets are met by similar fortuitous algebraic blessing. Change-point detection is one of these situations. In this case, the discrete parameter is the partition of the observations into contiguous temporal segments, and summing over the N -1 K-1 segmentations of time 1; N into K segments can actually be performed in O(KN 2 ) time. Our goal is not to dress an extensive list of algebraic tricks. But it turns out that the aforementioned ones can be combined to perform efficient detection of change-points in the dependence structure of a multivariate time-series, which was our second objective.

Chapter 1

This chapter is meant to supply the background underpinning Chapters 2, 3 and 4. The first section gives a detailed presentation of graphical models. It introduces Markov properties with respect to undirected graphs, as well as various extensions required for network inference. A few words are dedicated to directed graphical models for the sake of completeness, despite no further use hereabouts. A second section presents two algebraic results that are respectively used in Chapters 2 and 3 to perform summations over spanning trees and segmentations. This chapter can be read independently from the others. Some results are redundantly stated in Chapters 2 and 3 as they are based on submitted and accepted articles.

Chapter 2

Chapter 2 is concerned with network inference. In such problems, the dependence structure of the multivariate distribution behind the observations has to be explicitly taken into account, and that is exactly what graphical models are meant to do. Following on from the works of [START_REF] Meilă | Tractable bayesian learning of tree belief networks[END_REF], [START_REF] Kirshner | Learning with Tree-Averaged Densities and Distributions[END_REF] and [START_REF] Lin | Learning Sparse Markov Network Structure via Ensemble-of-Trees Models[END_REF], we make use of the hyper Markov property introduced by [START_REF] Dawid | Hyper Markov Laws in the Statistical Analysis of Decomposable Graphical Models[END_REF] to provide a full and formal framework for Bayesian inference in tree-structured graphical models. We consider a hierarchical model in which the top level is made of a spanning tree T , leading to a distribution π satisfying the Markov property with respect to this graph, and finally to observations X drawn according to π. Advantage is taken of the Matrix-Tree theorem presented in Section 1.2.1 and its generalisation to forests to show that the inference of this model can basically be performed in cubic time with respect to the number of vertices. In particular, the posterior probability for a given edge to appear in a random tree is derived in an exact manner, allowing for arbitrary prior edge appearance probability. We provide a new proof of a result formulated by [START_REF] Kirshner | Learning with Tree-Averaged Densities and Distributions[END_REF] on the complexity of posterior edge probability computation. A simulation study addresses the influence of the tree-assumption on the accuracy of structure inference for non-tree-structured graphical models. The last section presents an application to flow cytometry data on a cellular signalling network. This chapter is available as a preprint on arXiv (Schwaller et al., 2015). We also worked on a variation of the approach presented here in which π is integrated out through a BIC approximation, published in the french journal Revue d'intelligence artificielle (Schwaller & Robin, 2015).

Chapter 3

We build on the network inference approach develop in Chapter 2 to perform change-point detection in the dependence structure of a multivariate time-series. The global model in which the inference is performed is depicted below.

θ Y m K T
For a number of segments fixed to K, a segmentation m of 1; N and a series of K spanning trees T = (T k ) K k=1 are drawn. Then each segment of m is assigned a set of parameters θ depending on the structure given by the corresponding tree. Finally, on each segment, observations Y are drawn independently according to a distribution governed by θ. Using a result of [START_REF] Rigaill | Exact posterior distributions and model selection criteria for multiple change-point detection problems[END_REF] allowing efficient summation over the segmentations of 1; N under some factorability assumption, we are able to compute quantities such as the posterior probability of a change-point occurring at a given time or posterior edge probabilities over time in O(N 2 p 3 ) time. We also provide a way to assess whether the status of an edge (or of the whole graph) remains identical throughout the time-series or not when segmentation m is given. We benchmarked our approach against a simpler one, in which dependence structure is not explicitly taken into account, on synthetic data. Finally, we tackled two datasets respectively originating from cellular biology and neuroscience.

This chapter has been published as an article in the journal Statistics & Computing (Schwaller & Robin, 2016).

Chapter 4

Chapter 4 presents some extensions and perspectives to the work presented in Chapters 2 and 3.

The first section explains how covariates can be integrated in what has been previously presented. When observations are assumed to be normally distributed, adapting the model to take covariates into account is quite straightforward. We show conjugate priors for the multivariate multiple linear regression induce hyperdistributions that can be used as a basis to build a compatible family of hyperdistributions. We also derive an expression for the marginal likelihood of the observations on any subset of vertices. Whenever observations cannot be modelled by a multivariate normal distribution, we suggested an approximate approach based on copulas, that we used on microbial ecology data (Jakuschkin et al., 2016).

In the second section, we explain how to introduce temporal dependence in the segmentation model described in Chapter 3. Our goal is to lift the independent process assumption that cannot reasonably be made in many practical cases. Our suggestion is to use Temporal Independence Models (TIMs) introduced by [START_REF] Siracusa | Tractable Bayesian Inference of Time-Series Dependence Structure[END_REF]. Doing so forces us to leave the framework of undirected models to use their directed counter-parts. Under factorisation assumptions similar to the undirected independent case, we show that inference within TIM whose dependence structures are directed trees can be performed with cubic complexity with respect to the number of vertices.

Finally, we investigate prior specification on segmentations. We make out a non-exhaustive bestiary of priors fitting in our framework. We also describe a different kind of prior in which taking into account an initial guess on the segmentation or information coming from a previous study is straightforward. The aim of this chapter is to extract the main tools needed in the following chapters from the existing literature. In order to be as self-contained as possible, a reasonable amount of definition is included. This chapter will also serve as a reference concerning the notations used afterwards.

The first section is concerned with graphical models, which are at the center of the work presented here. As suggested by their explicit name, the idea behind these models is to make use of graphs to represent the conditional independence relationships satisfied by the distribution of a multivariate random variable. Each component of the random variable is Background pictured by a vertex, and edges are supposed to represent dependences between the different components. It is however of crucial importance to carefully define what properties of the distribution can be derived from the graph, as there is actually more than one way to decode a graph into conditional independence statements. Their being equivalent in many practical cases should not hide the fact that they are not in general. We mainly focus on undirected graphical models, stating the different associated Markov properties and how they relate to each other. This section is mostly based on the work of [START_REF] Dawid | Hyper Markov Laws in the Statistical Analysis of Decomposable Graphical Models[END_REF], [START_REF] Lauritzen | Graphical Models[END_REF] and [START_REF] Byrne | Hyper and Structural Markov Laws for Graphical Models[END_REF].

The second section details two algebraic results that motivated our committed stance on exact Bayesian inference in graphical models. Indeed, Bayesian inference typically relies on integrated quantities such as the marginal likelihood of the data, where parameters are integrated out. In graphical models, one of these parameters is the graph depicting the dependences. In this particular case, computing marginal likelihood requires to sum over a set of frightening cardinality. However, some subsets of graphs can be explored efficiently by making extensive use of algebra. Spanning trees form one of these sets, for which summation can be efficiently performed by relying on a result called the Matrix-Tree theorem [START_REF] Chaiken | A Combinatorial Proof of the All Minors Matrix Tree Theorem[END_REF]. To be able to put this result to good use for Bayesian inference, the prior distributions on the other parameters have to be chosen carefully, and this is were the (hyper) Markov properties given in the previous section come into play. The second result that we present is related to segmentation problems, where an other obviously discrete parameter is involved. Integrating a function over the set of segmentations of 1; N into K segments is tractable as soon as this function factorises over the segments. An algorithm based on dynamic programming principles can be used to perform the summation in O(KN 2 ) time [START_REF] Rigaill | Exact posterior distributions and model selection criteria for multiple change-point detection problems[END_REF]. Informally speaking, a graphical model is a probabilistic model whose conditional (in)dependence structure between random variables is given by a graph. This framework has received a fair amount of attention recently, but the ideas can be traced back as far as the end of the 19th century with J. Willard Gibbs. Indeed, one of the scientific areas that popularised graphical models is statistical physics. As an example, let us consider a simple model, named after the physicist Ernst Ising, that can be used to describe a large system of magnetic dipoles (or spins). Spins can be in one of the two states {±1}. They are spread on a graph (commonly a lattice) and can only interact with their neighbours. If σ = (σ 1 , σ 2 , ... ) is a state of the system giving assignments for all spins, the energy of σ and the associated Gibbsian distribution are respectively given by:

H(σ) . . = -J i∼j σ i σ j -H i σ i , p(σ) . . = 1 Z exp (-H(σ)) , σ =
where i ∼ j means that vertices i and j are neighbours in the lattice. The graphical aspect of this model is rather obvious, since the definition of H depends on the neighbourhood of each spin. This is an example of undirected graphical model. Such models are also called Background Markov random fields.

Graphical models also naturally arise for instance when designing hierarchical models with sequentially drawn variables. Let us consider a classical Bayesian framework where observations X are drawn according to a distribution with parameters θ, and where θ is itself drawn from a distribution with (hyper)parameters λ. This model can be depicted by a directed graph. θ λ X Here we have an example of directed graphical model. The idea is that the graph indicates a way to factorise the joint probability distribution of all variables as a product of conditional probability distributions. For this factorisation to be possible, the graph cannot have any directed cycle. It has to be a directed acyclic graph (DAG). These models are often referred to as Bayesian networks [START_REF] Jensen | Bayesian Networks and Decision Graphs[END_REF]. Other classical examples include Markov chains and hidden Markov models (HMM).

Our main concern will be with undirected graphical models. Nonetheless, as the models that we are especially interested in, namely tree-structured graphical models, can be equally represented in both formalism, we give a brief description of directed graphical models at the end of this section.

Graphs & Markov Properties

1.1.1.a Conditional Independence

In graphical models, the probabilistic notion that graphs are meant to represent is conditional (in)dependence. We give the most general definition of conditional independence for a triplet of random variables defined on their respective probability spaces, as stated by [START_REF] Lauritzen | Graphical Models[END_REF].

Definition 1.1 (Conditional independence). Let X, Y, Z be random variables with a joint distribution P . X is said to be conditionally independent of Y given Z, written

X |= Y |Z,
if, for any measurable set A in the sample space of X, there exists a version of P (A|Y, Z) that is a function of Z alone.

Whenever P admits a density with respect to some product measure, conditional independence can be more easily expressed.

Proposition 1.1. Let X, Y, Z be random variables with a joint distribution P that admits a density p with respect to some product measure. Then

X |= Y |Z ⇔ p(x, y|z) = p(x|z)p(y|z) a.s.
Here are some properties satisfied by the ternary relation X |= Y |Z.

Proposition 1.2. Let h be an arbitrary measurable function on the sample space of X.

Then, it holds that

Background (C1) if X |= Y |Z then Y |= X|Z; (C2) if X |= Y |Z and U = h(X), then U |= Y |Z; (C3) if X |= Y |Z and U = h(X), then X |= Y |(Z, U ); (C4) if X |= Y |Z and X |= W |(Y, Z), then X |= (W, Y )|Z.
Whenever the joint distribution of X, Y, Z admits a continuous and positive density, it also holds that

(C5) if X |= Y |Z and X |= Z|Y then X |= (Y, Z).

1.1.1.b Undirected Graphs

Let V = {1, ..., p}, p 2, and let P 2 (V ) denote the subsets of V of size 2.

Definition 1.2. For E ⊆ P 2 (V ), G = (V, E) is the undirected graph with vertices V and edges E. A graph G = (V, E) is said to be complete if E = P 2 (V ). For a subset A ⊂ V , the subgraph of G induced by A is defined as G A . . = (A, E A ) with E A . . = {(i, j) ∈ E|i ∈ A, j ∈ A}. Whenever G A is complete for A ⊆ V , A is said to be a clique of G.
For α ∈ V and G = (V, E), we define the following set of vertices:

boundary of α : bd(α) = {β ∈ V \ {α}|{α, β} ∈ E}, closure of α : cl(α) = bd(α) ∪ {α}.
Vertices in bd(α) are also called the neighbours of α.

Definition 1.3 (Path)

. Let G = (V, E) be an undirected graph and α, β be two distinct vertices in V . A path from α to β is a sequence α = γ 0 , ... , γ n = β, n 1, of distinct vertices such that, for all 1 i n, {γ i-1 , γ i } ∈ E.

A graph is connected if, ∀{α, β} ∈ P 2 (V ), there is a path between α and β.

Definition 1.4. A (spanning) tree is a connected graph on V with cliques of size at most 2. On undirected graphs, the core notion for graphical models is separation. 10

Definition 1.5 (Separation). Let A, B, C be subsets of V . C is said to separate A from B if any path from α ∈ A to β ∈ B intersects C. We let A |= G B|C denote the statement "C separates A from B in G".
Background (C2') if A |= G B|C and D ⊆ A, then D |= G B|C.
Properties (C1) to (C5) can in fact be seen as pure formal statements related to the notion of "irrelevance". Any ternary relation on the subsets of a finite set that satisfies properties (C1) to (C5), where h is replaced by set inclusion, is called a graphoid. If (C5) does not hold, it is only a semi-graphoid. The name graphoid goes back to [START_REF] Pearl | Graphoids: a graph-based logic for reasoning about relevance relations[END_REF] who noticed that these properties where in fact perfectly captured by graphs. The rationale behind graphical models is to use a graphoid defined by a graph to represent the graphoid (or semi-graphoid) induced by a probability distribution.

Example 1.1. A and B are separated by C.

A C B

Before actually connecting the notions of conditional independence and separation, we describe a particular class of undirected graphs, which are of particular interest because they can be recursively broken down into their maximal cliques.

Definition 1.6 (Decomposition). A pair

(A, B) of subsets of V is said to be a decomposition of G if V = A ∪ B, the subgraph induced by G on A ∩ B is complete and A ∩ B separates A from B.
If A and B are both proper subsets of V , the decomposition is said to be proper.

A graph G is said to be decomposable if it is either complete or if there exists a proper decomposition of G into two decomposable subgraphs.

Notice that the definition of a decomposable graph is recursive. From a conceptual point, the equivalent notion of chordality is easier to handle.

Definition 1.7 (Chordal graphs).

A graph G is chordal if all cycles of four or more vertices have a chord, which is an edge that is not part of the cycle but connects two vertices of the cycle.

Proposition 1.3 (Lauritzen, 1996, Prop. 2.5)

. (G is decomposable) ⇔ (G is chordal)
Example 1.2. The smallest example of non-decomposable graph is the cycle on four vertices. Whenever a graph is decomposable, one can build the set of its minimal complete separators. Let G be a decomposable graph. The set S of minimal separators for G is given by the following algorithm. We begin with an empty set of separators S = ∅. We consider a proper decomposition (A, B) of G such that A ∩ B is of minimal cardinality. If there is no Background such decomposition, it means that G is complete and we stop. Otherwise, we add A ∩ B to S and apply the same procedure to the subgraphs G A and G B . For any subset S ∈ S, we let ν(S) denote the number of times it appeared in the procedure. An other product of this algorithm is the set C containing the maximal cliques of G. This is often referred to as the junction tree algorithm. Indeed, the result of this algorithm can be represented by a factor graph, i.e. a bipartite graph whose vertices are indexed by C and S. We put an edge between C ∈ C and S ∈ S if and only if C ∩ S = S. The graph obtain this way is a tree, hence the name given to this algorithm (see Example 1.3). If d(S) denotes the degree of the node corresponding to S in the junction tree, we have that ν(S) = d(S) -1. We will see in Section 1.1.1.d that junction trees can be used to obtain an explicit form for the factorisation of a distribution. This definition of the junction tree of a decomposable graph is different from the one given in [START_REF] Lauritzen | Graphical Models[END_REF]. It is sometimes referred to as an 'Almond tree', after [START_REF] Almond | Optimality issues in constructing a Markov tree from graphical models[END_REF].

Example 1.3. A decomposable graph and its junction tree.

c Markov Properties

Let X = (X 1 , . . . , X p ) be a random vector taking values in a product space X = p i=1 X i . The set of probability distributions on X is denoted by Π(X ). For any A ⊂ V , we let X A denote (X α ) α∈A . The same notation is used for any x ∈ X . In the following, for any subsets A, B, C ⊂ V , we will use the short notation

Background

The Markov properties (P ), (L) and (G) are related as described in the proposition below.

Proposition 1.4 [START_REF] Pearl | Graphoids: a graph-based logic for reasoning about relevance relations[END_REF]. For any undirected graph G and any probability distribution P on X , it holds that (G) ⇒ (L) ⇒ (P ).

If P has a positive and continuous density with respect to a product measure, then

(G) ⇔ (L) ⇔ (P ).
Proof. We prove (G) ⇒ (L) ⇒ (P ) in the general case. The proof for (P ) ⇒ (G) in the case where P has a positive and continuous density with respect to a product measure is obtained through the Hammersley-Clifford theorem (Theorem 1.1 in Section 1.1.1.d).

(G) ⇒ (L): For α ∈ V , any path from α to a vertex in V \ cl(α) has to go through a neighbour of α. Therefore, bd(α) separates α from V \ cl(α) and (G) ⇒ (L).

(L) ⇒ (P ): We suppose that (L) is true. Let α and β be two distinct vertices of V such that {α, β} / ∈ E.

By (L), α |= V \ cl(α)|bd(α). As {α, β} / ∈ E, β ∈ V \ cl(α). By statement (C3), we obtain that α |= V \ cl(α)| bd(α) ∪ (V \ {cl(α) ∪ β} V \{α,β}
and (C2) gives that α |= β|V \ {α, β}.

Whenever statement (C5) holds, all Markov properties are equivalent. But when it does not hold, it is possible to find distributions that satisfy one but not the other. We give an example of such a distribution below.

Example 1.4. Let X, Y, Z be three random variables with values in {0, 1}. Let X = Y = Z and P (X = 0) = P (X = 1) = 1/2. The probability distribution of (X, Y, Z) satisfies the pairwise Markov property but not the local Markov property with respect to the graph given below.

X

Y Z Indeed, we have that X |= Y |Z and X |= Z|Y , as X = Y = Z thus leading to trivial statements. But conditioning on the neighbours of X, i.e. not conditioning at all, does not make it independent from Y and Z.

For the complete graph on V , it is not possible to find disjoint subsets A, B, S such that S separates A from B. Thus, any distribution is globally Markov with respect to the complete graph. More generally, if P is globally Markov with respect to a graph G = (V, E) and if G = (V, E ) is such that E ⊆ E , then P is globally Markov with respect to G . Whenever a graph G perfectly describes the conditional independence properties of a distribution P , we are talking about Markov faithfulness. Definition 1.9 (Markov faithfulness). A distribution P is said to be Markov faithful to a graph G if, for any triple (A, B, S) of disjoint subsets of V , it holds that

Background A |= G B|S ⇔ A |= B|S.
Markov faithfulness forbids the distribution to satisfy any conditional independence statement that is not encoded in the graph. If P is Markov faithful to a graph G, G is said to be a perfect map (or P-map) for P . If only the implication

A |= G B|S ⇒ A |= B|S.
holds, then G is said to be an independence map (or I-map). On the contrary, if we have that

A |= G B|S ⇐ A |= B|S.
then G is called a dependence map (or, as you might have guessed by now, a D-map). A Pmap is both a D-map and an I-map. Saying that P is Markov with respect to G is equivalent to saying that G is an I-map for P . There is just a change of focus in the terminology.

1.1.1.d Factorisation

Closely related to conditional independence is the property of factorisation, as can be seen in Proposition 1.1.

Definition 1.10 (Factorisation). A probability measure P on X is said to factorise over G (or to satisfy (F)) if it admits a density p with respect to some product measure on X of the form

p(x) = A∈A ψ A (x A ), ∀x ∈ X , (F )
where A are complete subsets of G or, equivalently, if

p(x) = C∈C ψC (x C ), ∀x ∈ X , (F )
where C are the maximal cliques of G.

The relationship between property (F ) and the different Markov properties is given by the following proposition.

Proposition 1.5 (Lauritzen, 1996, Prop. 3.8). For any undirected graph G and any probability distribution P on X , it holds that

(F ) ⇒ (G) ⇒ (L) ⇒ (P ).
If a distribution with density p is such that, ∀x ∈ X , p(x) > 0, we have seen that all Markov properties were equivalent. In this case, we also have that (P ) implies (F ).

Theorem 1.1 [START_REF] Hammersley | Learning bayesian networks: The combination of knowledge and statistical data[END_REF]. A probability distribution P with positive and continuous density p with respect to some product measure on X satisfies the pairwise Markov property with respect to a an undirected graph G if and only if it factorises according to G. Then

(F ) ⇔ (G) ⇔ (L) ⇔ (P ).
Proof. The proof, drawn from [START_REF] Lauritzen | Graphical Models[END_REF], is based on Möbius inversion lemma.

14

Background

Lemma (Möbius inversion). Let φ be a function defined on the subsets of V , taking values in an Abelian group. We let ζ be the function defined by

ζ(A) . . = B⊆A (-1) |A\B| φ(B), ∀A ⊆ V.
Then, it holds that

φ(A) = B⊆A ζ(B), ∀A ⊆ V.
We assume that (P ) is true. Let x * be an arbitrary but fixed element of X . For A ⊆ V , we define φ A and ζ A by

   φ A (x) = log p(x A , x * V \A ), ζ A (x) = B⊆A (-1) |A\B| φ B (x), ∀x ∈ X .
Both φ A and ζ A are in fact functions of x A . Möbius inversion lemma yields that

φ V (x) = log p(x) = A⊆V ζ A (x), ∀x ∈ X .
It remains to show that ζ A ≡ 0 as soon as A is not complete. Let A be a subset containing two distinct vertices α and β such that {α, β} / ∈ E. Subsets of A can be sorted according to their containing α and β, so that, if

C = A \ {α, β}, ζ A = B⊆C (-1) |C\B| φ B -φ B∪{α} -φ B∪{β} + φ B∪{α,β} .
Let D = V \ {α, β}. Then, using the pairwise Markov property, ∀x ∈ X ,

φ B∪{α,β} (x) -φ B∪{α} (x) = log p(x B , x α , x β , x * D\B ) p(x B , x α , x * β , x * D\B ) = log p(x α |x B , x * D\B )p(x β , x B , x * D\B ) p(x α |x B , x * D\B )p(x * β , x B , x * D\B ) = log p(x * α |x B , x * D\B )p(x β , x B , x * D\B ) p(x * α |x B , x * D\B )p(x * β , x B , x * D\B ) = log p(x B , x * α , x β , x * D\B ) p(x B , x * α , x * β , x * D\B ) = φ B∪{β} (x) -φ B (x),
therefore showing that ζ A ≡ 0. We have proved that (F ) is true.

In the general case, (F ) is stronger than (G), has shown by the example below.

Example 1.5. [START_REF] Moussouris | Gibbs and Markov Random Systems with Constraints[END_REF] Here is an example of a distribution satisfying (G) but not (F ). Consider the uniform distribution on the 8 (out of 16) configurations of {0, 1} 4 displayed below.

1 1 1 1 1 1 1 0 0 1 1 0 0 0 1 0 0 0 0 0 0 0 0 1 1 0 0 1 1 1 0 1 15

Background

To show that this distribution is globally Markov with respect to , one only has to show that two opposite vertices are independent conditionally on the other two. But when the values of two diagonally opposite vertices are fixed, so is the value of one of the other two vertices, the global Markov property is therefore satisfied. We now show that the density does not factorise. Let us assume that it does factorise. Then p(0, 0, 0, 0) = 1/8 = ψ 12 (0, 0)ψ 23 (0, 0)ψ 34 (0, 0)ψ 41 (0, 0), so these factors are all strictly positive. Similarly reasoning on all 8 possible configurations yields that all factors ψ are strictly positive. This contradicts the fact that 8 configurations have probability 0.

For decomposable graphs, (F ) and (G) are equivalent even if p is not positive, hence the fondness caused by decomposability.

Proposition 1.6 (Lauritzen, 1996, Prop. 3.19). Let G be a decomposable graph. Then it holds that

(F ) ⇔ (G).
Moreover, whenever G is decomposable, the junction tree algorithm described in Section 1.1.1.b can be used to obtain an explicit formula for the density of a distribution that factorises over G using marginal distributions over complete subsets. Indeed, if P is a distribution with density p that factorises over a decomposable graph G, and if S and C respectively denote the set of minimal separators and maximal cliques yielded by the junction tree algorithm, it can be shown that

p(x) S∈S p S (x S ) ν(S) = C∈C p C (x C ), ∀x ∈ X ,
where, for A ⊆ V , p A is the marginal density of p on X A . If p is positive, then

p(x) = C∈C p C (x C ) S∈S p S (x S ) ν(S) , ∀x ∈ X .
(1.1)

Example 1.6. On the graph given in Example 1.3, (1.1) would yield

p(x) = p(x 1 , x 2 , x 3 )p(x 2 , x 3 , x 4 )p(x 4 , x 5 )p(x 4 , x 6 ) p(x 2 , x 3 )p(x 4 ) 2 .
Whenever G = (V, E) is a tree, the factorisation is directly given by

p(x) = {i,j}∈E p ij (x i , x j ) i∈V p i (x i ) d(i)-1 = i∈V p i (x i ) {i,j}∈E p ij (x i , x j ) p i (x i )p j (x j ) , ∀x ∈ X ,
where d(i) stands for the degree of vertex i in G.

Background

1.1.1.e Graphical Models

In the previous sections, we have seen how graphs could be used to represent the conditional independence properties of a distribution. This connection between graphs and distributions can be used to define the notion of graphical model.

For any undirected graph G, we let M(G) denote the set of distributions on X that are globally Markov with respect to G. What we call a model is a family of distributions on X . Thus, a graphical model is a family of distributions sharing some Markov property with respect to a graph. The formal definition that we adopt here is the following.

Definition 1.11. A graphical model is given by a graph G describing its structure and a family of distributions F ⊂ M(G) globally Markov with respect to G.

If all distributions in F admit a positive density, one might use any of the Markov properties, as they are all equivalent. Alternatively, one might require the distributions in F to be Markov faithful to G. This definition is much more restrictive.

Network Inference in Graphical Models

Many scientific fields are concerned with retrieving an underlying network that is supposed to exist between entities of interest. This kind of problem can be given many names, but ours is network inference. In this situation, the graph and the distribution of the observations are themselves parameters of a hierarchical model that can be depicted as follows. In this model, π is assumed to be globally Markov with respect to G, and X is distributed according to π. The distribution π is drawn from a family F G of distributions that satisfy the global Markov property with respect to G. Thus, (G, F G ) is a graphical model, and network inference can be seen as a model selection problem between different graphical models.

There exists an abundant literature discussing network inference in a frequentist framework, maybe even vaster than its Bayesian equivalent. As we still decided to go the Bayesian way, this paragraph has no further claim than mentioning a very small fraction of this literature. Gaussian Graphical Models (GGMs) have been of particular interest because of the specific properties of normal distributions. Indeed, for such distributions, conditional independence properties can directly be read in the inverse covariance (or precision) matrix. Learning the dependence structure of the distribution therefore boils down to identifying the null entries in the precision matrix. This task can be seen as a set of variable selection problems, where the neighbours of each variable have to be selected individually. In a frequentist setting, sparse regularisation techniques such as the Lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso: A retrospective[END_REF] are well suited for this task. Various algorithms have been especially developed in this context, including 'neighbourhood selection' [START_REF] Meinshausen | High-dimensional graphs and variable selection with the Lasso[END_REF] and the 'Graphical Lasso' [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF]. Score-based approaches have also been considered for network inference. In the case of DAGs for instance, [START_REF] Hammersley | Learning bayesian networks: The combination of knowledge and statistical data[END_REF] introduced the BDe metric that was used in a greedy search algorithm called GES [START_REF] Chickering | Optimal structure identification with greedy search[END_REF].

Background

When considering Bayesian inference, a prior distribution has to be specified for G. Without any restriction on the support of this distribution, exact inference can only be contemplated as long as there are no more than thirty or so variables of interest [START_REF] Parviainen | Exact Structure Discovery in Bayesian Networks with Less Space[END_REF]. In this context, sampling approaches such as MCMC have been extensively used to get access to the posterior graph distribution, whether it be for undirected graphs [START_REF] Roverato | Hyper inverse wishart distribution for non-decomposable graphs and its application to bayesian inference for gaussian graphical models[END_REF][START_REF] Atay-Kayis | A Monte Carlo method to compute the marginal likelihood in non decomposable graphical Gaussian models[END_REF][START_REF] Green | Sampling decomposable graphs using a markov chain on junction trees[END_REF] or DAGs [START_REF] Madigan | Bayesian graphical models for discrete data[END_REF][START_REF] Friedman | Being Bayesian about network structure. A Bayesian approach to structure discovery in Bayesian networks[END_REF][START_REF] Niinimäki | Partial order mcmc for structure discovery in bayesian networks[END_REF]. Some interest has also been dedicated to specific subsets of graphs, such as spanning trees [START_REF] Chow | Approximating Discrete Probability Distributions with Dependence Trees[END_REF][START_REF] Meilă | Learning with Mixtures of Trees[END_REF][START_REF] Meilă | Tractable bayesian learning of tree belief networks[END_REF][START_REF] Kirshner | Learning with Tree-Averaged Densities and Distributions[END_REF][START_REF] Lin | Learning Sparse Markov Network Structure via Ensemble-of-Trees Models[END_REF]. In this particular case, the inference can be performed with cubic complexity with respect to the number of variables, as we shall see in the following chapter.

Defining a full and formal framework for Bayesian inference in tree-structured graphical models, we will find it usefull to extend Markov properties to the distributions of π and G. These extensions are respectively called hyper and structural Markov properties and are described in Sections 1.1.2.a and 1.1.2.c.

1.1.2.a Hyper Markov Properties for Hyperdistributions

Hyper Markov properties were introduced by [START_REF] Dawid | Hyper Markov Laws in the Statistical Analysis of Decomposable Graphical Models[END_REF] as an extension of Markov properties for distributions on sets of Markov distributions. In order to establish a clear distinction between the two levels of involved distributions, a distribution for π will be called a hyperdistribution. Whenever π is taken from a parametric family of distributions depending on a set of parameters θ ∈ Θ, a hyperdistribution ρ can be specified by giving a distribution on Θ. This distribution can itself be taken in a parametric family, in which case these parameters are traditionally called hyperparameters in a Bayesian setting. Hence the term hyperdistribution that we have adopted to refer to ρ. In order to remain as general as possible, all results are nonetheless given in a non-parametric framework.

Let G = (V, E) be a decomposable graph and let m G = (G, F G ) be a graphical model with structure G. We consider a hyperdistribution ρ defined on F G . For π ∈ F G and A, B ⊆ V , we let π A denote the marginal distribution obtained from π on the variables X A , and π B|A denote the collection of conditional distributions of X B |X A under π. We also let ρ A and ρ B|A respectively denote the marginal hyperdistribution induced by ρ on π A and the collection of hyperdistributions induced by ρ on π B|A . Definition 1.12. A hyperdistribution ρ is said to be (weak) hyper Markov with respect to G if, for any decomposition (A, B) of G, it holds that, under ρ,

π A |= π B |π A∩B . (1.2)
As stated by [START_REF] Dawid | Hyper Markov Laws in the Statistical Analysis of Decomposable Graphical Models[END_REF], condition (1.2) is equivalent to either

π A|B |= π B |π A∩B , π A|B |= π B|A |π A∩B .
The hyper Markov property can also be characterised through separation.

Theorem 1.2 (Byrne, 2011, Th. 1.3.1)

Background

We will actually need a stronger version of the hyper Markov property, for reasons that we will expose shortly after. Definition 1.13. A hyperdistribution ρ is said to be strong hyper Markov with respect to G if, for any decomposition (A, B) of G, it holds that, under ρ,

π A |= π B|A .
If G = (V, E) and G = (V, E ) are such that E ⊂ E and if ρ is weak hyper Markov with respect to G, then it is also hyper Markov with respect to G . This is not true for the strong hyper Markov property. But as the strong hyper Markov property implies the weak version, ρ being strong hyper Markov with respect to G implies that ρ is weak hyper Markov with respect to G . Proposition 1.7 (Dawid & Lauritzen, 1993, Prop. 5.1). Let ρ be a hyperdistribution that is (weak) hyper Markov with respect to a graph G. Then, for any decomposition (A, B) of G, it holds that

(X A , π A ) |= (X B , π B|A )|(X A∩B , π A∩B ).
If ρ is strong hyper Markov with respect to G, for any decomposition (A, B) of G, it further holds that

(X A , π A ) |= (X B , π B|A )|X A∩B .
We remind that ρ is a prior distribution for the distribution π of X. If we abusively let ρ either denote the distribution itself or its density with respect to some product measure η on Π(X ), the density of the marginal distribution for X can be written as p(x) . . = π(x)ρ(π)dη(π).

(1.3) Whenever ρ is strong hyper Markov, for any decomposition (A, B) of G, X A∩B is enough to isolate (X A , π A ) from (X B , π B|A ), so that it is possible to integrate over π, and the marginal distribution over X remains globally Markov with respect to G. This is not the case when ρ is only weak hyper Markov with respect to G.

Proposition 1.8 (Dawid & Lauritzen, 1993, Prop. 5.6). If ρ is strong hyper Markov with respect to G, the marginal distribution defined in (1.3) is globally Markov with respect to G.

If G is chosen to be decomposable and if ρ is strong hyper Markov with respect to G, Propositions 1.6 and 1.8 show that the marginal likelihood factorises over G. The integral given in (1.3) can then be computed on the maximal cliques C and minimal separators S of G as

p C (x C ) = π C (x)ρ C (π C )dη(π C ), ∀C ∈ C, (1.4) p S (x S ) = π S (x)ρ S (π S )dη(π S ), ∀S ∈ S, (1.5)
and the complete marginal likelihood is obtained through (1.1).

Background

In the model described in Figure 1.2, we are actually interested in defining a collection of hyperdistributions {ρ G } G∈G , where G is the support of a graph distribution ξ. In this situation, it might be desirable for ρ G and ρ G to be 'similar' whenever two graphs G and G of G are close. This can be achieved by building what is called a (hyper) compatible family of hyperdistributions. Let us consider a hyperdistribution on Π(X ) such that, for any A ⊆ V , under ρ,

π A |= π V \A|A . (1.6)
This means that ρ is strong hyper Markov with respect to the complete graph over V .

Proposition 1.9. (Dawid & Lauritzen, 1993, §6.2) For any decomposable graph G, there exists a unique hyperdistribution ρ G on M(G) that is strong hyper Markov with respect to G and such that, for every clique C of G,

ρ G C = ρ C .
If G is a graph family, {ρ G } G∈G is a (hyper) compatible family of strong hyper Markov hyperdistributions.

Example 1.7. We illustrate the hyper Markov property on the special case of normal distributions. We let N p denote the set of multivariate normal distributions on X = R p with mean vector 0 p . We also let P p denote the set of symmetric positive-definite real-valued matrices of size p × p. There is a one-to-one correspondence between P p and N p through the map

Υ p : P p -→ N p Λ -→ N (0 p , Λ -1 )
where N (0 p , Λ -1 ) stands for the multivariate normal distribution with mean vector 0 p and inverse covariance matrix (or precision matrix) Λ. For any undirected graph G, we also consider the set N G p of normal distributions that are globally Markov with respect to G. As normal distributions are positive, all Markov properties are equivalent and we can actually drop the "globally". The distributions in N G p can be characterised by a condition on their precision matrices.

Proposition 1.10. For a graph G = (V, E), we define

P G p = {Λ ∈ P p : ∀{i, j} / ∈ E, Λ ij = 0}.
Then, it holds that

N G p = Υ p (P G p ).
Proof. See for instance (Lauritzen, 1996, Prop. 5.2).

The Wishart distribution W p (α, Ψ) with α > p -1 degrees of freedom and scale matrix Ψ is defined on P p and has density

h(Λ) = |Ψ| α 2 2 αp 2 Γ p ( α 2 ) |Λ| α-p-1 2 e -1 2 tr(ΨΛ) , 20
Background where the multivariate gamma function Γ p is defined by

Γ p x 2 = π p(p-1) 4 p j=1 Γ x + 1 -j 2 .
Therefore, any Wishart distribution induces a hyperdistribution ρ on N p through Υ p :

ρ(π) = h(Υ -1 p (π)), ∀π ∈ N p .
In order to build a compatible family of strong hyper Markov hyperdistributions through Proposition 1.9, one has to prove that Condition (1.6) is satisfied for a hyperdistribution on N p induced by a Wishart distribution. For π ∈ N p with precision matrix Λ, classic results on normal distributions state that, for any partitioning (A, B) of V ,

π A = N (0 |A| , (Λ AA -Λ AB Λ -1 BB Λ BA ) -1 ), π B|X A =xa = N (-Λ -1 BB Λ BA x A , Λ -1 BB ).
where

[Λ AA , Λ AB , Λ BA , Λ BB ] is a block partitioning of Λ according to (A, B). Thus, show- ing that π A |= π B|A boils down to show that (Λ AA -Λ AB Λ -1 BB Λ BA ) is independent of {Λ BB , Λ AB }. Proposition 1.11. Let (A, B) be a partition of V , with |A| = a and let Λ ∼ W p (α, Ψ). We denote Λ AA•B the matrix Λ AA -Λ AB Λ -1 BB Λ BA . Then, it holds that Λ AA•B |= {Λ BB , Λ AB }.
Moreover,

Λ AA•B ∼ W a (α -p + a, Ψ AA ), (1.7)
where [Ψ AA , Ψ AB , Ψ BA , Ψ BB ] is a block partitioning of Ψ according to (A, B).

Proof. This is a particular case of Theorem 5 in [START_REF] Geiger | Parameter priors for directed acyclic graphical models and the characterization of several probability distributions[END_REF].

Therefore, any Wishart distribution W p (α, Ψ) can be used to build a compatible family of strong hyper Markov hyperdistributions, and (1.7) states that the clique hyperdistributions of this family are given by

ρ A (π A ) = h A (Υ -1 a (π A )), ∀A ⊆ V,
where h A stands for the density of W a (α -p + a, Ψ AA ). Exact formulas can be derived for the marginal likelihoods p C and p S given in (1.4) and (1.5) using classic results on conjugate priors.

1.1.2.b Meta Markov Models

Let ρ be a hyperdistribution on Π(X ). The support family of ρ is the subfamily of Π(X ) defined by

F ρ . . = {π ∈ Π(X ) : ρ(π) > 0}.
The support family of a (strong) hyper Markov hyperdistribution has some intrinsic Markov property of its own. Expressing this property requires to define variation independence.

Background

Definition 1.14. Let ϕ, ς and be functions on a common domain D. ϕ is said to be conditionally variation independent of ς given , written

ϕ ‡ ς| [D],
if, for all v ∈ ς(D), w ∈ (D), it holds that

ϕ (ς, ) -1 (v, w) = ϕ -1 (w) .
If is trivial, we write ϕ ‡ ς, meaning that, ∀v ∈ ς(D), ϕ ς -1 (v) = ϕ(D).

Basically, the statement "ϕ ‡ς| [D]" means that, whenever the value taken by is fixed, the range of values that can be attained by ϕ does not change with the value taken by ς. Conditional variation independence is a property related to a domain, just like conditional independence is related to a joint probability distribution. It can be seen as weaker than conditional independence, as the lack of variation independence proscribe the corresponding probabilistic independence statement.

Proposition 1.12 (Dawid & Lauritzen, 1993, Lemma 4.2). Properties (C1) to (C4) defined for conditional independence in Proposition 1.2 hold for conditional variation independence.

For A, B ⊆ V , we consider the following functions on Π(X ):

Ξ A : π -→ π A , Ξ B|A : π -→ π B|A .
Definition 1.15. Let G be an undirected graph. A model F ⊆ Π(X ) is said to be (weak) meta Markov with respect to G if, for any decomposition (A, B) of G it holds that

Ξ A ‡ Ξ B |Ξ A∩B [F].
A model F ⊆ Π(X ) is said to be strong meta Markov with respect to G if, for any decomposition (A, B) of G it holds that

Ξ A ‡ Ξ B|A [F].
Within a strong meta Markov model, any value taken by π A is logically compatible with any value taken π B|A . The following propositions describe two cases involving meta Markov models.

Proposition 1.13. Let G be a graph. The family M(G) made of all the distributions in Π(X ) that are Markov with respect to G is a meta Markov model with respect to G. Proposition 1.14. Let ρ be a hyperdistribution that is (strong) hyper Markov with respect to a graph G. Then the support family F ρ of ρ is a (strong) meta Markov model with respect to G.

Therefore, if F is not (strong) meta Markov with respect to a graph G, it will not be possible to find a hyperdistribution with support F that is (strong) hyper Markov with respect to G.

Example 1.8. The family of multivariate normal distributions on R p is a strong meta Markov model with respect to the complete graph on V (see Example 1.7). The family of multinomial distributions on 1; r p , r ∈ N * , is a strong meta Markov model with respect to the complete graph on V .

Background

1.1.2.c Structural Markov Property for Graph Distributions

In the model depicted in Figure 1.2, the graph itself is a parameter, for which we therefore need a prior distribution. Pursuing an approach similar to what [START_REF] Dawid | Hyper Markov Laws in the Statistical Analysis of Decomposable Graphical Models[END_REF] had done with hyperdistributions, [START_REF] Byrne | Hyper and Structural Markov Laws for Graphical Models[END_REF] proposed an extension of Markov properties for graph distributions.

We let U denote the set of undirected decomposable graphs on V . A covering pair of V is a pair of sets (A, B) such that A ∪ B = V . For any family of graphs G ⊆ U and for any covering pair (A, B), we define G(A, B) to be the set of graphs G ∈ G for which (A, B) is a decomposition. In particular, U(A, B) is the set of decomposable graphs for which (A, B) is a decomposition. We remind that, for G ∈ U(A, B), G A∩B is complete.

Definition 1.16. A graph distribution ξ is said to be structurally Markov if, for any covering pair (A, B) of V , it holds that

G A |= G B |{G ∈ U(A, B)}.
Structurally Markov graph distributions combine nicely with compatible families of strong hyper Markov hyperdistributions.

Proposition 1.15 (Byrne, 2011, Th. 4.4.4). Let ξ be a structurally Markov graph distribution and let {ρ G } G∈U be a compatible family of strong hyper Markov distributions. Then, for any covering pair (A, B) of V , it holds that

(π A , G A ) |= (π B|A , G B )|{G ∈ U(A, B)}.
Under these conditions, the posterior graph distribution is also structurally Markov.

Example 1.9. A graph distribution ξ on the set of trees, denoted T , can easily be specified through an edge weight matrix. Indeed, if we are given a symmetric matrix ω . . = (ω i,j ) 1 i,j p with non-negative entries, a distribution on T can be obtained through

ξ(T ) = 1 Z {i,j}∈E T ω i,j , ∀T = (V, E T ) ∈ T , (1.8)
where Z is a normalising constant defined by Z . . = T ∈T {i,j}∈E T ω i,j . In this kind of distributions, the weight of a tree is just the product of the weights of its edges. Thus, an edge shared by two different trees contributes in the same way in both trees. It can easily be shown that the distributions on T of the form given in (1.8) are structurally Markov.

These distributions form what [START_REF] Byrne | Structural markov graph laws for bayesian model uncertainty[END_REF] call the clique exponential family on T .

The meta Markov property that we defined for families of distributions (or models) can similarly be defined for graph families. For A ⊂ V , we consider the following function on U:

G A : G -→ G A .
Definition 1.17. A graph family G ⊆ U is said to be structurally meta Markov if, for any covering pair (A, B) of V , it holds that

G A ‡ G B [G(A, B)]. Background For a covering pair (A, B) of V , if G, G ∈ U(A, B) , we let G A ⊗ G B denote the one and only graph in U(A, B) such that (G A ⊗ G B ) A = G A and (G A ⊗ G B ) B = G B .
The structural meta Markov property can equivalently be defined as follows.

Theorem 1.3 (Byrne, 2011, Th. 4.3.1). A graph family G ⊆ U is structurally meta Markov if and only if, for any covering pair (A, B) of V and for G, G ∈ G(A, B), it holds that

G A ⊗ G B ∈ G(A, B).
Structural and structural meta Markov properties are closely related.

Theorem 1.4 (Byrne, 2011, Th. 4.3.2). The support of a structurally Markov graph distribution is a structurally meta Markov family.

Example 1.10. The family of spanning trees T is a structurally meta Markov family.

Directed Graphical Models & Markov Equivalence

Sections 1.1.1.b to 1.1.2.c described how graph theory on undirected graphs and conditional independence could be weaved together to produce undirected graphical models. The same results can be obtained for directed acyclic graphs (DAGs), with all the Markov properties defined for distributions, hyperdistributions or graph distributions having directed counterparts. We will not roll out the whole theory of directed graphical models. We will just define the core notion of separation on DAGs, as it is less straightforward than in the undirected case.

Definition 1.18 (Directed graph). A directed graph is a pair D = (V, E), where V is a set of vertices and the set of edges E is a subset of V × V of ordered pairs of distinct vertices. Notice that we do not allow edges from a vertex to itself. Definition 1.19 (Directed path). Let D = (V, E) be a directed graph and α, β be two vertices in V . A (directed) path from α to β is a sequence α = γ 0 , ... , γ n = β, n 1, of distinct vertices such that, for all 1 i n, (γ i-1 , γ i ) ∈ E. If there is a path from α to β in D, we say that α leads to β and write α → β. A directed cycle is a directed path from a vertex α to itself.

Definition 1.20 (DAG). A directed acyclic graph or DAG is a directed graph with no directed cycles.

The skeleton of a DAG D is the undirected graph obtained by forgetting the direction of the edges in D. For a DAG D = (V, E) and for α ∈ V , we define the following set of vertices:

parents of α: pa(α) = {β ∈ V |(β, α) ∈ E}, descendants of α: de(α) = {β ∈ V |α → β}, non-descendants of α: nd(α) = V \ (de(α) ∪ {α}).

Background

Markov properties

π Markov property with respect to G = (V, E) Similarly to the undirected case, weaker directed Markov properties can de defined on DAGs, but they are in fact equivalent to the global property without any further assumption on P . For more details, we refer the reader to [START_REF] Lauritzen | Graphical Models[END_REF]. Directed versions of the hyper Markov and structural Markov properties can also be easily derived [START_REF] Dawid | Hyper Markov Laws in the Statistical Analysis of Decomposable Graphical Models[END_REF][START_REF] Byrne | Hyper and Structural Markov Laws for Graphical Models[END_REF].

pairwise {α, β} / ∈ E ⇒ α |= β|V \ {α, β} local ∀α ∈ V , α |= V \ cl(α)|bd(α) global ∀A, B, S disjoint subsets of V , A |= G B|S ⇒ A |= B|S ρ Hyper Markov property with respect to decomposable graph G = (V, E) weak ∀(A, B) decomposition of G, π A |= π B |π A∩B strong ∀(A, B) decomposition of G, π A |= π B|A ξ Structural Markov property ∀(A, B) covering pair of G, G A |= G B |{G ∈ U(A, B)}

Meta Markov properties

Ξ A : π -→ π A Ξ A|B : π -→ π A|B F Meta Markov property with respect to decomposable graph G = (V, E) weak ∀(A, B) decomposition of G, Ξ A ‡ Ξ B |Ξ A∩B [F] strong ∀(A, B) decomposition of G, Ξ A ‡ Ξ B|A [F] G A : G -→ G A G Structural Meta Markov property ∀(A, B) covering pair of G, G A ‡ G B [G(A, B)]
Two undirected graphs induce the same graphoid if and only if they are equal. With DAGs, the situation is more complicated, as several DAGs can lead to the same graphoid. It is also interesting to determine in which cases a DAG and an undirected graph can describe the same graphoids (see Example 1.13).

Definition 1.23 (Markov equivalence). Let K and K be two graphs, each of them either an undirected graph or a DAG. K and K are said to be Markov equivalent if, for any triple (A, B, S) of disjoint subsets of V , it holds that

A |= K B|S ⇔ A |= K B|S.
Proposition 1.16 [START_REF] Frydenberg | The chain graph markov property[END_REF]. The first graph in Example 1.13 can be seen as the representative of the Markov equivalence class formed by itself and the three DAGs beside. More generally, partially directed graphs (PDG) were introduced to serve as representatives for Markov equivalence classes. As hinted by their names, these graphs are made of both directed and undirected edges. The PDG representing the Markov equivalence class of a DAG can be obtained by removing the directions of all edges that are not involved in an immorality. All the DAGs of a Markov equivalence class described by a PDG can be obtained by giving directions to undirected edges in a way that does not create additional immoralities. Consequently, the Markov equivalence class of an undirected spanning tree is made of all the DAGs that can be obtained by choosing a vertex as a root and directing all edges away from this root. The problems at hand in Chapters 2 and 3 are similar in the sense that they both require summations over large discrete sets. At first sight, these sets can seem quite hard to handle. In the first case, when we consider p vertices, we are dealing with p p-2 spanning trees. For p = 10, we are talking about a hundred million trees. When it comes to segmentations, the cardinality of the set M K made of all possible segmentations of 1; N into K segments also grows fast, as it is equal to N -1 K-1 ≈ (N /K) K . Nonetheless, in both cases, the summation can in fact be performed in polynomial time whenever the function to be summed up factorises in the right manner. The resulting algorithms heavily rely on algebraic results that we detail in the following sections.

Background

Summing over Trees

1.2.1.a Undirected Trees

We remind that V = {1, . . . , p} and that T denotes the set of spanning trees on V . Let f be a non-negative function defined on T . Our goal is to compute

Σ(f ) . . = T ∈T f (T ).
Let us assume that f can be written as

f (T ) = {i,j}∈E T ω i,j , ∀T = (V, E T ) ∈ T , (1.9)
where ω = (ω i,j ) 1 i,j p is a symmetric matrix with non-negative entries giving the weight of any possible edge. The diagonal of ω is set to zero. If the graph G ω = (V, E ω ), where

E ω = {{i, j} ∈ P 2 (V ) : ω i,j > 0},
is not connected, f ≡ 0. If this factorisation assumption is true, Σ(f ) is given by the following result.

Theorem 1.5 (Matrix-Tree theorem). Let ∆ be the Laplacian matrix associated to ω, whose general term is given by

∆ i,j = -ω i,j if i = j, k∈V ω kj if i = j.
Then all cofactors of ∆ are equal to Σ(f ).

If all edges are given a weight equal to one, this theorem yields the number of trees in T , which is p p-2 . This particular case is known as Cayley's formula. When the weight of an edge is set to either 0 or 1, it gives the number of spanning trees contained in the graph whose adjacency matrix is given by ω. This result is often referred to as Kirchhoff's theorem.

We now prove a series of lemma that will be put together in the proof of Theorem 1.5. We consider a function f of the form given in (1.9), defined by a weight matrix ω, and the associated Laplacian matrix ∆. We assume that the graph G ω = (V, E ω ) induced by the support of ω is connected.

Definition 1.24. The incidence matrix D G of a graph G = (V, E) is a p×|E| matrix, whose rows are indexed by V and columns are indexed by E. The l-th column of D G , corresponding to {i, j} ∈ E with i < j, is defined by

D G k,l =    1 if k = i, -1 if k = j, 0 otherwise.
For any matrix A, we let rk(A) denote the rank of A.

Lemma 1.1. A graph G is connected if and only if rk(D G ) = p -1.

Background

Proof. Suppose that G is connected. Each column having two non-zero entries, 1 and -1, summing all rows shows that rk(D G ) p -1. Let us consider a linear combination α of d p-1 rows of D G whose coefficients are all non-zero and let us say that this combination is equal to zero. Since G is connected, there exists an edge e j connecting a vertex corresponding to a chosen row to a vertex corresponding to a non-chosen one. There are only two non-zero entries in the j-th column of D G and only one is in the d chosen rows: the j-th entry of α cannot be equal to zero and that is absurd. So D G has rank p -1. Conversely, suppose that rk(D G ) = p -1 and that G is not connected, i.e. that G has a connected component of size d < p. By summing the rows corresponding to these d vertices, we would get a linear combination of size d equal to zero. That is not possible since rk(D G ) = p -1.

The Laplacian matrix associated to ω can be expressed using the incidence matrix D of the complete graph on V . If m = |P 2 (V )| = p(p -1)/2, we let Ω denote the diagonal matrix of size m × m giving the weights in ω in the order defined by D.

Lemma 1.2. ∆ = DΩD .

Proof.

[DΩD ] i,j is the inner product of the i-th row of DΩ with the j-th row of D. DΩ can be seen as a weighted incidence matrix for the complete graph. If i = j, the only nonzero entry in common between these two rows is in the column corresponding to the edge {i, j}. In this case, these two non-zero entries are either -ω i,j and 1 or ω i,j and -1, so that (DΩD ) i,j = -ω i,j . Similarly, if i = j, we got [DΩD ] i,i = k∈V ω ik .

Note that this result remains true whenever one consider the incidence matrix of the graph G ω associated to the support of ω instead of D. If Ω is the matrix obtained by removing the rows and columns corresponding to edges with weight zero from Ω, it holds that ∆ = (D Gω )Ω(D Gω ) .

(1.10)

We now give and prove two results on the incidence matrix that will be useful in the proof of the final result.

Lemma 1.3. Any square submatrix of D has determinant equal to 0 or 1 or -1.

Proof. The result can be proved recursively on the size of the submatrix. It is obviously true for submatrices of size 1. If 1 < k p, let us consider a square submatrix of D of size k. If all its entries are 0, its determinant is 0. Otherwise, two different cases can arise. All columns might have two non-zero entries, namely 1 and -1. Then we replace the first row by the sum of all rows and we get a null determinant. If not, there is one column with only one non-zero entry, be it 1 or -1, and we can develop our determinant according to this column, using our recursive hypothesis to conclude.

Lemma 1.4. Let U be a subset of P 2 (V ) with |U | = p-1. Let D U denote the (p-1)×(p-1)
submatrix of D consisting of the intersection of the p -1 columns corresponding to the edges in U and any set of p-1 rows of D. Then D U is invertible if and only if the graph G U = (V, U ) is a spanning tree.

Proof. Suppose that G U = (V, U ) is a spanning tree. Then the submatrix D U consists of p -1 rows of the incidence matrix

D G U of G U (D G U being a p × (p -1) submatrix of D itself). Since U is connected, the rank of D is p -1 and D U is invertible.

29

Background

Conversely suppose that D U is invertible. Then the incidence matrix D G U of G U has an invertible (p -1) × (p -1) submatrix and thus has rank p -1. It follows from Lemma 1.1 that G U is a connected graph on V , with p -1 edges. It is therefore a spanning tree.

Finally, we give the following result on the determinant of a product of two rectangular matrices of transpose shape.

Proposition 1.17 (Cauchy-Binet formula). Let n, m be two integers such that n m. Let A and B be two matrices of respective sizes n × m and m × n. For U ⊆ 1; m , we let A •U and B U • denote the matrices respectively obtained by retaining the columns indicated by U in A and the rows indicated by U in B. Then, it holds that

det(AB) = U ⊆ 1;m , |U |=n det(A •U ) det(B U • ).
Proof. See for instance (Broida & Williamson, 1989, §4.6).

We now have all the elements necessary for the proof of Theorem 1.5.

Proof of Theorem 1.5. We first show that all cofactors of ∆ are equal. Let C ∆ denote the cofactor matrix of ∆. Then ∆C ∆ = det(∆)I = 0 and every column of C ∆ is in the kernel of ∆. As G ω is connected, by Lemma 1.1, D Gω has rank p -1 and by (1.10), so has ∆. It follows that ∆ has a kernel of dimension 1. The vector (1, ..., 1) obviously belongs to that kernel. So each column of C ∆ has identical entries. But ∆ being symmetric, so is C ∆ and all cofactors are equal.

Let D denote the matrix obtained from D by removing the last row. Then det(DΩD ) is a cofactor of ∆ by Lemma 1.2 and this determinant can be expanded thanks to the Binet-Cauchy formula:

det(DΩD ) = U ⊂P2(V ), |U |=p-1 det([DΩ] U ) det(D U ) = U ⊂E, |U |=p-1 det(D U ) 2 {i,j}∈U ω i,j
where D U and [DΩ] U are respectively the square submatrices of D and DΩ whose p -1 columns correspond to the edges in U ⊂ P 2 (V ). Lemma 1. 3 and1

.4 yield that det(D U ) 2 is equal to 1 if G U = (V, U
) is a spanning tree on V and 0 otherwise, hence the result.

Therefore, computing Σ(f ) can be done in O(p 3 ) time as it only requires to compute the determinant of a (p -1) × (p -1) matrix.

Other sums that Σ(f ) might be of interest. For instance, in Chapter 2, we will actually need to sum f on the trees containing a particular edge {k, l}. The sum over the trees that do not contain edge {k, l} can be obtain by putting ω k,l = ω l,k to zero in ω and by applying the Matrix-Tree theorem to this new weight matrix, and the difference with Σ(f ) is the wanted sum.

Background

A more general version of the Matrix-Tree theorem can be given for graphs whose connected components are spanning trees on their respective sets of vertices. Such graphs include spanning trees, but also non-connected graphs. They are called forests, for a rather obvious reason.

Theorem 1.6 [START_REF] Chaiken | A Combinatorial Proof of the All Minors Matrix Tree Theorem[END_REF]. Let ∆ be the Laplacian of a weight matrix ω and let U ⊆ V . We let F U denote the set of forests on V with |U | connected components such that, for any two vertices u 1 , u 2 ∈ U , u 1 and u 2 are not in the same connected component. We also let ∆ U denote the matrix obtained from ∆ by removing the rows and columns corresponding to the vertices in U . Then it holds that

|∆ U | = F ∈F U {i,j}∈E F ω i,j .
Briefly speaking, U can be seen as a set of "roots" (even though the models are not directed) for the trees of the forests in F U . If U is taken equal to a single vertex, then the forests in F U only have one connected component which is a tree and we get Theorem 1.5.

Example 1.14. Three examples of forests in F {1,2,3} . 1 2 3 1 2 3 1 2 3
A straightforward result can be deduced from Theorem 1.6 on the summation over the trees containing a given edge {k, l}.

Corollary 1.6.1. Let {k, l} be an edge in P 2 (V ). Then it holds that

T ∈T , {k,l}∈E T {i,j}∈E T ω i,j = ω k,l |∆ {k,l} |.
Proof. There is a one to one correspondence between the trees that contain edge {k, l} and the forests in F {k,l} . Each forest in F {k,l} can uniquely be transformed in a tree by adding edge {k, l}. Reciprocally, a tree containing edge {k, l} induces a forest in F {k,l} by removing {k, l}. Therefore, using Theorem 1.6, we have that

T ∈T , {k,l}∈E T {i,j}∈E T ω i,j = ω k,l F ∈F {k,l} {i,j}∈E F ω i,j = ω k,l |∆ {k,l} |.

1.2.1.b Directed Trees

Both versions of the Matrix-Tree theorem given above can actually be stated for directed trees and forests. A directed tree is a DAG D whose skeleton is an undirected spanning trees Background and whose edges are oriented away from one root vertex. A directed forest is a DAG whose connected components are directed trees. We let -→ T and -→ F respectively denote the sets of directed trees and forests on V . The results that we are about to give will only be used in Section 4.2, when introducing temporal dependence in the modelling of multivariate time-series.

Let ω be a p × p matrix with entries in R + . The diagonal terms of ω are assumed to be zeros. For (i, j) ∈ V 2 , ω i,j gives the weight of the edge going from i to j. ω is not assumed to be symmetric. The Laplacian matrix -→ ∆ associated to ω is defined as in the undirected case by

-→ ∆ i,j = -ω i,j if i = j, k∈V ω kj if i = j.
The directed counter-part of Theorem 1.5 is given below.

Theorem 1.7 (Directed Matrix-Tree theorem). If we let C denote the cofactor matrix of -→ ∆ and -→ T r , r ∈ V , denote the set of directed trees rooted at r, then it holds that

-→ Z r (ω) . . = D∈ -→ Tr (i,j)∈E D ω i,j = C r,r , ∀r ∈ V.
If ω is taken to be symmetric, we find that all -→ Z r (ω) are equal, as stated by the undirected version of the theorem.

Theorem 1.6 also admits an adaptation to directed forests. Stating this result, which is actually the one proved by [START_REF] Chaiken | A Combinatorial Proof of the All Minors Matrix Tree Theorem[END_REF], requires some definitions. The signature of A ⊆ V is defined by

(A) . . = (-1) n(A) , n(A) . . = |{{i, j} ∈ (V \ A) × A : i < j}|. Let A, B be two subsets of V such that |A| = |B|. We let -→ F A,B denote the set of directed forests on V such that F ∈ -→ F A,B if and only if • F contains exactly |A| = |B| directed trees,
• each tree in F contains exactly one vertex in A and one vertex in B,

• each edge in F is oriented away from the vertex in B of the tree containing that edge. -→ F A,A is the set of forests with skeleton in F A , the undirected forests rooted in A, and in which each edge is oriented away from the vertex in A of the tree containing that edge.

Each forest F in -→ F A,B defines a bijective map π F from A to B so π F (j) = i if and only if i and j are in the same tree of F . A pair {i, i } ∈ A is an inversion of π F if i < i and π F (i) > π F (i ).
The number of inversions in π F is denoted by n(π F ) and the signature of π F is defined by

(π F ) . . = (-1) n(π F ) .
For every F ∈ -→ F A,A , π F is the identity and (π F ) = 1.

32

Background

Theorem 1.8 (All Minors Matrix-Tree theorem, [START_REF] Chaiken | A Combinatorial Proof of the All Minors Matrix Tree Theorem[END_REF]. Let -→ ∆ be the directed Laplacian matrix associated to a directed edge weight matrix β. Let A, B be two subsets of V such that |A| = |B|. We let -→ ∆ A,B denote the matrix obtained from -→ ∆ by removing the rows corresponding to the vertices in A and the columns corresponding to the vertices in B. For any DAG D = (V, E), we define

ω D . . = (i,j)∈E ω i,j .
Then, it holds that

| -→ ∆ A,B | = (A) (B) F ∈F A,B (π F )ω F .
Theorem 1.8 can be used to sum over the directed trees rooted in a vertex r and containing a given edge (i, j).

Corollary 1.8.1. Let r be a vertex in V and (u, v) be a directed edge such v = r. Then, it holds that v},{r,v} .

D∈ -→ Tr, (u,v)∈E D ω D = ω u,v | -→ ∆ {u,v},{r,v} | where | -→ ∆ {u,v},{r,v} | is the absolute value of the determinant of -→ ∆ {u,
Proof. Let A = {u, v} and B = {r, v}. Note that, as v = r and u = v, neither set can be reduced to a single vertex. Let D be a tree rooted in r and containing (u, v). Removing (u, v), we are left with two directed trees. Vertex v cannot be in the connected component of r. Otherwise, the last edge in the directed path from r to v and (u, v) would form a v-structure at vertex v. Therefore, the unique directed path between r and v goes through (u, v) and v and r are in two distinct connected components when (u, v) is removed. This also implies that in the connected component containing v, all edges are oriented away from v. The forest obtained by removing (u, v) from D thus belongs to -→ F A,B . Reversely, each forest in -→ F A,B can be converted to a tree rooted in r and containing edge (u, v) by adding (u, v). There is a one-to-one mapping between these two sets. Now, for all F ∈ -→ F A,B , note that π F is always such that π F (v) = v and π F (u) = r, so that (π F ) has the same value no matter F , denoted * . Using Theorem 1.8, we obtain that

D∈ -→ Tr, (u,v)∈E D ω D = ω u,v D∈ - → F A,B ω D = * (A) (B)ω u,v | -→ ∆ A,B | = ω u,v | -→ ∆ A,B | .

Summing over Segmentations

For 0 < i < j, we let i; j and i; j respectively denote the sets {i, . . . , j} and {i, . . . , j -1}.

Background

Definition 1.25. For 0 < i < j, a segmentation m of i; j is a partition of {i, . . . , j -1} into sets of consecutive elements, called segments. Thus, if m has K 0 segments, it can be written as

m = { τ k ; τ k+1 } K k=1 = {r k } K k=1 , with i = τ 0 < τ 1 < • • • < τ K-1 < τ K = j. {τ 1 , . . . , τ K-1 } are called the change-points of m.
For K 1, we let M K ( i; j ) denote the set made of all the segmentations of i; j into K segments. Let N 0. We let M K be a shorthand for M K ( 1; N + 1 ). Let f be a non-negative function defined on M K . Our goal is to compute

Σ(f ) . . = m∈M K f (m).
Let us assume that f can be written as

f (m) = i;j ∈m A i,j , ∀m ∈ M K , (1.11)
where A is a strictly upper triangular matrix of size (N + 1) × (N + 1). Matrix A basically gives the weight of all possible segments. Function f can therefore be extended to M k ( i; j ) for all k 0 and any segment i; j ⊆ 1; N + 1 .

Theorem 1.9 [START_REF] Rigaill | Exact posterior distributions and model selection criteria for multiple change-point detection problems[END_REF]. Let f be a function of the form given in (1.11) defined by a strictly upper triangular matrix A. Then, for k 1, and i; j ⊆ 1; N + 1 , it holds that

m∈M k ( i;j ) f (m) = [A k ] i,j .
The proof of Theorem 1.9 relies on the following lemma.

Lemma 1.5. Let A be a square matrix of size n × n, n 1. For k ∈ N * , we define f A,k as

f A,k (i, j) = t0=i, t k =j (t1,...,t k-1 )∈ 1;n k-1 k-1 s=0
A ts,ts+1 .

(1.12)

Then, it holds that

f A,k (i, j) = [A k ] i,j . Proof. f A,k (i, j) = [A k ] i,j holds for k = 1. Suppose that f A,k (i, j) = [A k ] i,j holds for k ∈ N * . For k + 1, we have that f A,k+1 (i, j) = t0=i, t k+1 =j (t1,...,t k )∈ 1;n k k s=0 A ts,ts+1 = n t=1 t0=i, t k =t (t1,...,t k-1 )∈ 1;n k-1 k-1 s=0 A ts,ts+1 = n t=1 f A,k (i, t)A t,j .

34

Background

Using our induction hypothesis and by definition of the matrix product, we obtain:

f A,k+1 (i, j) = n t=1 [A k ] i,t A t,j = [A k+1 ] i,j .
(1.13)

Proving Theorem 1.9 from Lemma 1.5 is straightforward.

Proof of Theorem 1.9. Whenever A is a strictly upper diagonal matrix, all the terms in f A,k (i, j) that do not correspond to a segmentation of i; j , i.e. for which we don't have

i = t 0 < t 1 < • • • < t k-1 < t k = j, are equal to zero. Indeed, if (t 0 , . . . , t k
) is not a segmentation of i; j , there exist s ∈ 0; k -1 such that t s t s+1 and A ts,ts+1 = 0 as A is strictly upper diagonal. Thus, we have that

f A,k (i, j) = m∈ i;j ts;ts+1 ∈m A ts,ts+1 = m∈M( i;j ) f (m),
and Lemma 1.5 gives the sought-after result.

Concerning complexity, (1.13) shows that, for i ∈ 1; N + 1 and K N , quantities in

{f A,k (i, j)} 1 i N 1 k K
can all be computed in O(KN 2 ) time as the i-th row of matrices A, A 2 , . . . , A K . The same quantities where j is fixed are obtained in O(KN 2 ) as the j-th column of matrices A, A 2 , . . . , A K .

2 Our angle on network inference is to assume that the structure of the underlying graphical model is a spanning tree. The chances are that the structure of interest do not belong to this particular class of graph. For this reason, instead of looking out for the most likely tree structure, we actually compute the posterior probability for any given edge to be borrowed by a random tree. This approach can somehow be justified by the fact that a graph can be close to a tree on a local scale, while having cycles or not being connected. Averaging over all spanning trees, we might therefore be able to derive meaningful information about local features, such as edges, while blurring out the strong global constraints imposed by individual trees. Having restrained our attention to spanning trees, the matrix containing the posterior probabilities of every possible edges can be computed in cubic time with respect to the number of vertices. This means that, from a complexity point of view, our approach could be used on problems involving up to a few hundreds variables in a reasonable amount of time. This chapter is available as a preprint on arXiv (Schwaller et al., 2015). Statistical models are getting more and more complex and can now involve very intricate dependency structures. Graphical models are both a natural and powerful way to depict such structures. Inferring a graphical model based on observed data is hence of great interest for many fields of applications. From a statistical point-of-view, considering the inference of a graphical model requires to consider the graphical model itself as a parameter. In a Bayesian context, it means that we have to define a full model and, more specifically, a prior distribution on graphical models, therefore on graphs themselves.

Regardless of whether we consider directed or undirected graphs, their sheer number make them difficult to deal with. Exact inference can only be contemplated as long as there are no more than thirty or so variables of interest [START_REF] Parviainen | Exact Structure Discovery in Bayesian Networks with Less Space[END_REF]. When exact inference is no longer tractable, sampling is used as a pragmatic alternative. Markov Chain Monte Carlo (MCMC) methods have for instance been used to sample from some sets of graphs, such as Directed Acyclic Graphs (DAGs) [START_REF] Madigan | Bayesian graphical models for discrete data[END_REF][START_REF] Friedman | Being Bayesian about network structure. A Bayesian approach to structure discovery in Bayesian networks[END_REF][START_REF] Niinimäki | Partial order mcmc for structure discovery in bayesian networks[END_REF] or decomposable graphs [START_REF] Green | Sampling decomposable graphs using a markov chain on junction trees[END_REF]. The decomposability assumption for undirected graphical models, also called Markov random fields, is commonly made in the literature, although some interest has been devoted to the less easy to handle non-decomposable graphs [START_REF] Roverato | Hyper inverse wishart distribution for non-decomposable graphs and its application to bayesian inference for gaussian graphical models[END_REF][START_REF] Atay-Kayis | A Monte Carlo method to compute the marginal likelihood in non decomposable graphical Gaussian models[END_REF]. The sampling schemes developed in the aforementioned papers are often subject to standard issues related to MCMC sampling in high-dimensional spaces, namely slow mixing and difficulty to get to the stationary distribution.

One way to bypass theses hurdles is to further restrict the exploratory space so as to make exact inference tractable. When a subset of graphs is considered, it sometimes becomes possible to get access to the full posterior distribution on graphs. The obvious drawback of this approach is that the "true" graph might not belong to this subset. In this case, computing a maximum a posteriori (MAP) estimate would for instance yield a systematically wrong answer. But usually, such methods are not intended to assess the global structure all at once but to separately assess a collection of local features of the graph (typically, edges). The idea is that the inference of such features is less affected by the restriction than the global structure. In that perspective, trees have been of particular interest as a subset of both decomposable graphs and DAGs [START_REF] Chow | Approximating Discrete Probability Distributions with Dependence Trees[END_REF][START_REF] Meilă | Learning with Mixtures of Trees[END_REF][START_REF] Meilă | Tractable bayesian learning of tree belief networks[END_REF][START_REF] Kirshner | Learning with Tree-Averaged Densities and Distributions[END_REF][START_REF] Lin | Learning Sparse Markov Network Structure via Ensemble-of-Trees Models[END_REF][START_REF] Burger | Disentangling direct from indirect co-evolution of residues in protein alignments[END_REF].

Our first contribution is to provide a well-defined fully Bayesian framework for graphical model inference based on trees. We use the work of [START_REF] Dawid | Hyper Markov Laws in the Statistical Analysis of Decomposable Graphical Models[END_REF] on hyper Markov laws to define priors on tree parameters and distributions that can easily be marginalized over. This framework spares us from requiring likelihood equivalence between Markov equivalent directed tree models, like [START_REF] Meilă | Tractable bayesian learning of tree belief networks[END_REF] did building on the work of [START_REF] Hammersley | Learning bayesian networks: The combination of knowledge and statistical data[END_REF]. We also point out that it fits within the recent work of [START_REF] Byrne | Structural markov graph laws for bayesian model uncertainty[END_REF] on structurally Markov graph distributions. We then go through a series of typical models befitting this framework, namely tree-structured copulas [START_REF] Kirshner | Learning with Tree-Averaged Densities and Distributions[END_REF], multinomial distributions [START_REF] Meilă | Tractable bayesian learning of tree belief networks[END_REF] and Gaussian distributions. Bayesian inference in this framework requires integration over the set of trees, that can be carried out exactly and efficiently thanks to an algebraic result called the Matrix-Tree theorem.

Our second contribution focuses on edge inference. When [START_REF] Meilă | Tractable bayesian learning of tree belief networks[END_REF] and [START_REF] Kirshner | Learning with Tree-Averaged Densities and Distributions[END_REF] were interested in the joint distribution of the observations, we are interested in the inference of the dependence structure. To this purpose, we are not concerned with the inference of the parameters but we need to account for the uncertainty of their estimates. The Bayesian construction we propose provides a natural framework to achieve this. We derive the exact posterior probability of any given edge in an exact manner, allowing for an arbitrary prior edge appearance probability.

Most works on tree-structured graphical model inference rely on the aforementioned Matrix-Tree theorem. As noticed by [START_REF] Kirshner | Learning with Tree-Averaged Densities and Distributions[END_REF], the computation of posterior probabilities for all the edges in this setting can be achieved with cubic complexity with respect to the number of variables. We provide a new proof of this result relying on a generalization of the Matrix-Tree theorem to forests.

Our last contribution is a simulation study which addresses the influence of the tree assumption on the accuracy of structure inference for non-tree-structured graphical models. We demonstrate that, as long as edge inference is concerned, the computational efficiency following from this assumption can be obtained at a limited cost.

An R-language package saturnin implementing the approach presented here is available from the Comprehensive R Archive Network.

In Section 2.2, we provide some background on graphical models and Markov properties before writing down the full model in which the inference is performed. Priors for tree structures and distributions are defined in Section 2.3. Section 2.4 deals with the inference of the model. Integrations with respect to distributions and structures are respectively discussed in Sections 2.4.1 and 2.4.2. The simulation study and its results are described in Section 2.5. An application to flow cytometry data is presented in Section 2.6.

Background

Markov Properties & Graphical Models

Let V = {1, ..., p} and X = (X 1 , ..., X p ) be a random vector taking values in a product space X = p i=1 X i . The set of distributions on X is denoted by F. For any subset A of V , X A denotes the subvector of X corresponding to the components in A. Let P 2 (V ) denote the subsets of V of size 2. For E ⊆ P 2 (V ), G = (V, E) is the undirected graph with vertices V and edges E. In the following, the notations of [START_REF] Dawid | Hyper Markov Laws in the Statistical Analysis of Decomposable Graphical Models[END_REF] will be used. We refer the reader to the appendix of their article for a quick introduction to graph terminology and graphical models, or to [START_REF] Lauritzen | Graphical Models[END_REF] for a more detailed overview.

A pair (A, B) of subsets of V is said to be a decomposition of G if V = A ∪ B, the subgraph induced by G on A ∩ B is complete and A ∩ B separates A from B. If A and B are both proper subsets of V , the decomposition is said to be proper. Here we restrain our attention to decomposable graphs, namely graphs that are either complete or for which there exists a proper decomposition into two decomposable subgraphs.

Definition 2.1. A distribution π ∈ F is said to be Markov with respect to (w.r.t.) a decomposable graph G if, for any decomposition

(A, B) of G, under π, X A |= X B |X A∩B .
Proposition 2.1. [START_REF] Hammersley | Learning bayesian networks: The combination of knowledge and statistical data[END_REF]) Let π ∈ F . If π is a positive distribution (for all x ∈ X , π(x) > 0), being Markov w.r.t. a decomposable graph G is equivalent to the existence of a factorisation of π on the (maximal) cliques of G.

Here we will focus on distributions that are Markov w.r.t. to connected graphs without any cycles. Such graphs are called spanning trees and their maximal cliques are of size 2. Thus, a positive distribution that is Markov w.r.t. a tree T = (V, E T ) can be factorised on the edges of the tree, using the marginal distributions of order 1 and 2:

∀x ∈ X , π(x) = i∈V π i (x i ) {i,j}∈E T π ij (x i , x j ) π i (x i )π j (x j )
.

Such distributions will be called tree distributions in the following.

Definition 2.2. A graphical model m G . . = (G, F G ) is given by a decomposable graph G and a family of distributions F G ⊆ F that are Markov w.r.t. G.
Let m G = (G, F G ) be a graphical model. To avoid any confusion, distributions on a set of distributions will be called hyperdistributions. If π ∈ F G and ρ is a hyperdistribution on F G , for any A, B ⊆ V , π A denotes the marginal distribution obtained from π on the variables X A and π B|A the collection of conditional distributions of X B |X A under π. We also denote ρ A the marginal hyperdistribution induced by ρ on π A and ρ B|A the collection of hyperdistributions induced by ρ on π B|A .

Definition 2.3. ρ is said to be strong hyper Markov w.r.t. G if, for any decomposition

(A, B) of G, under ρ, π A |= π B|A .
Such hyperdistributions will be useful to define priors on distribution spaces.

Model for Bayesian Inference of Graphical Models Based on Trees

Let T denote the set of spanning trees on V . For any tree T ∈ T , we consider a graphical model m T = (T, F T ) with a family of positive distributions F T ⊆ F Markov w.r.t. T . Here we consider a Bayesian framework. We therefore need to define prior distributions for T and for π conditional on T . This is dealt with in Section 2.3. The full Bayesian model consists in first drawing a random tree T * , then a distribution π in F T and finally X according to π (Figure 2.1). Defining a prior on tree distributions could be especially troublesome since it needs to be defined for every graphical model m T . The idea is to require these hyperdistributions to be strong hyper Markov w.r.t. to their trees, so that they can be built from local hyperdistributions defined on the edges and chosen once and for all trees. This choice of prior and the fact that we only consider tree structures for the graphical models make the inference of the graph in our model tractable in an exact manner, thanks to the Matrix-Tree theorem. Restraining the explored set of graphs to the spanning trees obviously helps to make the inference easier to perform, but it still leaves us with a super-exponential number, p p-2 , of graphs. Nonetheless, a suitable choice of priors on tree structures and parameters leads to a tractable situation. [START_REF] Meilă | Tractable bayesian learning of tree belief networks[END_REF] defined what they call decomposable priors under which parameters can be dealt with at the edge level. The integration over the set of trees can then be performed exactly thanks to algebra. We will use strong hyper Markov hyperdistributions [START_REF] Dawid | Hyper Markov Laws in the Statistical Analysis of Decomposable Graphical Models[END_REF] to define our prior but the idea is basically the same. Let D = (x (1) , ..., x (n) ) be an independent sample of size n 1 drawn from X. Our goal is to define a prior distribution on (T, θ) such that the posterior distribution on trees ξ(•|D) factorises over the edges, i.e.

Priors on

ξ(T |D) = 1 Z {i,j}∈E T ω i,j , ∀T ∈ T , (2.1)
where ω = (ω i,j ) (i,j)∈V 2 is a symmetric matrix with non-negative values and

Z = T ∈T {i,j}∈E T ω i,j (2.2)
is a normalising constant. Both ω and Z obviously depend on the data D but we drop the dependence in the notations for sake of clarity.

Prior on Tree Structures

Let β = (β i,j ) (i,j)∈V 2 be a symmetric matrix with non-negative values such that its support graph G β = (V, E β ), where E β = {{i, j} ∈ P 2 (V ), β i,j > 0}, is connected. We consider a prior distribution ξ on T that factorises over the edges,

ξ(T ) = 1 Z 0 {i,j}∈E T β i,j , ∀T ∈ T . (2.3)
The assumption about β is here to serve as a guarantee that β induces a proper distribution on trees; ξ can typically be taken as a uniform on T .

These distributions belong to the family of structurally Markov graph distributions described by [START_REF] Byrne | Structural markov graph laws for bayesian model uncertainty[END_REF] (see Section 2.3.3). (2.4)

Prior on Tree Distributions

For every T ∈ T , we have to define a prior distribution on F T such that the marginal likelihood p(D|T ) can also be factorised on the edges. [START_REF] Meilă | Tractable bayesian learning of tree belief networks[END_REF] built their prior on multinomial tree distributions around three main assumptions, namely likelihood equivalence, parameter independence and parameter modularity. The first assumption requires that the prior treats all possible para-metrizations consistent with a given tree T (be it directed or undirected) as indistinguishable. As we only consider undirected parametrizations in our construction, we shall not need this assumption here. As for the parameter independence assumption, it can be broken down into local and global independences [START_REF] Spiegelhalter | Sequential updating of conditional probabilities on directed graphical structures[END_REF]. Strong hyper Markov hyperdistributions satisfy global independence but not necessarily local independence. The latter is in fact not needed to get the desired factorisation property for the marginal likelihood. Finally, the parameter modularity assumption is ensured by the construction of a compatible family of strong hyper Markov hyperdistributions.

Let T be a tree and ρ T a strong hyper Markov hyperdistribution on F T . Such hyperdistributions have an interesting property regarding the marginal likelihood p(D|T ).

Proposition 2.2. (Dawid & Lauritzen, 1993, Prop. 5.6 This means that the marginal likelihood can be factorised on the edges of T . For i ∈ V , let D i = {x 1 i , ..., x n i } be the observed data restricted to X i . The integral given in (2.4) can then be rewritten as

p(D|T ) = p(D|π)p(π|T )dπ = π(D)ρ T (π)dπ = i∈V p(D i |T ) {i,j}∈E T p(D i , D j |T ) p(D i |T )p(D j |T ) (2.5)
where, for all (i, j)

∈ V 2 , p(D i , D j |T ) = π ij (D i , D j )ρ T ij (π ij )dπ ij ; (2.6) p(D i |T ) = π i (D i )ρ T i (π i )dπ i . (2.7)
The calculation of these integrals will be addressed in Section 2.4.1.

We now explain how to choose ρ T for all T so that the distributions of {π ij } {i,j}∈P2(V ) do not depend on T . Let us consider a general hyperdistribution ρ on F such that, for any A ⊆ V and under ρ,

π A |= π V \A|A . (2.8)
This means that ρ is strong hyper Markov w.r.t. the complete graph over V .

Proposition 2.3. (Dawid & Lauritzen, 1993, §6.2) For any tree T ∈ T , there exists a unique hyperdistribution ρ T on F T that is strong hyper Markov w.r.t. T and such that, for every edge {i, j} ∈ E T ,

ρ T ij = ρ ij .
(2.9)

{ρ T } T ∈T is said to be a (hyper) compatible family of strong hyper Markov hyperdistributions.

Proposition 2.3 guarantees that all ρ T are strong hyper Markov w.r.t. T . By Proposition 2.2, for all T ∈ T , the marginal likelihood under ρ T is Markov w.r.t. T . Moreover, the compatibility of the family {ρ T } T ∈T makes the dependence on T in the local marginal distributions given in (2.6) and (2.7) irrelevant. They can be computed once and for all for every {i, j} ∈ P 2 (V ). This choice of priors for the distributions assures that (2.1) is satisfied with

ω i,j = β i,j p(D i , D j ) p(D i )p(D j ) , ∀(i, j) ∈ V 2 .
(2.10)

The model under which we are now working is fully described in Figure 2.2.

Proposition 2.3 shows that we do not need to have access to the full basis hyperdistribution to specify a compatible family of strong hyper Markov hyperdistributions. It is indeed enough to provide a consistent family of pairwise hyperdistributions {ρ ij } P2(V ) , where the consistency property must be understood in the sense that two hyperdistributions involving a common vertex should induce the same marginal hyperdistribution on this vertex. This is automatically satisfied when {ρ ij } P2(V ) is obtained from a fully specified hyperdistribution ρ. In order to obtain strong hyper Markov hyperdistributions when combining these pairwise hyperdistributions, we shall additionally require that, for all i, j ∈ V , π i|j |= π j under ρ ij (Dawid & Lauritzen, 1993, Prop. 3.16), meaning that ρ ij is strong hyper Markov w.r.t. the graph on {i, j} where vertices i and j are connected.

Structural Markov Property and Structurally Meta Markov Families

The purpose of this section is to show how the model that we have described so far is related to the structural Markov property defined by [START_REF] Byrne | Structural markov graph laws for bayesian model uncertainty[END_REF]. Indeed, trees have specific algebraic properties that will be taken advantage of in Section 2.4 for the inference of the model, but the model itself can be extended to other subsets of decomposable graphs. [START_REF] Byrne | Structural markov graph laws for bayesian model uncertainty[END_REF] defined an extension of the (hyper) Markov properties described in [START_REF] Dawid | Hyper Markov Laws in the Statistical Analysis of Decomposable Graphical Models[END_REF] to undirected decomposable graphs (and to directed acyclic graphs, but this will not be discussed here) called the structural Markov property.

Let U be the set of undirected decomposable graphs on V . A pair of subsets (A, B) of V is called a covering pair if A ∪ B = V . For any family of graphs G ⊆ U and for any covering pair (A, B), we define G(A, B) to be the set of graphs G ∈ G for which (A, B) is a decomposition. Definition 2.4. A distribution ξ for G ∈ U is said to be structurally Markov if for any covering pair (A, B) such that ξ(U(A, B)) > 0,

T * ∼ ξ; π ∼ ρ T ; X ∼ π.
G A |= G B |{G ∈ U(A, B)} (2.11) under ξ.
The support graph families of structurally Markov graph distributions have structural properties. They satisfy the so called structural meta Markov property.

Definition 2.5. Let G be a family of undirected decomposable graphs on V . Then G is structurally meta Mar-kov if for any covering pair (A, B), the set

{G A |G ∈ G(A, B), G B = J} is the same for all J ∈ {G B |G ∈ G(A, B)}.
The set of trees T is an example of such a family (Byrne & Dawid, 2015, Ex. 3.1) and the distributions that we considered in Section 2.3.1 are structurally Markov.

These graph distributions naturally interact with (strong) hyper Markov hyperdistributions and Markov distributions when they are chosen carefully. Compatible hyperdistribution families as described in Proposition 2.3 conjugate nicely with graph distributions factorised on the edges, so that all hyperdistribution updates can be performed down to the edge level. But compatibility can be defined for any structurally meta Markov family G (Byrne & Dawid, 2015, Definition 3.4). Then, the update can be performed locally on

C G = G∈G C G where C G denotes the cliques of graph G.
We finish this section by laying stress upon the fact that, among structurally meta Markov graph families, trees are of particular computational interest given their algebraic properties. One of the main difficulties in assessing graph distributions is to compute normalisation constants, but closed-form expressions can be derived for these constants in the case of trees (see Section 2.4.2). Different inference tasks can be performed on graphical models. One might be interested in estimating the emission distribution of X. [START_REF] Chow | Approximating Discrete Probability Distributions with Dependence Trees[END_REF] gave an algorithm to get the tree distribution maximizing the likelihood of discrete multivariate data in the frequentist equivalent of the model described in the previous section. It can easily be adapted to MAP estimation in a full Bayesian framework [START_REF] Meilă | Learning with Mixtures of Trees[END_REF]. It is also possible to look at the posterior predictive distribution [START_REF] Meilă | Tractable bayesian learning of tree belief networks[END_REF],

Inference in

p(x|D) = T ∈T p(x|T )ξ(T |D).
In some other situations, the structure of dependence between the variables, that is the graph G, might be the only object of interest. [START_REF] Lin | Learning Sparse Markov Network Structure via Ensemble-of-Trees Models[END_REF] were for instance interested in the probability of an edge appearing in a tree. They looked out for the matrix β maximizing the likelihood of the data under a mixture of all possible tree models, where the probability of a tree model is defined just as in (2.3). The parameters of the models are estimated with plug-in estimators. The distribution on trees cannot be called a prior in the traditional sense but the likeness to the model that we have described is obvious.

Here we are also interested in the probability for edges to appear in a tree, but in a full Bayesian framework. Formally, we would like to compute P ({k, l} ∈ E T * |D, ξ) for any edge {k, l},

P ({k, l} ∈ E T * |D, ξ) = T ∈T , E T {k,l} ξ(T |D).
(2.12)

The previous section shows that achieving this requires two things. First, we have to get access to ω by computing local marginal likelihoods, which amounts to integrating w.r.t. π (Section 2.4.1). Then comes in the integration over the set of trees, that can be performed exactly thanks to an algebra result called the Matrix-Tree theorem (Section 2.4.2).

Integration with Respect to π

Thanks to the strong hyper Markov property required for the hyperdistributions, the integration on π can be performed locally and the compatibility ensures that these local integrated quantities can be passed from one model to another whenever they are needed. Thus, the integrations are always made on sets of bivariate distributions, with p(p+1) 2 of them to be computed. The small dimension of each of the involved problems makes it possible to consider numerical or Monte Carlo integration. We begin by describing a framework based on tree-structured copulas where it might be needed, depending on the choice of local copulas. We then present two models where the local integrated likelihood terms can even be computed exactly thanks to conjugacy.

2.4.1.a Tree-Structured Copulas

Let us assume that X = [0, 1] p . If we make the assumption that the marginal distribution of each variable is uniform, the joint distribution for X is called a copula. Here we are interested in a subset of these distributions called the tree-structured copulas [START_REF] Kirshner | Learning with Tree-Averaged Densities and Distributions[END_REF]. We denote by U the uniform distribution on [0, 1]. Let us assume that, for all i ∈ V , X i ∼ U. We are basically considering a copula model where the marginal data distributions have been dealt with in a relevant manner, independently from our model. For any i ∈ V , the marginal hyperdistribution ρ i for π i is then a Dirac distribution concentrated on U, denoted δ U . Defining a compatible family of hyperdistributions requires that we consider pairwise hyperdistributions whose marginals are concentrated on U. Such hyperdistributions are in fact defined on bivariate copulas.

As an example, we consider the particular class of Archimedean copulas [START_REF] Nelsen | An Introduction to Copulas[END_REF]. Such copulas have simple expressions for their cdf. Let ψ : [0, 1] → R + ∪ {∞} be a continuous, strictly decreasing function such that ψ(1) = 0. Its pseudo-inverse ψ [-1] : R + ∪ {∞} → [0, 1] is the continuous function defined by

∀t ∈ R + ∪ {∞}, ψ [-1] (t) = ψ -1 (t) if 0 t ψ(0), 0 otherwise. Let us remark that if ψ(0) = ∞, ψ [-1] = ψ -1 .
The cdf of the Archimedean copula generated by ψ is given by

C ψ (x i , x j ) = ψ [-1] (ψ(x i ) + ψ(x j )).
ψ is said to be a generator of the copula C ψ . There is an extensive list of commonly used families of generators, many of them being governed by one or more parameters. Once again, we refer the reader to [START_REF] Nelsen | An Introduction to Copulas[END_REF] for a detailed list of such generators. We can mention the well-known Gumbel copulas for instance, whose generator and inverse generator are given by ) regulating the strength of the dependence.

ψ θ (x) = (-log(x)) θ ∀x ∈ [0, 1], ψ -1 θ (t) = exp(-t 1/θ ) ∀t ∈ R + ∪ {∞}, with θ ∈ [1, ∞
Let {i, j} be a given edge. If we consider an identifiable parametric family of Archimedean copulas {C θ } θ∈Θ , Θ ⊆ R, defined by parametric generators {ψ θ } θ∈Θ , there is a one-to-one mapping Υ between θ and the distributions π ij on (X i , X j ). A pairwise hyperdistribution ρ ij for π ij is then easily defined by any distribution κ for θ through the identity

ρ ij (π ij ) = κ Υ -1 (π ij )
and the integrated pairwise distribution p(x i , x j ) is given by

p(x i , x j ) = Θ ∂ 2 C θ ∂x i ∂x j (x i , x j )κ(θ)dθ, ∀(x i , x j ) ∈ [0, 1] 2 .
(2.13) Such a family of pairwise hyperdistributions is bound to be consistent since all marginals are equal to δ U . Morever, the global hyperdistributions that we obtain from this family are strong hyper Markov since it holds that, for i, j ∈ V , π i|j |= π j under ρ ij . The integrals given in (2.13) shall be computed exactly or through numerical integration depending on the choice of the copula family. This choice needs not be the same for all the edges. In the case of the Gumbel copula, a numerical or Monte Carlo integration is required. Obviously, bivariate Gaussian copulas would also be a valid choice. The pairwise hyperdistributions could then be specified through Wishart distributions for the precision matrices of the copulas, just like in the full Gaussian case described in Section 2.4.1.c.

2.4.1.b Multinomial Distributions

We now consider the case where all X i are discrete, taking their values in finite spaces X i of size r i respectively. Let X be the Cartesian product of spaces X i . A distribution for X is given by a probability vector θ in

Θ = θ ∈ [0; 1] |X | x∈X θ(x) = 1 . (2.14)
Θ is the set of multinomial distributions on X . It happens that the conjugate Dirichlet distribution is satisfying the condition given in (2.8) that is necessary to build a compatible family of strong hyper Markov hyperdistributions. Let λ = (λ(x)) x∈X be a family of positive numbers indexed by X . We denote D(λ) the Dirichlet distribution for θ ∈ Θ, with density

f (θ|λ) ∝ x∈X θ(x) λ(x)-1 .
Proposition 2.4. (Dawid & Lauritzen, 1993, Lemma 7.2) Let θ ∼ D(λ). Then for all

A ⊆ V and B = V \ A, i. θ A ∼ D(λ A ); ii. θ A |= θ B|A ; with λ A (x A ) = y,y A =x A λ(y) for all x A ∈ X A .
All these properties result from the fact that, if {Y k } K k=1 are independent random variables distributed as Γ(λ k , θ) respectively and

V = K k=1 Y k , (Y 1 /V, ..., Y K /V ) ∼ D(λ).
(ii) assures that any λ gives rise to a hyperdistribution ρ on the multinomial family of distributions from which we can build a family of compatible strong hyper Markov hyperdistributions. (i) states that the marginal hyperdistributions are also Dirichlet distributed. The conjugacy can then be used locally to compute ω. These hyperdistributions were referred to as hyper-Dirichlet laws in (Dawid & Lauritzen, 1993, §7.2.2).

As mentioned in Section 2.3.2, specifying a full set of hyperparameters λ is in fact not necessary to define the family of hyperdistributions {ρ T } T ∈T . We only need a consistent family of {λ i,j } (i,j)∈V 2 , in the sense that, for (i, j, k) ∈ V 3 , λ i,j and λ ik should induce the same λ i . A possibility is to set the prior hyperparameters on the edges thanks to an equivalent sample size N ,

λ i,j . . = N /r i r j , λ i . . = N /r i . (2.15)
If all X i are of equal size r, a possibility is to set N = r 2 /2 so that all λ i,j are equal to 1/2 to mimic Jeffreys priors for the bivariate distributions on the edges. However, this choice will not induce global Jeffreys priors, which are not hyper-Dirichlet hyperdistributions [START_REF] York | Bayesian methods for estimating the size of a closed population[END_REF]. For an edge {i, j}, we denote λ i,j the updated hyperparameters for the edge {i, j}:

λ i,j ( , ) = λ i,j ( , ) + n k=1 δ x k i , δ x k j , ∀( , ) ∈ X i × X j ,
where δ x, = 1 if x = and 0 otherwise.

The matrix ω defined in (2.10) is then given by [START_REF] Meilă | Tractable bayesian learning of tree belief networks[END_REF])

ω i,j =β i,j ∈Xi Γ(λ i ( )) Γ(λ i ( )) ∈Xj Γ(λ j ( )) Γ(λ j ( )) ( , )∈Xi×Xj Γ(λ i,j ( , )) Γ(λ i,j ( , )) (2.16)
where Γ denotes the gamma function. If R = max i∈V r i , computing ω requires O(np 2 R 2 ) operations [START_REF] Meilă | Tractable bayesian learning of tree belief networks[END_REF].

Let us finish this section by a remark on parameter independence. The following property of the Dirichlet distribution can be added to Proposition 2.4 even though it is not used here.

Proposition 2.5. (Dawid & Lauritzen, 1993, Lemma 7

.2) Let θ ∼ D(λ). Then for all A ⊆ V and B = V \ A, θ B|A (•|x A ) are all independent and distributed as D(λ B|A (.|x A )) with λ B|A (x B |x A ) = λ(x) for all x ∈ X (up to a rearrangement of the components of x).
Thus, although not required here, the local independence assumption made by [START_REF] Meilă | Tractable bayesian learning of tree belief networks[END_REF] is in fact satisfied. In the multinomial case, [START_REF] Geiger | A Characterization of the Dirichlet Distribution Through Gloabl and Local Parameter Independence[END_REF] even showed that, together with likelihood equivalence, global parameter independence and parameter modularity, the local parameter independence assumption constrains the prior to be locally Dirichlet distributed.

2.4.1.c Gaussian Distributions

Whenever X is real-valued, one might work under the assumption that X is Gaussiandistributed with mean µ and inverse covariance matrix Λ. The conjugate normal-Wishart distribution is then a natural choice of prior for (µ, Λ). The normal-Wishart distribution is denoted by nW(ν, λ, α, Ψ) and is hierarchically defined by

Λ ∼ W(α, Ψ), µ|Λ ∼ N (ν, (λΛ) -1 ),
where W(α, Ψ) is the Wishart distribution with α > p -1 degrees of freedom and positivedefinite parametric matrix Ψ. [START_REF] Geiger | Parameter priors for directed acyclic graphical models and the characterization of several probability distributions[END_REF] showed that the normal-Wishart distribution satisfy the parameter independence property given in (2.8). They further proved that this property coerces the distribution to be normal-Wishart when p 3. It can thus be used to build a compatible family of strong hyper Markov hyperdistributions. Moreover, for any partitioning (A, B) of V ,

X A ∼ N (µ A , Λ A -Λ AB Λ -1 B Λ T AB -1 ) and (µ A , Λ A -Λ AB Λ -1 B Λ T AB ) is also normal-Wishart-distributed with parameters (ν A , λ, α -p + l, Ψ A -Ψ AB Ψ -1 B Ψ T AB )
where all indices are understood as partitioning of the corresponding vectors and matrices according to (A, B).

The pairwise marginal likelihoods can then be computed by updating the hyperparameters of the basis hyperdistribution to (ν , λ , α , Ψ ) thanks to classical Bayesian updating formulae. The locally updated hyperparameters are then derived from the globally updated ones and

p(D i , D j ) ∝ |Ψ {i,j} | α-p+2 2 |Ψ {i,j} | α -p+2 2 , p(D i ) ∝ |Ψ i | α-p+1 2 |Ψ i | α -p+1 2 , (2.17)
where, for a matrix M and i, j ∈ V , M {i,j} denotes the submatrix of size 2 corresponding to vertices i and j. This result is given in the work of [START_REF] Kuipers | Addendum on the scoring of gaussian directed acyclic graphical models[END_REF] as a correction to the erroneous result stated in [START_REF] Geiger | Parameter priors for directed acyclic graphical models and the characterization of several probability distributions[END_REF].

The compatible hyperdistributions built on (µ, Λ) are called hyper-normal-Wishart distributions. One can notice that Λ -1 follows a hyper-inverse-Wishart distribution (Dawid & Lauritzen, 1993, §7.3.2).

Integration with Respect to T

We assume that we have knowledge of ω. Consequently, we know ξ(•|D) up to the normalising constant Z. For an edge {k, l}, gaining access to P ({k, l} ∈ E T * |D, ξ) means being able to sum the posterior tree distribution over the trees that possess the edge {k, l}. Because we are only considering trees, this is tractable thanks to the Matrix-Tree theorem.

Let ω = (ω i,j ) (i,j)∈V 2 be a symmetric weight matrix such that, for all i ∈ V , ω i,i = 0, the off-diagonal terms being non-negative. The weight of a graph G = (V, E G ) is defined as the product of the weights of its edges,

ω G . . = {i,j}∈E G ω i,j .
The Laplacian ∆ = (∆ i,j ) (i,j)∈V 2 of ω is given by

∆ i,j =    -ω i,j if i = j, j ω i,j if i = j.
For U ⊆ V , we defined ∆ U as the matrix obtained from ∆ by removing the rows and columns corresponding to U .

Theorem 2.1 [START_REF] Chaiken | A Combinatorial Proof of the All Minors Matrix Tree Theorem[END_REF]. Let ∆ be the Laplacian of a weight matrix ω. Then all minors |∆ {u} |, u ∈ V , are equal and the following identity holds:

|∆ {u} | = T ∈T ω T .
(2.18)

We directly get the normalising constant of ξ(T |D) from this result.

There is a more general version of this theorem concerning graphs whose connected components are spanning trees on their respective sets of vertices. Such graphs are called forests.

Theorem 2.2 (All Minors Matrix-Tree theorem, [START_REF] Chaiken | A Combinatorial Proof of the All Minors Matrix Tree Theorem[END_REF]. Let ∆ be the Laplacian of a weight matrix ω and U ⊆ V . Let F U be the set of forests on V with |U | connected components such that, for any two vertices u 1 , u 2 ∈ U , u 1 and u 2 are not in the same connected component. Then

|∆ U | = F ∈F U ω F .
(2.19)

Briefly speaking, U can be seen as a set of "roots" (even though the models are not directed) for the trees of the forests in F U . If U is taken equal to a single vertex, then the forests in F U only have one connected component which is a tree and we get Theorem 2.1. This theorem will be used in the proof of the following result that was first stated by [START_REF] Kirshner | Learning with Tree-Averaged Densities and Distributions[END_REF].

Theorem 2.3 [START_REF] Kirshner | Learning with Tree-Averaged Densities and Distributions[END_REF]. Let ω be defined as in (2.10) and ∆ be the associated Laplacian. Let u be a vertex in V . We define Q = ∆ {u} -1 . Then, for all {k, l} ∈ P 2 (V ),

P ({k, l} ∈ E T * |D, ξ) =    ω k,l (Q k,k + Q l,l -2Q k,l ) if k = u, l = u, ω k,u Q k,k if l = u, ω u,l Q l,l if k = u.
(2.20)

A proof of this result is provided in the extended version of [START_REF] Kirshner | Learning with Tree-Averaged Densities and Distributions[END_REF] available online. We provide a shorter version relying on the generalized version of the Matrix-Tree theorem given above (see Eq. 2.21).

Proof. Let {k, l} be an edge in P 2 (V ). Let Z, Z + kl and Z - kl respectively denote the sums of ω T over the sets T , {T ∈ T : {k, l} ∈ E T } and {T ∈ T : {k, l} ∈ E T }. It is immediate to see that Z = Z + kl + Z - kl . Lemma 3 of [START_REF] Meilă | Tractable bayesian learning of tree belief networks[END_REF] states that

∂Z ∂ω k,l = M k,l |∆ {u} | = M k,l Z
where M is a symmetric matrix with 0 diagonal defined by

M i,j =    (Q i,i + Q j,j -2Q i,j ) if i = u, j = u, Q i,i if j = u, Q j,j if i = u.
It is easy to see that Z - kl can be obtained by applying Theorem 2.1 to a weight matrix equal to ω except for the terms ω k,l and ω lk that are set to 0. This means that Z - kl does not depend on ω k,l and

∂Z ∂ω k,l = ∂Z + kl ∂ω k,l .
By Theorem 2.2,

Z + kl = ω k,l F ∈F {k,l} ω F = ω k,l |∆ {k,l} |.
(2.21)

|∆ {k,l} | does not depend on ω k,l since the only terms of ∆ that depend on ω k,l are ∆ k,l , ∆ lk , ∆ k,k , ∆ l,l and these terms are all withdrawn in ∆ {k,l} . Therefore,

|∆ {k,l} | = ∂Z + kl ∂ω k,l = ∂Z ∂ω k,l = M k,l Z. (2.22)
Combining (2.21) and (2.22) with the fact that P ({k, l} ∈ E T |D, ξ) = Z + kl /Z, we get the claimed result.

Theorem 2.3 shows that all posterior probabilities for the edges can be computed at once by inverting a matrix of size p -1, amounting to a complexity of O(p 3 ).

In a Bayesian framework, the posterior entropy gives insight about the concentration of the posterior distribution, which is for instance of particular interest when a MAP approach is considered. The computation of this quantity is not always straightforward, but here, it can be obtained at small cost once posterior probabilities for the edges have been computed.

Proposition 2.6. The entropy of the posterior distribution on trees ξ(•|D) can be computed with complexity O(p 3 ).

Proof. We show that the entropy has a simple expression depending on Z and (P ({k, l} ∈ E T * |D, ξ)) {k,l}∈P2(V ) which can both be computed with complexity O(p 3 ) through Theorems 2.1 & 2.3. Indeed,

H(ξ(•|D)) = - T ∈T ξ(T |D) log (ξ(T |D)) = T ∈T 1 Z {i,j}∈E T ω i,j   log(Z) - {k,l}∈E T log(ω k,l )   = log(Z) - {k,l}∈P2(V ) log(ω k,l ) Z T {k,l} {i,j}∈E T ω i,j = log(Z) - {k,l}∈P2(V ) log(ω k,l )P ({k, l} ∈ E T * |D, ξ).

Controlling prior edge probability

If the distribution on trees is not strongly peaked, the prior probability for an edge to appear in a random tree can be quite small. For instance, the uniform distribution on T leads to any edge appearing with probability 2/p. Indeed, no edge is favoured and each tree borrows p -1 of the p(p -1)/2 possible edges. We consider an edge {k, l} ∈ P 2 (V ) and the event E kl . . = {T : {k, l} ∈ E T }. We let p 0 kl and p kl respectively denote the prior and posterior probabilities of event E kl . These probabilities are obtained through Theorem 2.3.

In a decision perspective, it might be desirable to allow some control on the prior probability of E kl . To this aim, we use a binary random variable kl ∼ B(λ kl ) explicitly controlling the status of edge {k, l} in the random tree:

p(T | kl , ξ) = ξ(T |E kl ) if kl = 1 ξ(T |E kl ) if kl = 0 .
In particular, the choice λ kl = 1/2 takes us back to a non-informative prior configuration regarding E kl . We obtain the model represented in Figure 2.3 in which the fully marginal likelihood can be written as

p(D) = λ kl p(D|E kl ) + (1 -λ kl )p(D|E kl ).
We are now interested in the posterior distribution of kl .

Proposition 2.7. Proof.

P ( kl = 1|D) = λ kl p kl p 0 kl • λ kl p kl p 0 kl + (1 -λ kl ) 1 -p kl 1 -p 0 kl -1 kl T * π X
P ( kl = 1|D) = p(D| kl = 1)P ( kl = 1) p(D) = λ kl p(D|E kl ) p(D) = λ kl p(D|E kl ) • λ kl p(D|E kl ) + (1 -λ kl )p(D|E kl ) -1 = λ kl p kl p 0 kl • λ kl p kl p 0 kl + (1 -λ kl ) 1 -p kl 1 -p 0 kl -1
The computation of P ( kl = 1|D) for all edges can be achieved in O(p 2 ) time from the posterior edge probability matrix {p kl } {k,l}∈P2(V ) . We can notice that P ( kl = 1|D) is a strictly increasing function of p kl . When the initial prior on trees ξ is uniform and all λ kl are taken equal, the order induced on the edges by {P ( kl = 1|D)} {k,l}∈P2(V ) is identical to the order induced by the posterior edge probability matrix. The ROC and PR curves that are commonly used to assess network inference accuracy therefore remain unchanged. In this section, we use synthetic data to meet a twofold objective. On one hand, the aim of this study is to show that there is an advantage in averaging over trees rather than considering a single MAP estimate. On the other hand, we show that assuming a tree structure is not substantially more detrimental to the accuracy of the inference of non-treestructured graphical models than assuming a DAG structure. To do so, we compare our method with another fully Bayesian inference method carried out on DAGs, described by [START_REF] Niinimäki | Partial order mcmc for structure discovery in bayesian networks[END_REF] and implemented in the BEANDisco software. Computations for our approach were performed with the R package saturnin.

Simulation Scheme

We have chosen four typical networks with p = 25 vertices, namely a tree and three Erdös-Rényi random graphs drawn with probabilities of connection p c = 2/p, 4/p and 8/p. These graphs are shown in Figure 2.4. Data sets are then simulated according to Gaussian graphical models. For all four adjacency matrices A, we used the Laplacian matrix of A augmented of 1 on the diagonal as precision matrix Λ A . This construction ensures that Λ A is non-singular. Independent samples were drawn according to N (0, Λ -1 A ) and discretized into r = 3 bins. For n = 25, 50, 75, 100 and 200, we generated 100 data sets of size n.

We then considered the Multinomial/Dirichlet framework described in Section 2.4.1.b, setting the prior on trees ξ to the uniform and the equivalent prior sample size N to r 2 /2 = 4.5 (see Eq. (2.15)). For each data set, we computed • the MAP tree structure in our model thanks to a Maximal Spanning tree algorithm applied to ω;

• the matrix of posterior edge probabilities P ({k, l} ∈ E T * |D) in our model. For all the edges, the prior appearance probability was brought back to q (kl) 0

= 1/2 (see Section 2.4.3); • an estimation of the matrix of posterior edge appearance probabilities in a random DAG obtained by MCMC sampling [START_REF] Niinimäki | Partial order mcmc for structure discovery in bayesian networks[END_REF]. We refer the reader to this paper for details on the prior distribution on DAGs. We ran the code provided by the authors with default parameters. The direction of the edges of the sampled DAGs was not taken into account to get empirical frequencies for all undirected edges.

We also considered the normal-Wishart/Gaussian framework described in Section 2.4.1.c on the undiscretized data.

The accuracy of the inference was evaluated against the true adjacency matrix according to the yielded outputs. In the case of the MAP estimate, we calculated the True and False Positives Rates (TPR, FPR) between the best tree and the true graph. These rates are constrained by the fact that a spanning trees on p vertices has exactly p -1 edges. For the (estimated) posterior edge appearance probability matrices, ROC and PR curves against the true adjacency matrix are plotted and summarized by the area under the curves.

Results

COMPARISON WITH MAP. Figure 2.5 simultaneously represents the (TPR, FPR) scores and the ROC curves obtained for the MAP estimate and the tree posterior edge appearance probability matrix respectively. It makes sense to plot both results on the same graph since a ROC curve is just a succession of (TPR, FPR) points computed as more and more edges are selected, going from the most to the least likely. When p -1 edges are selected, both methods behave similarly. So, if there is external evidence that the true graph is in fact a tree, a MAP approach could be considered but using posterior edge probabilities would do as well. Nonetheless, when the true graph is not a tree, the MAP approach is penalized by its lack of flexibility. Computing posterior appearance probabilities for the edges allows to retain an arbitrary number of edges. The balance between selectivity and sensibility achieved by the MAP approach can obviously be improved by selecting more edges. An other argument in favor of considering the whole posterior distribution on trees instead of the MAP is presented in Figure 2.6. For all four simulation scenarios, posterior tree distributions are not really peaked around their modes, especially for small samples, with the second most probable tree always being very close to the MAP. Moreover, the entropy of the posterior distribution on trees behaves similarly across all scenarios. We could have q q q q q q q q q q q q q q q q q q q q q q q q q Sample Size 25 Tree q q q q q q q q q q q q q q q q q q q q q q q q q q Sample Size 25 Erdös-Rényi, pc = 2/p q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Sample Size 25 Erdös-Rényi, pc = 4/p q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Sample Size 25 expected a singular behaviour in the tree scenario, but it is only slightly observed for n = 200.

Sample

Sample

Sample

Sample

INFLUENCE OF THE TREE ASSUMPTION.

We now study the influence of the tree assumption on the accuracy of structure inference when the true graphical model is not tree-structured. With this end in view, we consider a similar model where DAGs are drawn instead of trees and use the posterior edge appearance probabilities yielded by this model as gold standard, as it achieves the same goal in terms of Bayesian inference within a larger class of graphs.

Results are given in Figure 2.7. Both algorithms seem to perform equally well in all four situations. The accuracy of the inference expectedly increases with sample size. The results we get here indicate that the posterior probabilities for the edges to belong to a random tree can be relevant even when the true network is not a tree, with no clear evidence in favor of considering an inference within the broader class of DAG structures.

RUNNING TIME. We conclude this section on synthetic data by mentioning running times (Table 2.1). While retaining similar accuracy to the algorithm based on MCMC sampling in q q q q q q q q q q q q q q q -20 -10 25 100 200

Sample Size log of MAP Tree Probability q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1 2 3 25 100 200

Sample Size Ratio to Second Best Tree Probability q q q q q q q q q q q q q q q 0.2 (p -2) log(p)).

q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.4 0.5 0.6 0.7 0.8 0.9

1.0 25 50 75 100 200 ROC q q q q q q q q q q q q q 0.0 Sample Size PR Tree q q q q q q q q q q q q q q q q q 0.4 0.5 0.6 0.7 0.8 0.9

1.0 25 50 75 100 200 ROC q q q q q q q q q q q q q q q q 0.0 Sample Size

PR

Erdös-Rényi, pc = 2/p q q q q q q q q q q q 0.4 0.5 0.6 0.7 0.8 0.9

1.0 25 50 75 100 200 ROC q q q q q q q q q q q q q q 0.0 Sample Size

PR

Erdös-Rényi, pc = 4/p q q q q q q q q q q q 0.4 0.5 0.6 0.7 0.8 0.9

1.0 25 50 75 100 200 ROC q q q q q q q q 0.0 the space of DAGs that we used as a point of comparison, our algorithm runs significantly faster than the MCMC sampling ran with default parameters, especially for large networks.

Of course, the accuracy of the MCMC sampling approach could be improved by augmenting the number of samples at the cost of even longer running times, but we have not observed any evidence going that way. network inference techniques. They are related to the Raf cellular signalling network, which is involved in many different processes, including the regulation of cellular proliferation in human immune cells. The activation levels of the 11 proteins and phospholipids that are part of this pathway can be measured by flow cytometry. The generally accepted structure of the Raf pathway is given in Figure 2.8, but the true structure of this network, despite considerable experimental and theoretical efforts, may be more subtle. The undirected skeleton of this network will, however, be used as the gold standard network in our study.

Application to

Data

In flow cytometry experiments, cells are suspended in a stream of fluid and go through a laser beam one at a time. Different parameters are then measured on each cell by recovering the light that is reemitted by diffusion or fluorescence. We are interested in the activation levels (also called phosphorylation levels) of the involved proteins and phospholipids. Such experiments typically produce samples of several thousands observations. Since all biological network inference problems are not met by such a profusion of data, [START_REF] Werhli | Comparative evaluation of reverse engineering gene regulatory networks with relevance networks, graphical gaussian models and bayesian networks[END_REF] sampled down 5 samples with 100 data points from the data provided by [START_REF] Sachs | Causal protein-signaling networks derived from multiparameter single-cell data[END_REF].

We discretized each sample into r=3 bins and performed the inference on each of them with our algorithm (Tree) and the MCMC sampling in DAGs algorithm (DAG), just like in the previous section. The accuracy of the inference was once again assessed by the area under the ROC and PR curves, averaged on all 5 samples.

Results

The results of the inference are reported in Table 2.2. The DAG approach performs globally better than our inference within trees. These results qualify those of the previous section. Nonetheless, we would like to make the following points. While not being as accurate, our approach still provides good results and might in fact be more adapted to bigger problems where MCMC sampling can hardly be contemplated. Moreover, unlike the simulation study, the gold standard network against which the accuracy of the inference is assessed here, shown in Figure 2.8, is not perfectly known and may still differ quite considerably from the truth. The difference observed between the two approaches is small enough that it could be reversed if the status of an edge was to be changed in the gold standard network.

Figure 2.9 gives a graphical representation of the results obtained on one of the five data sets, offering a more detailed overview. We note that the gold standard network as defined here has 20 edges. The two likeliest trees in the posterior tree distribution are given in Figure 2.9a. Both trees have 9 true positives out of the p -1 = 10 edges they respectively selected. As expected, most of these edges also have strong posterior probabilities (Figure 2.9b). When the prior probabilities of all edges is brought back to 1/2, we get 13 edges with posterior probabilities strictly greater than 1/2, among which the same true positives as in the MAP estimate. More generally, one could consider using the histogram of posterior probabilities to empirically find a more appropriate cut-off.

We did not represent the empirical edge frequencies obtained for DAGs since prior appearance probabilities could not be easily accounted for in this case, thus making direct comparison with posterior edge probabilities in trees impossible.

As a conclusion, these results lead us to believe that it might be preferable to favour inference using DAGs for small problems. When that is no longer possible in a reasonable DAG Tree ROC 0.767 (0.068) 0.729 (0.047) PR 0.725 (0.070) 0.690 (0.051) We consider the problem of change-point detection in multivariate time-series. The multivariate distribution of the observations is supposed to follow a graphical model, whose graph and parameters are affected by abrupt changes throughout time. Building on the approach develop in Chapter 2, we assume that all graphical models are tree-structured, and make the best use of this assumption to compute a variety of quantities of interest in polynomial time with respect to the number of vertices and the length of the series. This chapter has been published as an article in the journal Statistics & Computing (Schwaller & Robin, 2016). We are interested in time-series data where several variables are observed throughout time.

An assumption often made in multivariate settings is that there exists an underlying network describing the dependences between the different variables. When modelling time-series data, one is faced with a choice: shall this network be considered stationary or not? Taking the example of genomic data, it might for instance be unrealistic to consider that the network describing how a pool of genes regulate each other remains identical throughout time. This network might slowly evolve, or undergo abrupt changes leading to the initialisation of new morphological development stages in the organism of interest. Here, we focus our interest on the second scenario.

The inference of the dependence structure ruling a multivariate time-series was first performed under the assumption that this structure was stationary (e.g. [START_REF] Friedman | Learning the structure of dynamic probabilistic networks[END_REF][START_REF] Murphy | Modelling gene expression data using dynamic bayesian networks[END_REF]). Non-stationarity has then been addressed in a variety of ways. Classical Dynamic Bayesian Networks (DBNs) can for instance be adapted to allow the directed graph (or Bayesian Network) describing the interactions between two consecutive time-points to change, leading to so-called switching DBNs [START_REF] Robinson | Learning non-stationary dynamic Bayesian networks[END_REF][START_REF] Lèbre | Statistical inference of the time-varying structure of gene-regulation networks[END_REF][START_REF] Grzegorczyk | Improvements in the reconstruction of time-varying gene regulatory networks: dynamic programming and regularization by information sharing among genes[END_REF]. Some models alternatively suppose that the heterogeneity is the result of parameters changing smoothly with time [START_REF] Zhou | Time varying undirected graphs[END_REF][START_REF] Kolar | Estimating time-varying networks[END_REF]. This is especially appropriate for Gaussian graphical models where the graph structure can directly be read in the non-zero terms of the precision (or inverse-covariance) matrix, therefore enabling smooth transitions within the otherwise discrete space of graphs. Hidden Markov Models (HMM) have also been used to account for heterogeneity in multivariate time-series [START_REF] Fox | Nonparametric Bayesian learning of switching linear dynamical systems[END_REF][START_REF] Barber | Graphical Models for Time-Series[END_REF]. In the aforementioned models, the inference can rarely be performed exactly, and often relies on sampling techniques such as Markov Chain Monte Carlo (MCMC).

The model that we consider here belongs to the class of product partition models (PPM) [START_REF] Barry | Product Partition Models for Change Point Problems[END_REF]. We assume that the observed data {y t } t=1,...,N are a realisation of a process {Y t } t=1,...,N where, for 1 t T , Y t is a random vector with dimension p 2. If m is a segmentation of {1, . . . , T } with change-points

1 = τ 0 < τ 1 < • • • < τ K-1 < τ K = N , the model has the general form Y t ∼ π(G r , θ r ), if t ∈ r and r = τ i , τ i+1 ,
where G r and θ r respectively stand for the graph describing the dependence structure and the distribution parameters on segment r. The parameters (G r , θ r ) are assumed to be independent between segments. This model is illustrated in Figure 3.1. We are interested in retrieving all change-points at the same time, therefore performing off-line detection. It has been shown that both off-line [START_REF] Fearnhead | Exact and efficient Bayesian inference for multiple changepoint problems[END_REF][START_REF] Rigaill | Exact posterior distributions and model selection criteria for multiple change-point detection problems[END_REF] and on-line detection [START_REF] Fearnhead | On-line inference for multiple changepoint problems[END_REF][START_REF] Caron | On-line changepoint detection and parameter estimation with application to genomic data[END_REF] of change-points can be performed exactly and efficiently in this model thanks to dynamic programming. [START_REF] Xuan | Modeling changing dependency structure in multivariate time series[END_REF] explicitly consider this framework in a multivariate Gaussian setting. They estimate a set of possible structures for their model by performing regularized estimation of the precision matrix on arbitrary overlapping time segments. This graph family is then taken as a starting point in an iterative procedure where the segmentation and the graph family are sequentially updated to get the best segmentation and graph series.

OUR CONTRIBUTION From a Bayesian point of view, the problem at hand raises an interesting and quite typical problem as both continuous and discrete parameter are involved in the model. Indeed, the location and scale parameters or, more specifically, the means and (conditional) covariances associated with each segments are continuous but the location of the change-points and the structure of the graphical model within each segments are not. Denoting θ the set of continuous parameters, Q the set of discrete parameters and y the observed data, Bayesian inference will typically rely on integrals such as the marginal likelihood of the data, that is

p(y) = Q∈Q p(Q) θ∈Θ p(y|θ, Q)p(θ|Q)dθ.
In many situations, the use of conjugate priors allows us to compute the integral with respect to θ in an exact manner. Still, the summation over all possible values for the discrete parameter Q is often intractable due to combinatorial complexity. One aim of this article is to remind that the algebraic properties of the space Q can sometimes help to actually achieve this summation in an exact manner, so that a fully exact Bayesian inference can be carried out.

We show that exact and efficient Bayesian inference can be performed in a multivariate product partition model within the class of undirected graphs called spanning trees. These structures are connected graphs, with no cycles. When p nodes are considered, we are left with p p-2 possible spanning trees, but exact inference remains tractable by using algebraic properties pertaining to this set of graphs. On each independent temporal segment, we place ourselves in the framework developed by Schwaller et al. (2015), in which the likelihood of a segment s; t , defined by p(y s;t ) . . = T ∈T p(y s;t |θ, T )p(θ|T )dθ, where T stands for the set of spanning trees, can be computed efficiently. We provide explicit and exact formulas for quantities of interest such as the posterior distribution of changepoints or posterior edge probabilities over time. We also provide a way to assess whether the status of an edge (or of the whole graph) remains identical throughout the time-series or not when the partition is given.

OUTLINE In Section 3.2, we provide some background on graphical models and product partition models. In particular, we give a more detailed presentation of the results of [START_REF] Rigaill | Exact posterior distributions and model selection criteria for multiple change-point detection problems[END_REF] on dynamic programming used for change-point detection problems. We also introduce tree-structured graphical models. The model and its properties are presented in Section 3.3. Section 3.4 enumerates a list of quantities of interest that can be computed in this model, while Section 3.5 deals with edge and graph status comparison, when the segmentation is known. Sections 3.6 and 3.7 respectively present the simulation study and the applications to both biological and neuroscience data. In this section we introduce two models involving a discrete parameter, for which exact integration over this parameter is possible.

Product Partition Models

Let Y = {Y t } t=1,...,N be an independent random process and let y be a realisation of this process. For any time interval r, we let Y r . . = {Y t } t∈r denote the observations for t ∈ r. PPMs as described in [START_REF] Barry | Product Partition Models for Change Point Problems[END_REF] work under the assumption that the observations can be divided in independent adjacent segments. Thus, if m is a partition of 1; N , the likelihood of y conditioned on m can be written as

p(y|m) = r∈m p(y r |r), p(y r |r) = t∈r p(y t |θ r ) p(θ r )dθ r ,
where θ r is a set of parameters giving the distribution of Y t for t ∈ r. For the sake of clarity, we let p(y r ) denote p(y r |r) in the following.

For K 1, we let M K denote the set made of the partitions of 1; N into K segments. The cardinality of this set is N -1 K-1 . More generally, we let M K ( s; t ) denote the partitions of any interval s; t into K segments. In order to get the marginal likelihood of y conditionally on K, one has to integrate out both m and θ = {θ r } r∈m :

p(y|K) = m∈M K p(m) r∈m p(y r ) = m∈M K p(m) r∈m t∈r p(y t |θ r ) p(θ r )dθ r .
If the distribution of m, conditional on K, factorises over the segments with an expression of the form

p(m|K) = 1 C K (a) r∈m a r , (3.1)
where a r are non-negative weights assigned to all segments and C K (a) = m∈M K r∈m a r is a normalising constant, these integrations can be performed separately. [START_REF] Rigaill | Exact posterior distributions and model selection criteria for multiple change-point detection problems[END_REF] introduced a matrix containing the weighted likelihood of all possible segments, whose general term is given by

A s,t = a s;t • p(y s;t ) if 1 s < t N + 1, 0 otherwise. (3.2)
This matrix can be used in an algorithm designed according to a dynamic programming principle to perform the integration on M K efficiently.

Proposition 3.1 [START_REF] Rigaill | Exact posterior distributions and model selection criteria for multiple change-point detection problems[END_REF].

[A K ] s,t = m∈M K ( s;t ) r∈m a r • p(y r )
where A k denotes the k-th power of matrix A and A k s,t its general term. Moreover,

A K . . = {[A k ] 1,t , [A k ] t,n+1 } 1 k K 2 t N can be computed in O(KN 2 ) time.
In particular,

[A K ] 1,n+1 = C K (a) • p(y|K).
Several quantities of interest share the same form: from A K , Rigaill et al. ( 2012) also derived exact formulas for the posterior probability of a change-point to occur at time t or for the posterior probability that a given segment r belongs to m (see Section 3.4.1). Classical Bayesian selection criteria for K are also given. One can notice that C K (a) can be recovered by applying Proposition 3.1 not to matrix A but to a matrix defined similarly from a. For the uniform distribution on M K , i.e. a r ≡ 1, [START_REF] Fearnhead | Exact and efficient Bayesian inference for multiple changepoint problems[END_REF] worked under a slightly different model where m is not chosen conditionally on K but is instead drawn sequentially by specifying the probability mass function for the time between two successive change-points. They presented a filtering recursion to compute the marginal likelihood of the observations under their model where the integrations over parameters and segmentations are also uncoupled. [START_REF] Fearnhead | On-line inference for multiple changepoint problems[END_REF] showed that on-line and exact inference is also tractable in this model.

we get C K (a) = N -1 K-1 .

Tree-structured Graphical Models

In a multivariate setting, graphical models are used to describe complex dependence structures between the involved variables. A graphical model is given by a graph, either directed or not, and a family of distributions satisfying some Markov property with respect to this graph. We concentrate our attention on undirected graphical models, also called Markov random fields. We refer the reader to [START_REF] Lauritzen | Graphical Models[END_REF] for a complete overview on the subject. Let V = {1, ..., p} and Y = (Y 1 , ..., Y p ) be a random vector taking values in a product space X = p i=1 X i . We consider the set T of connected undirected graphs with no cycles. These graphs are called spanning trees. For T ∈ T , we let E T denote the edges of T .

We consider a hierarchical model where one successively draws a tree T in T , the parameters θ of a distribution that factorises according to T , and finally a random vector Y according to this distribution. The marginal likelihood of the observations under this model, where both θ and T are integrated out, is given by p(y) = T ∈T p(T ) p(y|T, θ)p(θ|T )dθ.

These integrations can be performed exactly and efficiently by choosing the right priors on T and θ [START_REF] Meilă | Tractable bayesian learning of tree belief networks[END_REF]Schwaller et al., 2015). The distribution on trees is taken to be factorised on the edges,

p(T ) = 1 Z(b) {i,j}∈E T b i,j , (3.3)
where b i,j are non-negative edge weights and

Z(b) . . = T ∈T {i,j}∈E T b i,j (3.4)
is a normalizing constant. The prior on θ has to be specified for all trees in T . The idea is to require each of these priors to factorise on the edges and to specify a prior on θ ij once and for all trees, θ ij designating the parameters governing the marginal distribution of (Y i , Y j ). These priors must be chosen coherently, in the sense that, for all i, j, k ∈ V , the priors on θ ik and θ jk should induce the same prior on θ k . Some local Markov property is also needed. Schwaller et al. (2015) especially detailed three frameworks in which this can be achieved, namely multinomial distributions with Dirichlet priors, Gaussian distributions with normal-Wishart priors and copulas. We elaborate a little more on the particular case of Gaussian graphical models (GGMs). In a multivariate Gaussian setting, θ = (µ, Λ) where µ and Λ respectively stand for the mean vector and precision matrix of the distribution. A classical result on GGMs states that if the (i, j)-th term of the precision matrix is equal to zero, there is no edge between nodes i and j. Thus, the support of p(θ|T ) is the set of sparse positive definite matrices whose non-zero terms are given by the adjacency matrix of T . The distribution of θ|T can be defined for all trees at once by using a general normal-Wishart distribution defined on all positive-definite matrices (Schwaller et al., 2015, Sec. 4.1.3). Marginal distributions of this normal-Wishart distributions are used to build distributions for {θ|T } T ∈T . When p(θ|T ) is carefully chosen, the integration on θ can be performed independently from the integration on T and p(y|T ) factorises on the edges of T :

p(y|T ) = i∈V p(y i ) {i,j}∈E T p(y i , y j ) p(y i )p(y j )
where p(y i , y j ) = p(y i , y j |θ ij )dθ ij , (3.5)

p(y i ) = p(y i |θ i )dθ i .
Computing {p(y|T )} T ∈T only requires p(p+1)/2 computations of low-dimensional integrals, where p is the dimension of the model. As both p(T ) and p(y|T ) factorise on the edges, integrating the likelihood over T can be performed in O(p 3 ) time.

Proposition 3.2. The marginal likelihood is given by

p(y) = Z(ω) Z(b) • i∈V p(y i )
where Z(•) is defined as in (3.4) and ω is the posterior edge weight matrix whose general term is given by ω i,j . . = b i,j p(y i , y j ) p(y i )p(y j ) .

(3.6)

Moreover, p(y) is obtained in O(p 3 ) time from b and ω.

Proof.

p(y) = T ∈T p(y|T )p(T ) = 1 Z(b) i∈V p(y i ) T ∈T {i,j}∈E T b i,j p(y i , y j ) p(y i )p(y j ) = Z(ω) Z(b) • i∈V p(y i ),
with ω as defined above. As Z(•) can be computed in O(p 3 ) time using the Matrix-Tree theorem, we get the announced complexity.

The posterior probability for an edge to belong to T , P ({i, j} ∈ E T |y), can also be obtained for all edges at once in O(p 3 ) time (Schwaller et al., 2015, Th. 3). Sections 3.2.1 and 3.2.2 presented two models in which Bayesian inference requires us to integrate out a fundamentally discrete parameter (either the segmentation m or the spanning tree T ) and other (usually continuous) parameters θ. In both situations, these integrations can be performed exactly and efficiently by uncoupling the problems. The integration over θ is performed "locally" and the results are stored to be used in an algorithm that heavily relies on algebra to integrate over the discrete parameter. This is made possible by a careful choice of priors for both parameters. Our aim is to show that these algebraic tricks can be combined to perform exact Bayesian inference of multiple change-points in the dependence structure of multivariate time-series.

Model

Model

It is assumed that the observed data y = {y t } N t=1 are a realisation of a multivariate random process Y = {Y t } N t=1 of dimension p 2. For 1 t N , Y t = (Y t 1 , ..., Y t p ) is a multivariate random variable of dimension p taking values in a product space X = p i=1 X i . We model Y by a PPM where, at each time-point, observations Y t are modelled by a tree-structured graphical model. If m is a segmentation with K segments, we let T = {T k } K k=1 and θ = {θ r } r∈m respectively denote the tree structures and parameters for each segment. For r ∈ m, we also let κ(r|m) denote the position of r in m. We use index r rather than k as, in the following, all possible segments r will have to be considered per se, whatever the segmentation they belong to. Our model can then be written as follows:

p(m|K) = 1 C K (a) r∈m a r , p(T|K) = K k=1 p(T k ) = 1 Z(b) K K k=1 {i,j}∈E T k b i,j , p(θ|m, T) = r∈m p(θ r |T κ(r|m) ), p(y|m, θ, T) = r∈m t∈r p(y t |T κ(r|m) , θ r ).
For r ∈ m, {Y t } t∈r are independent and identically distributed with structure T κ(r|m) and parameters θ r . The priors on m and each of T k are respectively taken of the form given in (3.1) and (3.3) through segment weights a and edge weights b. The distribution of θ r |{T κ(r|m) = T } is assumed to factorise over the edges of T , coherently between all spanning trees T ∈ T , as described in Section 3.2.2. A graphical representation of this model is given in Figure 3.2.

Factorisation Properties

In the model that we have described, the marginal likelihood of the observation, conditionally on K, is given by p(y|K) = m∈M K T∈T K p(y, m, θ, T|K)dθ.

(3.7)

Integrating out the discrete parameters (m, T) requires to sum over a set of cardinality

|M K | • |T K | = N -1 K -1 • p K(p-2) ≈ N p p-2 K K .
Nonetheless, the joint distribution of (m, θ, T), conditionally on K, factorises at different levels and integration can therefore be performed by combining the results given in Section 3.2. Proposition 3.3. The marginal likelihood p(y|K) can be computed in O(max(K, p 3 )N 2 ) time, where p and N respectively stand for the dimension of the model and the length of the series, from the posterior edge weight matrices computed on all possible segments r, whose general terms are given by ω (r) i,j

. . = b i,j p(y r i , y r j ) p(y r i )p(y r j )

.

(3.8) p(y r i , y r j ) and p(y r i ) are local integrals on θ computed on edges and vertices as defined in (3.5).

Proof. For any segmentation m ∈ M K of 1; N into K segments, {(T κ(r|m) , θ r )} r∈m are independent, so that p(y, m|K) can be written as

p(y, m|K) = 1 C K (a) r∈m a r p(y r ),
where p(y r ) stands for the locally integrated likelihood of y r on segment r, p(y r ) = Computing matrix A requires to integrate the likelihood over tree structure T and parameters θ for all possible segments r ⊆ 1; N . On each segment, we fall back to the tree-structured model described in Section 3.2.2 and the integrated likelihood can be expressed using the local terms computed on vertices and edges that were defined in (3.5). Indeed, for r ⊂ 1; N , p(y r ) is obtained through Proposition 3.2 applied to ω (r) (defined in (3.8)):

p(y r ) = Z(ω (r) ) Z(b) • i∈V p(y r i ).
where we remind that Z(•) is the function giving the normalising constant of a tree distribution. As a consequence, A is computed in O(p 3 N 2 ) time from the posterior edge weight matrices {ω (r) } r⊆ 1;N , hence the total complexity.

The components of the matrices ω (r) result from the integration over θ, which can be made separately and locally thanks to the assumptions made on its prior distribution in Section 3.3.1. This integration comes down to remove node θ in the global graphical model displayed in Figure 3.2. Marginal likelihood is only one of many quantities than might be of concern in this model. Yet, once matrix A has been calculated, other quantities of interest with respect to our model can be obtained at low cost. The next section provides a non-exhaustive list of such quantities.

Proposition 3.4.

p(K|y) ∝ p(K)[A K ] 1,N +1 C K (a) .
Proof. Bayes' rule states that p(K|y) ∝ p(K)p(y|K) and by Proposition 3.1, p(y|K) =

[A K ] 1,N +1 /C K (a).
The best segmentation a posteriori can also be recovered efficiently by using matrix A in the Segment Neighbourhood Search algorithm [START_REF] Auger | Algorithms for the optimal identification of segment neighborhoods[END_REF]. Thus, if one's interest lies in retrieving the number of segments K, two estimators can be considered

K1 = arg max K p(K|y), K2 = K(arg max m p(m|y)).
where K(m) stands for the number of segments in m.

Posterior Edge Probability

For any segment r ⊆ 1; N , it is possible to compute the posterior edge probabilities corresponding to segment r:

P ({i, j} ∈ E T |y r ), ∀{i, j} ∈ P 2 (V ),
where T is a random tree distributed as T 1 , . . . , T K . Whenever m is unknown, the segmentation can be integrated out to obtain instant posterior edge probabilities at any given time t. Conditionally on K, the instant posterior appearance probability of edge {i, j} at time t can be written as

p K ij (t) . . = m∈M K p(m|y, K)P ({i, j} ∈ E T κ(t|m) |y, m),
where κ(t|m) gives the position of the segment containing t in m.

Proposition 3.5. The instant posterior probability of edge {i, j} at time t is given by

p K ij (t) = r t S K (r)P ({i, j} ∈ E T |y r ). (3.10) {p K ij (t)} 1 i,j p 1 t N can be computed in O max(K, p 3 )N 2 time from A and {ω (r) } r⊆ 1;N .
Proof. This formula is similar to the one giving the posterior mean of the signal in [START_REF] Rigaill | Exact posterior distributions and model selection criteria for multiple change-point detection problems[END_REF]. If r ∈ m and t ∈ r, then P ({i, j} ∈ E T κ(t|m) |y, m) = P ({i, j} ∈ E T |y r ), hence the result. {S K (r)} r∈ 1;N is obtained with complexity O(KN 2 ) and {P ({i, j} ∈ E T |y r )} r∈ 1;N with complexity O(p 3 N 2 ), and that gives an upper bound on total complexity.

One could be interested in computing the posterior probability for an edge to keep the same status throughout time when m is integrated out, given K. Nonetheless, it would require to integrate on subsets of M K ⊗ T K that are in direct contradiction with the factorability assumption, making the results that we have presented so far useless. Indeed, we would effectively be introducing dependency between segments, thus breaking up the factorability of p(y, m) with respect to r ∈ m. In this situation, Proposition 3.1 can no longer be used. A drastic workaround is to work under a fixed segmentation instead of integrating out m, and this is what we do in the following section. We now turn to the specific case where m is known and has K segments (r 1 , . . . , r K ). This situation is far less general than the framework we considered until now. Still, it corresponds to some practical situations where segment comparison is interesting and for which further exact inference can be carried out.

Edge

Edge Status Comparison

Let i, j be two distinct nodes in V . We are interested in computing the posterior probability of the subsets of T K defined by

E + ij = {T = (T 1 , . . . , T K )|∀k ∈ 1, K , {i, j} ∈ E T k }, E - ij = {T = (T 1 , . . . , T K )|∀k ∈ 1, K , {i, j} / ∈ E T k }, E ij = E + ij ∪ E - ij ,
that respectively correspond to the situations where edge {i, j} is always present, always absent, or has the same status in all trees. If T belongs to E ij = T K \ E ij , it means that there exists two segments in which {i, j} does not have the same status. We let (q - 0 , q 0 , q + 0 ) respectively denote the prior probabilities of E - ij , E ij and E + ij . These probabilities can be written as

q - 0 = K k=1 P ({i, j} / ∈ E T k ) = P ({i, j} / ∈ E T ) K , q + 0 = P ({i, j} ∈ E T ) K , q 0 = 1 -q - 0 -q + 0 ,
where T is a tree distributed as T 1 , . . . , T K , and are obtained for all edges at once in O(p 3 ) time by computing the prior edge probability matrix (P ({i, j} / ∈ E T )) 1 i j p . Posterior probabilities (q -, q, q + ) for E - ij , E ij and E + ij can be computed similarly but one posterior edge probability matrix has to be calculated per segment:

q -= K k=1 P ({i, j} / ∈ E T k |y r k ), q + = K k=1 P ({i, j} ∈ E T k |y r k ), q = 1 -q --q + ,
However, if the prior distribution on trees is not strongly peaked, as events E + ij and E - ij only account for a relatively small number of tree series in T K , q - 0 and q + 0 (as well as q - and q + ) will always be very small. To allow some control on the prior probabilities of these events, we use a random variable ij taking values {-1; 0; 1} with probabilities (λ -, λ, λ + ) and explicitly controlling the status of edge {i, j} in all trees:

p(T| ij ) =    p(T|E + ij ) if ij = 1, p(T|E ij ) if ij = 0, p(T|E - ij ) if ij = -1
. We obtain the model described in Figure 3.3, in which Proposition 3.6. The vector of posterior probabilities for ij is proportional to λ -q - q - 0 , λ q q0 , λ + q + q + 0 .

p(y) = λ + p(y|E + ij ) + λ -p(y|E - ij ) + λp(y|E ij ).
Proof. We have that

p( ij = 1|y) = λ + p(y|E + ij ) p(y) = λ + p(y|E + ij ) λ + p(y|E + ij ) + λ -p(y|E - ij ) + λp(y|E ij ) = λ + q + q + 0 λ + q + q + 0 + λ -q - q - 0 + λ q q0 .
We reason similarly with p( ij = -1|y) to get the result.

Structure Comparison

The same reasoning can be applied for the global event

E = {T = (T 1 , . . . , T K )|∃T ∈ T , ∀k ∈ 1, K , T k = T },
which corresponds to a constant dependency structure across all segments (we remind that, in this section, the segments are known a priori), with possible changes for the parameters. The prior probability of E is given by

q 0 . . = P (E) = 1 Z(b) K T ∈T {i,j}∈E T b K i,j = Z(b K ) Z(b) K ,
where b K stands for the element-wise K-th power of matrix b. On each segment r k , the posterior distribution on trees factorises as

p(T k |y r k ) = 1 Z(ω (k) ) {i,j}∈T k ω (k) i,j ,
and the posterior probability of E is therefore given by

q . . = T ∈T K k=1 p(T |y r k ) = Z k ω (k) k Z(ω (k) )
where denotes the element-wise matrix product. Just as in the edge status comparison, we let a binary variable ∼ B(π) control the prior probability of E, with p(T| = 1) = p(T|E), and derive a similar formula for the posterior distribution of . Proposition 3.7. |y ∼ B(π * ) with π * . . = π q q 0 π q q 0 +(1-π) 1-q 1-q 0 . Proof. Similar to Proposition 3.6. Our approach was especially concerned with explicitly modelling the structure of the graphical model within each segment, but a simpler model could be considered in which the structure remains implicit. In a Gaussian setting, that would mean that the precision matrix governing the distribution on a given segment would be drawn without any zero-term constraints. One goal of this simulation study is to show how both models (with and without structure constraints) comparatively behave when one is interested in retrieving the number of segments or the location of the change-points.

Simulations

Another concern addressed by these simulations is the cost of the tree assumption when the true model is not tree-structured. How well can the number of segments, the changepoints or even the structures be recovered when the true networks are not trees?

Simulation Scheme

For this study, we generated time-series of size N = 70, 140 and 210. We choose segmentations with four segments of lengths 3 7 N , 1 7 N , 2 7 N and 1 7 N such that the relative length of each segment is kept identical through all sample sizes. The number of variables was fixed to p = 10. To give an idea of the sizes of the discrete sets we are working with, for N = 210, the cardinalities of the segmentation and tree sets are respectively |M 4 | ≈ 1.5 • 10 6 and |T | = 10 8 , so the size of the space to be explored is ≈ 1.5 • 10 38 . We built three structure scenarios by sampling structures from the uniform distribution on spanning trees, or from an Erdös-Rényi random graph distribution with connection probability p C = 2/p or 4/p. Thus, for each scenario, we got a series {∆ r } r∈M N of adjacency matrices describing the structure of the graphical model on all segments. The observations on a segment r were then drawn according to a multivariate Gaussian distribution with mean vector zero and precision matrix Λ r equal to the Laplacian matrix of ∆ r augmented of 1 on the diagonal, rescaled so that each variable as unit variance. For each sample size and structure series, 100 datasets were generated.

As described in the introduction of this section, the inference was then performed in the two following models. The first one is the full precision matrix model, without any structure constraint, and is given by

{Λ r } r∈m i.i.d., Λ r ∼ W(α, Ψ), (3.11) {Y t } N t=1 independent, Y t ∼ N (0 p , Λ -1 r ), ∀t ∈ r.
where W(α, Ψ) stands for the Wishart distribution with α degrees of freedom and scale matrix Ψ. The second one is the corresponding model with tree-structure assumption, as described in Section 3.3.1, and given by

{T k } K k=1 i.i.d., T k ∼ U (T ), {Λ r } r∈m independent, Λ r ∼ hW(α, Ψ, T κ(r|m) ), (3.12) {Y t } N t=1 independent, Y t ∼ N (0 p , Λ -1 r ), ∀t ∈ r,
where we let hW(α, Ψ, T ) denote the hyper-Wishart distribution based on W(α, Ψ) and with structure T (Schwaller et al., 2015). In both cases, we set α = p + 10 and Ψ = (α -p -1) • I p , where I p stands for the identity matrix of size p. The distribution of m|K is set to the uniform on M K and K follows a Poisson distribution with parameter γ = 4, truncated to 1; 10 . We emphasize the fact that, when the tree-structured model is considered, the series of precision matrices {Λ r } r∈m used to generate the data only belongs to the support of the law in the first structure scenario. The graphs drawn from the Erdös-Rényi distributions are not trees and therefore cannot induce precision matrices in the support of a tree-structured hyper-Wishart distribution. On the contrary, the full model obviously allows such precision matrices.

Finally, for the sake of clarity, we limited our study to centered data and null mean models, but one could allow the mean to vary between segments by using a (hyper) normal-Wishart distribution for (µ r , Λ r ), where µ r stands for the mean on segment r.

Results

We plotted the posterior probability of a change-point intervening at time t, integrated over K, as a function of t in the tree-structured and full models (Figure 3.4). In both cases, Full Tree Full Tree Full Tree q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.5 0.6 0.7 0.8 0.9

••• • •• • • ••• •••
••• • • • •••• •• • •••
• • • •• • •• •• • •• •••••
1.0 0 20 40 60 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.5 0.6 0.7 0.8 0.9

1.0 0 50 100 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.5 0.6 0.7 0.8 0.9

1.0 0 50 100 150 200 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.5 0.6 0.7 0.8 0.9

1.0 0 20 40 60 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.5 0.6 0.7 0.8 0.9

1.0 0 50 100 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.5 0.6 0.7 0.8 0.9

1.0 0 50 100 150 200 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.5 0.6 0.7 0.8 0.9

1.0 0 20 40 60 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.5 0.6 0.7 0.8 0.9

1.0 0 50 100 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.5 with respect to the true adjacency matrix at time t. We set K to the true number of segments (K = 4). The curve represents the mean value obtained from the 100 samples and the ribbon gives the standard deviation. change-points are hardly retrieved in the high-density Erdös-Rényi scenario, the inference performing better in the other two low-density scenarios. The standard deviations across samples are lower for the tree-structured model than for the full model. We can also observe a smoother behaviour with respect to time in the tree-structured model. Results on simulations with a greater number of segments (not shown) confirmed these observations. As expected, the shortest segments are hardly detected when the length of the series is small. These results seem to show that, when one is interested in retrieved change-point locations, the tree-structured model that we have presented can be considered in non-tree scenarios without any meaningful drop in performances.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 0.
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 0.
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 0.
NUMBER OF SEGMENTS For each sample, we computed K = arg max K p(K|y) and K( m) = K(arg max m p(m|y)). The results are given in Figure 3.5. In the full model, the number of segments selected by K and K( m) varies a lot across samples and is usually higher than in the tree model. In the tree-structured model, both K and K( m) tend to slightly underestimate the number of segments, especially in the highly-connected Erdös-Rényi scenario. They also display a more stable behaviour in the tree model. On small samples, K( m) seems to achieve better stability.

POSTERIOR EDGE PROBABILITY For t ∈ 1; N , we computed the posterior edge probability matrix defined in (3.10) for K = 4. Figure 3.6 shows the area under the ROC curve of this matrix against the true adjacency matrix at time t. In all scenarios, the structure is better retrieved on long segments. A drop in the accuracy is systematically observed near true change-points. While presenting lower accuracy compared to the other two scenarios, the structure inference in the highly connected scenario still provides meaningful results.

EDGE STATUS COMPARISON

The posterior probability for an edge to keep the same status throughout time was computed for all edges as explained in Section 3.5. We set the prior probability to change status at λ = 0.5 and the prior probabilities to be always present or absent to λ + = λ -= 0.25. We expected edges changing status during the time-series to be given low posterior probabilities. For small samples and across all scenarios, the posterior probability to have the same status remains close to the prior probability 0.5 for all edges. When samples grow bigger, a small contrast sets up according to the edges effectively changing status or not. We nonetheless observe a large variability across samples and edges, that could be explained by the fact that some configurations are harder to detect than others. An edge only present on a small segment might for instance be considered absent through the whole series.

Drosophila Life Cycle Microarray Data

The life-cycle of Drosophila melanogaster is punctuated by four main stages of morphological development: embryo, larva, pupa and adult. The expression levels of 4028 genes of wildtype Drosophila were measured by [START_REF] Arbeitman | Gene expression during the life cycle of Drosophila melanogaster[END_REF] at 67 time-points throughout their life-cycle. We have here restricted our attention to eleven genes involved in wing muscle development and previously studied by [START_REF] Zhao | Inferring gene regulatory networks from time series data using the minimum description length principle[END_REF] and [START_REF] Dondelinger | Non-homogeneous dynamic Bayesian networks with Bayesian regularization for inferring gene regulatory networks with gradually timevarying structure[END_REF].

The expectation was that our approach would find change-points corresponding to the four different stages of development observed for Drosophila melanogaster.

We used the normal-Wishart version of the model described in the simulation study. When using the naive prior parameters given in Section 3.6, we obtained poor results (Figure 3.8.a), probably because of the small number of time-points. We noticed that the results could be improved by using data-driven prior specification. We centered the data and set the prior scale matrix φ of the normal-Wishart distribution with α = p + 10 degrees of freedom to φ = (α -p -1) • Σ y where Σ y stands for the sample covariance matrix. By doing this, the normal-Wishart distribution that we get has expectation (0 p , Σ y ). We then obtained the results given in Figure 3.8.b. For this prior, we looked closer to the results for K = arg max K p(K|y) = 5 segments, i.e. one more than the number of development stages. The best segmentation m5 with 5 segments has change-points at positions (19,32,41,53). The posterior probability of observing a change-point at these locations is quite high (Figure 3.8.b). The larva stage is almost exactly recovered, with a shift of one position for the end of the segment. The embryo stage is divided into two segments and the separation between pupa and adult states is missed, the last segment including both adulthood and part of the pulpa stage. These results are nonetheless encouraging. For each segment r of m5 , we computed the posterior edge probability matrix given by (P ({i, j} ∈ E T |y r )) 1 i,j p . On each segment, the prior probability for an edge to appear was set to 0.5 with an approach similar to what was done in Section 3.5. We give a graphical representation of the results in Figure 3.9. In the first segment, fewer edges have large posterior probabilities. However, this higher contrast in probabilities might just be a consequence of this segment being larger than the others.

Finally we compared our results with those obtained by [START_REF] Dondelinger | Non-homogeneous dynamic Bayesian networks with Bayesian regularization for inferring gene regulatory networks with gradually timevarying structure[END_REF] on the same dataset. As for the probability of change-point along time, the results we give in Figure 3.8.b are very similar to those displayed in Figure 12 of this reference. The comparison in terms of inferred networks is more complex as the networks they displayed correspond to the expected stages (embryo, larva, pupa and adult) and not to the one they actually inferred. We found good concordances between the network they inferred for the embryo stage and those that we obtained on segments [1-18] and [19][20][21][22][23][24][25][26][27][28][29][30][31] (both in the embryo stage). We also found similarities at the larva stage (which is close to our inferred [32-40] segment). .9 -Graphical representation of posterior edge probability matrix for each segment of the best segmentation with 5 segments. The width of an edge is proportional to its posterior probability. Edges with probability higher than 0.5 are coloured in blue. Edges with probability lower than 0.2 were not represented.

Embryo

Functional MRI Data

Functional magnetic resonance imaging (fMRI) is commonly used in neuroscience to study the neural basis of perception, cognition, and emotion by detecting changes associated with blood flow. This second application focuses on fMRI data collected by [START_REF] Cribben | Dynamic connectivity regression: Determining state-related changes in brain connectivity[END_REF]. We give a brief description of the experiment but we refer the reader to their article for a more detailed description. Twenty participants were submitted to an anxiety-inducing experiment. Before scanning, participants were told that they would have two minutes to prepare a speech on a subject given to them during scanning. Afterwards, they would have to give their speech in front of expert judges, but they had a "small chance" not to be selected. The subject of the speech was given after two minutes of recording. After two minutes of preparation, participants were told that they would not have to give the speech. The recording continued for two minutes afterwards. A series of 215 images at two-second intervals were acquired during the experiment. [START_REF] Cribben | Dynamic connectivity regression: Determining state-related changes in brain connectivity[END_REF] preprocessed the data and determined five regions of interest (ROIs) in the brain on which the signals were averaged. Thus, we have p = 5 and N = 215, for U = 20 participants. We standardised the data across all participants.

Each participant can be analysed individually by using the same approach as in the previous application. To analyse all participants together, we make the assumption that the dependence structure between the different ROIs of the brain is the same across participants, while being allowed to vary throughout time. Nonetheless, on a given temporal segment, therefore for a given structure, parameters are independently drawn for each participant, so that the likelihood on a segment r can be written as

p(y r ) = T ∈T U u=1 t∈r p(y t,u |θ u )p(θ u |T )dθ u (3.13)
where y t,u stands for the vector of observations at time t for participant u. The distribution p(θ u |T ) and p(y t,u |θ u ) are respectively taken to be normal-Wishart and Gaussian distributions, as in the individual model. In practice, when we tried to perform the inference of the joint model, we were faced with numerical issues, occurring at different levels. The summation over trees was problematic for some segments, especially the largest one. Indeed, we are summing very small quantities and the product over participants in p(y|T ) brings us to deal with quantities of the order of machine precision. Moreover, while searching for the best segmentation can be achieved through log(A) = [log(A s,t )] 1 s,t N +1 , integrating over segmentations requires the actual computation of matrix A. Thus, the exponentiation of the segment log-likelihood matrix leads to other numerical issues. During the dark red interval, the subject of the speech was revealed to participants, who prepared their speech during the light red interval. This preparation is ended by a statement saying that they would not have to give the speech.

1 0.5 0.2 q q q q q 1 2 3 4 5 1-19 q q q q q 1 2 3 4 5 20-38 q q q q q 1 2 3 4 5 39-62 q q q q q 1 2 3 4 5 63-72 q q q q q 1 2 3 4 5

73-111

q q q q q 1 2 3 4 5 112-132 q q q q q 1 2 3 4 5 133-150 q q q q q 1 2 3 4 5 151-187 q q q q q 1 2 3 4 5

188-215

Figure 3.11 -Graphical representation of posterior edge probability matrix for each segment of the best segmentation with K = 9 segments on fMRI data with non-tempered likelihood. The width of an edge is proportional to its posterior probability. Edges with probability higher than 0.5 are coloured in blue.

Our pragmatic answer to these issues was to considered a tempered version of the likelihood given in (3.13):

p * α (y r ) = T ∈T U u=1 t∈r p(y t,u |θ u )p(θ u |T )dθ u 1/α , (3.14)
with α > 1. Tempering the likelihood does not change the mode of the posterior distribution on m, if the matrix a giving prior segment weights is tempered similarly. By doing this, we are actually reducing the effective sample size: a very big α would yield a posterior distribution on m close to the prior.

Figures 3.10a and 3.10b sum up the results obtained participant per participant for change-point location. They vary a lot across participants, as shown by the five given examples, as well as the mean and standard deviation curves. The left panel of Figure 3.10c shows the posterior probability of observing a change-point when participants are jointly considered with a tree-structured model or with a non-structured model, with likelihood tempered at α = U /2 = 10 and α = U = 20. For both values of α, the profiles are quite similar, with an expected more peaked behaviour for α = 10. The strongest peak is observed during the announcement of the speech topic. The right panel of Figure 3.10c gives the posterior distribution of K for both models and for different values of α. We observe flatter distributions for the full model, with a mode at 11 segments. In the tree-structured model, 9 segments are selected. For this value of K, we looked at the best segmentation and computed the posterior edge probability matrices for its segments. A graphical representation of the results is given in Figure 3.11. [START_REF] Cribben | Dynamic connectivity regression: Determining state-related changes in brain connectivity[END_REF] retrieved 8 segments from these data. There is no clear correspondence between our segmentation and theirs. A remark that can nonetheless be made is that, in our case, each change-point is associated with a clear change in the topology of the network. These structure changes are less obvious in [START_REF] Cribben | Dynamic connectivity regression: Determining state-related changes in brain connectivity[END_REF]. This might be a consequence of our model explicitly modelling the structure, thus encouraging change-points to mark out abrupt changes in structure rather than in parameters. In this paper, we showed how exact Bayesian inference could be achieved for change-points in the structure of a multivariate time-series with careful specification of prior distributions. Essentially, prior distributions have to factorise over both segments and edges. For the sake of clarity, we assumed that, within a segment r, observations Y t were independent conditionally on T and θ. While convenient and leading to comfortable formulas, this independence assumption is hardly realistic in many applied situations, including those that we have considered here. Yet, time dependency could be considered within segments, as long as p(y r |T ) still factorises over the edges of T . One could for instance consider using the work of [START_REF] Siracusa | Tractable Bayesian Inference of Time-Series Dependence Structure[END_REF] to achieve this. Trees would then be used to model the dependences between two consecutive times instead of instantaneous dependences.

Discussion

The framework that we have described is also convenient for Bayesian model comparison. When one is faced with an alternative in modelling, Bayes factors between two models are easily obtained, as fully marginal likelihood can be computed exactly and efficiently. For instance, the question of whether changes should be allowed in the mean of a Gaussian distribution or not can be addressed by computing p(y) in both cases and by looking at their ratio. This is by no mean specific to our approach, but exact computation makes it completely straightforward.

The exactness of the inference also creates a comfortable framework to precisely study the effect of the prior distribution on segmentations. Once again, as the inference does not rely on stochastic integration, the impact of prior specification could be evaluated at low cost and in an exact manner.

We finish this discussion by mentioning numerical issues. As explained in Section 3.7, when the number of observations increases, we have to deal with elementary probabilities that differ from several order of magnitudes. Because the summations over the huge spaces of both segmentations and trees are carried out in an exact manner, these quantities have to be added to each others, resulting in numerical errors. Obviously, no naive implementation would work and some of these errors can be avoided with careful and skilful programming. At this stage, this is still not sufficient and the likelihood tempering approach that we propose is not satisfying. Further numerical improvements could be considered such as the systematic ordering of the terms when computing a determinant in a recursive way.

The R code used in the simulations and the applications is available from the authors upon request. A package will soon be available from the Comprehensive R Archive Network. In the previous chapters, we used distinctive features pertaining to the sets of spanning trees and segmentations, and leading to tractable inference in various models. A blatant similarity between these sets is that they are both made of elements that are intrinsically modular. A spanning tree is a collection of edges, just as a segmentation is collection of segments. Modularity is actually one of the key features allowing for an efficient browse through these large collections. The algorithms used to sum over spanning trees and segmentations assume that the function to be summed up factorises according to these basic units, and express the global constraints on the collections of units in a smart algebraic fashion. A spanning tree on p vertices is a collection of p -1 edges, with the additional constraint that the resulting graph is connected. This connectivity assumption can be encoded as a condition on the incidence matrix corresponding to the collection of edges, and through this property, the summation can be performed by computing a single determinant of size p -1. Similarly, in the algorithm developed by [START_REF] Rigaill | Exact posterior distributions and model selection criteria for multiple change-point detection problems[END_REF] for segmentations, a crucial ingredient is that the matrix giving the weights of all segments is a strictly upper triangular matrix. The K-th power of this matrix is used to sum over M K and all irrelevant terms in the summation are set to zero thanks to this property.

All of this is easy enough to notice when we are presented with both examples, but is not quite helpful in finding other situations for which such algebraic tricks can be conjured. Nonetheless, it can further be noticed that maximisation and summation problems are closely linked in this matter, in the sense that they are often jointly tractable or not. Indeed, the maximum (or minimum) spanning tree problem can be solved through Kruskal's or Prim's algorithms, and classic dynamic programming can be used to obtain the best segmentation. This has to do with the fact that max-sum and sum-product problems can often be seen as two particular instances of a more generic problem expressed with abstract operators. This link could be useful to dig out other interesting situations in which exact inference is tractable.

As our goal was to perform Bayesian inference on models involving spanning trees and segmentations, we had to consider distributions on these sets. Their modularity obviously lead to the exponential family with a binary vector indexed by all possible units (either edges or segments) as sufficient statistics. The normalising constant, or partition functions, of these distributions could then be computed through the aforementioned algorithms, as the factorisation assumptions are obviously satisfied. Of course, this is what motivated our considering these particular sets in the first place. For many distributions in the exponential family, computing the normalising constant is not possible in a reasonable amount of time. There is for instance no explicit form for this constant in the case of decomposable graphs.

Models were then built so as to preserve this factorisation property for posterior distributions. If Q either denotes a tree or a segmentation, the distribution π of the observations conditionally on Q was assumed to factorise over the elements of Q. So was the prior distribution on π (or on its parameters θ). Under these assumptions, expressed through (hyper) Markov properties, the integration over π or θ could be performed locally and separately from the summation on Q. We merely assumed that these low-dimensional integrals could be computed through conjugacy and explicit formulas, for the sake of exact inference. If one has to stray from conjugate priors, numerical or stochastic integration can be considered, as well as approximation techniques such as Laplace's approximation (see Schwaller & Robin, 2015).

Among the various assumptions that we have made in the models, some were fundamental for the proper functioning of the approach, and some were merely convenient, leaving room for some interesting tweaks. This chapter is meant to provide some extensions and perspectives to the work presented before.

The first section explains how covariates can be taken into account in what has been previously presented. When observations are assumed to be normally distributed, adapting the model is straightforward. We show that conjugate priors for the multivariate multiple linear regression induce hyperdistributions that can be used as a basis to build a compatible family of hyperdistributions. We also derive an expression for the marginal likelihood of the observations on any subset of vertices. Whenever observations cannot be modelled by a multivariate normal distribution, we suggested a heuristic approach based on copulas, that we illustrate on microbial ecology data.

In the second section, we consider a different set of graphs for which exact structure inference is possible, namely the set of directed trees. This set inherits many of the interesting properties of its undirected counter-part, as all directed trees can be obtained from an undirected tree by choosing a root. We use these graphs to introduce temporal dependence between observations. The independence assumption is indeed hardly realistic in many applied situations. Our suggestion is to use the Temporal Independence Models (TIMs) introduced by [START_REF] Siracusa | Tractable Bayesian Inference of Time-Series Dependence Structure[END_REF], for which structure inference can also be performed with cubic complexity with respect to the number of vertices, as long as possible structures are restricted to directed trees. More than introducing temporal dependence, we actually exhibit a slightly different example of discrete set on which summation is possible through algebra.

Finally, we investigate prior specification on segmentations. We make out a non-exhaustive bestiary of priors fitting in our framework. We also describe a different kind of prior in which taking into account an initial guess on the segmentation, or information coming from a previous study, is straightforward. These prior distributions still factorise over segments, but the contribution of a segment to a given segmentation depends on its position within the segmentation. Previously, the position of a segment had no influence on its contribution. The algorithm used to perform the inference has to be adapted to this new kind of prior distribution. In some situations, observations of a set of variables come with observations of associated variables that, while not being of direct interest, have to be taken into account in the analysis. The latter are often called covariates. We were drawn to such problems through a collaboration with Corinne Vacher (UMR Biogeco, Biodiversity, Genes & Communities, Université de Bordeaux) on microbial ecology data.

Plan-inhabiting micro-organisms can interact with each other forming complex interaction networks. These interactions can either be direct (predation, parasitism, mutualism, etc) or mediated by the environment. The first kind of interactions are especially interesting when trying to understand how a new microbial species might integrate in an existing ecosystem. The microbial ecologists we worked with were interested in better understanding powdery mildew, one of the most common tree diseases in European forests caused by a fungus called Erysiphe alphitoides. To that end, they measured the abundances of the microbial species that could be found on 120 leaves taken from three different oak trees. For each leaf, many covariates were also measured, including the distance of the leaf to the base of the branch, to the trunk or to the ground. Our goal was to infer the network giving the direct interactions between species while taking these covariates into account.

We show that, for multivariate normal distributions, including covariates in the framework developed in Chapters 2 and 3 is straightforward. Nonetheless, this is basically made possible by specific properties of normal distributions. In the general case, the hyper Markov property for the prior on the distribution of the observations is hard to retain while including covariates. For the count data described above, we used an approach based on copulas as a pragmatical solution. The results were given in an article published in the journal Microbial Ecology (Jakuschkin et al., 2016).

Multivariate Multiple Linear Regression

We consider a setting where we have N observations of p real-valued variables of interest and k real-valued covariates. These observation are denoted by (Y

(n) , X (n) ) ∈ R p × R k ,
1 n N and gathered in matrices Y and X: 1) . . . 1) . . .

Y =    Y (
Y (N )    , X =    X (
X (n)    .
For U ⊂ V , we let Y U and B U respectively denote matrices Y and B restricted to the columns corresponding to the variables in U . Our variables of interest are now the univariate series given by column vectors Y i , i ∈ V , for which we only have one observation, along with the associated covariates. The sample space we are working on is

X = i∈V X i = i∈V R N .
A distribution π on X can be seen as a distribution on R N p , which is obviously isomorphic to X . Similarly, if π is a distribution on X and A is a subset of V of size a, the marginal distribution π A is basically a distribution on R ap .

Let us consider the regression model given by

Y = XB + E, vec(E) ∼ N (0, Λ -1 ⊗ I N ),
where "vec" is the linear transformation stacking the columns of a matrix on top of one another, and ⊗ stands for the Kronecker product of two matrices. Matrix B is a k × p real-valued matrix giving the regression coefficients of the model. Matrix Λ is the precision matrix of size p × p describing the interactions between the different series. This model can be rewritten as vec(Y)|X, B, Λ ∼ N (vec(XB), Λ -1 ⊗ I N ).

(4.1)

We remind vec(XB) = (I p ⊗ X)vec(B).

Any prior distribution on (B, Λ) induces a hyperdistribution ρ on the distribution π(•|X) of Y. Classic results on conjugate priors for regression models suggest this choice of prior distributions:

Λ ∼ W(α, Ψ), (4.2) vec(B)|Λ ∼ N (vec(B 0 ), Λ -1 ⊗ K -1 0 ).
where we remind that W(α, Ψ) is the Wishart distribution with α degrees of freedom and parametric matrix Ψ. In the following, we use the shorthands

β = vec(B), β A = vec(B A ), ∀A ⊆ V, β 0 = vec(B 0 ), β 0 A = vec((B 0 ) A ), ∀A ⊆ V, = vec(E).

4.1.1.a Hyper Markov Property

We now show that any distribution for (B, Λ) of the form given in (4.2) induces a hyperdistribution that is strong hyper Markov with respect to the complete graph, and can therefore be use to build a compatible family of strong hyper Markov hyperdistributions.

Proposition 4.1. Let R ⊂ V and S = V \ R, of respective sizes r and s = p -r. Then

vec(Y R )|X, β, Λ ∼ N ((I r ⊗ X)β R , (Λ R•S ) -1 ⊗ I N ), (4.3) vec(Y S )|Y R = y r , X, β, Λ ∼ N (µ S•R , Λ -1 SS ⊗ I N ), (4.4)
where [Λ RR , Λ SR , Λ RS , Λ SS ] is a partitioning of Λ according to (R, S) and where

Λ R•S = Λ RR -Λ RS Λ -1 SS Λ SR , µ S•R = (I s ⊗ X)β S -(Λ -1 SS Λ SR ⊗ I N ) (vec(y R ) -(I r ⊗ X)β R ) . Proof. Let R ⊂ V and S = V \ R
of size r and s = p -r respectively. Without loss of generality, we suppose that β and Λ can be written as

β = β R β S ; Λ = Λ RR Λ RS Λ SR Λ SS .
The mean vector of vec(Y) can be expressed using β R and β S :

vec(XB) = (I p ⊗ X)β = I r ⊗ X 0 0 I s ⊗ X β R β S = (I r ⊗ X)β R (I s ⊗ X)β S .
If Σ . . = Λ -1 and considering a partitioning of Σ similar to Λ, the covariance and precision matrices of vec(Y) are respectively given by

Λ -1 ⊗ I N -1 = Λ ⊗ I N = Λ RR ⊗ I N Λ RS ⊗ I N Λ SR ⊗ I N Λ SS ⊗ I N ; Λ -1 ⊗ I N = Σ ⊗ I N = Σ RR ⊗ I N Σ RS ⊗ I N Σ SR ⊗ I N Σ SS ⊗ I N .
and Σ -1 RR = Λ R•S . Classic results on the marginal and conditional distributions of a multivariate normal distribution conclude this proof. Proposition 4.2. Let R ⊂ V and S = V \ R. Then

(β R , Λ R•S ) |= β S + Λ -1 SS Λ SR ⊗ I k β R , Λ SS , Λ RS . (4.5)
Proof. The probability density function of (β, Λ) is given by

f (β, Λ) ∝ |Λ| (α-p+k-1)/2 exp - 1 2 (β -β 0 ) (Λ ⊗ K 0 )(β -β 0 ) exp - 1 2 tr (ΨΛ) . (4.6)
Let us define

β S•R . . = β S + Λ -1 SS Λ SR ⊗ I k β R . β 0 S•R
is defined similarly. We want to make the following change of variables in (4.6):

(β R , β S , Λ RR , Λ SS , Λ RS ) -→ (β R , β S•R , Λ R•S , Λ SS , Λ RS )
The Jacobian matrix of this change of variables can be written as [START_REF] Massam | On Transformations and Determinants of Wishart Variables on Symmetric Cones 1[END_REF] showed that |Z| = 1 . So we only need to express (4.6) as a function of

β R β S•R (Λ R•S , Λ SS , Λ RS ) β R I rk Λ -1 SS Λ SR ⊗ I k 0 β S 0 I sk 0 (Λ RR , Λ SS , Λ RS ) 0 * Z = J. (4.7) So |J| = |Z|
(β R , β S•R , Λ R•S , Λ SS , Λ RS ). f (β, Λ) ∝ |Λ R•S | (α-p+k-1)/2 exp - 1 2 (β R -β 0 R ) (Λ R•S ⊗ K 0 )(β R -β 0 R ) exp - 1 2 tr (Ψ RR Λ R•S ) • |Λ SS | (α-p+k)/2 exp - 1 2 (β S•R -β 0 S•R ) (Λ SS ⊗ K 0 )(β S•R -β 0 S•R ) (4.8) • exp - 1 2 tr Ψ RR Λ RS Λ -1 SS Λ SR + 2tr (Ψ SR Λ RS ) + tr (Ψ SS Λ SS )
Hence the independence of (β R , Λ R•S ) and (β S•R , Λ SS , Λ RS ). We also get the distribution of (β R , Λ R•S ):

Λ R•S ∼ W(α + k -r, Ψ RR ), (4.9) β|Λ R•S ∼ N (β 0 R , Λ -1 R•S ⊗ K -1 0 ).
Proposition 4.3. The hyper distribution ρ on π(•|X) induced by the joint distribution of (Λ, β) is strong hyper Markov with respect to the complete graph on V .

Proof. We have to show that, for R ⊆ V and

S = V \ R, it holds that π R (•|X) |= π S|R (•|X). (4.10)
These distributions are given by Proposition 4.1 and it is enough to notice that

µ S•R = (I s ⊗ X) β S + (Λ -1 SS Λ SR ⊗ I k )β R -(Λ -1 SS Λ SR ⊗ I k )vec(y R ) = (I s ⊗ X)β S•R -(Λ -1 SS Λ SR ⊗ I k )vec(y R )
to get (4.10) from Proposition 4.2.

We can therefore put Proposition 1.9 to good use and build a compatible family of strong hyper Markov hyperdistributions from the hyperdistribution ρ induced by the distribution on (B, Λ) given in (4.2).

4.1.1.b Marginal Likelihood

Now that we have our compatible family of strong hyper Markov hyperdistributions, we have to compute the marginal likelihood of the observations on any subset of V . For A ⊆ V , it can be written as

p(Y A |X) = π A (Y A |X)ρ A (π A )dν A (π A ).
We are in fact mainly interested in subsets of size 1 and 2 but we derive the formula for any subset.

Proposition 4.4. The marginal likelihood of Y is given by

p(Y|X) = 1 π pN 2 Γ p α+N 2 Γ p α 2 |K 0 | p 2 |K N | p 2 |Ψ| α 2 |Ψ (N ) | α+N 2
where

Ψ (N ) . . = Ψ + B 0 K 0 B 0 + B N K N B N + Y Y, K N . . = K 0 + X X, (4.11) B N . . = K -1 N (X Y + K 0 B 0 ).
Proof. We have that p(Y, Λ, B|X)

= |Λ| α-p+k+N -1 2 |K 0 | p 2 |Ψ| α 2 2 αp 2 (2π) p(N +k) 2 Γ p α 2 • exp - 1 2 tr (Λ [(Y -XB) (Y -XB) + (B -B 0 ) K 0 (B -B 0 ) + Ψ]) .
If K N and B N are defined as in 4.11, then

(Y -XB) (Y -XB) + (B -B 0 ) K 0 (B -B 0 ) = (B -B N ) K N (B -B N ) + B N K N B N + B 0 K 0 B 0 + Y Y.
In the right-hand side of the equation above, only the first term depends on B and the corresponding term in p(Y, Λ, B|X) is a Gaussian kernel. We therefore get that

p(Y, Λ|X) = |Λ| α-p+N -1 2 |Ψ| α 2 2 αp 2 (2π) pN 2 Γ p α 2 |K 0 | p 2 |K N | p 2 exp - 1 2 tr (Λ [B N K N B N + B 0 K 0 B 0 + Y Y + Ψ]) . = |Λ| α-p+N -1 2 |Ψ| α 2 2 αp 2 (2π) pN 2 Γ p α 2 |K 0 | p 2 |K N | p 2
exp -1 2 tr ΛΨ (N ) .

In the equation above, we recognise that

|Λ| α-p+N -1 2 exp - 1 2 tr ΛΨ (N )
is the unnormalised density of a Wishart distribution with α + N degrees of freedom and parametric matrix Ψ (N ) , hence the sought-after result.

Proposition 4.5. For A ⊆ V of size a, the marginal likelihood of Y A is given by

p(Y A |X) = 1 π aN 2 Γ a α-p+a+N 2 Γ a α-p+a 2 |K 0 | a 2 |K N | a 2 |Ψ AA | α-p+a 2 |Ψ (N ) AA | α-p+a+N 2 .
Proof. It follows from (4.9) and Proposition 4.4.

We remind that network inference in tree-structured graphical models relies on a posterior edge weight matrix given by

ω i,j = β i,j p(Y i , Y j ) p(Y i )p(Y j ) , ∀{i, j} ∈ P 2 (V ),
where β is a prior edge weight matrix. From Proposition 4.5, we get that

p(Y i , Y j |X) p(Y i |X)p(Y j |X) ∝ |Ψ (N ) i,i | α-p+1+N 2 |Ψ i,i | α-p+1 2 |Ψ (N ) j,j | α-p+1+N 2 |Ψ j,j | α-p+1 2 |Ψ {i,j}{i,j} | α-p+2 2 |Ψ (N ) {i,j}{i,j} | α-p+2+N 2
.

For network inference, knowing p(Y i , Y j |X)/p(Y i |X)p(Y j |X) up to a constant identical for all edges is enough. For segmentation however, the remaining terms have to be included as they depend on the number of observations and the covariates, which vary depending on the segment that is considered.

Copulas: a Pragmatical Approach

We now come back to the reason we considered covariates in the first place, that is to say the microbial ecology data gathered on oak trees that we described at the beginning of this chapter. These data were obtained by high-throughput sequencing of the DNA that could be found on the collected leaves. For each of the N leaves, the abundance of a microbial species was given by a number of reads. We are therefore dealing with multivariate count data. If we consider that all covariates are real-valued for the sake of simplicity, we are now in a situation where (Y (n) , X (n) ) ∈ N p × R k , 1 n N . Y and X are defined as in the previous section.

4.1.2.a A Model for Multivariate Count Data based on Copulas

Our first idea to model multivariate count data with covariates was to rely on univariate generalized linear models (GLMs) weaved together with a copula. We wanted to remain as close as possible to the Gaussian case, and we hoped to get similar results when using Gaussian copulas. But we will see that the strong hyper Markov property that we showed for multivariate multiple regression does not hold in this particular model.

For i ∈ V , we model the abundance of species i on all leaves conditionally on X as follows:

Y (n) i ∼ P(λ n i ), log(λ n i ) = X (n) B i , {Y (n)
i } 1 n N independent |X. (4.12)

We could also have modelled Y i |X with a negative binomial GLM. The dependence between the different Y i is then handled by a Gaussian copula. Let Φ -1 denote the inverse cumulative density of the univariate standard normal distribution, and ψ Λ p the probability density function of the multivariate normal distribution with mean vector zero and partial correlation matrix Λ. The probability density function of the Gaussian copula C Λ p with partial correlation matrix Λ is given by The distributions of the form given in (4.13) do not form a strong meta Markov model with respect to the complete graph on V . We will therefore not be able to build a compatible family of strong hyper Markov hyperdistributions for this model. Indeed, let A ⊆ V and B = V \ A. Evaluating p(y B |y A , Λ, B, X) requires more information about what happens on A than just y A . We actually need to know P(•|λ n i ), i ∈ A, 1 n N . This means that we cannot mimic what we have done in the multiple linear regression framework.

c Λ p (u) = ψ Λ p (Φ -1 (u 1 ), . . . , Φ -1 (u p )), ∀u ∈ [0; 1] p .

4.1.2.b Pragmatical Approach

All our attempts to model multivariate count data while taking covariates into account ended up leading to models that were not strong meta Markov. We therefore decided to use the copula-based model described above, while dropping the Bayesian paradigm for the univariate part of the model involving GLMs. In so doing, we betrayed our committed stance on fully exact Bayesian inference, but this moral lapse brought us back to the simple framework of tree-structured copulas described in Section 2.4.1.a.

The approach splits up in two steps. The first one deals with covariates by estimating the regression coefficients in the univariate GLMs described in (4.12). For all i ∈ V , we obtain an estimation Bi of B i . The effect of covariates is then eliminated by computing the Pearson residuals associated to each fitted model. For Poisson GLMs, they are given by

R (n) i = Y (n) i - λ(n) i λ(n) i , λ(n) i = exp(X (n) Bi ), ∀i ∈ V, ∀n ∈ 1; N .
We let Fi denote the empirical cdf of {R

i } N n=1 and we define

Ỹ(n) i . . = Fi (R (n) i ), ∀i ∈ V, ∀n ∈ 1; N .
Having (somehow) dealt with covariates in the first step, we come back to Bayesian inference. We consider that { Ỹ(n) } N n=1 is an independent sample drawn from a distribution on [0; 1] p with uniform marginal distributions. The second step is a direct application of Section 2.4.1.a. The model that we used is the following:

T ∼ ξ, (4.14) Λ i,j |T ∼    δ 1 if i = j, U([-1; 1]) if {i, j} ∈ E T , δ 0 otherwise, {Λ i,j } 1 i,j p independent |T, Ỹ|Λ ∼ C Λ p .
In this model, network inference can be performed as described in Chapter 2. The results obtained for Ỹ are interpreted as a proxy for the desired results on Y|X.

The main flaw of this pragmatical approach is that no uncertainty is taken into account at the regression level, as residuals of the GLMs are plugged in the tree-copula model through their empirical cdfs. We are trudging troubled waters mixing both frequentist and Bayesian paradigms. But what is lost in clarity and exactness is partially made up for in flexibility. Indeed, both the univariate GLMs and the copula can be tailored to the specific needs pertaining to a given modelling problem. If need be, each GLM can be individually adapted, and so can the bivariate copulas describing the pairwise interactions in Ŷ.

4.1.2.c Results for Microbial Ecology Data

We briefly present the results obtained on the microbial ecology data through the copula approach described above. Data were collected from 120 leaves taken from three oaks trees presenting different observed susceptibility to powdery mildew. DNA was extracted from these leaves and sequenced. After standard preprocessing, for each leaf, we were given the number of reads per species or, more accurately, operational taxonomic unit (OTU). We were also given the distance of the leaf to the base of the branch, to the trunk and to the ground, as well as the orientation (either south-west or north-east) of the branch.

A subset of the most abundant OTUs was used for the analysis. These OTUs represented at least 0.5% of the sequences for at least one of the three oak trees. We retained 114 OTUs, among which 48 were bacterial and 66 were fungal (including the pathogen responsible for powdery mildew, Erysiphe alphitoides). In most cases, the number of reads for each OTU was fitted with a negative binomial GLM with environmental variables as predictors and the total number of reads per sample (log-transformed) as an offset. Poisson GLMs were used for a few OTUs. The analysis was performed for each tree separately, as well as for all trees together with an additional covariate indicating the origin of the leaf.

The network inference model was chosen as presented in (4.14), with ξ put to the uniform distribution on T . Posterior edge probabilities were computed and the prior appearance probability of each edge was set to 0.5 as described in Section 2.4.3. For each edge, the sign of the interaction was determined by looking at the maximum a posteriori estimate for the partial correlation of the corresponding bivariate copula.

Figure 4.1 shows the posterior probability matrices obtained when all trees are jointly considered and for each individual tree. The former is more contrasted as a result of the grouped analysis working with more data. The block-structure that can more or less be seen on all four matrices indicates that intra-kingdom interactions (fungus to fungus or bacteria to bacteria) are more frequent than inter-kingdom interactions. Of particular interest were the interactions of Erysiphe alphitoides with other OTUs on the susceptible tree, as the main goal of the study was to get a better understanding of the microbial mechanisms behind powdery mildew. Truncating the posterior edge probability matrix at the level set for prior appearance probability, namely 0.5, we found that E. alphitoides was involved in 26 interactions, among which 13 with fungal OTUs (6 positive, 7 negative) and 13 with bacterial OTUs (6 positive, 7 negative). These interactions are represented in Figure 4.2. All but one of these interactions were previously considered likely by microbial ecologists. In the models described in Chapters 2 and 3, observations are assumed to be independent given the parameters. As discussed in Section 3.8, this assumption is hardly realistic in many practical cases. Let us take the example of the fMRI data treated in Section 3.7.2. Such data are typically recorded with high temporal resolution. Here, an image was taken every two seconds, but it is not unusual to sample at higher frequencies. With observations this close to one another in time, it would desirable for the state of the brain of a participant at time t to directly influence on the state of its brain at time t+1. In order to be entirely satisfactory, the dependence model should preserve some factorisation property with respect to edges, so that we can keep using the Matrix-Tree theorem to efficiently integrate on structures. Directed graphical models prove to be more suited to this situation than their undirected counter-parts. [START_REF] Siracusa | Tractable Bayesian Inference of Time-Series Dependence Structure[END_REF] introduced Temporal Interaction Models as a way to include structured temporal dependence in the modelling of multivariate time-series. When the time dependence structure of these models is taken to be a directed tree, a directed version of the Matrix-Tree theorem can be used to perform calculation in the same manner as in the undirected case.

When considering the segmentation model, the factorability assumption of [START_REF] Rigaill | Exact posterior distributions and model selection criteria for multiple change-point detection problems[END_REF] is not violated as long as there is no temporal dependence between segments. Within segments however, there is no problem introducing dependence if we are able to compute the marginal likelihood.

Temporal Interaction Models

Let y = {y t } N t=1 be a realisation of a multivariate random process Y = {Y t } N t=1 of dimension p 2. For 1 t N , Y t = (Y t 1 , ..., Y t p ) is a multivariate random variable taking values in a product of measurable spaces X = p i=1 X i . We let D denote the set of DAGs on V . For D ∈ D and i ∈ V , we let pa(i, D) denote the parents of vertex i in D. When there is no ambiguity on the the DAG, we drop D from the notation.

Directed graphical models are expressed in terms of conditional distributions rather than marginal distributions. A conditional distribution is a collection of distributions indexed by the values that can be taken by the conditioning variables. Such a collection is called a Markov kernel. Definition 4.1. Let (X, A) and (Z, B) be two measurable spaces. A Markov kernel from X to Z is a map κ : X × B -→ [0; 1] such that,

• for all B ∈ B, the map x → κ(x, B) is A-measurable,

• for all x ∈ X, the map B → κ(x, B) is a probability measure on ((Z, B).

We begin by describing Temporal Interaction Models with unitary lag or TIM(1). Such models are given by a DAG D and a set of kernels κ = {κ i|pa(i,D) } i∈V , where, for i ∈ V , κ i|pa(i,D) is a kernel from X pa(i,D)∪{i} to X i . We also need a set of initial conditions y 0 ∈ X . The likelihood of y is then given by p(y|D, κ, y 0 ) = N t=1 i∈V κ i|pa(i) (y t-1 pa(i)∪{i} , y t i ).

The kernels in κ describe the distribution of Y t |Y t-1 and this distribution factorises according to D, so that Y t i depends on Y t-1 i and Y t-1 pa(i) . In a parametric framework, these kernels would just be a set of parameters {θ i|pa(i) }, where θ i|pa(i) specifies the distribution of

Y t i |Y t-1 i , Y t-1 pa(i)
. We refer to D as the dependence structure of the model.

Temporal Interaction Models with an arbitrary lag , or TIM( ) in short, are obtained similarly, but Y t is assumed to depend on (Y t-l , . . . , Y t-1 ). For i ∈ V , κ i|pa(i) is now a kernel from (X pa(i)∪{i} ) to X i . We also need to provide extended initial conditions (y -+1 , . . . , y 0 ). The graphical representation of an example of TIM(1) and TIM(2) sharing the same dependence structure is given in Figure 4.3.

Y t-1 Y t Y t+1 Y t-1 Y t Y t+1 Y t-1 2 Y t-1 1 Y t-1 3 Y t 1 Y t 2 Y t 3 Y t+1 1 Y t+1 2 Y t+1 3 Y t-1 2 Y t-1 1 Y t-1 3 Y t 1 Y t 2 Y t 3 Y t+1 1 Y t+1 2 Y t+1 3 Y 2 Y 1 Y 3 TIM(1) TIM(2)

Bayesian Inference of Tree-structured TIM

We limit our attention to TIMs whose dependence structures are directed trees, i.e. connected DAGs in which each vertex has at most one parent. We remind that -→ T denotes the set of such DAGs. Performing Bayesian inference in these models requires to define prior distributions on D and κ. If these priors factorise properly, the marginal likelihood of the observations, integrated on both D and κ, can be obtained in polynomial time.

4.2.2.a Prior Distribution on D

As D ∈ -→ T , for i ∈ V , pa(i) is either empty or reduced to a single vertex. A distribution on -→ T can therefore be given by a p × p matrix β through

p(D) = 1 -→ Z (β) i∈V β pa(i),i , ∀D ∈ -→ T ,
where, ∀i ∈ V , β ∅,i = β i,i . Contrary to the undirected case, β is not necessarily symmetric. The normalising constant -→ Z (β) is given by

-→ Z (β) = D∈ - → T i∈V β pa(i,D),i .
The sum in -→ Z (β) involves p p-1 terms. Indeed, we have seen that there are p p-2 undirected spanning trees on p vertices, and each of these trees can be oriented in p different ways by choosing a root and directing all edges away from this root. As stated by [START_REF] Siracusa | Tractable Bayesian Inference of Time-Series Dependence Structure[END_REF], the summation over the trees rooted at a vertex r can be performed in O(p 3 ) time thanks to the directed version of the Matrix-Tree theorem. We define the Laplacian matrix associated to a directed edge weight matrix β as the matrix whose general term is given by

-→ ∆ i,j (β) = -β i,j if i = j, p l=1 β l,j -β i,i if i = j.
If we let -→ T r , r ∈ V , denote the set of directed trees rooted at r, Theorem 1.7 directly yields that

-→ Z r (β) . . = D∈ -→ Tr (i,j)∈E D β i,j = | -→ ∆ {r} (β)|, ∀r ∈ V,
where we remind that -→ ∆ {r} (β) is the matrix obtained from -→ ∆(β) by removing row and column r. Using this formula to compute

-→ Z (β) yields a time complexity of O(p 4 ), as -→ Z (β) = r β r,r -→ Z r (β).
It can in fact be obtained in O(p 3 ) time.

Proposition 4.6 (Koo et al., 2007, Prop. 1). Let -→ ∆(β) be the Laplacian matrix associated to a directed edge weight matrix β. We denote -→ ∆ * (β) the matrix obtained by replacing the first row in

-→ ∆(β) by [β 1,1 , . . . , β p,p ]. Then it holds that -→ Z (β) = | -→ ∆ * (β)|.

4.2.2.b Prior Distribution on κ

The distribution on κ, conditionally on D, is taken so that the kernels in {κ i|pa(i,D) } i∈V are mutually independent. The distribution for each individual kernel is also chosen once and for all dependence structures. In a parametric framework, [START_REF] Siracusa | Tractable Bayesian Inference of Time-Series Dependence Structure[END_REF] refer to these properties as parameter independence and parameter modularity. Under these assumptions, the distribution of κ conditionally on D is specified for all D ∈ -→ T through a collection of distributions for individual kernels {κ i|j } (i,j)∈V 2 where κ i|j is a kernel from (X {i,j} ) , being the lag of the model, to X i . By convention, κ i|i is the kernel associated to pa(i) = ∅. In a Gaussian setting, these kernels are linear regressions models with Y t i as response variable and (Y t-1 i , Y t-1 j ) as predictor variables, and the prior distributions on the regression coefficients and the covariance matrix of the noise are drawn independently between the different models. Assume that, for (i, j) ∈ V 2 the distribution of κ i|j admits a density ς i|j with respect to some measure ν i|j on the kernels from (X {i,j} ) to X i . Then the density of the distribution of κ|K with respect to the appropriate product measure is given by

ς(κ|D) = i∈V ς i|pa(i) (κ i|pa(i,D) ).
When the prior on κ is chosen of this form, the posterior distribution on D continues to factorises on directed edges and can be expressed as

p(D|y) = 1 -→ Z (ω) i∈V ω pa(i),i , ∀D ∈ -→ T ,
where, by convention, ω ∅,i = ω i,i , and

∀(i, j) ∈ V 2 , ω i,j = β i,j N t=1
κ j|i (y t- i,j , . . . , y t-1 i,j , y t j ) ς j|i (κ j|i )dν j|i (κ j|i ).

Matrix ω contains the posterior weights associated to every possible directed edges that can be borrowed by a tree in -→ T .

The kernels in κ define a distribution for (Y 1 , . . . , Y p ), so that ς induces a hyperdistribution on the set of distributions on X N = i∈V X N i . In Chapter 1, we have not defined hyper Markov properties for DAGs, but in this case, parameter independence is actually an other name for the strong hyper Markov property with respect to a DAG. Definition 4.2. Let D be a DAG on V . A hyperdistribution ρ over a random distribution π is said to be strong hyper Markov with respect to D, if it holds that {π v|pa(v) } v∈V are mutually independent under ρ.

4.2.2.c Computing Posterior Edge Probabilities

We now give a proof a theorem, similar to Theorem 2.3, on the computation of posterior edge probabilities in the case of directed trees. This result is mentioned in [START_REF] Koo | Structured prediction models via the matrix-tree theorem[END_REF]. In this section, we drop ω from the notations and

-→ ∆, -→ ∆ * , -→ Z are used as shorthands for -→ ∆(ω), -→ ∆ * (ω), -→ Z (ω).
We begin by proving the following lemma on the derivatives of -→ Z .

Lemma 4.1. Let -→ M be the matrix whose general term is given by

-→ M k,l =    ( -→ ∆ * ) -1 1,l if k = l, ( -→ ∆ * ) -1 l,l -(1 -δ 1k ) ( -→ ∆ * ) -1 k,l if l = k.
Then, it holds that, for {k, l} ∈

V 2 , ∂ -→ Z ∂ω k,l = -→ M k,l • | -→ ∆ * |. Proof. By Proposition 4.6, -→ Z = | -→ ∆ * |. We assume that -→ Z = 0, so that -→ ∆ * is non-singular and ∂| -→ ∆ * | ∂ -→ ∆ * i,j = | -→ ∆ * | • ( -→ ∆ * ) -1 i,j , ∀(i, j) ∈ V 2 . For k = l and k = 1, the only elements of -→ ∆ * that depend on ω k,l are -→ ∆ * k,l = -ω k,l and -→ ∆ {r} l,l = u =l ω u,l , so that ∂ -→ Z ∂ω k,l = ∂| -→ ∆ * | ∂ω k,l = (i,j)∈V 2 ∂ -→ ∆ * ∂ -→ ∆ * i,j ∂ -→ ∆ * i,j ∂ω k,l = | -→ ∆ * | ( -→ ∆ * ) -1 l,l -( -→ ∆ * ) -1 k,l .
The result is similarly obtained for k = 1, k = l and k = l.

Theorem 4.1. Posterior edge probabilities can be computed for all edges at once in O(p 3 ) time from the posterior edge weight matrix ω as

P ({k, l} ∈ E D |y) = ω k,l • -→ M k,l ,
where -→ M is defined as in Lemma 4.1

Proof. This proof follows the outline of the proof written for Theorem 2.3 in the case of undirected trees. Let (k, l) be a directed edge. We define

-→ Z + . . = D∈ - → T (k,l)∈E D i∈V ω pa(i,D),i , -→ Z -. . = D∈ - → T (k,l) ∈E D i∈V ω pa(i,D),i , so that -→ Z = -→ Z + + -→ Z -and P ({k, l} ∈ E D |y) = -→ Z + / -→ Z . Using Lemma 4.1, as -→ Z -does not depend on ω k,l , we have that ∂ -→ Z + ∂ω k,l = ∂ -→ Z ∂ω k,l = -→ M k,l • | -→ ∆ * |. (4.15)
Using Corollary 1.8.1, we also get that

-→ Z + = r∈v ω r,r D∈ -→ Tr (k,l)∈E D (i,j)∈E D ω i,j = ω k,l r∈v ω r,r | -→ ∆ {k,l},{r,l} | = ω k,l • L k,l (
∂ -→ Z + ∂ω k,l = L k,l .
(4.17)

Combining (4.15), (4.16) and (4.17) yields

-→ Z + = ω k,l • -→ M k,l • | -→ ∆ * |, P ({k, l} ∈ E D |y) = ω k,l • -→ M k,l . (4.18)
The complexity of the computation is dominated by the complexity of the inversion of -→ ∆ * required to compute -→ M , and all posterior edge probabilities can therefore be obtained in O(p 3 ) time. In Chapter 3, we computed, among other things, the posterior probability of a change-point occurring at a given time t. Assuming that the number of segments is equal to K, we are actually evaluating the posterior probability associated to the subset of M K made of the segmentations having a change-point at time t. This subset was denoted by B K (t) in Chapter 3. The posterior probability given to this subset obviously depends on the prior distribution that was set on M K . We remind that we limited our attention to the distributions on M K that factorise over the segments. Such distributions are given through a upper triangular matrix a with non-negative terms by

Prior

p(m) = 1 C K (a) r∈m a r , ∀m ∈ M K . (4.19)
Under such prior distributions, the prior probability of B K (t) can be computed exactly as described in Section 3.4.1 by replacing A by a. A uniform prior on M K can be obtained by setting all terms above the diagonal in a to 1. But matrix a can be used to enforce additional constraints on the allowed segmentations. For instance, segments with length smaller than L can be forbidden by putting the terms in a = (a s,t ) such that t = s + l, 1 l L, to zero. Similarly, it is possible to enforce that all segments have an even length by settings all terms corresponding to odd segments to zero. In a, each superdiagonal corresponds to segments of a given length, so that designing a matrix allowing for a list of given lengths is straightforward. In the aforementioned examples, whenever all non-zero terms are set to one, a induces a uniform distribution on its support. But we are by no means limited to binary values. Segmentations with segments of similar lengths can for instance be favoured a priori by setting

a s,t = |t -s| d , 1 s < t N + 1, d > 0. (4.20)
The greater d is taken, the more we tip the scale in favour of homogeneous segmentations in terms of segment length. Similar results can be obtained when power functions are replaced with any strictly increasing non-negative function on R + . On the contrary, heterogeneous segmentations are favoured when a s,t is a decreasing function of the length of segment s; t .

Obviously it is perfectly conceivable to combine (4.20) with hard constraints to obtain complex prior distributions. Figure 4.5 illustrates the different examples that we have given on the segmentations of 1; 20 .

A Different Kind of Prior

Now imagine that we are given a segmentation m, either from previous analyses or from expert elicitation, in which we believe to some extent. Our wish would be to specify a prior distribution on m that concentrates around this initial guess. This is possible to some extent with the prior distributions we have considered so far. We consider a different kind of prior for m in which it is straightforward, as well as the consequences of this change on the inference.

For where f (•|d + 1, . . . , d + 1) stands for the density of the Dirichlet distribution on

L K = ∈ [0; 1] K K k=1 k = 1
with parameters (d + 1, . . . , d + 1). In a general Dirichlet distribution, the parameters are not assumed to be identical. In our context, that would mean that the position of a segment within a segmentation influences its contribution to the probability of the segmentation, with a distribution of the form

p(m) ∝ K k=1 k (m) d k . (4.21) If (d 1 , . . . , d K ) = (d 1 + 1, . . . , d K + 1), the distribution of ( 1 (m), . . . , K (m)) is basically a discretised version of D(d 1 , . . . , d K ).
Indeed, it can be seen as a distribution on L K whose density f N is proportional to

f N ( 1 , . . . , K |d 1 , . . . , d K ) ∝ f N 1 N , . . . , N K N d 1 , . . . , d K , ∀ ∈ L K . (4.22) Now, if ∼ D(d 1 , . . . , d K )
, the mean and mode of are respectively given by

E[ k ] = d k K k=1 d k , max k = d k K k=1 d k , ∀k ∈ 1; K .
Therefore, if we are given a prior guess m on the segmentation, the vector d can be taken proportional to ( m), so that the Dirichlet distribution underpinning the distribution of m has an expected value corresponding to m.

Using such prior distributions for m in the general model described in Chapter 3 requires some adjustments. We now have to consider one weighted segment likelihood matrix A per position within the segmentation. The matrix linked with the k-th segment is denoted by A (k) , and its general term is given by

A (k) s,t = a (k) s,t • p(y s;t ) if s < t, 0 otherwise, a (k) s,t = |t -s| d k if s < t, 0 otherwise, ∀k ∈ 1; K .
For 1 k 1 k 2 K, we let a (k1:k2) . . = k2 k=k1 a (k) and A (k1:k2) . . = k2 k=k1 A (k) . The latter are meant to replace the powers of matrix A involved in the different results of Section 3.4. The normalising constant in (4.22) is equal to [a (1:K) ] 1,N +1 and the marginal likelihood of the observations is therefore given by

p(y|K) = [A (1:K) ] 1,N +1 [a (1:K) ] 1,N +1 .
Similarly, the posterior probabilities associated to are now given by

B K,k (t) = {m ∈ M K |τ k = t}, B K (t) = K k=1 B K,k (t), S K,k ( s; t ) = {m ∈ M K |r k = s; t }, S K ( s; t ) = K k=1 S K,k ( s; t ), d = (3D, 7D) d = (3D, 4D ,
B K,k (t) = A (1:k) 1,t A (k+1:K) t,N +1 A (1:K) 1,N +1 , B K (t) = K-1 k=1 B K,k (t), S K,k ( s; t ) = A (1:k-1) 1,s A (k) s,t A (k+1:K-) t,N +1 A (1:K) 1,N +1 , S K ( s; t ) = K k=1 S K,k ( s; t ).
As these probabilities can all be computed from the first line of matrices {A (1:k) } K k=1 and the last column of matrices {A (k:K) } K k=1 , we retain a complexity of O(KN 2 ). Prior probabilities are likewise obtained through matrices {a (k1:k2) } 1 k1 k2 K (see Figure 4.6 for some examples).

Transferring Knowledge from One Dataset to Another

We begin this section with a practical example. Genome annotation is a segmentation problem in which one tries to delimit the regions of a DNA sequence that are transcribed to produce RNA, also called genes, from intergenic regions. Different technologies are currently available to measure gene expression along a genome. DNA microarrays have been used for almost twenty years. Roughly speaking, these microarrays consist of a collection of short single-strand DNA sequences, or probes, fixed and arranged on a solid support. Each probe is represented a great number of time. The RNA sequences contained in a biological sample of interest are retro-transcribed into single-strand DNA sequences and marked with a fluorescent compound. These DNA sequences are put on the chip and a property of DNA called hybridisation brings them together with probe sequences if their nucleotidic sequences are matching. After washing off unhybridised sequences, the fluorescence associated to each group of probes is measured and used as a proxy for gene expression level at the locus of the probe. The number of probes is usually small with respect to the length of the sequence of interest. This technology therefore produces low-resolution continuous measurements of gene expression.

More recent technologies are based on high-throughput sequencing, also called nextgeneration sequencing or NGS. They are referred to RNA-seq methods. In such methods, the RNA sequences of a biological sample are cut into small pieces and sequenced. After aligning all these sequences on a reference genome, we get a profile indicating the number of times each nucleotide has been seen in an RNA sequence, so that gene expression is given by count data, at a much higher resolution than microarray data.

Although RNA-seq methods are progressively overshadowing microarrays, a lot of data are available from times preceding the advent of NGS. Imagine a situation in which we have have both microarray and RNA-seq data. Ideally, the microarray data could be used to provide some kind of prior information for the analysis of the finer RNA-seq data. This is actually possible in the framework that we have described.

Let y 1 , y 2 be realisations of two random processes Y 1 , Y 2 with N 1 , N 2 time-points respectively. Each random process is modelled with a segmentation model, for which we assume that the number of segments is known to be K:

p(m 1 |d 1 ) ∝ K k=1 k (m) d 1 k , p(m 2 |d 2 ) ∝ K k=1 k (m) d 2 k , p(θ 1 |m 1 , α 1 ) = r∈m 1 p(θ 1 r |α 1 ), p(θ 2 |m 2 , α 2 ) = r∈m 2 p(θ 2 r |α 2 ), p(y 1 |m 1 , θ 1 ) = r∈m 1 p(y 1,r ), p(y 2 |m 2 , θ 2 ) = r∈m 2 p(y 2,r ).
The marginal likelihood of y 1 is given by

p(y 1 |d 1 , α 1 ) ∝ m∈M K ( 1;N1 ) K k=1 k (m) d 1 k p(y 1,r |α 1 ).
Let d1 be the value of d 1 maximising this expression. We are therefore dealing with the first dataset in a frequentist paradigm. Taking d 2 = d1 , we then use what we have learned on the distribution of m 1 to specify a prior distribution on m 2 and perform classical Bayesian inference on the second dataset. A schematic representation of this procedure is given in Figure 4.7. The key element is that a given Dirichlet distribution can be discretised at any level N ∈ N to obtain a distribution on M K ( 1; N ). This procedure is thus allowed by the specific form that we have chosen for the prior distribution of m. It would not be possible with a prior of the form p(m) = r∈m a r .

α 1 d 1 m 1 Y 1 θ 1 θ 2 Y 2 α 2 Frequentist inference Bayesian inference } d1 d 2 m 2
If we were to perform both steps of the procedure on the same dataset, we would be considering what is called an empirical Bayes approach. Such methods are straying from standard Bayesian inference, as information is drawn from the data before actually performing the inference. There is no proper theoretical basis behind such approaches, and all qualms raised by the facts that data are used twice are well founded.

We finish this section with a toy example on synthetic data in the context of gene annotation. We took datasets of respective sizes N 1 = 100 and N 2 = 500. In both cases, the true segmentation was chosen with segment lengths proportional to (0.3, 0.2, 0.1, 0.4). Datasets were simulated as follows: µ = (20, 10, 20, 10), σ = 2.5, λ = (5, 4, 5, 4),

Y 1,t ∼ N (µ k , σ 2 ), ∀t ∈ r k , Y 2,t ∼ P(λ k ), ∀t ∈ r k .

En Bref

This chapter, written in French, is a standalone and substantial abstract of the remainder of the dissertation.

Que l'on parle de réseaux de régulation de gènes, de réseaux informatiques ou encore de réseaux sociaux, il est indéniable que le concept de réseau est central dans de nombreux domaines. On peut définir un réseau comme un système d'entités inter-connectées. Dans les problèmes qui nous intéressent, le réseau n'est pas connu. Il ne s'agit donc pas d'étudier les propriétés d'un réseau que l'on observe, mais plutôt de travailler sur un réseau latent. Dans un premier temps, nous nous intéressons à ce qui est communément appelé l'inférence de réseaux, dans lequel le but est de reconstruire ce réseau latent à partir d'observations individuelles sur chacune des entités impliquées. Lorsque ces observations sont organisées en série temporelle, il est souvent peu réaliste de supposer que le réseau est stationnaire. Nous concentrons notre attention sur le cas où le réseau subit un certain nombre de brusques changements au cours du temps. La détermination de ces points de rupture est un problème de segmentation.

Le formalisme naturel dans lequel traiter les problèmes impliquant des réseaux est celui des modèles graphiques. Un graphe est un ensemble de sommets reliés par des arêtes, qui peuvent être dirigées ou non. Dans un modèle graphique, les relations de dépendance et d'indépendance conditionnelles vérifiées par une distribution sont représentées à l'aide d'un graphe. Ces propriétés abstraites deviennent ainsi plus facilement appréhendables. Dans le cadre des modèles graphiques, l'inférence de réseaux consiste donc en l'inférence du graphe décrivant le modèle graphique. Le graphe est alors lui-même un paramètre du modèle général.

Nous avons décidé de considérer les problèmes d'inférence de réseaux et de segmentation d'un point de vue bayésien. Il est donc nécessaire de définir des distributions a priori pour les différents paramètres impliqués, qu'ils soient à valeurs discrètes, tels que le graphe ou la segmentation, ou à valeurs continues. En choisissant ces distributions avec soin, il est possible d'effectuer une inférence bayésienne exacte dans le cadre des modèles graphiques en un temps polynomial par rapport au nombre de variables. Cette complexité est en particulier permise par une restriction sur l'espace des graphes explorés. En effet, une sommation efficace sur l'ensemble des graphes non-dirigés, connectés et acycliques, aussi appelés arbres couvrants, est permise par un résultat algébrique connu sous le nom de théorème arbre-matrice. Il s'agit d'une voie différente des méthodes basées sur l'échantillonnage, telles que les approches de type Markov Chain Monte Carlo (MCMC), qui sont couramment utilisées dans le contexte de l'inférence de réseaux. L'échantillonnage au sein de l'ensemble des graphes est dans notre cas remplacée par une exploration exhaustive d'un sous-ensemble de graphes. Le premier chapitre fournit le contexte dans lequel s'inscrivent les chapitres suivants. La première partie est consacrée aux modèles graphiques, plus particulièrement aux modèles graphiques non-dirigés. Elle se fonde principalement sur les travaux de [START_REF] Dawid | Hyper Markov Laws in the Statistical Analysis of Decomposable Graphical Models[END_REF], [START_REF] Lauritzen | Graphical Models[END_REF][START_REF] Byrne | Hyper and Structural Markov Laws for Graphical Models[END_REF]. Les modèles dirigés sont évoqués bien que n'étant pas au centre des travaux présentés ici. La seconde présente deux résultats d'algèbre utiles pour la suite.

Modèles Graphiques

Formellement, un graphe non-dirigé est donné par un ensemble de sommets V , que l'on prend ici égal à {1, . . . , p}, et par un ensemble d'arêtes E contenu dans l'ensemble des parties de V de taille 2, dénoté P 2 (V ). On note G = (V, E). Pour tout graphe G, on note également E G les arêtes de G. Un chemin γ entre deux sommets distincts α, β d'un graphe G est une suite α = γ 0 , . . . , γ n = β, n 1, telle que, pour tout

1 i n, {γ i-1 , γ i } ∈ E G .
Définition (Séparation). Soit G un graphe non-dirigé et (A, B, S) un triplet de sousensembles de V . L'ensemble S sépare A de B si tout chemin joignant un sommet de A à un sommet de B dans G intersecte S.

La notion de séparation est fondamentale pour établir le lien entre graphe et indépendance conditionnelle. Soit X = (X i ) i∈V un vecteur aléatoire à valeurs dans un espace produit X = i∈V X i . Pour tout sous-ensemble A de V , X A dénote (X i ) i∈A .

Propriété (Markov globale). Soit G un graphe. La distribution de X satisfait la propriété de Markov (globale) par rapport à G si, pour tout triplet (A, B, S) de sous-ensembles disjoints de V , on a

S sépare A de B ⇒ X A |= X B |X S .
Ainsi, lorsque qu'une distribution satisfait la propriété de Markov par rapport à un graphe G, ses propriétés d'indépendance conditionnelle sont représentées par le graphe. Il n'y a cependant qu'une implication entre séparation et indépendance conditionnelle, et toute distribution est par exemple Markov par rapport au graphe complet sur V . Pour tout graphe G, M(G) dénote l'ensemble des distributions sur X qui sont Markov par rapport à G. Définition (Modèle graphique). Un modèle graphique est un couple formé d'un graphe G, appelé structure du modèle, et d'une famille de distributions F telle que F ⊆ M(G).

Les notions de propriété de Markov et de modèle graphique peuvent être adaptées au cas des graphes dirigés acycliques, plus connus sous l'acronyme anglais de DAG.

L'inférence de la structure d'un modèle graphique dans un cadre bayésien suppose que la distribution des observations ainsi que le graphe donnant la structure du modèle soient euxmêmes des objets aléatoires. [START_REF] Dawid | Hyper Markov Laws in the Statistical Analysis of Decomposable Graphical Models[END_REF] et [START_REF] Byrne | Hyper and Structural Markov Laws for Graphical Models[END_REF] ont respectivement étendu la propriété de Markov aux distributions sur les distributions et sur les graphes. Ces extensions permettent de donner un cadre théorique complet à l'inférence de réseaux.

Algèbre et Algorithmique

L'ensemble des arbres couvrants sur V , dénoté T , est de cardinal p p-2 , bien inférieur au cardinal de l'ensemble des graphes non-dirigés qui est de 2 p(p-1)/2 , tout en restant conséquent. Cependant, sous une hypothèse de factorisation sur les arêtes, il est possible d'intégrer une fonction définie sur T en un temps polynomial.

Théorème (arbre-matrice). Soit ω = (ω i,j ) 1 i,j p une matrice symétrique à valeurs positives, de diagonale nulle. Soit ∆ la matrice laplacienne associée à ω, donnée par

∆ i,j = -ω i,j if i = j,
k∈V ω k,j if i = j.

Alors tous les cofacteurs de ∆ sont égaux à

Σ . . = T ∈T {i,j}∈E T ω i,j .
Ce résultat est attribué, sous différentes formes, à Kirchhoff (1847) ou [START_REF] Cayley | A theorem on trees[END_REF]. Le calcul de Σ est donc possible en O(p 3 ) en tant que le déterminant d'une matrice de taille p -1. Une version généralisée de ce résultat [START_REF] Chaiken | A Combinatorial Proof of the All Minors Matrix Tree Theorem[END_REF] permet de relâcher la contrainte de connexité sur les graphes considérés, et ainsi de passer des arbres aux forêts. Le théorème arbre-matrice et sa version généralisée sont utilisés au Chapitre 2 dans le calcul des probabilités d'apparition d'arêtes a posteriori. Nous fournissons une nouvelle preuve d'un résultat du à [START_REF] Kirshner | Learning with Tree-Averaged Densities and Distributions[END_REF].

Le second résultat énoncé dans cette partie concerne les partitions de 1; N , N 1, en sous-ensembles de la forme i; j , appelés segments. De telles partitions sont appelées segmentations. L'ensemble des segmentations d'un segment i; j en un nombre K de soussegments est noté M K ( i; j ). Le résultat suivant permet d'intégrer une fonction sur les segmentations factorisant sur les segments.

Théorème [START_REF] Rigaill | Exact posterior distributions and model selection criteria for multiple change-point detection problems[END_REF]. Soit f une fonction définie par f (m) = s;t ∈m A s,t , ∀1 K N, ∀ 1 i < j N + 1, ∀m ∈ M K ( i; j ), où A est une matrice triangulaire supérieure à valeurs positives, de diagonale nulle et de taille N + 1. Alors

F K i,j . . = m∈M K ( i;j ) f (m) = [A K ] i,j .
À partir de ce résultat, il est possible de calculer de nombreuses quantités dans les modèles de détection de rupture. Ce chapitre s'intéresse à l'inférence de réseaux. Comme annoncé dans l'introduction, l'inférence est réalisée dans un cadre bayésien. L'espace des graphes explorés est restreint aux arbres couvrants afin de permettre une inférence exacte. [START_REF] Chow | Approximating Discrete Probability Distributions with Dependence Trees[END_REF] sont parmi les premiers à avoir étudié l'inférence de distributions dont la structure de dépendance est arborescente. Depuis, de nombreux travaux se sont penchés sur le sujet. On citera notamment ceux de Meilă & Jordan (2001), [START_REF] Meilă | Tractable bayesian learning of tree belief networks[END_REF], [START_REF] Kirshner | Learning with Tree-Averaged Densities and Distributions[END_REF][START_REF] Lin | Learning Sparse Markov Network Structure via Ensemble-of-Trees Models[END_REF]. Cependant, la connexité et l'absence de cycles sont des contraintes fortes inhérentes aux arbres couvrants, qui ont le défaut d'être peu réalistes dans de nombreux cas d'application. Localement, de nombreux réseaux ont cependant une structure arborescente, tout en ne satisfaisant pas ces contraintes globales. C'est pourquoi, au lieu de chercher à apprendre le meilleur arbre, nous nous sommes intéressés à la probabilité d'apparition de structures locales, telles que les arêtes, dans un arbre aléatoire. En moyennant sur l'espace des arbres couvrants, il est ainsi possible d'effacer en partie les contraintes individuelles de chaque arbre.

La contribution majeure de ce chapitre est la construction d'un cadre bayésien complet pour l'inférence de modèles graphiques à structure d'arbre. Ce cadre est fondé sur les travaux de [START_REF] Dawid | Hyper Markov Laws in the Statistical Analysis of Decomposable Graphical Models[END_REF] généralisant la propriété de Markov aux distributions sur des espaces de distributions elles-mêmes Markov par rapport à un graphe. De telles distributions sont appelées hyperdistributions pour plus de clarté. La propriété hyper Makov forte sur les hyperdistributions permet de définir une distribution a priori sur la distribution des observations (ou sur la loi des paramètres de la distribution dans un cadre paramétrique) respectant la structure du modèle graphique. Lorsque ce choix de distribution est combiné à une distribution a priori sur les arbres couvrants qui factorise sur les arêtes, il est en fait possible d'intégrer à la fois sur la distribution des observations et sur la structure de manière exacte et efficace.

Nous fournissons également une nouvelle preuve d'un résultat énoncé par Kirshner (2007) sur le calcul des probabilités d'apparition d'arêtes a posteriori. Un calcul direct à partir du théorème arbre-matrice conduit à une complexité en O(p 5 ) pour le calcul de ces probabilités, mais on montre ici qu'il est en fait possible de se ramener à une complexité en O(p 3 ).

Une étude simulatoire permet de quantifier empiriquement l'impact de l'hypothèse arborescente sur l'inférence de réseaux dans le cas où le réseau à récupérer n'est pas un arbre. Nous observons des performances similaires à celles obtenues par une inférence sur l'espace des graphes dirigés acycliques. Une application sur des données de cytométrie en flux est également présentée, pour laquelle on remarque un léger avantage des graphes dirigés acycliques sur les arbres.

Un paquet R implémentant l'approche développée dans ce chapitre est disponible sur le CRAN sous le nom de saturnin.

Modèle

Reprenant les notations du chapitre précédent, X = (X i ) i∈V est un vecteur aléatoire à valeurs dans un espace produit X = i∈V X i . Pour tout arbre T ∈ T , on considère le modèle graphique m T = (T, M(T )), où l'on rappelle que M(T ) désigne l'ensemble des distributions sur X qui sont Markov par rapport à T . Le modèle global consiste à tirer un arbre T dans T selon une certaine loi ξ. Une distribution π est ensuite tirée dans le modèle graphique m T selon une hyperdistribution ρ T définie sur M(T ). Enfin, X est tiré selon π.

Une représentation graphique de ce modèle est donnée ci-dessous. Une spécification complète de ce modèle nécessite de donner une distribution a priori sur les arbres ξ et une collection d'hyperdistribution {ρ T } T ∈T comme distribution a priori de π|T . Cette collection peut en fait être définie, sous certaines conditions, en donnant les hyperdistributions marginales sur la loi de (X i , X j ) pour tout {i, j} ∈ P 2 (V ).

Distributions a priori sur T et π [START_REF] Dawid | Hyper Markov Laws in the Statistical Analysis of Decomposable Graphical Models[END_REF].

Pour tout couple (A, B) de sous-ensembles de V , on note π A la distribution marginale obtenue à partir de π sur les variables X A et π B|A la collection de distributions conditionnelles de X B |X A sous π. On note également ρ T A l'hyperdistribution induite par ρ T sur π A et ρ T B|A la collection d'hyperdistributions induites par ρ T sur π B|A . Définition. ρ T est dite fortement hyper Markov par rapport à T si, pour tous A, B ⊂ V tels que A∩B est complet et sépare A de B, {π B|A , π A|B , π A∩B } sont mutuellement indépendants sous ρ T .

De telles hyperdistributions ont des propriétés intéressantes vis-à-vis de p(•|T ).

Proposition 1 [START_REF] Dawid | Hyper Markov Laws in the Statistical Analysis of Decomposable Graphical Models[END_REF]. (5.2)

Il faut maintenant choisir toutes les hyperdistributions ρ T de telle manière qu'elles soient hyper Markov par rapport à leur arbre respectif et que les intégrales locales utilisées en (5.2) ne dépendent pas de T .

Proposition 5.1 [START_REF] Dawid | Hyper Markov Laws in the Statistical Analysis of Decomposable Graphical Models[END_REF]. Soit ρ une hyperdistribution générale sur F telle que, pour tout A ⊂ V , π A |= π V \A|A sous ρ. Alors, pour tout arbre T ∈ T , il existe une unique hyperdistribution ρ T sur F T fortement hyper Markov par rapport à T et telle que, pour toute arête {i, j} ∈ E T , ρ T ij = ρ ij .

(5.3)

{ρ T } T ∈T est une famille compatible d'hyperdistributions.

On requiert en fait que l'hyperdistribution ρ servant de base à la famille compatible soit fortement hyper Markov par rapport au graphe complet. Si {ρ T } T ∈T est construite ainsi, la Proposition 5.1 garantit que ρ T est fortement hyper Markov par rapport à T pour tout T ∈ T et que les vraisemblances locales intégrées p(x i , x j |T ) et p(x i |T ) ne dépendent en fait pas de T . Ces termes peuvent être calculés une fois pour toutes et utilisés pour tous les arbres. Ainsi, on a ξ(T |x) = 1 Z(ω) {i,j}∈E T ω i,j , ∀T ∈ T, ω i,j = β i,j p(x i , x j ) p(x i )p(x j ) , ∀{i, j} ∈ P 2 (V ).

Dans la pratique, les intégrales locales permettant d'obtenir ω peuvent être calculées simplement dans un cadre conjugué, mais peuvent aussi être obtenues par des méthodes d'intégration numérique ou stochastique. Trois cas particuliers sont détaillés, à savoir le cas où X suit une distribution multinomiale, une distribution gaussienne ou une copule. Une utilisation directe du théorème arbre-matrice fournit Z(ω). La seconde somme peut être calculée de manière similaire en appliquant le théorème à une matrice de poids ω (kl) obtenue à partir de ω en mettant le poids ω k,l = ω l,k à zéro. Effectué de cette manière, le calcul de {P ({k, l} ∈ E T |x)} {k,l}∈P2(V ) est de complexité O(p 5 ). Ces probabilités peuvent en fait être calculées en temps cubique.

Probabilité d'apparition d'arête

Théorème [START_REF] Kirshner | Learning with Tree-Averaged Densities and Distributions[END_REF]. Soit u un sommet de V . On note ∆ {u} la matrice obtenue à partir de la matrice laplacienne ∆ associée à ω en enlevant la ligne et le colonne correspondant à u et Q = (∆ {u} ) -1 . Soit M la matrice dont le terme général est donné par

M k,l =    (Q k,k + Q l,l -2Q k,l ) if k = u, l = u, Q k,k if l = u, Q l,l if k = u.
Alors, pour toute arête {k, l} ∈ P 2 (V ), P ({k, l} ∈ E T |x) = ω k,l • M k,l .

Ainsi, le calcul de {P ({k, l} ∈ E T |x)} {k,l}∈P2(V ) revient au calcul de l'inverse d'une matrice de taille p -1, de complexité O(p 3 ). La question à laquelle se propose de répondre ce chapitre est celle de la détection de brusques changements, aussi appelés ruptures, dans la structure de dépendance d'une série temporelle multivariée. On parle aussi de problème de segmentation. En se fondant sur le cadre développé au chapitre précédent, le modèle graphique décrivant la structure de dépendance à un instant donné est supposé arborescent. En tirant partie de cette hypothèse et en utilisant le résultat de Rigaill et al. (2012)

Modèle

On suppose que les observations {y t } N t=1 sont une réalisation d'un processus aléatoire multivarié {Y t } N t=1 de dimension p 2. Pour 1 t N , Y t = (Y t 1 , ..., Y t p ) est un vecteur aléatoire à valeurs dans un espace produit X = p i=1 X i . Pour tout intervalle temporel r ⊆ 1; N , Y r . . = {Y t } t∈r dénote le processus restreint à r. Le modèle suppose qu'il existe une partition m de 1; N en intervalles r telle que, sur chaque intervalle r de m, Y r est indépendant et identiquement distribué selon un modèle graphique à structure d'arbre. Ce modèle appartient à la famille des modèles à partition décrite par [START_REF] Barry | Product Partition Models for Change Point Problems[END_REF]. Si m possède K segments r 1 , . . . , r K , on note respectivement T = {T k } K k=1 et θ = {θ r } r∈m les arbres et paramètres donnant le modèle graphique de chaque segment. Pour r ∈ m, κ(r|m) indique la position du segment r dans la segmentation m. Les observations {Y t } t∈r sont indépendantes, identiquement distribuées et suivent un modèle graphique de structure T κ(r|m) et de paramètres θ r . La structure et les paramètres relatifs à chaque segment sont également indépendants et identiquement distribués, selon un modèle de la forme décrite dans le Chapitre 2. La distribution de chacun des arbres de T est donnée à travers une matrice de poids d'arêtes b. La distribution a priori sur m est choisie comme factorisant sur les segments, conformément à l'hypothèse de factorisation requise par [START_REF] Rigaill | Exact posterior distributions and model selection criteria for multiple change-point detection problems[END_REF]. Elle est donnée par une matrice triangulaire supérieure a de taille N + 1, dans laquelle a s,t donne le poids du segment s; t . Une description complète du modèle est donnée en Figure 5.2. de covariables peut être effectuée à peu de frais. Nous montrons par exemple que le cadre de la régression linéaire multiple multivariée permet de conserver les propriétés de Markov nécessaires au bon fonctionnement des résultats algébriques utilisés dans l'inférence.

L'exemple ayant motivé notre intérêt pour la question des covariables concernait des données de comptage. Le passage au modèle linéaire généralisé n'est dans ce cas pas sans frais, et empêche de poursuivre dans la même direction. C'est pourquoi nous proposons une approche pragmatique fondée sur les copules. L'abondance de chacune des espèces est modélisée individuellement par un modèle linéaire généralisé utilisant les covariables comme variables explicatives. Les paramètres de ces modèles sont ajustés de manière classique. Les résidus de Pearson correspondant à chaque espèce sont calculés et ramenés à l'intervalle [0; 1] grâce à leur fonction de répartition empirique. L'inférence du réseau sous-jacent est ensuite effectuée sur ces fractiles en utilisant des distributions définies sur [0; 1] p dont les marginales sont uniformes, aussi appelées copules. Cette approche, bien qu'ayant le mérite d'être entièrement générique, n'est bien évidemment qu'un pis-aller, dans la mesure où l'incertitude sur les paramètres de régression n'est pas du tout prise en compte.

Dans la seconde section, nous expliquons comment introduire de la dépendance temporelle dans le modèle de segmentation décrit au Chapitre 3. En effet, l'hypothèse d'indépendance faite dans ce modèle est souvent peu vraisemblable en pratique. En ce qui concerne l'intégration sur l'ensemble des segmentations, la seule indépendance nécessaire est celle des observations entre les segments, afin que la loi a posteriori sur m factorise sur les segments. De ce point de vue là, rien n'empêche d'introduire de la dépendance temporelle au sein des segments. Notre suggestion est d'utiliser les modèles d'indépendance temporelle (abrégé par TIM, d'après l'acronyme anglais correspondant) introduit par [START_REF] Siracusa | Tractable Bayesian Inference of Time-Series Dependence Structure[END_REF] afin de gérer la dépendance temporelle intra-segment. Ce sont des modèles graphiques dirigés, dans lesquels la dépendance d'un instant au suivant est décrite par un graphe dirigé. Si pa(i) désigne l'ensemble des parents du sommet i dans le graphe décrivant la dépendance, alors Y t i dépend de Y t-l pa(i) , . . . , Y t-1 pa(i) , où l est la latence du modèle. Un exemple est donné Figure 5.3. Si l'on suppose que la structure de dépendance est un arbre, il est en fait possible d'utiliser une version dirigée du théorème arbre-matrice pour effectuer l'inférence. Les hypothèse d'indépendance sur les paramètres doivent également être adaptées à cette nouvelle configuration.

Enfin, la dernière partie se penche un peu plus sur la question de la distribution a priori sur la segmentation m du modèle donné Figure 5.2. Une liste non-exhaustive de contraintes pouvant être encodées dans la matrice a des poids de segments est donnée. Il est par exemple possible d'interdire les segments de longueur inférieure à une certaine taille, ou encore de Titre : Inférence bayésienne exacte dans les modèles graphiques : inférence de réseaux à structure arborescente et segmentation Title: Exact Bayesian Inference in Graphical Models: Tree-structured Network Inference and Segmentation

Y t-1 2 Y t-1 1 Y t-1 3 Y t 1 Y t 2 Y t 3 Y t+1 1 Y t+1 2 Y t+1 3 Y 2 Y 1 Y 3
Keywords: Bayesian inference, graphical models, Matrix-Tree theorem, network inference, segmentation, spanning tree.

Abstract: In this dissertation we investigate the problem of network inference. The statistical framework tailored to this task is that of graphical models, in which the (in)dependence relationships satisfied by a multivariate distribution are represented through a graph.

We consider the problem from a Bayesian perspective and focus on a subset of graphs making structure inference possible in an exact and efficient manner, namely spanning trees. Indeed, the integration of a function defined on spanning trees can be performed with cubic complexity with respect to number of variables under some factorisation assumption on the edges, in spite of the super-exponential cardinality of this set. A careful choice of prior distributions on both graphs and distribution parameters allows to use this result for network inference in tree-structured graphical models, for which we provide a complete and formal framework. We also consider the situation in which observations are organised in a multivariate time-series. We assume that the underlying graph describing the dependence structure of the distribution is affected by an unknown number of abrupt changes throughout time. Our goal is then to retrieve the number and locations of these changepoints, therefore dealing with a segmentation problem.

Using spanning trees and assuming that segments are independent from one another, we show that this can be achieved with polynomial complexity with respect to both the number of variables and the length of the series.

Université Paris-Saclay Espace Technologique / Immeuble Discovery Route de l'Orme aux Merisiers RD 128 / 91190 Saint-Aubin, France

"

  Nine worlds I knew, the nine in the tree With mighty roots beneath the mold." -The poetic Edda iii 0 Introduction Foreword . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 -

 1 Figure 1 -Log-cardinality of the set of undirected graphs and spanning trees as a function of the number of vertices.

  Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Graphical Models & Markov Equivalence 24 1.2 Algebra & Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 1

1. 1

 1 Graphical Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 . 1 -

 11 Figure 1.1 -An example of tree.

Figure 1 . 2 -

 12 Figure 1.2 -Hierarchical model for network inference.

c

  does not d-separate a and b. f d-separates a and b. {c, f } d-separates a and b. For a DAG D, we let A |= D B|C denote the statement "C d-separates A from B in D". Just like |= G for undirected graphs G, the ternary relation |= D is a graphoid. It can therefore be used to represent conditional independence. Definition 1.22 (Directed global Markov property). Let D = (V, E) be a DAG. A probability distribution P on X is said to obey the (directed) global Markov property relative to D, if for any triple (A, B, S) of disjoint subsets of V such that S d-separates A from B in D, it holds that A |= B|S.

2. 1

 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  & Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2 . 1 -

 21 Figure 2.1 -Global hierarchical model.

  Tree Structures & Distributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  As Bayes' rule states that ξ(T |D) ∝ ξ(T )p(D|T ), we are now interested in the marginal likelihood of the data under a tree model m T , p(D|T ) = F T p(D|π)p(π|T )dπ.

  Figure 2.2 -Compatible strong hyper Markov tree model.

  Tree Graphical Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2 . 3 -

 23 Figure 2.3 -Model with variable kl explicitly controlling the status of edge {k, l} in T * .

2. 5

 5 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2 . 4 -

 24 Figure 2.4 -Gold standard networks in the simulation study.

Figure 2 . 5 -

 25 Figure 2.5 -ROC curves for the posterior edge probabilities and (TPR, FPR) scores for the MAP estimate on data sets of size 25, 50 and 200 (from top to bottom). For the ROC curves, the mean curve is plotted in bold line. The color of a (TPR, FPR) point expresses its frequency within the 100 samples.

Figure 2 . 6 -

 26 Figure 2.6 -Posterior probability of the MAP tree, ratio to the posterior probability of the second best tree and entropy of the posterior tree distribution (normalised by the entropy of the uniform distribution on T , i.e. (p -2) log(p)).

Figure 2 . 7 -

 27 Figure2.7 -Area under the ROC (top) & PR (bottom) curves computed for the output of our approach and of the MCMC sampling algorithm in the DAGs on the multinomial samples and for the output of our approach in the Gaussian setting on the undiscretized data.

  This section presents an application of our approach to flow cytometry data. They have been collected by[START_REF] Sachs | Causal protein-signaling networks derived from multiparameter single-cell data[END_REF] and were used by[START_REF] Werhli | Comparative evaluation of reverse engineering gene regulatory networks with relevance networks, graphical gaussian models and bayesian networks[END_REF] in a

  Figure 2.8 -Raf pathway.

  Most likely (left) and second most likely (right) trees in the posterior distribution on trees. Posterior probabilities for the edges in the tree model (with change of prior probability to q

Figure 2 . 9 -

 29 Figure 2.9 -Graphical representation of the results obtained on one of the five data sets. The edges of the golden standard network are colored in blue.

3. 1

 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3 . 1 -

 31 Figure 3.1 -Illustration of the change-point detection problem in the tree structure of a graphical model.

3. 2

 2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  & Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3 . 2 -

 32 Figure 3.2 -Global graphical model.

  Status & Structure Comparisons . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3 . 3 -

 33 Figure 3.3 -Model for edge status comparison.

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3 . 4 -

 34 Figure 3.4 -Posterior probability of observing a change-point for the tree-structured model (blue) and for the full model (red). The curve represents the mean value obtained from the 100 samples and the ribbon gives the standard deviation.

Figure 3 . 5 -

 35 Figure 3.5 -Boxplot of K = arg max K p(K|y) and K( m) = K(arg max m p(m|y)) against sample size N for the full model (Full) and the tree-structured model (Tree).

Figure 3 . 6 -

 36 Figure 3.6 -Area under the ROC curve computed for the posterior edge probability matrix p K ij (t) p i,j=1

Figure 3 . 7 -

 37 Figure3.7 -Boxplot of the posterior probability for an edge to have the same status throughout the time-series. Edges were separated according to their true status (either identical in all graphs or not). Each boxplot aggregates the results for all edges with a given status and all datasets.

Figure 3 . 8 -

 38 Figure 3.8 -Posterior probability of a change-point occurring at time t as a function of time integrated on K (left) and posterior distribution for K (right) for the full (Full) and tree-structured (Tree) models.

  Posterior change-point probability for five participants with the tree-structured model. Mean and standard deviation of posterior change-point probability across participants with the tree-structured model. Posterior change-points probability (left) and posterior distribution on K (right) when participants are jointly considered, with likelihood tempered at a level α for the full (Full) and tree-structured (Tree) models.

Figure 3 .

 3 Figure3.10 -Change-point location for the fMRI data. During the dark red interval, the subject of the speech was revealed to participants, who prepared their speech during the light red interval. This preparation is ended by a statement saying that they would not have to give the speech.

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.1.1 Multivariate Multiple Linear Regression 83 4.1.2 Copulas: a Pragmatical Approach 88 4.2 Temporal Dependence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 4.2.1 Temporal Interaction Models 92 4.2.2 Bayesian Inference of Tree-structured TIM 93 4.3 Prior Distribution on Segmentations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 4.3.1 A Different Kind of Prior 97 4.3.2 Transferring Knowledge from One Dataset to Another 100

4. 1

 1 Covariates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  If we let P(•|λ) denote the cdf of the Poisson distribution with parameter λ, the distribution of Y can be written as

4. 2

 2 Temporal Dependence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 4 . 1 -

 41 Figure 4.1 -Posterior edge probability matrix.

Figure 4 . 2 -

 42 Figure 4.2 -Interactions of Erysiphe alphitoides (Ea) with other OTUs on the susceptible oak tree. The names of the fungal OTUs that could be assigned to species level are indicated.

Figure 4 . 3 -

 43 Figure 4.3 -Graphical representation of a TIM(1) and TIM(2) model with tree dependence structure. The temporal and dependence structures are respectively depicted in blue and red.

Figure 4 . 4 -

 44 Figure 4.4 -Hierarchical model for structure inference in TIM.

  {k,l},{r,l} | stands for the absolute value of the determinant of the matrix obtained from -→ ∆ by removing rows k and l and columns r and l. As | -→ ∆ {k,l},{r,l} | , r ∈ V , does not depend on ω k,l , neither does L k,l , and

  Distribution on Segmentations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 4 . 5 -

 45 Figure 4.5 -Prior probability of observing a change-point as a function of time, for different choices of prior matrices a (represented above each set of curves) and for a number of segments K ∈ 2; 5 . Probabilities are normalised by the number of change-points of the corresponding segmentations, i.e. K -1.

DFigure 4 . 6 -

 46 Figure 4.6 -Prior probability of observing a change-point as a function of time, for different values of d.Probabilities are normalised by |d| -1.

Figure 4 . 7 -Figure 4 . 8 -

 4748 Figure 4.7 -Schematic representation of the knowledge transfer procedure.

Chapitre 1 .

 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Chapitre 2 .

 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  La distribution a priori sur T est prise de la formeξ(T ) = 1 Z(β) {i,j}∈E T β ij , ∀T ∈ T , (5.1) avec β une matrice symétrique à valeurs positives, de diagonale nulle. En prenant ξ sous cette forme, la constante de normalisation Z(β) = T ∈T {i,j}∈E T β ij peut être calculée grâce au théorème arbre-matrice. L'idée est maintenant de choisir {ρ T } T ∈T de telle sorte à conserver cette factorisation pour la distribution a posteriori ξ(•|x). Par la formule de Bayes, on obtient que ξ(T |x) ∝ ξ(T )p(x|T ), où p(x|T ) est la vraisemblance marginale des observations donnée par p(x|T ) = π(x)ρ T (π)dπ, ∀T ∈ T , ∀x ∈ X . Pour un arbre donné T , ρ T doit donc être choisie de telle manière que p(x|T ) reste factorisable sur les arêtes. Ce choix doit de plus être effectué de manière cohérente entre les différents arbres afin d'obtenir une écriture similaire à (5.1) pour ξ(•|x). Ceci est effectué en construisant une famille compatible d'hyperdistributions fortement hyper Markov (Dawid

  Si ρ T est fortement hyper Markov par rapport à T , alors p(•|T ) est Markov par rapport à T . Cela signifie que p(x|T ) peut être factorisée sur les arêtes de T .p(x|T ) = i∈V p(x i |T ) {i,j}∈E T p(x i , x j |T ) p(x i |T )p(x j |T ), ∀x ∈ X .

  La probabilité d'apparition d'une arête {k, l} a posteriori est donnée parP ({k, l} ∈ E T |x) = T ∈T , E T {k,l} ξ(T |x) = 1 -1 Z(ω) T ∈T , E T {k,l} {i,j}∈E T ω i,j .

Figure 5 . 1 -

 51 Figure 5.1 -Illustration du problème de détection de ruptures dans la structure d'un modèle graphique.

Figure 5 . 3 -

 53 Figure 5.3 -Représentation graphique d'un modèle d'indépendance temporelle de latence 1 à structure de dépendance arborescente.

  The ternary relation |= G satisfies Properties (C1) to (C4) given for conditional independence in Proposition 1.2, when h is replaced by set inclusion. Property (C5) holds whenever A, B, C are disjoint subsets of V . Let us illustrate this by rewriting (C2) for |=

G :

Table 1 .1 -Handbook

 1 of Markov properties.

	25

  Two DAGs are Markov equivalent if and only if they have the same skeleton and the same immoralities. A DAG D and an undirected graph G are Markov equivalent if and only if D has no immoralities and its skeleton is G.

	Example 1.13. All these graphs are Markov equivalent:
	but not Markov equivalent to:

  Model for Bayesian Inference of Graphical Models Based on Trees 40 2.3 Priors on Tree Structures & Distributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 Inference in Tree Graphical Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Application to Cytometry Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	2.3.1 Prior on Tree Structures	41
	2.3.2 Prior on Tree Distributions	41
	2.3.3 Structural Markov Property and Structurally Meta Markov Families	43
	2.4 44
	2.4.1 Integration with Respect to π	45
	2.4.2 Integration with Respect to T	48
	2.4.3 Controlling prior edge probability	51
	2.5 52
	2.5.1 Simulation Scheme	52
	2.5.2 Results	53
	2.6 55
	2.6.1 Data	56
	2.6.2 Results	56

.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 2.2 Background & Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.2.1 Markov Properties & Graphical Models 39 2.2.2

  ) If ρ T is strong hyper Markov w.r.t. T , then the marginal likelihood p(D|T ) is Markov w.r.t. to T .

Table 2 . 1 -

 21 Average running times for different network sizes with our method (Tree) and the MCMC approach on DAGs (DAG MCMC) on data sets of size n = 100.

Table 2 .

 2

2 -Inference results on flow cytometry data. Area under the ROC and PR curves for different discretization levels (standard deviation).

  énoncé au Chapitre 1, nous montrons qu'il est possible de calculer un certain nombre de quantités d'intérêt vis-a-vis du problème posé, telles que la probabilité a posteriori d'observer un point de rupture à un temps donné, ou encore la loi a priori du nombre de segments, en un temps polynomial par rapport à la longueur de la série et au nombre de variables. Lorsque la segmentation est donnée, nous fournissons également un moyen d'évaluer la probabilité a posteriori qu'une arête (ou que le graphe dans son intégralité) conserve le même statut au cours du temps.

  Mots clefs : arbre couvrant, inférence bayésienne, inférence de réseaux, modèles graphiques, segmentation, théorème arbre-matrice.Résumé : Cette thèse porte sur l'inférence de réseaux. Le cadre statistique naturel à ce genre de problèmes est celui des modèles graphiques, dans lesquels les relations de dépendance et d'indépendance conditionnelles vérifiées par une distribution multivariée sont représentées à l'aide d'un graphe. Il s'agit alors d'apprendre la structure du modèle à partir d'observations portant sur les sommets. Nous considérons le problème d'un point de vue bayésien. Nous avons également décidé de nous concentrer sur un sous-ensemble de graphes permettant d'effectuer l'inférence de manière exacte et efficace, à savoir celui des arbres couvrants. Il est en effet possible d'intégrer une fonction définie sur les arbres couvrants en un temps cubique par rapport au nombre de variables à la condition que cette fonction factorise selon les arêtes, et ce malgré le cardinal super-exponentiel de cet ensemble. En choisissant les distributions a priori sur la structure et les paramètres du modèle de ma-nière appropriée, il est possible de tirer parti de ce résultat pour l'inférence de modèles graphiques arborescents. Nous proposons un cadre formel complet pour cette approche. Nous nous intéressons également au cas où les observations sont organisées en série temporelle. En faisant l'hypothèse que la structure du modèle graphique latent subit un certain nombre de brusques changements, le but est alors de retrouver le nombre et la position de ces points de rupture. Il s'agit donc d'un problème de segmentation. Sous certaines hypothèses de factorisation, l'exploration exhaustive de l'ensemble des segmentations est permise et, combinée aux résultats sur les arbres couvrants, permet d'obtenir, entre autres, la distribution a posteriori des points de ruptures en un temps polynomial à la fois par rapport au nombre de variables et à la longueur de la série.

For the inference, the prior distributions on µ, σ and λ were chosen as in[START_REF] Rigaill | Exact posterior distributions and model selection criteria for multiple change-point detection problems[END_REF]. On the first dataset, the prior distribution on segmentations was set to the uniform. On the second one, it was either set to the uniform or using d1 , which was estimated through numerical optimisation. Figure4.8 shows the posterior probability of observing a change-point for the second dataset on 50 repetitions and for both prior distribution on segmentations. As expected, on this very simple example, we notice that there is an interest in transferring knowledge from the first dataset to the second. By doing so, we are able to detect the smallest of the four segments that is otherwise overlooked.

Remerciements

Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

We thank Ivor Cribben (Alberta School of Business, Canada) for kindly providing the fMRI data. We also thank Sarah Ouadah (AgroParisTech, INRA, Paris, France) for fruitful discussions.

Change-point Location

For m ∈ M K , we let 1 = τ 0 < τ 1 < • • • < τ K = N denote the change-points of m and, for 1 k K, we let r k = τ k-1 ; τ k denote its k-th segment. In this section we are interested in computing the posterior probabilities of the following subsets of M K ,

Subsets B K (t) and S K ( s; t ) are respectively the set of segmentations having a change-point at time t and the set of segmentations containing segment s; t . We let B K,k (t), B K (t), S K,k ( s; t ) and S K ( s; t ) denote the respective posterior probabilities of these subsets. [START_REF] Rigaill | Exact posterior distributions and model selection criteria for multiple change-point detection problems[END_REF] showed that, with the convention that A 0 t1,t2 = 1 for all 1 t 1 < t 2 N + 1, these probabilities could be expressed as

{B K,k (t)} N t=1 provides the exact posterior distribution of the k-th change-point when m has K segments. Posterior segment probabilities {S K ( s; t )} 1 s<t N +1 will be useful in the following.

Once {B K (t)} K 2 is computed, the posterior probability B(t) of a change-point occurring at time t integrated on K is obtained as

The computation of the posterior distribution on K is addressed below.

Number of Segments

The posterior distribution on K can also be derived from Proposition 3.1.

Inférence

Dans ce modèle, la vraisemblance marginale des observations, conditionnellement à K, est donnée par

(5.4) L'intégration sur les paramètres discrets nécessite de sommer sur un ensemble de taille

. Ce calcul peut cependant être effectué en un temps polynomial en utilisant le théorème de sommation sur les segmentations énoncé au Chapitre 1 sur une matrice de terme général A s,t = a s,t • p(y s,t ) si s < t, 0 sinon, où p(y s,t ) est la vraisemblance du segment s, t , intégrée sur l'arbre et les paramètres.

Ces vraisemblances sont obtenues grâce au théorème arbre-matrice. Une matrice de poids d'arêtes a posteriori ω s,t est calculée par segment. On a alors

Le résultat de Rigaill et al. (2012) indique alors que la vraisemblance marginale donnée en (5.4) est donnée par

. Cette vraisemblance marginale permet en particulier de calculer la loi a posteriori du nombre de segments K, puisque p(K|y) ∝ p(K)p(y|K). Il est ainsi possible d'obtenir un estimateur du maximum a posteriori pour K.

D'autres quantités peuvent être calculées à partir des puissances de la matrice A. La probabilité a posteriori B K,k (t) que le k-ème des K -1 points de rupture d'une segmentation à K segments intervienne à l'instant t peut par exemple être exprimée comme

La probabilité qu'un segment donné apparaisse dans m peut être obtenue de manière similaire. Une formule exacte est également donnée pour la probabilité d'apparition a posteriori au cours du temps.

Il n'est en revanche pas possible de calculer de la même manière la probabilité qu'une arête conserve le même statut tout au long de la série temporelle lorsque m est aléatoire. En effet, ce calcul nécessite d'intégrer sur des sous-ensembles de segmentations en contradiction directe avec l'hypothèse de factorisation sur les segments. Ce calcul devient néanmoins possible dès que la segmentation est fixée.

Simulations et applications

Une étude simulatoire a été effectuée afin d'étudier le comportement de notre méthode, notamment lorsque l'hypothèse arborescente n'est pas vérifiée par les graphes servant à générer les données. Nous nous sommes placés dans le cadre classique des modèles graphiques gaussiens. Nous avons comparé notre modèle à celui obtenu en n'imposant aucune structure sur la matrice de précision. Les résultats semblent indiquer que l'hypothèse arborescente pénalise très peu notre approche en termes de segmentation quand la densité des réseaux reste faible. Dans tous les cas, l'inférence semble plus stable dans le modèle que nous avons décrit, en comparaison avec le modèle non-structuré, tout en permettant l'inférence de la structure.

Nous avons également appliqué notre approche à des données d'expression de gènes récupérées au cours du cycle de vie de la drosophile [START_REF] Arbeitman | Gene expression during the life cycle of Drosophila melanogaster[END_REF] et concernant onze gènes impliqués dans le développement des muscles des ailes. Les résultats obtenus semblent cohérents avec les différents stades de la morphogénèse observés chez la drosophile. Une seconde application à des données d'imagerie par résonance magnétique fonctionnelle [START_REF] Cribben | Dynamic connectivity regression: Determining state-related changes in brain connectivity[END_REF] Dans ce chapitre, nous tentons de dégager les conditions génériques sous lesquelles les techniques d'inférence bayésiennes exactes présentées dans ce manuscrit peuvent s'appliquer. Nous présentons ensuite diverses extensions et perspectives aux travaux des deux chapitres précédents. Celles-ci se distinguent en trois parties. La première section s'intéresse à l'intégration de covariables dans les problèmes d'inférence de réseaux et de segmentation. La motivation vient d'un jeu de données issu de l'écologie microbienne où l'abondance d'un certain nombre d'espèces a été mesurée sur des feuilles d'arbres. Ces données sont accompagnées de plusieurs covariables, telles que l'arbre dont proviennent les feuilles, ou leur position dans la canopée. Il semble pertinent d'intégrer cette information à la procédure d'inférence de réseaux, dans la mesure où ce sont généralement les interactions directes entre espèces que les écologues cherchent à mettre en évidence. Dans le cas où les observations peuvent être modélisées par une loi normale multivariée, l'intégration