Graphical models for mixed data with categorical latent variables

Graphical Models for Mixed Data with Categorical Latent Variables A thesis submitted in fulfillment of the requirements for the Bachelor’s degree in Mathematics Bachelor’s degree in Data Science and Engineering Author: Luis Sierra Muntané Supervisors: Piotr Zwiernik (University of Toronto) Juan José Rué Perna (UPC) May 2023 Acknowledgments To the many wonderful people at the DoSS of UofT, and especially to my supervisor, P. Zwiernik thanks to whom I got a glimpse of what life in academia is like, and for making the cold Canadian winter feel like home. To my family and closest friends, for bearing with my numerous peaks and valleys these last 4 years. 1 Abstract This thesis aims to provide an overview of probabilistic graphical models in the context of other widely used machine learning methods, and how these methods can be formalised using conditional independence models and algebraic statistics. By comparing the extremely popular Gaussian Mixture Models and Neural Networks as generative graphical models, we are able to propose a graphical model for mixed data (discrete and continuous components) that provides a solid theoretical background and a way to analyse the Gaussian-Bernoulli Restricted Boltzmann Machine. This is used to model variables with a Gaussian conditional distribution, with discrete latent variables. On top of that, this thesis then goes into learning and sampling procedures as well as using techniques from algebraic statistics to further depict the expressibility of the model, making use of independence and mixture models for cumulant semi-algebraic varieties. Keywords: Graphical Models, Statistical Learning, Gaussian Mixtures, Markov Random Fields, CG Distributions, Restricted Boltzmann Machine, Gaussian-Bernoulli RBM MSC: 46A32, 62H22, 62R01, 62R07 Resum Aquesta tesi pretén proporcionar una visió general dels models gràfics probabilístics en el context d’altres mètodes d’aprenentatge automàtic àmpliament utilitzats, i com aquests mètodes es poden formalitzar utilitzant models de independència condicionada i estadística algebraica. A base de comparar les Mixtures Gaussianes amb les Xarxes Neuronals interpretades com a models generatius, proposem un model gràfic per a dades mixtes (variables discretes i contítnues) que proporciona una base teòrica sòlida i una manera d’analitzar la Màquina de Boltzmann Restringida Gaussiana-Bernoulli. Això s’utilitza per modelar variables amb una distribució gaussiana condicionada, amb variables latents discretes. A més a més, aquesta tesi es centra en els procediments d’aprenentatge i mostreig, així com en l’ús de tècniques d’estadística algebraica per a descriure la expressivitat del model, fent servir models d’independència i de mixtura per a les varietats semi-algebraiques dels cumulants. Mots clau: Models Gràfics, Aprenentatge Estadístic, Mixtures Gaussianes, Camps Aleatoris de Markov, Distribucions CG, Màquina de Boltzmann Restrictiva, RBM Gaussiana-Bernoulli MSC: 46A32, 62H22, 62R01, 62R07 2 Resumen Esta tesis pretende proporcionar una visión general de los modelos gráficos probabilísticos en el contexto de otros métodos de aprendizaje automático ampliamente utilizados, y cómo estos métodos pueden ser formalizados utilizando modelos de independencia condicional y estadística algebraica. Al comparar los populares modelos de Mezclas Gaussianas con las redes neuronales vistas como modelos generativos, proponemos un modelo gráfico para datos mixtos (variables discretas y continuas) que sirve de base teórica sólida a la vez que permite analizar la Máquina de Boltzmann Restringida Gaussiana-Bernoulli. Esto se utiliza para modelar variables con una distribución gaussiana condicionada, con variables latentes discretas. Además, esta tesis se explaya sobre los procedimientos de aprendizaje y muestreo del modelo, así como en el uso de técnicas de estadística algebraica para describir la expresividad del mismo, utilizando modelos de independencia y mezcla para las variedades semi-algebraicas de los cumulantes. Palabras clave: Modelos Gráficos, Mezclas Gaussianas, Campos Aleatorios de Markov, Distribuciones CG, Máquina de Boltzmann Restrictiva, RBM Gaussiana-Bernoulli MSC: 46A32, 62H22, 62R01, 62R07 3 Contents 1 Introduction 5 2 Literature Review 8 2.1 Gaussian Mixture Models ................................ 8 2.2 Probabilistic PCA ..................................... 10 2.3 Product of Experts .................................... 11 2.4 Neural Networks Overview ................................ 12 3 Graphical Models 15 3.1 Bayesian Networks .................................... 15 3.2 Markov Random Fields .................................. 20 3.3 Restricted Boltzmann Machines ............................. 26 3.4 CG Distribution for Mixed Data ............................. 28 4 Algebraic Statistics 32 4.1 Graphical Models ..................................... 32 4.2 Mixture Models ...................................... 35 5 Gaussian-Bernoulli RBM 38 5.1 Model Definition and Parameters ............................ 38 5.2 Latent Tree Model .................................... 41 5.3 Learning and Sampling Procedures ........................... 45 5.3.1 Gibbs Sampling .................................. 45 5.3.2 Expectation-Maximisation Algorithm ...................... 45 5.4 Expressibility ....................................... 46 5.4.1 Approximation Properties ............................ 47 5.4.2 Model Cumulants ................................. 48 6 Conclusion 53 4 1 Introduction In a measure space (Ω,A, ν)where Ωis our sample space, Ais a σ-algebra on Ω, and νis a dominating measure, which for us will always be the Lebesgue measure, a statistical model Mis a collection of measures {Pθ}θ∈Θdominated by ν, where in this work, the index set Θwill be taken to be a finite dimensional subset of Rd, in other words, a parametric model. We define the density as the Radon-Nikodym derivative of the measure with respect to the dominating measure ν pθ(x) = p(x;θ) = dPθ dν(x) When encountering a statistical learning task, the aim is to find a statistical model that can approximate the true population distribution, whatever that may be. Examples of statistical models may be the Bernoulli distribution for tossing a coin, a Poisson distribution to model the frequency of car accidents etc. Therefore, we need some way to compare statistical models to measure their effectiveness at their chosen task. The canonical way to do so is using the maximum likelihood estimation, especially for parametric models, due some important theoretical results, but as a way to start with the geometric spirit of this report, we may define the maximum likelihood as a sideline of a distance consideration. Definition 1.1. Let Mbe a set of probability measures. A divergence refers to a function D(·||·) : M×M → R∪{∞} satisfying (i)D(p||q)≥0 (ii)D(p||q) = 0,only if p=q Note that symmetry is not one of the requirements, and so it is often useful to consider the dual divergence D∗(p||q) = D(q||p). Now for P, Q two probability measures on (X,A), the Kullback-Leibler Divergence, also referred to as relative entropy is given by DKL(P||Q) = ZX log dP dQdP, if P≪Q(1.1) which is clearly not symmetric. If Q≪Pthen DKL(P||Q) = ∞. When p, q are the densities of P, Q respectively, this is equivalent to the more common expression DKL(p||q) = ZX plog p qdν 5 Suppose we observe a sample of ndata points given by x1, . . . , xn, then the log-likelihood function can be written as ˆ θMLE = arg max θ∈Θ{ℓ(θ)}= arg max θ∈Θ(n X i=1 log pθ(xi)) Assuming the true density is p(·)(mind you it need not be in M) then for a fixed θ∈Θ, assuming we obtain our datapoints xifrom a r.v. Xi∼pthen the expected normalized log-likelihood is given by E1 nℓ(θ)=E"1 n n X i=1 log pθ(Xi)#=ZX p(x) log pθ(x)dν(x)(1.2) Now by computing the KL Divergence of pand pθwe find DKL(p||pθ) = ZX p(x) log p(x) pθ(x)dν(x) =ZX p(x) log p(x)dν(x)−ZX p(x) log pθ(x)dν(x) =−H(p)−ZX p(x) log pθ(x)dν(x) where we have used Hfor Shannon’s entropy. Comparing this previous equation with that in 1.2 we can interpret the maximum likelihood as the minimisation of DKL(p||pθ)“in the sense of expectation", because by assuming constant ground truth p, it is found by minimising the second term from the previous equation, which is just the log-likelihood. In this way, the maximum likelihood estimation can be seen as having an information geometry interpretation, and in this way we can think of finding the best possible distribution inside a given family of models for extracting the most information from a specific real-world task. The MLE can also already serve as an example of the discrepancy between models for discrete and for continuous variables, since the likelihood function for discrete data is well defined in simple terms, whereas for continuous variable models the likelihood loses probabilistic meaning. Models for mixed variables, discrete and continuous pose great theoretical difficulties as they fail to combine nicely in many instances like oil and water. This work will be dedicated at exploring how to deal with this problem and what simplifications can be made to give rise to meaningful models for such a combination. In recent years, the tremendous progress and success in machine learning applications has come about mostly through cruder engineering work, and the statistical knowledge of the models used in ML has lagged behind. Maximum likelihood has been replaced by the mean squared 6 error (MSE), and empirical risk measures based on a more pragmatic and experimental heuristics (albeit some equivalences have been found). A current concern of modern ML models is their lack of interpretability, and with the new architectures coming out for image classification, large language models and other industry requirements, this statistical understanding is straggling even more. This phenomenon is probably a consequence of the very difference between statistics/ mathematics and machine learning, since the former is interested in obtaining theoretical results and extracting inferences from data, the latter has a particular goal centered in applications. Such a contrast is akin to that of physics and engineering, where theory and application in a same field can lead to vastly different cultures in a community. The aim here is that of recovering some of the lost ground from the mathematical statistics perspective, starting with this project, where an exploration into the probabilistic foundations of neural networks and other related models will be given. On top of that, an analysis will be made of a network-like model involving independent Gaussian variables which has been seldom been looked at from an algebraic or statistical perspective, seeing how we can use various tools in graphical models, algebraic statistics and tensors to examine its properties. Each of these areas individually are huge; graphical models as a way to analyse conditional dependencies between variables are an active area of study with applications in many different domains, and the invention of algebraic statistics as a way to bridge the power and generality of algebraic geometry into statistical challenges has proved to be very fruitful in the past decades. We hope that looking at their intersection in this work can lead to some interesting new ways to look at the problem of data modelling and have future contributions in data science. To motivate this intersection and analysis, we will start with a look at several well-known models for data analysis and feature extraction which will contain many of the key ideas we will use later on. As such, this work is meant to serve as an overview of how mathematics, statistics and machine learning can fit together, provide some interesting results about why some things work the way they do, and try to leverage their respective strengths for each other’s benefit: new model ideas, better understanding of their properties and better results in their subsequent applications. 7 2 Literature Review 2.1 Gaussian Mixture Models What can be done when a dataset that could reasonable come from a normal distribution does not fit the model? This question came to Karl Pearson in the 19th century when exploring a dataset on crab dimensions [1]. He observed that the distribution was not symmetric as a normal distribution should be, to a larger extent than what could have happened by random chance. His solution to the problem was creative as well as founded: assume that the crabs came from two different populations, normally distributed, and the observations observed were the result of the sum of two independent normal distributions. He then fitted the distribution functions using methods that today would seem extremely inefficient, by solving a system of polynomial equations in the first two moments of each distribution (the means and the variances) as explained in [2] which relates this technique to algebraic statistics, which will be explored later on. However, the main rationale of assuming a distribution consisting of an additive combination of other distributions was born, and is known as a mixture distribution. As such, the idea would be encapsulated by the following equation, where we assume our overall distribution comes from the combination of Kseparate, independent Gaussians Xi∼ N (µi,Σi). p(x) = K X k=1 πiϕk(x), K X k=1 πk= 1,0≤πk≤1(2.1) Where the ϕiare the probability density functions for each variable Xi. In modern notation, the actual object of study is known as a mixture distribution, defined as follows. Definition 2.1. Amixture distribution consists of a family of probability distributions FX|T=t with distribution function FX|T=t(x)where tis not a fixed parameter but comes from a family of distributions FT, in such a way that we have FX=FX|T=t∧FTthe mixture, where X|T=thas distribution function FX|T=t(x). Proposition 2.2. The distribution of FXis given by FX(x) = ETFX|T=t Proof. Directly from the definition, taking FX(x) = ZR fX(y)dy =ZRZR f(X,T )(y, t)dtdy =ZRZR fT(t)fX|T=t(y, t)dtdy and now using Fubini we get ZR fT(t)ZR fX|T=t(y, t)dydt =ZR fT(t)FX|T=t(x)dt =ET[FX|T=t] 8 3 Graphical Models In this chapter we will see a natural way to model random vectors with dependence relationships between their components, which generalise the factorisation of the previously seen models of Gaussian mixtures 2.1, Probabilistic PCA 2.2, factor analysis, and beyond, using graph theory. Before starting with any definitions, a short note on conditional independence in terms of the covariance. Definition 3.1. The covariance between two random variables X, Y is given by Cov(X, Y ) = E[(X−E[X]) (Y−E[Y])] For a set of random variables X1, . . . , Xnwe may define the covariance matrix as Σ = Cov(X) = {Cov(Xi, Xj)}ij. When any two variables are independent, then their covariance is equal to 0 since in particular they will be mean independent, but the reverse implication is only true whenever X follows a multivariate normal distribution. We may also define the precision matrix as K= Σ−1, which exists whenever the covariance matrix is full rank. By direct computation using a diagonal change of basis we obtain that zeros in the i, j entry of the precision matrix correspond to the variables Xi, Xjbeing conditionally uncorrelated given all other variables. Kij = 0 ⇐⇒ ρXi,Xj|X\{i,j}= 0 For this reason, this matrix is quite relevant in the handling of graphical models. 3.1 Bayesian Networks Graphical Models have been used in a variety of different contexts for modelling, most prominently in statistical physics and phylogenetics, and as such we can find some of the terminology spilling in from these fields. In essence, their aim is to model the dependency relationships between several attributes, modelled in turn by random variables. For a graph G= (V, E)with vertices Vand edges Ecomprising of pairs of vertices, we can distinguish between directed and undirected graphs, according to whether the elements of Eare ordered or unordered pairs, respectively. Accordingly, for a graphical model, we take Vto be a set of random variables. When the graph is directed, the graphical model is referred to as a Bayesian Network while when the graph is undirected, the standard term in the literature is a Markov Random Field. We will start off with the case where Gis a directed graph. In order for the dependence relations 15 to be defined correctly, this underlying graph Gmust not contain any cycles, as will shall see in 3.2, that is, it has to be a Directed Acyclic Graph. But first, we shall see how the graph is used to encode the dependencies between the variables. When the graph contains an edge (a, b)∈E, we say ais a parent of band they are often referred to as the cause and the effect, respectively. This nomenclature points towards the usefulness of these models in the study of causality and general causal inference, (see for example [18]). Let v∈V, we will denote as pa(v), the set of parents of v, given by the vertices wsuch that (w, v)∈E. The set de(v)of descendants of vis given by the vertices wfor which there exists a directed path v→win G, and we may also define the set nd(v)of non-descendants as V\({v}∪de(v)). The local Markov property of conditional independence states that Xv⊥⊥ Xnd(v)\pa(v)|Xpa(v), v ∈V(3.1) That is, given the parent nodes of a vertex v, the variable Xvis independent from all others. From it, we obtain the key property used to define graphical models, that of their factorisation. Using this notion for conditional independence, we may recursively factorise the joint density of the random vector X= (X1, . . . , Xn), where V={Xi}n i=1, as: p(x) = n Y i=1 p(xi|pa(xi)) (3.2) Where it is understood that p(x) := p(X=x). This is often known as the chain rule of Bayesian networks, following the chain rule of conditional probability and Markov chains. This can be intuitively seen as a design in which the effects are being factorised in terms of their causes, in a way where cause and effect are related sequentially, and two variables are independent if there is no path between them, or if we know the value of a variable found inside the path connecting them. To illustrate this point is the following illustrative example of a Bayesian Network. Example 3.1.1. Consider the graphical model associated to the following graph over the variables A, B, C, D. 16 A B C D Where we would have the factorization p(a, b, c, d) = p(a|b, c)p(c|b)p(d|b)p(b) Now, we are ready to see why in order for Bayesian Networks to be well defined, the underlying graph must be acyclic. Proposition 3.2. The underlying graph in a Bayesian network must be a DAG. Proof. Let A, B, C be three binary variables with dependency graph given by the following directed graph, and where the value of each parent completely determines the value of each child, meaning P(A=x|B=x)=1, as depicted in the probability tables below. A B C a0a1 b01 0 b10 1 b0b1 c01 0 c10 1 c0c1 a01 0 a10 1 In this way, if we calculate the sum of all events (which should be 1) using the Bayesian Network factorisation from 3.2 yields X A×B×C p(A=a, B =b, C =c) = X A×B×C p(B|A)p(A|C)p(C|B)=2>1 17 The problem that arises is due to the fact that it is not possible to guarantee that the distribution is correctly normalised in this case, whereas when we have a DAG, we can factorise the distribution starting from a root node (which exists due to there being no cycles) and can easily see how the factorisation produces a distribution with total measure 1. This definition is deceptively simple, and we can find, even with small examples, already some interesting phenomena. Example 3.1.3. Consider three r.v.s X, Y, Z in a v-structure: X→Y←Z. Then we have that p(x, y, z) = p(x)p(z)p(y|x, z). We can observe that X⊥⊥ Zsince their marginal distributions factor completely, but when we condition on Ywe obtain p(x, z|y) = p(x)p(z)p(y|x, z) p(y) And so despite the marginal independence of Xand Z, they are not independent when we condition on Y! In some sense, we can think of two parents being independent but when we condition on a child, whatever is not explained by one of the parents must be explained by the other one, so the independence is broken. One of the most well-known general classifiers is that of the Naive Bayes model, introduced as a simple and fast to train classification scheme that is often taught as a way introduce more complex classifiers later on. Example 3.1.4. (Naive Bayes Classifier) Consider a collection of so called "features", that is, some vector of random variables X= (X1, X2, . . . , Xm)and a variable Zrepresenting to which of k classes an observation may belong. In order to predict for an observation of data x= (x1, . . . , xm) to which class zit belongs, Bayes’ rule is applied in order to compute the posterior probability of the class: P(z|x) = P(z)P(x|z) P(x)=P(z)Qm i=1 P(xi|z) P(x) Where now we would estimate P(xi|z)according to the observed frequencies from a given dataset with respect to the class of the datapoint, and using Bayes’ decision rule we would assign the label with the highest posterior probability: z= arg max z{P(z) n Y i=1 P(xi|z)} And thus the inference of the parameters of the model, which would consist on the class frequencies, is performed in a simple and direct manner. This ease of calculation comes at a cost of 18 course: such an independence assumption is a rather strong constraint on the data, and so this "naivety" of the model which comes in assuming the conditional independence of features (variables) of the observations given the class Z, is a rather optimistic way of modelling data. As such, this could be represented as the following Bayesian Network. Z X1X2X3Xm The node for the variable Zis coloured in black to indicate that this is a latent variable, meaning that it is not observed, as with the population in a Gaussian Mixture. In this way, Naive Bayes can be thought of as a generative graphical model with the simple dependence structure Xi⊥⊥ Xj|Z. Compare such a simple model to the task of inference or learning in a general Bayesian Network with several observable and latent variables. When learning from complete data, one can go about maximising the log-likelihood directly with respect to the parameters θ, or with some other improved numerical approximations, but in any case the approach is straightforward for a sample of size N. ˆ θ= arg max θ∈Θ(N X i=1 ln p(xi|θ))(3.3) Where the expression for p(xi|θ)will be computable whenever the Bayesian Network structure is known. We will not be looking at validating the structure of a Bayesian network, or any graphical model in this analysis, as we will take Gas known. Note that, so far, no assumptions are made on the distribution of each of the individual variables in the graphical model, and this general structure is agnostic to such a choice. We will impose concrete distributions later on, but for now, assume each variable in the network has a known distribution. In fact, for the factor analysis model, a key assumption is that the variables are jointly distributed as a multivariable normal, and the graphical model is a Bayesian Network. A classical maximum likelihood scheme would work to fit the parameters once we assume said distribution, but for models with latent variables, despite the simplicity of the Naive Bayes model 19 where there was a single latent for which we had known prior distributions, in general we run into the problem that we do not know the full expression of the likelihood function, and thus it is not possible to maximise the expression in 3.3 directly. As such, the main approach to avoid this pitfall is that of using the Expectation Maximisation Algorithm, introduced in the late 1970s, which follows the same ideas for the MLE as in the introduction, but averaging out the latent variables. Definition 3.5. (Expectation-Maximization) The EM algorithm is an iterative procedure used to find the maximum likelihood estimation when there are latent variables involved. Let V= (Z, X) where Zare latent and Xare observables and define, for a current set of parameters θ(t)the function Q(θ, θ(t))as the partial likelihood defined below. E-Step: Q(θ, θ(t)) := EZ∼p(·|X,θ(t))[ln p(X, Z|θ)] =X Z ln p(X, Z|θ)p(Z|X, θ(t)) M-Step: θ(t+1) := arg max θ{Q(θ, θ(t))} These two steps are performed iteratively until a given precision is satisfied. For the M-Step, assuming some regularity for the distribution, typical gradient ascent schemes are often used [19]. Due to the procedure being monotone increasing in the estimated parameters, it produces a consistent estimator (convergent in probability) whenever the likelihood has a single local maximum, which for practical purposes is quite a strong assumption. As always, when working with a non-convex likelihood function, using multiple starting points is the practical solution to avoid missing the global maximum, without going into concerns about the degree of non-convexity [17]. 3.2 Markov Random Fields On the other hand, we have the case where the graph Gis undirected, which can be thought of as having both directional edges for each pair of connected variables. As before, we take G= (V, E) with Va set of random variables, but this change in Ecauses the graphical models to be very different, given that the training/inference procedures and whole description change substantially. They often take the name in the literature of Markov Random Fields and many of their uses came 20 from statistical mechanics, as a way to model the energy of particles set in a lattice structure, where the energy in a node depends on the energy in the neighbouring nodes. As before, we may define conditional independence via graph separation, only now the factorisation we defined in 3.2 doesn’t make sense. This time, the factorisation of the joint distribution is given in terms of potential functions, again borrowing nomenclature from physics. In fact, the way in which we need to preserve the notion of locality (which previously was just path-connectedness in a DAG) in the factorisation in order for the conditional independence to be well defined, will be done by factorising the joint density into functions of the correct subsets of V. Consider two nodes Xi, Xj∈Vwhere {Xi, Xj}/∈E. We should require that these two be independent given the rest of the variables, since they will be separated in G p(xi, xj|x\{i,j}) = p(xi|x\{i,j})p(xj|x\{i,j}) This restriction is known as the pairwise Markov property, which we will name (P), but a stronger condition, the local Markov property or (L) of the graph is that of requiring that only the variables that "surround" a node xibe required for conditional independence with any other node, that is, if we define bd(v)to be nodes that are connected to v∈V, then (L) states: v⊥⊥ V\{bd(v)∪{v}}| bd(v) This only refers to a single variable, and so the strongest condition, known as the global Markov property or (G) states that: XA⊥⊥ XB|XCwhenever Cseparates A and B in G(3.4) In many texts such as [20], the previous condition is also written abusing the notation by using A, B, C ⊂Vand saying: A⊥⊥ B|Cwhenever Cseparates A and B in G Remark 3.1. In the previous equations for condition (G), we only mentioned the reverse implication instead of stating it as an if and only if condition. The reason for this is that the implication from left to right states that the graph Gis faithful to the distribution, and would be concern when trying to find a graph that could model a given distribution. In practice this is something that would be evaluated according to conditions from the data. This technical detail will not be relevant to our study since we take the graph Gas given, but it is nevertheless an important condition for the general study of graphical models. 21 Proposition 3.1. It is easy to see that for any graph Gand probability measure Pon Xwe have that (G) =⇒(L) =⇒(P) Proof. Pretty straightforward, just notice that bd(v)is a separating set between vand V\{bd(v)∪ {v}}. Then, for the second implication, by non-adjacency we have that w∈V\{bd(v)∪{v}} and we may observe that bd(v)∪((V\{bd(v)∪{v}})\{w}) = V\{v, w}. The reverse implications require quite a bit more work, and an extra assumption, in the form of the property that for disjoint sets A, B, C, D we have that Proposition 3.2. Suppose the density over the graph Gis strictly positive, then A⊥⊥ B|(C∪D)and A⊥⊥ C|(B∪D) =⇒A⊥⊥ (B∪C)|D(3.5) This condition makes intuitive sense but is not for free; some degenerate graphs may not satisfy it, and it is equivalent on some condition on the dual graph of Gthat we will not go into. Fortunately, any probability measure with continuous and positive density will satisfy this, and this is already a more than strong enough assumption [20]. Proof. The two conditional independences imply the density factorisations pA∪B|C∪D(xA,xB|xC,xD) = pA|C∪D(xA|xC,xD)pB|C∪D(xB|xC,xD)(3.6) pA∪C|B∪D(xA,xC|xB,xD) = pA|B∪D(xA|xB,xD)pC|B∪D(xC|xB,xD)(3.7) Multiplying 3.6 by pC∪D(xC,xD)and 3.7 by pB∪D(xB,xD)we obtain the full density, and so we can equate the results on the right hand side, and dividing by the (positive) density pB∪C∪D(xB,xC,xD) we obtain pA|C∪D(xA|xC,xD) = pA|B∪D(xA|xB,xD) Since the left-hand-side does not depend on xB, we have that pA|C∪D(xA|xC,xD) = pA|D(xA|xD) and so if we take the product we had for 3.6 and condition on xDwe get: pA∪B∪C|D(xA,xB,xC|xD) = pA|D(xA|xD)pB∪C|D(xB,xC|xD) Which implies the independence result we wanted. Theorem 2. (Pearl and Paz) If a probability distribution on Xis such that 3.5 holds for all disjoint sets A, B, C, D, then (G)⇐⇒ (L)⇐⇒ (P) 22 Proof. It is sufficient to show that (P) =⇒(G). To do so, we will show that for sets A, B, S where Sseparates the non-empty sets A, B we have that condition (G)holds. By reverse strong induction on the size of S, let n=|S|and consider n=|V|−2, then, since both Aand Bconsist of a single vertex, (P)is trivially enough for the conditional independence. Now assume (G)holds for all sets of more than nelements and take n < |V|−2. Take V=A⊔B⊔Sand so either of A, B has more than one element. Without loss of generality, assume Ahas more than one element. Take α∈Aand note S∪{α}separates A\{α}and B, as well as S∪(A\{α})separates {α}and B. Thus, from the induction hypothesis we have the two conditions A\{α} ⊥⊥ B|S∪{α}and {α} ⊥⊥ B|S∪(A\{α}) =⇒A⊥⊥ B|S where the implication follows from 3.5. We can of course have that A∪B∪S⊂V, in which case we can simply take α∈V\(A∪B∪S)and then observe that S∪{α}separates A, B, and either S∪Aseparates {α}, B or S∪Bseparates {α}, A from which we can apply 3.5 to yield our desired result. These so-termed Markov conditions will force us to factorise p(x)into a product of the distribution over the maximal cliques of G, which are the maximal complete induced subgraphs of G. In other words, if we let C(G)be the set of maximal cliques in G, then we can factorise the joint density as p(x) = 1 ZY C∈C(G) ψC(xC), Z =ZXY C∈C(G) ψC(xC)(3.8) Where Zis the normalising function with a similar role as that in the product of experts model we had in 2.3, and the functions ψCare known as the potential functions, which are non-negative in order for the density to be non-negative. The fact that this factorisation is consistent with the conditional independence condition in 3.4 is provided by the Hammersley-Clifford theorem. We will need a small combinatorial lemma to do so. Lemma 3.1. (Möbius Inversion) Let Vbe a finite set and let Ψ,Φbe functions defined on P(V) taking values in an abelian group. Then, the following two statements are equivalent: (1) for all a⊆V: Ψ(a) = Pb⊆aΦ(b) (2) for all a⊆V: Φ(a) = Pb⊆a(−1)|a\b|Ψ(b) Proof. Both directions are proven analogously, using a kind of sophisticated inclusion-exclusion 23 idea. For (2) =⇒(1), take X b⊆a Φ(b) = X b⊆aX c⊆b (−1)|b\c|Ψ(c) =X c⊆a Ψ(c)(X b:c⊆b⊆a (−1)|b\c|) =X c⊆a Ψ(c)  X h⊆a\c (−1)|h|   the last sum is equal to zero except for the case where a\c=∅, since for a non-empty set there is the same amount of sets of even degree as of odd degree. Now we are ready to prove our factorisation theorem. Theorem 3. (Hammersley-Clifford) A probability distribution Pwith positive and continuous density pwith respect to a product measure satisfies (P)with respect to a graph Gif and only if it factorizes according to G. Proof. We of course only need to prove the reverse implication, since the factorisation already implies conditional independence and thus the (P) property. To do so, first rewrite 3.8 (thanks to the positivity of the density) as log p(x) = X a∈C(G) ϕa(xa)(3.9) where ϕa(xa) = log ψa(xa). We are only assuming Pis pairwise Markov. Let us fix an element x∗∈ X of the sample space and define for all a⊆Vthe functions Ha(x) := log p(xa, xaC∗), and also ϕa(x) = X b⊆a (−1)|a\b|Hb(x) we write these as a functions of xbut of course each one only depends on xa. Now, by lemma 3.1 applied to ϕawe have that log p(x) = HV(x) = X a⊆V ϕa(x) We now need to show that ϕa(x)≡0when adoes not induce a clique in V. To do so, let us take α, β ∈asuch that α∼ βand consider the following decomposition of ϕausing inclusion-exclusion: ϕa(x) = X b⊆c (−1)|c\b(Hb−Hb∪{α}−Hb∪{β}+Hb∪{α,β})(3.10) 24 moment parameters g(h) = log p(h)−n 2log 2π+1 2log det K(h)−1 2G(i)TK(h)−1G(h)(3.22) = log p(h)−n 2log 2π−1 2log det Σ(h)−1−1 2µ(h)TΣ−1µ(h)(3.23) In this way we have found explicit ways to change the parameters of the distribution, from what Lauritzen calls the canonical characteristics (g, G, K)(due to them actually being equal to the canonical parameters when written in exponential form) and the moment characteristics (p, µ, K−1) If we specify further the nature of the categorical variables, and looking to transition into the model family we will look at in later chapters, we can find a more concrete expression for p(h), such as, for instance, when we have mbinary variables H1, . . . , Hm(as happens when we have a directed graphical model to represent the conditional independence of our variables). In this case we can compute the probability of a particular state has: p(h) = Pr(Hi=hi:i∈[m]) = Y i∈[m] bihi(1 −bi)1−hi(3.24) Where it is understood that Hi∼Bern(bi)so we have transformed a categorical variable on m states into mindependent binary random variables. This binarisation is very apt for calculations and can be thought of as being reminiscent to a sort of one-hot encoding to display the variables in a more convenient manner, but allowing more than one state to be active at once, so there are actually 2mpossible configurations. For us, this configuration is a way to reparametrise our latent space that will be useful when defining Gaussian-Bernoulli RBM models. 31 4 Algebraic Statistics In the previous chapter, we defined the notion of a graphical model and looked at the two main types in terms of the type of edges of the underlying graph. In this section, we will see how these models can be defined in an algebraic way, and what can we gain from such an interpretation. At a high level, the gist is that of using tools from algebraic geometry to attain results in statistics [27]. This idea came from the seminal paper from Sturmfels and Diaconis [28], which defined the concept of a Markov basis and marked the birth of algebraic statistics as a separate field by using commutative algebra to sample from a conditional distribution. 4.1 Graphical Models We will start by providing some definitions and then move on to how this algebraic setting translates the graphical models we defined in the previous chapter. The point of all of this is that we can express graphical models algebraically in terms of the independence ideals that are generated by the conditional independence equations. Most of the notation and style is thanks to [27]. Definition 4.1. Ak-simplex is a k-dimensional polytope given by the convex hull of k+ 1 (independent) vertices u0, . . . , ukas in C=(λ0u0+···+λkuk k X i=0 λi= 1, λi≥0)(4.1) Since our context is that of probability theory, all of our vertices will be given by the canonical basis, as in ∆k=λ∈Rk+1 :λ0+···+λk, λi≥0(4.2) And here we don’t need to specify the use of barymetric coordinates. These simplices are useful because they are a geometrical representation of all the probability distributions on a discrete set of k+ 1 elements. We may think of a distribution as a point pi1...ik∈∆k−1. In the case of the binary distribution, nrandom variables take joint states in {0,1}n, the probability of each state ai1...in∈ {0,1}ncan be described as pi1...in= Pr(H=a) = Pr({Hj=ij:j= 1,...n}). Therefore, we have that (pi1...in)i1...in, ij∈ {0,1}is a 2×n) . . .×2tensor storing the probability of each state, and the elements lie in the simplex ∆2n−1. Observation 4.2. Following from the previous definition of pi1...init is standard to define the 32 marginal probabilities with the notation pi+:= X j pij = Pr(X1=i)(4.3) Then, we can think that two variables are independent when the joint probability tensors factor as pij =pi+p+j, for all i∈[r1], j ∈[r2]. What does independence mean algebraically? The discrete case will be simpler to illustrate. Take X= (X1, . . . , Xm)to be discrete r.v.s and let [rj]be the set of values taken by Xj, that is, Im Xj= [rj], and so Xtakes values in R=Qm j=1 [rj]. Then, we shall see that conditional independence constraints correspond to a system of quadratic polynomial equations. Proposition 4.3. The conditional independence statement XA⊥⊥ XB|XCholds if and only if piA,iB,iC,+·pjA,jB,iC,+−piA,jB,iC,+·pjA,iB,iC,+= 0 (4.4) for all iA, jA∈ RA, iB, jB∈ RB, iC∈ RC. Observation 4.4. Two random variables are marginally independent, X1⊥⊥ X2, if and only if the r1×r2matrix (pij)has rank 1. Definition 4.5. The conditional independence ideal IA⊥⊥B|Cis generated by all quadratic polynomials in 4.4, and for two random variables X1, X2, it would be represented as IX1⊥⊥X2=⟨pi1i2pj1j2−pi1j2pj1i2|i1, j1∈[r1], i2, j2∈[r2]⟩(4.5) To provide a taste of the usefulness of this theory, [29] provides us with this result using fully algebraic methods. Proposition 4.6. The Restricted Boltzmann Machine with 2 hidden units and 3 observables RBM2,3is described in the interior of the ∆7simplex by the union of the algebraic sets {p000p011 ≥p001p010, p100p111 ≥p101p110} And the other 5permutations of indices. These sets represent the conditional independence between variables, which we know from 3.1. Note that the model is over-parameterised. In this case, the model has 11 parameters but only lives in ∆7. 33 The main result of [29] is concerning the tree model with a 3-state latent categorical variable and 3 binary leaves: Theorem 4. M3,3=RBM2,3with M3,3=RBM2,3on the interior of the simplex. The proof of this result involves purely algebraic methods, and is a nice way to illustrate how this theory can be useful to find model equivalences. We will be using to find to describe the expressibility of a concrete model, but if we are able to find a model with the same algebraic description, then we can show it is equivalent to ours. This was for categorical distributions, where it is easy to represent the entirety of a distribution as a single point. For absolutely continuous distributions, whose domains are uncountable subsets of Rn, this would not work. However, a huge benefit of modelling our continuous variables using Gaussian distributions is that it is actually possible to describe one with only two values, by using a classical result from Marcinkiewicz: Theorem 5. (Marcinkiewicz 1935) Xfollows a multivariate Gaussian iff KX(t)is a polynomial. We have not yet defined what KX(t)means, but we will do so in 5.1. The point is that a Gaussian distribution is completely determined by its mean vector and covariance matrix, and all its higher order cumulants are 0. Proposition 4.7. Let X∼ N(µ, Σ). The conditional independence statement XA⊥⊥ XB|XCif and only if the covariance matrix ΣA∪C,B∪C= Cov(XA∪C, XB∪C)has rank |C|. Proof. For this proof we will need to calculate the conditional distribution of a Gaussian, which is a well known result1: XA∪B|XC=xc∼ N µA∪B+ ΣA∪B,CΣ−1 C,C(xC−µC),ΣA∪B,A∪B−ΣA∪B,CΣ−1 C,CΣC,A∪B Thanks to the variables following a multivariate Gaussian distribution, the conditional independence statement will hold if and only if the A, B components (Arows and Bcolumns or vice-versa) of the covariance matrix are equal to zero, that is, whenever ΣA∪B,A∪B−ΣA∪B,CΣ−1 C,CΣC,A∪B= ΣA,B −ΣA,CΣ−1 C,CΣC,B = 0 Since ΣC,C is invertible of rank |C|due to being a covariance matrix, and the previous expression is the Schur Complement of the matrix ΣA∪C,B∪C= ΣA,B ΣA,C ΣC,B ΣC,C! 1thanks to stackexchange for the full derivation 34 Hence, we have that rank(ΣA∪C,B∪C) = |C|. Therefore, since the condition of having a particular rank kcorresponds to a set of equations being equal to 0 (the k+ 1 minors of the matrix) then this conditional independence condition may also be expressed as an algebraic constraint. As before, this allows us to define an independence ideal like in 4.5, which Sullivant calls The Gaussian conditional independence ideal JA⊥⊥B|C: JA⊥⊥B|C=⟨(|C|+ 1) −minors of ΣA∪C,B∪C⟩(4.6) And for a collection of conditional independence statements over our variables given by C= {XA1⊥⊥ XB1|XC1, XA2⊥⊥ XB2|XC2. . . }we may define the Gaussian conditional independence ideal as JC=JXA1⊥⊥XB1|XC1+JXA2⊥⊥XB2|XC2+. . . (4.7) As we have seen, we may use the independence ideals to turn the conditional independence constraints of a given graphical model into algebraic ideals, such the conditions (P)and (G)of Markov Random Fields using either 4.5 for discrete models or 4.6 for Gaussian graphical models. The challenge is that these conditions are easily defined when all the variables are discrete, as with the RBM, or when they are all Gaussian, where we can look at zeros in the precision matrix (equivalently vanishing minors of the covariance matrix). For a mixed model, the task would be quite hard, but since our particular mixed model will be such that the latent variables will be conditionally binary and the observables will be conditionally Gaussian, we are able to sidestep a lot of the complexity using mixture models. 4.2 Mixture Models Even though we have not gone into much depth on the issue of hidden variables, simply stating that they were not observed and so they were not available to calculate the MLE, they suppose a challenge because standard asymptotic theory does not apply to them. The solution is for us, since we will be taking discrete hidden variables, is to model them using mixture models. This is because given a state of the latent variables, we will have a different observable distribution, and the marginal distribution will be a mixture of the different observables. Definition 4.1. Let V, W be two parameter sets, their mixture model Mixt (V, W)is given by Mixt (V, W) = {λv + (1 −λ)w:v∈V, w ∈W, λ ∈[0,1]} 35 If you recall the factor analysis model discussed in the literature review, it could be modelled as a Gaussian graphical model with a bipartite structure, and with some of the conditional independence conditions that we had for the RBM and will want for our conditional Gaussian model in the next section. As such, it is a good example to illustrate what we are looking for. The factor analysis model with mlatent variables (H1, . . . , Hm)and nobservables (X1, . . . , Xn) denoted Fm,n, satisfies the conditional independence constraints X1⊥⊥ X2⊥⊥ . . . ⊥⊥ Xn|(H1, . . . , Hm) Where the variables follow a joint multivariate Gaussian distribution. Note that as a graphical model, this would be given by a DAG with the edges pointing towards the visible variables. Thanks to the conditional Gaussian assumption, we can parametrise the model as the space of vectors µ∈Rnfor the mean and covariance matrix Σlying in the cone: Fm,n =Ψ + ΛΛT∈Rn×n: Ψ ≻0diagonal,Λ∈Rn×m(4.8) Which is a semi-algebraic set. To see this, note that the Cov X H!= Σ Λ ΛTΦ! Which is an (n+m)×(n+m)matrix, and by 4.6 we have that the conditional independence assumptions of the model translate to the vanishing of the m+ 1 minors of the covariance matrix containing all the terms for the hidden variables, since those are the variables we are conditioning on. Assuming i=jwe obtain: det σij Λi∗ ΛT j∗Φ!= det (Φ) ·σij −Λi∗Φ−1ΛT j∗= 0 By the Laplace expansion on the first entry. Hence, since det (Φ) = 0 it will be invertible and we can define the Schur complement Ψ=Σ−ΛΦ−1ΛTand now solving for the covariance we get Σ = Ψ + ΛΦ−1ΛT, and taking Φ≽0then we can use its Cholesky decomposition and we get the algebraic description for the covariance in equation 4.8. And we have an algebraic description of the Factor Analysis model. This is only valid for a joint Gaussian distribution of all variables, which is not the case for the model we will be looking at, only the observable will follow a conditional Gaussian. For our mixture model, let’s assume a hidden categorical variable Hwith state space [s], and for each j∈[s], we have a distribution of the observable random variable X, given by Pj. Let 36 πj=P(Y=j)so that π∈∆s−1. In this way, we’ve made the marginal distribution of Xa mixture distribution. If we let p(j)∈ P, then we may write our distribution as Mixts(P) = (s X j=1 πjp(j):π∈∆s−1, p(i)∈ P)(4.9) In the following chapter, we will be taking p(j)be conditionally independent Gaussians and analyse the model’s properties. 37 5 Gaussian-Bernoulli RBM Now that we have defined graphical models and seen some of their properties as related to their factorisation, inference and general computations, as well as seen how to model mixed variables through the Conditional Gaussian, we are ready to look at a class of models that have seldom been studied in statistical circles. Their use has been limited to the machine learning community for feature extraction and classification tasks in high-impact articles such as [30] or [31] and so results about them are more centered around their performance with respect to specific tasks or data sets, rather than having results on their theoretical/statistical foundations. Therefore, in this main section we will provide a motivated and sound definition for the model, properly define it in the context of probabilistic graphical models, we will examine its statistical properties, and refrain from going into the machine learning discussions concerning hyperparameter selection or adequateness for a specific learning task, for which there already is extensive literature. Building upon the Restricted Boltzmann Machine from 3.13, we will now consider the case where the observable variables are continuous and normally distributed, and we will consider the effect of different categorical latent structures, which will allow us model continuous data while keeping the simple discrete latent structure as before. In this sense, our model is suited for mixed data, but with the requirement that the discrete variables be latent while the observables be continuous (which we will simplify as Gaussian). For this reason, and as seen in the literature review, this model has a Gaussian mixture structure when marginalising over the visible variables but at the same time will retain the product structure that eases training 5.1 Model Definition and Parameters As previously mentioned, our model is based on the bipartite structure of the Restricted Boltzmann Machine, with mconditionally binary latent variables H= (H1, . . . , Hm)and nconditionally Gaussian observable variables Y= (Y1, . . . , Yn). Yi⊥⊥ Yj|H=h Hi⊥⊥ Hj|Y=y(5.1) Definition 5.1. The Gaussian-Bernoulli graphical model with mlatent variables and nobservables, denoted GRBMm,n is defined by the conditional independence rules in 5.1 and the conditional 38 distributions Yi|H=h∼ N (µ(h),Σ(h)) Hi|Y=y∼Bern (b(y)) We will find the explicit expressions for the moment parameters b, µ, Σlater on. In fact, the parameters of the model will be written following convention as θ= (a,b,Σ, W), where Σ(h)=Σ so there will no dependence on the latent state to simplify the model, and since we will assume conditional independence by the graphical model structure, this matrix will also be diagonal. Then, awill affect the mean of the Gaussians and Wwill be the interaction term as in 3.1. Y1Y2Y3Y4Y5Y6Y7 X1X2X3X4X5 Figure 3: The graphical representation of GRBM7,5with 7 latent and 5 observable variables. We may leverage theorem 3along with the fact that the maximal cliques in a bipartite graph are just the edges to directly construct the density of the GRBM, which will be very similar to the RBM. p(h,y;θ) = 1 Z(θ)exp{−E(h,y)}, θ = (W, Σ,a,b) E(h,y) = −1 2(y−a)TΣ−1(y−a)+hTb+yTWΣ−1h (5.2) Note the term highlighted in blue is quadratic in yand will provide us with the conditional Gaussian distribution of the observable variables. Proposition 5.2. The density in 5.2 indeed yields the properties promised in 5.1. Proof. To check the conditional distributions, we may compute directly from 5.2, where the integral is rather direct as a Gaussian integral, we later complete the square, and the only problem is the 39 somewhat hairy notation. p(y|h) = p(h,y) Ryp(h,y) dy =exp−1 2(y−a)TΣ−1(y−a) + hTb+yTWΣ−1h Ryexp−1 2(y−a)TΣ−1(y−a) + hTb+yTWΣ−1hdy =Qn i=1 expn−(yi−ai)2 2σi+bTh+yi σiPm j=1 Wijhjo Qn i=1 σi√2πexp1 2Pm j=1 Wijhj2+bTh+bi σiPm j=1 Wijhj = n Y i=1 1 σi√2πexp −1 2σi2 yi−ai−σi m X j=1 Wijhj!2  Which we can see is the density of a multivariate Gaussian with diagonal covariance matrix. For the conditional density of a single latent with respect to the observables: p(hj= 1|y) = Phk=jp(hj= 1,hk=j,y) Phexp{−E(h,y)} = expnPn i=1 yi σiWij +bjoPhk=jexp{−E(hj= 0,hk=j,y)} Phk=jexp{−E(hj= 0,h,y)}+Phk=jexp{−E(hj= 1,h,y)} =1 1 + expn−Pn i=1 xi σiWij +bjo Corollary 5.1. We can now see the explicit expression for the parameters of the conditional distributions: Y|H=h∼ N (a+ ΣWh,Σ) Hj|Y=y∼Bern S n X i=1 yi σi wij +bj!! Where S(x)denotes the sigmoid function, which again makes an appearance, of a similar nature as in the RBM case. We can also infer the conditional independence between variables so that the 40 5.4.1 Approximation Properties Recall we had a Universal Approximation result in theorem 1for neural networks with arbitrary width and height. We are interested in whether there is a similar result for our GRBMs. Intuitively, we could have a positive result given that the marginal distribution of the GRBM is given by a Gaussian mixture, and Gaussian Mixtures Models are universal approximators, as mentioned in page 65 of [38]. That being said, the GRBM marginal is unlike the general Gaussian mixture since in our case we have a mixture of 2mindependent Gaussians, and the means of each component of the mixture are structured in a parallelepiped, as explained in 5.3. For this reason, as already mentioned in the definition, modelling with these GRBMs has its complications due to the restrictions in the model structure [32]. In fact, the universal approximation of GRBMs is still an open problem, but when replacing the sigmoid activation with a ReLU, it was shown in 2022 [39] that the model becomes a universal approximator. On top of that, the paper also showed that the stacking of two latent layers as if increasing the depth of a neural network, which is often termed the Gaussian-Bernoulli Deep Belief Network or GB-DBN, increases the representational power with respect to our GRBM enough to be able to prove a universal approximation result. A depiction of this extension with three hidden layers and directed conditionals, with the undirected RBM component on top, is provided in figure 6.2 Figure 6: Illustration of the GB-DBN extension of the GRBM, from [39] 2In case one wishes to look further, this type of graph with mixed edges is known as a chain graph. 47 For results on the approximation rate, from the machine learning perspective we have the same article from 2022 [39] providing the convergence rate of the GB-DBN model. In this extension of the model, they provide the result of O(ϵ−2)hidden units per layer for a double latent layer, meaning that, module a constant value, we would need to quadruple the the number of hidden units in each of the two layers to half the approximation error. 5.4.2 Model Cumulants Now, moving onto the second question, we shall be looking at the possible moments that we can attain with our model. That being said, we wont be looking exactly at the moments but at the cumulants of our model, defined as follows. Definition 5.1. The nth-order cumulant κnof a random variable Xis given by the nth term of the power series expansion of the cumulant generating function, which is defined in terms of the moment-generating function as KX(t) = log MX(t) = log E[etX ](5.7) Thus, the nth order cumulant κnis given by Kn(0) in K(t) = X n=1 κn n!=κ1t+κ2 t2 2+κ3 t3 6+···+κk tk k!+. . . (5.8) Remark 5.1. The cumulant generating function of a joint distribution X= (X1, . . . , Xk)is given by KX(t1, . . . , tk) = log E[exp X j tjXj!] The reason for us to prefer working with cumulants instead of moments is due them having some very desirable statistical and combinatorial properties. Proposition 5.2. Some basic properties of cumulants are as follows: (i) For n > 1and constant c,κn(c+X) = κn(X)(translation invariant). (ii) For constant c,κn(cX) = cnκ(X)(homogeneous of degree n). (iii) For X1, . . . , Xmindependent r.v.s we have κn(X1+···+Xm) = κn(X1) + ···+κ(Xm) The first and second order cumulant of Xare, conveniently, the expectation E[X]and the covariance Cov(X, X). On top of that, we had Marcinkiewicz’s result from 5, concerning the 48 cumulants of the Gaussian distribution. Another strong reason to prefer cumulants is a somewhat obscure result by Brillinger generalising the law of total covariance known as the Law of Total Cumulance [40]: Theorem 6. (Law of total cumulance) κ(X1, . . . , Xk) = X π (κ(κ(Xi:i∈B|Y) : B∈π)) Where Yis the random variable we are conditioning on and πare partition blocks of {1, . . . , k}. Using this result, we may easily find the higher order cumulants of any set of observables (without repetitions) from our model thanks to the independence between them when conditioning on the latent variables. For instance, consider the covariance between any two variables Yi, Yj, i =j from our model: Cov(Yi, Yj) = E(Cov(Yi, Yj|H)) + Cov (E(Yi|H),E(Yj|H)) =(((((((( ( E(Cov(Yi, Yj|H)) + Cov (E(Yi|H),E(Yj|H)) =E(E(Yi|H)E(Yj|H)) −E(E(Yi|H)) E(E(Yj|H)) =E" X h µ(i) h 1 {H=h}! X h µ(j) h 1 {H=h}!#−E(Yi)E(Yj) (5.9) Where we used µ(i) hfor the conditional mean of Yigiven the latent state H=h, and 1 Afor an indicator function of the event A. We can see from the previous expression how we could characterise the covariance of our model in terms of a multinomial product For a general n-th order cumulant we would have: κ(Yi1, . . . , Yir) = κ(κ(Yi1|H), . . . , κ(Yir|H)) E[Yi|H] = X h µ(i) h 1 {H=h}(5.10) Unlike the tree case where we had the simple linear expressions for the expectation and the covariance, as in equations 5.4 and 5.5, now these will not be linear. That being said, there is still little to say about the expectation of the model, as we can trivially attain any expectation we wish by moving the aparameters for the H=0, and the others will be given by the edges of the parallelepiped given by Wh. With the law of total expectation, we may describe the mean of any observable Yias: E[Yi] = E[E[Yi|H]] = X h µ(i) hP(H=h)(5.11) 49 As the combination of the means for each of the 2mlatent states. Concerning the second cumulant, the covariance, things were not so trivial. As with the tree case, we will want to look at the normalised covariance (correlation) to make things more visual. Recall that for three variables, the correlation matrix is given by: R=   1x y x1z y z 1    Whose determinant is therefore det(R) = 1 −x2−y2−z2+ 2xyz and to the set of points (x, y, z)where this is positive, we may call it the correlation variety. We saw that with the tree model there were many different correlations that weren’t attainable, but with the GRBM this is no longer the case. In figure 7we have the results of some simulations using uniformly generated parameters, with one of the three correlation values on each axis, thanks to the fact that we can use Gibbs Sampling to generate samples of our distribution, where we increased the number of hidden variables. Even though the images represent a noisy version of the actual correlation space, we can clearly appreciate how increasing the number of hidden variables increases the space of correlations attainable by our model, increasing the volume spanned by the samples correlations produced, and we conjecture that the space described by our model can be made to be uniformly and arbitrarily close to that of the equation det(R)=1−x2−y2−z2+ 2xyz, |x|,|y|,|z| ≤ 1. For the code used to generate the samples in Python3, check out this repository. 50 (a) Tree Case (Single H) (b) |H|= 10 (c) |H|= 50 (d) |H|= 100 Figure 7: Simulated correlations of a GRBM model with 3 observable variables and different numbers of latent variables Moving on to higher order cumulants, the parametrisation in 5.11 we found to be rather cumbersome, and so there was a tensor notation that we found was more illuminating; in many instances it is more useful to talk about the means with another tensor basis, noting that H= (H1, . . . , Hm) can be re-written using Hi= (1, Hi), and then H=H1⊗H2⊗···⊗Hm:R2×···×2−→ R And using this notation we can think the conditional mean as consisting of a base term corresponding to H1=··· =Hm= 0, and some offsets that are added when certain latent variables are "active" giving us a tensor ai1i2···im, ij∈0,1: A(i)∈R2×···×2, A(i)⋆ H =E[Yi|H] 51 Where we define the ⋆operation on tensors so as to combine them to give what we describe in the following expression of our conditional expectation: E[Yi|H] = a0...0+a10...0H1+. . . +a0...01Hm+···+a1...1H1···Hm (5.12) Unlike the tree case, where the conditional mean was a linear function of the latent variables, here we have the conditional mean is a polynomial function of the latent variables, and so we need to incorporate higher order tensors. In this way, the arbitrary-order cumulants we had in 5.10 could be algebraically manipulated using the properties of cumulants to get: κ(Yi1, . . . , Yir) = κA(i1)⋆ H, . . . , A(ir)⋆ H= (A(i1), . . . , A(ir))⋆ κ(H, . . . , H) We already saw for the covariance how this ended up being a multinomial product of two terms, and in general we will have a multinomial product of rterms for the r-th order cumulant. And so our problem of calculating these cumulants reduces to calculating multinomial products and the the cumulants of binary random vectors, which to the best of our knowledge is a combinatorial open problem. Therefore, in order to progress in this direction, combinatorial tools must be developed and applied in this realm, and the cumulants of this conditional Gaussian model with fall with them. 52 6 Conclusion As a conclusion, this work has been an overview of probabilistic models with recent developments and applications, how to formalise them as statistical models and, hopefully, an interesting way to motivate why they work the way they do. The main constructions of Restricted Boltzmann Machines and our GRBM model have been provided, and an review of their strengths and weaknesses with respect to concerns of identifiability, parameter estimation, expressibility, and relation to other known models has been covered. The parametrisation provided for the higher order cumulants of the GRBM model serve as a possible starting point for further research, and demonstrate once again the power of tensor algebra to advance our knowledge of statistics. In addition, we gave background on the use of the somewhat mysterious sigmoid function, usually given in machine learning courses as a simple non-linear transformation of the outputs, and how it ties together the conditional density in harmonium structures. There are also several open problems concerning these data modelling challenges which touch upon purer areas of mathematics such as algebraic geometry and especially combinatorics, and further work in this direction would help to better understand the backbone of the machine learning models used in the cutting-edge technologies and state-of-the-art algorithms that dominate the computer vision, natural language processing and reinforcement learning paradigms. Some work has been done on general exponential family distributions over the nodes in a harmonium, but again this work, while recent, mostly comes from the machine learning community and so the results are more practical than theoretical in nature, so this is definitely another possible future extension. What the author also finds really fascinating about this area is how easily you can find yourself down deep rabbit holes that go in completely different directions. For instance, when looking at graphical models one can ask questions about sampling theory and MCMC, or inquire about how the graph properties determining the conditional independence structure can be used to model causality, or even probe into the depths of exponential families and their rich information-theoretic properties. That without even mentioning the vast jungle of numerical methods and operations research algorithms used to train learning algorithms. At every junction there is a possibility for new explorations, and despite the difficulties this posed to the author’s focus, this was an invaluable introduction into such a plethora of topics for further study. 53 References [1] Karl Pearson and Olaus Magnus Friedrich Erdmann Henrici. “III. Contributions to the mathematical theory of evolution”. In: Philosophical Transactions of the Royal Society of London. (A.) 185 (1894), pp. 71–110. doi:10 . 1098 / rsta . 1894 . 0003. eprint: https : //royalsocietypublishing.org/doi/pdf/10. 1098 / rsta . 1894 .0003.url:https: //royalsocietypublishing.org/doi/abs/10.1098/rsta.1894.0003. [2] Carlos Améndola, Marta Casanellas, and Luis David García-Puente. “Tapas of Algebraic Statistics”. In: Notices of the American Mathematical Society (2018). [3] Christopher M. Bishop. Pattern recognition and machine learning. Information Science and Statistics. Springer, New York, 2006. isbn: 0387310738. doi:10.1007/978-0-387-45528-0. [4] R. Fletcher. Practical Methods of Optimization. A Wiley-Interscience publication. Wiley, 2000. isbn: 9780471494638. [5] Jorge Nocedal and Stephen J Wright. Numerical optimization. Springer, 1999. [6] Michael JD Powell. “Radial basis functions for multivariable interpolation: a review.” In: Algorithms for the Approximation of Functions and Data. (1985). [7] Carlos Améndola, Alexander Engström, and Christian Haase. “Maximum number of modes of Gaussian mixtures”. In: Information and Inference: A Journal of the IMA 9.3 (June 2019), pp. 587–600. doi:10.1093/imaiai/iaz013.url:https://doi.org/10.1093/imaiai/ iaz013. [8] Karl Pearson F.R.S. “LIII. On lines and planes of closest fit to systems of points in space”. In: Philosophical Magazine Series 1 2 (1901), pp. 559–572. [9] Michael E. Tipping and Christopher M. Bishop. “Probabilistic Principal Component Analysis”. In: Journal of the Royal Statistical Society. Series B (Statistical Methodology) 61.3 (1999), pp. 611–622. issn: 13697412, 14679868. url:http://www.jstor.org/stable/ 2680726. [10] A.T. Basilevsky. Statistical Factor Analysis and Related Methods: Theory and Applications. Wiley Series in Probability and Statistics. Wiley, 1994. isbn: 9780471570820. [11] Geoffrey E Hinton. “Training products of experts by minimizing contrastive divergence”. In: Neural computation 14.8 (2002), pp. 1771–1800. [12] Iqbal Sarker. “Deep Learning: A Comprehensive Overview on Techniques, Taxonomy, Applications and Research Directions”. In: SN Computer Science 2 (Aug. 2021). doi:10.1007/ s42979-021-00815-1. 54 [13] Kai Fong Ernest Chong. “A closer look at the approximation capabilities of neural networks”. In: International Conference on Learning Representations. 2020. [14] Kurt Hornik, Maxwell B. Stinchcombe, and Halbert L. White. “Multilayer feedforward networks are universal approximators”. In: Neural Networks 2 (1989), pp. 359–366. [15] Moshe Leshno et al. “Multilayer feedforward networks with a nonpolynomial activation function can approximate any function”. In: Neural Networks 6.6 (1993), pp. 861–867. issn: 0893-6080. doi:https://doi.org/10.1016/S0893-6080(05)80131-5. [16] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton. “ImageNet Classification with Deep Convolutional Neural Networks”. In: Advances in Neural Information Processing Systems. Ed. by F. Pereira et al. Vol. 25. Curran Associates, Inc., 2012. [17] Yann LeCun. “Who is afraid of non-convex loss functions?” NIPS. 2007. url:https://cs. nyu.edu/~yann/talks/lecun-20071207-nonconvex.pdf. [18] Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Elements of Causal Inference: Foundations and Learning Algorithms. Adaptive Computation and Machine Learning. Cambridge, MA: MIT Press, 2017. isbn: 978-0-262-03731-0. url:https://mitpress.mit.edu/ books/elements-causal-inference. [19] A. P. Dempster, N. M. Laird, and D. B. Rubin. “Maximum Likelihood from Incomplete Data via the EM Algorithm”. In: Journal of the Royal Statistical Society. Series B (Methodological) 39.1 (1977), pp. 1–38. issn: 00359246. [20] S.L. Lauritzen. Graphical Models. Oxford Statistical Science Series. Clarendon Press, 1996. isbn: 9780191591228. [21] R.J. Baxter. Exactly Solved Models in Statistical Mechanics. Dover books on physics. Dover Publications, 2007. isbn: 9780486462714. [22] Paul Smolensky. “Information processing in dynamical systems: Foundations of harmony theory”. In: Parallel Distributed Process 1 (Jan. 1986). [23] Yoav Freund and David Haussler. “Unsupervised Learning of Distributions of Binary Vectors Using 2-Layer Networks”. In: NIPS. 1991. [24] Miguel Á. Carreira-Perpiñán and Geoffrey Hinton. “On Contrastive Divergence Learning”. In: Proceedings of the Tenth International Workshop on Artificial Intelligence and Statistics. Ed. by Robert G. Cowell and Zoubin Ghahramani. Vol. R5. Proceedings of Machine Learning Research. Reissued by PMLR on 30 March 2021. PMLR, Jan. 2005, pp. 33–40. 55 [25] Guido Montúfar. “Restricted Boltzmann Machines: Introduction and Review”. In: ArXiv abs/1806.07066 (2018). [26] S. L. Lauritzen and N. Wermuth. “Graphical Models for Associations between Variables, some of which are Qualitative and some Quantitative”. In: The Annals of Statistics 17 (1989), pp. 31–57. doi:10 . 1214 / aos / 1176347003.url:https : / / doi . org / 10 . 1214 / aos / 1176347003. [27] Mathias Drton, Bernd Sturmfels, and Seth Sullivant. Lectures on Algebraic Statistics. Vol. 39. Oberwolfach Seminars. Springer, 2009. doi:10.1007/978-3-7643-8905-5. [28] Persi Diaconis and Bernd Sturmfels. “Algebraic algorithms for sampling from conditional distributions”. In: The Annals of Statistics 26.1 (1998), pp. 363–397. doi:10.1214/aos/ 1030563990.url:https://doi.org/10.1214/aos/1030563990. [29] Anna Seigal and Guido Montufar. “Mixtures and products in two graphical models”. In: (2017). doi:10.48550/ARXIV.1709.05276.url:https://arxiv.org/abs/1709.05276. [30] Alex Krizhevsky, Geoffrey Hinton, et al. “Learning multiple layers of features from tiny images”. In: (2009). [31] Jan Melchior, Nan Wang, and Laurenz Wiskott. “Gaussian-binary restricted Boltzmann machines for modeling natural image statistics”. In: PLOS ONE 12.2 (Feb. 2017), pp. 1–24. doi:10.1371/journal.pone.0171015. [32] Oswin Krause et al. “Approximation properties of DBNs with binary hidden units and realvalued visible units”. In: Proceedings of the 30th International Conference on Machine Learning. Ed. by Sanjoy Dasgupta and David McAllester. Vol. 28. Proceedings of Machine Learning Research 1. Atlanta, Georgia, USA: PMLR, June 2013, pp. 419–426. [33] Piotr Zwiernik. “Latent tree models”. In: Handbook of graphical models. Chapman & Hall/CRC Handbooks of Modern Statistical Methods. CRC Press, 2019, pp. 265–288. url:https: //books.google.ca/books?id=g7J5DwAAQBAJ. [34] M. Maathuis et al. Handbook of Graphical Models. Chapman & Hall/CRC Handbooks of Modern Statistical Methods. CRC Press, 2018. isbn: 9780429874239. [35] Paul Erdos and Alfréd Rényi. “Asymmetric graphs”. In: Acta Math. Acad. Sci. Hungar 14.295-315 (1963), p. 3. [36] “Full reconstruction of Markov models on evolutionary trees: Identifiability and consistency”. In: Mathematical Biosciences 137.1 (1996), pp. 51–73. issn: 0025-5564. doi:https://doi. org/10.1016/S0025-5564(96)00075-2. 56