Hierarchical Multinomial-Dirichlet model for the estimation of conditional probability tables

I Introduction

A Bayesian network is a probabilistic model constituted by a directed acyclic graph (DAG) and a set of conditional probability tables (CPTs), one for each node. The CPT of node $X$

contains the conditional probability distributions of

X

given each possible configuration of its parents. Usually all variables are discrete and the conditional distributions are estimated adopting a Multinomial-Dirichlet model, where the Dirichlet prior is characterised by the vector of hyper-parameters

α

. Yet, Bayesian estimation of multinomials is sensitive to the choice of

α

and inappropriate values cause the estimator to perform poorly [hausser2009entropy]. Mixtures of Dirichlet distributions have been recommended both in statistics [casella2009assessing, good1987robustness]

and in machine learning

[nemenman2002entropy] in order to obtain more robust estimates. Yet, mixtures of Dirichlet distributions are computationally expensive; this prevents them from being widely adopted. Another difficulty encountered in CPT estimation is the presence of rare events. Assuming that all variables have cardinality

k

and that the number of parents is

q

, we need to estimate

k^{q}

conditional distributions, one for each joint configuration of the parent variables. Frequently one or more of such configurations are rarely observed in the data, making their estimation challenging.

We propose to estimate the conditional distributions by adopting a novel approach, based on a hierarchical Multinomial-Dirichlet model. This model has two main characteristics. First, the prior of each conditional distribution is constituted by a mixture of Dirichlet distributions with parameter $α$ ; the mixing is attained by treating $α$

as a random variable with its own prior and posterior distribution. By estimating from data the posterior distribution of

α

, we need not to fix its value a priori. Instead we give more weight to the values of

α

that are more likely given the data. Secondly, the model is hierarchical since it assumes that the conditional distributions within the same CPT (but referring to different joint configurations of the parents) are drawn from the same mixture. The hierarchical model jointly estimates all the conditional distributions of the same CPT, called columns of the CPT. The joint estimates generate information flow between different columns of the CPT; thus the hierarchical model exploits the parameters learned for data-rich columns to improve the estimates of the parameters of data-poor columns. This is called borrowing statistical strength [murphy2012machine, Sec 6.3.3.2] and it is well-known within the literature of hierarchical models [gelman2013bayesian]. Also the literature of Bayesian networks acknowledges [koller2009probabilistic, Sec.17.5.4] that introducing dependencies between columns of the same CPT could improve the estimates, especially when dealing with sparse data. However, as far as we known, this work is the first practical application of joint estimation of the columns of the same CPT.

To tackle the problem of computational complexity we adopt a variational inference approach. Namely, we compute a factorised approximation of the posterior distribution that is highly efficient. Variational inference appears particularly well suited for hierarchical models; for instance the inference of Latent Dirichlet Allocation [LDA] is based on variational Bayes. By extensive experiments, we show that our novel approach considerably improves parameter estimation compared to the traditional approaches based on Multinomial-Dirichlet model. The experiments show large gains especially when dealing with small samples, while with large samples the effect of the prior vanishes as expected.

The paper is organised as follows. Section II introduces the novel hierarchical model. Section III provides an analytical study of the resulting estimator, proving that the novel hierarchical approach provides lower estimation error than the traditional approaches, under some mild assumptions on the generative model. Section IV presents some simulation studies showing that, given the same network structure, hierarchical estimation yields both a better fit of the joint distribution and a consistent improvement in classification performance, with respect to the traditional estimation under parameter independence. Section V reports some concluding remarks.

Ii Estimation Under Multinomial-Dirichlet Model

We want to induce a Bayesian network over the set of random variables $X = {X_{1}, \dots, X_{I}}$ . We assume that each variable $X_{i} \in X$ is discrete and has $r_{i}$ possible values in the set $X_{i}$ . The parents of $X_{i}$ are denoted by ${P a}_{i}$ and they have $q_{i}$ possible joint states collected in the set $P_{i}$ .
We denote by $θ_{x | p a}$ the probability of $X_{i}$ being in state $x \in X_{i}$ when its parent set is in state $p a \in P_{i}$ , i.e., $θ_{x | p a} = p (X_{i} = x | {P a}_{i} = p a) > 0$ . We denote by $θ_{X_{i} | p a}$ the parameters of the conditional distribution of $X_{i}$ given ${P a}_{i} = p a$ . A common assumption [koller2009probabilistic, Sec.17] is that $θ_{X_{i} | p a}$ is generated from a Dirichlet distribution with known parameters. The collection of the conditional probability distributions $θ_{X_{i}} = (θ_{X_{i} | {p a}_{1}}, \dots, θ_{X_{i} | {p a}_{q_{i}}})$ constitutes the conditional probability table (CPT) of $X_{i}$ . Each vector of type $θ_{X_{i} | p a}$ , with $p a \in P_{i}$ , is a column of the CPT.

The assumption of local parameter independence [koller2009probabilistic, Sec. 17] allows to estimate each parameter vector $θ_{X_{i} | p a}$ independently of the other parameter vectors. The assumed generative model, $\forall i \in 1, \dots, I$ , is:

	$p (θ_{X_{i} \| p a})$	$= D i r (s α)$	$p a \in P_{i},$
	$p (X_{i} \| {P a}_{i} = p a, θ_{X_{i} \| p a})$	$= C a t (θ_{X_{i} \| p a})$	$p a \in P_{i},$

where $s \in R$ denotes the prior strength, also called equivalent sample size, and $α \in R^{r_{i}}$ is a parameter vector such that $\sum_{x \in X} α_{x} = 1$ . The most common choice is to set $α_{x} = 1 / r_{i}$ and $s = 1 / q_{i}$ , which is called BDeu prior [koller2009probabilistic, Sec.17].

If there are no missing values in the data set $D$ , the posterior expected value of $θ_{x_{|} p a}$ is [gelman2013bayesian]:

E [θ_{x_{|} p a}] = \frac{n_{x, p a} + s α_{x}}{n_{p a} + s},

where $n_{x, p a}$ is the number of observations in $D$ characterised by $X_{i} = x$ and ${P a}_{i} = p a$ , while $n_{p a} = \sum_{x \in X_{i}} n_{x, p a}$ .

Iii Hierarchical Model

The proposed hierarchical model estimates the conditional probability tables by removing the local independence assumption. In order to simplify the notation we present the model on a node $X$ with a single parent $Y$ . $X$ has $r$ states in the set $X$ , while $Y$ has $q$ states in the set $Y$ . Lastly, we denote by $n_{x y}$ the number of observations with $X = x$ and $Y = y$ and by $n_{y} = \sum_{x \in X} n_{x y}$ the number of observations with $Y = y$ , where $x \in X$ and $y \in Y$ .

As described in Section II, $θ_{X_{i} | p a}$ , for $i = 1, \dots, I$ and $p a \in P_{i}$ , are usually assumed to be independently drawn from a Dirichlet distribution, with known parameter $α$ . On the contrary, the hierarchical model treats $α$ as a hidden random vector, thus making different columns of the CPT dependent. Specifically, we assume $α$ to be drawn from a higher-level Dirichlet distribution with hyper-parameter $α_{0}$ .

We assume that $(x_{k}, y_{k})$ for $k = 1, \dots, n$ are $n$ $i . i . d .$ observations from the hierarchical Multinomial-Dirichlet model:

$p (α \| α_{0})$	$= Dirichlet (α_{0})$
$p (θ_{X \| y} \| s, α)$	$= Dirichlet (s α)$	$y \in Y$	(1)
$p (X \| Y = y, θ_{X \| y})$	$= Cat (θ_{X \| y})$	$y \in Y$

where $s \in R$ is the equivalent sample size, and $α_{0} \in R^{r}$ is a vector of hyper-parameters.

Iii-a Posterior moments for $θ_{X | y}$

We now study the hierarchical model, deriving an analytical expression for the posterior average of $θ_{x | y}$ , which is the element $x$ of vector $θ_{X | y}$ , and for the posterior covariance between $θ_{x | y}$ and $θ_{x^{'} | y^{'}}$ . To keep notation simple, in the following we will not write explicitly the conditioning with respect to the fixed parameters $s$ and $α_{0}$ . We introduce the notation $E^{D} [\cdot] = E [\cdot | D]$ to represent the posterior average and ${C ov}^{D} (\cdot, \cdot) = C ov (\cdot, \cdot | D)$ to represent the posterior covariance.

Definition 1.

We define the pointwise estimator ${^θ}_{x | y}$ for the parameter $θ_{x | y}$ as its posterior average, i.e., ${^θ}_{x | y} = E^{D} [θ_{x | y}]$ , and the pointwise estimator ${^α}_{x}$ for the element $x$ of the parameter vector $α$ as its posterior average, i.e., ${^α}_{x} = E^{D} [α_{x}] .$

Theorem 1.

Under model (1), the posterior average of $θ_{x | y}$ is

{^θ}_{x | y} = E^{D} [θ_{x | y}] = \frac{n_{x y} + s {^α}_{x}}{n_{y} + s},

(2)

while the posterior covariance between $θ_{x | y}$ and $θ_{x^{'} | y^{'}}$ is

{C ov}^{D} (θ_{x | y}, θ_{x^{'} | y^{'}}) = δ_{y y^{'}} \frac{{^θ}_{x | y} δ_{x x^{'}} - {^θ}_{x | y} {^θ}_{x^{'} | y}}{n_{y} + s + 1} + \frac{s^{2} {C ov}^{D} (α_{x}, α_{x^{'}})}{C_{y y^{'}}},

where $C_{y y^{'}}$ is defined as

C_{y y^{'}} = {\begin{matrix} (n_{y} + s) (n_{y^{'}} + s) & if% y \neq y^{'} (n_{y} + s) (n_{y} + s + 1) & if y = y^{'} . \end{matrix}

The posterior average and posterior covariance of $α$ cannot be computed analytically. Some results concerning their expression and numerical computation, together with the complete proof of Theorem 1, are detailed in Appendix.

Notice that the pointwise estimator ${^θ}_{x | y}$ is a mixture of traditional Bayesian estimators obtained under (non-hierarchical) Multinomial-Dirichlet models with $α$ fixed, i.e, $\frac{n_{x y} + s α_{x}}{n_{y} + s}$ . Indeed, thanks to the linearity in $α_{x}$ , we obtain

{^θ}_{x | y} = \frac{n_{x y} + s {^α}_{x}}{n_{y} + s} = \int \frac{n_{x y} + s α_{x}}{n_{y} + s} p (α | D) d α .

This mixture gives more weight to the values of $α$ that are more likely given the observations.

Iii-B Properties of the estimator ${^θ}_{X | y}$

We study now the mean-squared error (MSE) of ${^θ}_{x | y}$ and we compare it to the MSE of other traditional estimators.
In order to study the MSE of ${^θ}_{x | y}$ we need to assume the generative model

	$p (θ_{X \| y} \| s, ~ α)$	$= Dirichlet (s ~ α)$	$y \in Y,$
	$p (X \| Y = y, θ_{X \| y})$	$= Cat (θ_{X \| y})$	$y \in Y,$		(3)

where $s$ and $~ α$ are the true underlying parameters. Moreover, since $θ_{X | y}$ is a random vector, we define the MSE for an estimator ${¯ θ}_{x | y}$ of the single component $θ_{x | y}$ as

MSE ({¯ θ}_{x | y}) = E_{θ} [E_{n} [{({¯ θ}_{x | y} - θ_{x | y})}_{x | y}^{2}]],

(4)

where $E_{θ} [\cdot]$ and $E_{n} [\cdot]$ represent respectively the expected value with respect to $θ_{x | y}$ and $n_{x y}$ , and the MSE for the estimator ${¯ θ}_{X | y}$ of the vector $θ_{X | y}$ as $MSE ({¯ θ}_{X | y}) = \sum_{x \in X} MSE ({¯ θ}_{x | y}) .$

Notice that the generative model (3) is the traditional, thus non-hierarchical, Multinomial-Dirichlet model, which implies parameter independence. Hence, the traditional Bayesian estimator satisfies exactly the assumptions of this model. The Bayesian estimator is usually adopted by assuming $~ α$

to be fixed to the values of a uniform distribution on

X

, i.e.,

α^{B} = 1 / r \cdot 1_{1 \times r}

, see e.g. [gelman2013bayesian]. However, since in general

~ α \neq 1 / r \cdot 1_{1 \times r}

, the traditional Bayesian approach generates biased estimates in small samples. On the contrary, the novel hierarchical approach estimates the unknown parameter vector

~ α

basing on its posterior distribution. For this reason the proposed approach can provide estimates that are closer to the true underlying parameters, with a particular advantage in small samples, with respect to other traditional approaches.

In order to study the MSE of different estimators, we first consider an ideal shrinkage estimator

θ_{X | y}^{*} = {~ ω}_{y} θ_{X | y}^{ML} + (1 - {~ ω}_{y}) ~ α,

(5)

where ${~ ω}_{y} \in (0, 1)$ and $θ_{x | y}^{ML} = \frac{n_{x y}}{n_{y}}$ is the maximum-likelihood (ML) estimator, obtained estimating from the observations each vector $θ_{X | y}$ independently of other vectors. This convex combination shrinks the ML estimator towards the true underlying parameter $~ α$ . Setting ${~ ω}_{y} = \frac{n_{y}}{n_{y} + s}$ , the ideal estimator corresponds to a Bayesian estimator with known parameter $~ α$ , i.e., $θ_{x | y}^{*} = \frac{n_{x y} + {~ α}_{x}}{n_{y} + s}$ . However, since $~ α$ represents the true underlying parameter that is usually unknown, the ideal estimator (5) is unfeasible. Yet it is useful as it allows to study the MSE.

The main result concerns the comparison, in terms of MSE, of the ideal estimator with respect to the traditional (non hierarchical) Bayesian estimator, which estimates ${~ α}_{x}$ by means of the uniform distribution $1 / r$ , i.e., $θ_{x | y}^{B} = \frac{n_{x y} + s / r}{n_{y} + s}$ , see e.g. [gelman2013bayesian]. The traditional Bayesian estimator can be written as

θ_{X | y}^{B} = ω_{y} θ_{X | y}^{ML} + (1 - ω_{y}) \frac{1}{r},

(6)

where $ω_{y} = \frac{n_{y}}{n_{y} + s}$ .

Theorem 2.

Under the assumption that the true generative model is (3), the MSE for the ideal estimator is

MSE (θ_{x | y}^{*}) = ({~ ω}_{y}^{2} \frac{s}{} n_{y} + (1 - {~ ω}_{y})^{2}) \frac{{~ α}_{x} - {~ α}_{x}^{2}}{s + 1},

while the MSE for the traditional Bayesian estimator is

MSE (θ_{x | y}^{B}) = MSE (θ_{x | y}^{*}) + (1 - ω_{y}^{2}) {({~ α}_{x} - \frac{1}{r})}_{x}^{2} .

If ${~ ω}_{y} = ω_{y}$ , $MSE (θ_{X | y}^{*}) \leq MSE (θ_{X | y}^{B})$ .

The proof is reported in Appendix.

Since in general ${~ α}_{x} \neq \frac{1}{r}$ , the second term in (2) is positive and the ideal estimator achieves smaller MSE with respect to traditional Bayesian estimator. To improve the estimates of the traditional Bayesian model in terms of MSE, we propose to act exactly on the second term of (2). Specifically, we can achieve this purpose by estimating the parameter vector $~ α$ from data, instead of considering it fixed.

The proposed hierarchical estimator defined in (2) has the same structure of the ideal estimator (5):

{^θ}_{X | y} = ω_{y} θ_{X | y}^{ML} + (1 - ω_{y})^α,

where $ω_{y} = \frac{n_{y}}{n_{y} + s}$ . This convex combination of $θ_{X | y}^{ML}$ and $^α$ shrinks the ML estimator towards the posterior average of $α$ , with a strength that is inversely proportional to $n_{y}$ .

Contrary to the traditional Bayesian estimator (6), the hierarchical estimator provides an estimate of $~ α$ that converges to the true underlying parameter as $n$ increases. Indeed, it is well known that the posterior average $E^{D} [α]$ converges to the true underlying parameter $~ α$ as $n$ goes to infinity, i.e., $^α \to ~ α$ as $n \to + \infty$ , [gelman2013bayesian]. As a consequence, ${^θ}_{x | y}$ converges to $θ_{x | y}^{*}$ and the MSE of ${^θ}_{x | y}$ converges to the MSE of $θ_{x | y}^{*}$ , i.e., $MSE ({^θ}_{x | y}) \to MSE (θ_{x | y}^{*})$ as $n \to + \infty$ .

In the finite sample assumption $MSE ({^θ}_{x | y})$ differs from $MSE (θ_{x | y}^{*})$ , since ${^θ}_{x | y}$ includes an estimator of $α$ . Since in the hierarchical model we cannot compute this quantity analytically, we verify by simulation that the hierarchical estimator provides good performances in terms of MSE with respect to the traditional Bayesian estimators.

In conclusion, as we will show in the numerical experiments, the hierarchical model can achieve a smaller MSE than the traditional Bayesian estimators, in spite of the unfavourable conditions of the generative model (3). This gain is obtained thanks to the estimation of the parameter $~ α$ , rather than considering it fixed as in the traditional approaches. In more general conditions with respect to (3), the true generating model could not satisfy parameter independence and the MSE gain of the hierarchical approach would further increase.

Iv Experiments

In the experiments we compute the proposed hierarchical estimator using variational Bayes inference in R by means of the rstan package [stan2]

. The variational Bayes estimates are practically equivalent to those yielded by Markov Chain Monte Carlo (MCMC), though being the variational inference much more efficient than MCMC (less than a second for estimating a CPT compared to a couple of minutes for the MCMC method). In the following we report the results obtained via variational inference. The code is available at

http://ipg.idsia.ch/software.php?id=139.

Fig. 1: Boxplots of MSE difference between the Bayesian ( $s = r$ ) and the hierarchical estimator in test 1 (light blue) and test 2 (green) with different dimension of the conditioning set ( $q = 2, 4, 6, 8$ ). Positive values favour the hierarchical model.

Iv-a MSE analysis

In the first study we assess the performances of the hierarchical estimator in terms of MSE. We consider two different settings, in which we generate observations from model (3), where $~ α$ is fixed. In the first setting (test 1) we sample $~ α$ from a Dirichlet distribution with parameter $1_{1 \times r}$ , while in the second setting (test 2) we sample it from a Dirichlet distribution with parameter $10^{6} \cdot 1_{1 \times r}$ . Under test 2 the parameters of the sampling distribution for $~ α$ are very large and equal to each other, implying ${~ α}_{x} ≃ 1 / r$ , $\forall x \in X$ . For this reason, test 2 is the ideal setting for the traditional Bayesian estimator, while test 1 is the ideal setting for the hierarchical estimator.

In both test 1 and test 2 we consider all the possible combinations of $r$ (the number of states of $X$ ) and $q$ (the number of conditioning states), with $r \in {2, 4, 6, 8}$ and $q \in {2, 4, 6, 8}$ . For each combination of $r$ and $q$ , for both test 1 and test 2, we generate data sets with size $n \in {20, 40, 80, 160, 320, 640}$ . We repeat the data sampling and the estimation procedure 10 times for each combination of $r$ , $q$ and $n$ . Then we compare the estimates yielded by the hierarchical estimator, with $s = r$ and $α_{0} = 1_{1 \times r}$ , and by the traditional Bayesian estimator assuming parameter independence, with $s = r$ . We compare the performance of different estimators by computing the difference in terms of average MSE, defined as $MSE ({¯ θ}_{X | Y}) = \sum_{y \in Y} \sum_{x \in X} MSE ({¯ θ}_{x | y}) / r q$ for an estimator ${¯ θ}_{X | Y}$ . In every repetition of the experiment ${1, \dots, 10}$ , we then compute $MSE (θ_{X | Y}^{B}) - MSE ({^θ}_{X | Y})$ and we represent it in Figure 1.

The results show that the hierarchical estimator mostly provides better or equivalent results in comparison to $θ_{X | Y}^{B}$ , especially for small $n$ and/or large $q$ . In a Bayesian network it is usual to have large values of $q$ , since $q$ represents the cardinality of the parents’ joint states set. In test 1 (light blue boxplots) the advantage of the hierarchical estimator over the Bayesian one is generally large, as expected. The advantage of the hierarchical model steadily increases as $q$ increases, becoming relevant for $q = 6$ or $q = 8$ . For large $n$ the gap between the two estimators vanishes, although it is more persistent when dealing with large $q$ . Interestingly, in test 2 (green boxplots), the traditional Bayesian estimator is just slightly better than the hierarchical one, even though the former is derived exactly from the true generative model. The traditional estimator has a small advantage only for $q = 4$ and small values of $n$ , and this advantage quickly decreases if either $q$ or $n$ increase.

Iv-B Joint distribution fitting

In the second study we assess the performance of the hierarchical estimator in the recovery of the joint distribution of a given Bayesian network.

We consider 5 data sets from UCI Machine Learning Repository: Adult, Letter Recognition, Nursery, Pen-Based Recognition of Handwritten Digits and Spambase. We discretise all numerical variables into five equal-frequency bins and we consider only instances without missing values. For each dataset we first learn, from all the available data, the associated directed acyclic graph (DAG) by means of a hill-climbing greedy search, as implemented in the bnlearn package [scutari2010learning]. We then keep such structure as fixed for all the experiments referring to the same data set, since our focus is not structural learning. Then, for each data set and for each $n \in {20, 40, 80, 160,$ $320, 640, 1280}$ we repeat 10 times the procedure of 1) sampling $n$ observations from the data set and 2) estimating the CPTs. We perform estimation using the proposed hierarchical approach, with $s = r$ and $α_{0} = 1_{1 \times r}$ , and the traditional BDeu prior (Bayesian estimation under parameter independence) with $s = 1$ and $s = 10$ . The choice of $s = 1$ is the most commonly adopted in practice, while $s = 10$ is the default value proposed by the bnlearn package. Conversely, we did not offer any choice to the hierarchical model. Indeed, we set the smoothing factor $s$ in the proposed model to the number of states of the child variable, which has the same order of magnitude of the smoothing factors used in the traditional Bayesian approach. In spite of the more limited choice for the parameter $s$ , the hierarchical estimator consistently outperforms the traditional Bayesian estimator, regardless whether the latter adopts a smoothing factor 1 or 10.

We then measure the log-likelihood of all the instances included in the test set, where the test contains all the instances not present in training set. We report in the top panels of Figure 2 the difference between the log-likelihood of the hierarchical approach and the log-likelihood of Bayesian estimation under local parameter independence, i.e., the log-likelihood ratio. The log-likelihood ratio approximates the ratio between the Kullback-Leibler (KL) divergences of the two models. The KL of a given model measures the distance between the estimated and the true underlying joint distribution.

The log-likelihood ratios obtained in the experiments are extremely large on small sample sizes, being larger than one thousand on all data sets (Figure 2, top panels). This shows the huge gain delivered by the hierarchical approach when dealing with small data sets. This happens regardless of the equivalent sample size used to perform Bayesian estimation under parameter independence: we note however that in general $s = 10$ yields better results than $s = 1$ . The lowest gains are obtained on the data set Nursery; the reason is that the DAG of Nursery has the lowest number of parents per node (0.9 on average, compared to about twice as much for the other DAGs). Thus, this data set is the less challenging from the parameter estimation viewpoint, with respect to the others. We point out that significant likelihood ratios are obtained in general even for large samples, even though they are not apparent from the figure due to the scale. For instance, for $n = 320$ the log-likelihood ratios range from 50 (Nursery) to 85000 (Letter).

Fig. 2: Boxplots of the logarithm of the likelihood ratio (top panels), accuracy gain (central panels) and area under the ROC gain (bottom panels) obtained comparing the hierarchical method with respect to BDeu ( $s = 1$ in orange and $s = 10$ in red) for the five machine learning datasets analysed. Positive values favour the hierarchical model.

Iv-C Classification

In the third study we assess the performance of the hierarchical estimator in terms of classification. We consider the same datasets of the previous experiment, discretised in the same way.

For each dataset we first learn the Tree-Augmented Naive Bayes (TAN) structure by means of the

bnlearn package. The networks are estimated on the basis of all the available samples and are kept fixed for all the experiments referring to the same data set, since our focus is not structural learning. Then, for each dataset and for each

n \in {20, 40, 80, 160,

320, 640, 1280}

, we sample

n

observations. We then estimate the CPTs of the Bayesian network from the sampled data by means of the hierarchical estimator (

s = r

and

α_{0} = 1_{1 \times r}

) and the traditional Bayesian estimators obtained under a BDeu prior (

s = 1

and

s = 10

). We repeat the sampling and the estimating steps 10 times for each value of

n

and each data set. We then classify each instance of the test set, which contains 1000 instances sampled uniformly from all the instances not included in the training set. We assess the classification performance by measuring accuracy and area under the ROC (ROC AUC) of the classifiers. In the central panels of Figure

2 we report the difference in accuracy between the hierarchical estimator and the traditional Bayesian ones, while in the bottom panels of the same figure we report the difference in ROC AUC between the same classifiers.

The area under the ROC is a more sensitive indicator for the correctness of the estimated posterior probabilities with respect to accuracy. According to Figure

2 (bottom panels), the hierarchical approach yields consistently higher ROC AUC than both the BDeu classifiers. The increase of ROC AUC in small samples (

n = 20

n = 40

) ranges between 2 and 20 points compared to both the BDeu priors. As

n

increases this gain in ROC AUC tends to vanish. However, for

n > 320

the gain in ROC AUC for the datasets Adult and Letter ranges between 1 and 5 points.

Figure 2 (central panels) shows also improvements in accuracy, even if this indicator is less sensitive to the estimated posterior probability than the area under ROC. Indeed, in computing the accuracy, the most probable class is compared to the actual class, without paying further attention to its posterior probability. In small samples ( $n = 20$ , $n = 40$ ) there is an average increase of accuracy of about 5 points compared to the BDeu prior with $s = 10$ and of about 10 points compared to the BDeu prior with $s = 1$ . The accuracy improvements tends to decrease as $n$ increases; yet on both Adult and Letter data sets an accuracy improvement of about 1-2 points is shown also for $n = 1280$ .

V Conclusions

We have presented a novel approach for estimating the conditional probability tables by relaxing the local independence assumption. Given the same network structure, the novel approach yields a consistently better fit to the joint distribution than the traditional Bayesian estimation under parameter independence; it also improves classification performance. Moreover, the introduction of variational inference makes the proposed method competitive in terms of computational cost with respect to the traditional Bayesian estimation.

In order to prove Theorem 1

, we first need to derive some results concerning the posterior moments for the vector

α

, whose general element

α_{x}

is associated to state

x \in X

for the variable

X

. Given a dataset

D

, the

k

-th posterior moment for the element

α_{x}

E [α_{x}^{k} ∣ ∣ D] = \int α_{x}^{k} p (α | D) d α .

The following proposition states a general result for computing any posterior moment of $α$ , whose general expression is $E [\prod_{x \in X} α_{x}^{k_{x}} | D]$ , where $k_{x} \in N$ represents the power of element $α_{x}$ .

Lemma 1.

Under the assumptions of model (2), the posterior average of the quantity $\prod_{x \in X} α_{x}^{k_{x}}$ , with $k_{x} \in N$ $\forall x \in X$ , is

E^{D} [\prod x \in X α_{x}^{k_{x}}] = γ \int \prod x \in X \prod \begin{matrix} y \in Y s.t. n_{x y} > 0 \end{matrix} n_{x y} \prod ν = 1 (s α_{x} + ν - 1) α_{x}^{{~ k}_{x}} d α,

(7)

where ${~ k}_{x} = α_{0, x} + k_{x} - 1$ and $γ$ is a proportionality constant such that

(8)

The element $x^{'}$ of the posterior average vector $E^{D} [α]$ is

{^α}_{x^{'}} = γ \int \prod x \in X \prod \begin{matrix} y \in Y s.t. n_{x y} > 0 \end{matrix} n_{x y} \prod ν = 1 (s α_{x} + ν - 1) α_{x}^{δ_{x, x^{'}}} d α,

(9)

where $δ_{x, x^{'}}$ is a Kronecker delta.
The element $(x^{'}, x^{''})$ of the posterior covariance matrix ${C ov}^{D} (α)$ is

{C ov}^{D} (α_{x^{'}}, α_{x^{''}}) = E^{D} [α_{x^{'}} α_{x^{''}}] - {^α}_{x^{'}} {^α}_{x^{''}},

(10)

where

Both integrals in (7) and (8) are multiple integrals computed with respect to the $r$ elements of vector $α$ , such that $\sum_{x \in X} α_{x} = 1$ . The space of integration is thus the standard $r$ -simplex.

Proof of Lemma 1.

Under the assumptions of model (2), the joint posterior density of $α, θ_{X | y_{1}}, \dots, θ_{X | y_{q}}$ is

	$p (α, θ_{X \| y_{1}}, \dots, θ_{X \| y_{q}} \| D)$
	$\propto \frac{Γ (s)}{\prod_{x \in X} Γ (s α_{x})} \prod y \in Y \prod x \in X (θ_{x \| y})^{n_{x y} + s α_{x} - 1} α_{x}^{α_{0, x} - 1} .$

Marginalising $p (α, θ_{X | y_{1}}, \dots, θ_{X | y_{q}} | D)$ with respect to $θ_{X | y_{1}}, \dots, θ_{X | y_{q}}$ , we obtain the marginal posterior density for $α$ , i.e.,

(11)

Thanks to the well-known property of the Gamma function

Γ (α + m) = m \prod ν = 1 (α + ν - 1) \cdot Γ (α), % for m \geq 1

we can write the posterior marginal density (11) as

p (α | D) \propto \prod x \in X \prod \begin{matrix} y \in Y s.t. n_{x y} > 0 \end{matrix} n_{x y} \prod ν = 1 (s α_{x} + ν - 1) α_{x}^{α_{0, x} - 1} .

(12)

The proportionality constant of the posterior marginal density is obtained by integrating the right term in (12) with respect to the $r$ elements of $α$ , such that $\sum_{x \in X} α_{x} = 1$ . The resulting proportionality constant is thus

γ = {⎛ ⎜ ⎜ ⎜ ⎝ \int \prod x \in X \prod \begin{matrix} y \in Y s.t. n_{x y} > 0 \end{matrix} n_{x y} \prod ν = 1 (s α_{x} + ν - 1) α_{x}^{α_{0, x} - 1} d α ⎞ ⎟ ⎟ ⎟ ⎠}^{- 1} .

The posterior average for the quantity $\prod_{x \in X} α_{x}^{k_{x}}$ can be derived directly from the posterior marginal density of $α$ as

E^{D} [\prod x \in X α_{x}^{k_{x}}] = γ \int \prod x \in X \prod \begin{matrix} y \in Y s.t. n_{x y} > 0 \end{matrix} n_{x y} \prod ν = 1 (s α_{x} + ν - 1) α_{x}^{{~ k}_{x}} d α,

where ${~ k}_{x} = α_{0, x} + k_{x} - 1$ . In the special case of $α_{0, x} = 1$ , $\forall x \in X$ , we have ${~ k}_{x} = k_{x}$ .

The posterior average for $α_{x^{'}}$ is obtained directly from (7), by choosing $k_{x} = δ_{x = x^{'}}$ , i.e., $k_{x} = 1$ for $x = x^{'}$ and $k_{x} = 0$ for $\forall x \neq x^{'}$ , while the posterior average for the product of $α_{x^{'}}$ and $α_{x^{''}}$ is obtained directly from (7), by choosing $k_{x} = δ_{x = x^{'}} + δ_{x = x^{''}}$ , i.e., $k_{x} = 1$ for $x \in {x^{'}, x^{''}}$ and $k_{x} = 0$ for $∀x\nin{x′,x′′}$ . ∎

Proof of Theorem 1.

Given $α$ , the marginal posterior density for $θ_{X | y}$ is a Dirichlet distribution with parameters $s α + n_{y}$ , where $n_{y} = (n_{x_{1} y}, \dots, n_{x_{r} y})^{'}$ . It is thus easy to compute

E [θ_{x | y} ∣ ∣ α, D] = E^{α, D} [θ_{x | y}] = \frac{n_{x y} + s α_{x}}{n_{y} + s},

where $x \in X$ , $y \in Y$ , and $E^{α, D} [\cdot] = E [\cdot | α, D]$ .

The posterior average of $θ_{x | y}$ can thus be computed by means of the law of total expectation as

E^{D} [θ_{x | y}] = E^{D} [E^{α, D} [θ_{x | y}]] = \frac{n_{x y} + s E^{D} [α_{x}]}{n_{y} + s} .

The posterior average of $α_{x}$ is obtained directly from Lemma 1.

In order to compute the posterior covariance between $θ_{x | y}$ and $θ_{x^{'} | y^{'}}$ we can use the law of total covariance, i.e.,

	${C ov}^{D} (θ_{x \| y}, θ_{x^{'} \| y^{'}}) =$	${C ov}^{D} (E^{α, D} [θ_{x \| y}], E^{α, D} [θ_{x^{'} \| y^{'}}]) +$
		$+ E^{s, D} [{C ov}^{α, D} (θ_{x \| y}, θ_{x^{'} \| y^{'}})] .$

The first quantity is:

	${C ov}^{D} (E^{α, D} [θ_{x \| y}], E^{α, D} [θ_{x^{'} \| y^{'}}])$

If $y^{'} = y$ , the second quantity is:

	$E^{D} [{C ov}^{α, D} (θ_{x \| y}, θ_{x^{'} \| y})]$
	$= E^{D} [\frac{(n_{x y} + s α_{x}) ((n_{y} + s) δ_{x x^{'}} - (n_{x^{'} y} + s α_{x^{'}}))}{(n_{y} + s)^{2} (n_{y} + s + 1)}]$
	$= \frac{(n_{x y} + s E^{D} [α_{x}]) δ_{x x^{'}}}{(n_{y} + s) (n_{y} + s + 1)} - \frac{E^{D} [(n_{x y} + s α_{x}) (n_{x^{'} y} + s α_{x^{'}})]}{(n_{y} + s)^{2} (n_{y} + s + 1)}$
	$= \frac{{^θ}_{x \| y} δ_{x x^{'}} - {^θ}_{x \| y} {^θ}_{x^{'} \| y^{'}}}{n_{y} + s + 1} - \frac{s^{2} (E^{D} [α_{x} α_{x^{'}}] - E^{D} [α_{x}] E^{D} [α_{x^{'}}])}{(n_{y} + s)^{2} (n_{y} + s + 1)}$
	$= \frac{{^θ}_{x \| y} δ_{x x^{'}} - {^θ}_{x \| y} {^θ}_{x^{'} \| y^{'}}}{n_{y} + s + 1} - \frac{s^{2} {C % ov}^{D} (α_{x}, α_{x^{'}})}{(n_{y} + s)^{2} (n_{y} + s + 1)} .$

Otherwise, if $y^{'} \neq y$ , $E^{D} [{C ov}^{α, D} (θ_{x | y}, θ_{x^{'} | y^{'}})] = 0$ , since $θ_{X | y} t o 0.0 p t $ ⊥ $ ⊥ θ_{X | y^{'}}$ given $α$ .

Exploiting the law of total covariance, we obtain

	${C ov}^{D} (θ_{x \| y}, θ_{x \| y^{'}}) = \frac{s^{2} {C ov}^{D} (α_{x}, α_{x^{'}})}{(n_{y} + s) (n_{y^{'}} + s)} +$
	$+ δ_{y y^{'}} ⎛ ⎝ \frac{{^θ}_{x \| y} δ_{x x^{'}} - {^θ}_{x \| y} {^θ}_{x^{'} \| y}}{n_{y} + s + 1} - \frac{s^{2} {C ov}^{D} (α_{x}, α_{x^{'}})}{(n_{y} + s)^{2} (n_{y} + s + 1)} ⎞ ⎠ .$

The posterior covariance between $α_{x}$ and $α_{x^{'}}$ is obtained directly from Lemma 1. ∎

Proof of Theorem 2.

Exploiting the linearity of the ideal estimator we obtain that

θ_{x | y}^{*} - θ_{x | y} = {~ ω}_{y} (θ_{x | y}^{ML} - θ_{x | y}) + (1 - {~ ω}_{y}) ({~ α}_{x} - θ_{x | y}) .

The MSE of $θ_{x | y}^{*}$ is thus

	$MSE (θ_{x \| y}^{*}) = {~ ω}_{y}^{2} MSE (θ_{x \| y}^{ML}) + (1 - {~ ω}_{y})^{2} MSE ({~ α}_{x}) +$
	$+ {~ ω}_{y} (1 - {~ ω}_{y}) E_{θ} [E_{n} [(θ_{x \| y}^{ML} - θ_{x \| y}) ({~ α}_{x} - θ_{x \| y})]] .$

Using the definition of MSE (4) and the assumptions of model (3), i.e., $E_{n} [n_{x y}] = n_{y} θ_{x | y}$ , $E_{θ} [θ_{x | y}] = {~ α}_{x}$ and ${V ar}_{θ} (θ_{x | y}) = \frac{{~ α}_{x} - {~ α}_{x}^{2}}{s + 1}$ , we obtain

	$MSE (θ_{x \| y}^{ML}) = E_{θ} [E_{n} [{(\frac{n_{x y}}{n_{y}} - θ_{x \| y})}_{x \| y}^{2}]]$
	$= E_{θ} [\frac{1}{n_{y}^{2}} {V % ar}_{n} (n_{x y})] = \frac{n_{y}}{n_{y}^{2}} E_{θ} [θ_{x \| y} (1 - θ_{x \| y})]$
	$= \frac{1}{n_{y}} (E_{θ} [θ_{x \| y}] - E_{θ} {[θ_{x \| y}]}_{x \| y}^{2} - {V ar%}_{θ} (θ_{x \| y}))$
	$= \frac{1}{n_{y}} ({~ α}_{x} - {~ α}_{x}^{2} - \frac{{~ α}_{x} - {~ α}_{x}^{2}}{s + 1}) = \frac{s}{n_{y}} \frac{{~ α}_{x} - {~ α}_{x}^{2}}{s + 1} .$		(13)

This quantity corresponds to the first term of $MSE (θ_{x | y}^{*})$ . The second term is obtained as

	$MSE ({~ α}_{x})$
		$= {V ar}_{θ} (θ_{x \| y}) = \frac{{~ α}_{x} - {~ α}_{x}^{2}}{s + 1},$

since ${~ α}_{x} - θ_{x | y}$ is independent of $n_{x y}$ and $E_{θ} [θ_{x | y}] = {~ α}_{x}$ . The last term , since ${~ α}_{x} - θ_{x | y}$ is independent of $n_{x y}$ and $E_{θ} [θ_{x | y}^{ML}] = θ_{x | y}$ .
The MSE for the ideal estimator is thus

	$MSE (θ_{x \| y}^{*})$	$= {~ ω}_{y}^{2} \frac{s}{n_{y}} \frac{{~ α}_{x} - {~ α}_{x}^{2}}{s + 1} + (1 - {~ ω}_{y})^{2} \frac{{~ α}_{x} - {~ α}_{x}^{2}}{s + 1}$
		$= ({~ ω}_{y}^{2} \frac{s}{n_{y}} + (1 - {~ ω}_{y})^{2}) \frac{{~ α}_{x} - {~ α}_{x}^{2}}{s + 1} .$

If ${~ ω}_{y} = \frac{n_{y}}{n_{y} + s}$ and $s > 0$ , $\forall x \in X$ ,

	$MSE (θ_{x \| y}^{*}) = (\frac{n_{y}^{2}}{(n_{y} + s)^{2}} \frac{s}{n_{y}} + \frac{s^{2}}{(n_{y} + s)^{2}}) \frac{{~ α}_{x} - {~ α}_{x}^{2}}{s + 1}$
	$= \frac{s}{n_{y} + s} \frac{{~ α}_{x} - {~ α}_{x}^{2}}{s + 1} < \frac{s}{n_{y}}$

	$E^{D} [{C ov}^{α, D} (θ_{x \| y}, θ_{x^{'} \| y})]$
	$= E^{D} [\frac{(n_{x y} + s α_{x}) ((n_{y} + s) δ_{x x^{'}} - (n_{x^{'} y} + s α_{x^{'}}))}{(n_{y} + s)^{2} (n_{y} + s + 1)}]$
	$= \frac{(n_{x y} + s E^{D} [α_{x}]) δ_{x x^{'}}}{(n_{y} + s) (n_{y} + s + 1)} - \frac{E^{D} [(n_{x y} + s α_{x}) (n_{x^{'} y} + s α_{x^{'}})]}{(n_{y} + s)^{2} (n_{y} + s + 1)}$
	$= \frac{{^θ}_{x \| y} δ_{x x^{'}} - {^θ}_{x \| y} {^θ}_{x^{'} \| y^{'}}}{n_{y} + s + 1} - \frac{s^{2} (E^{D} [α_{x} α_{x^{'}}] - E^{D} [α_{x}] E^{D} [α_{x^{'}}])}{(n_{y} + s)^{2} (n_{y} + s + 1)}$
	$= \frac{{^θ}_{x \| y} δ_{x x^{'}} - {^θ}_{x \| y} {^θ}_{x^{'} \| y^{'}}}{n_{y} + s + 1} - \frac{s^{2} {C % ov}^{D} (α_{x}, α_{x^{'}})}{(n_{y} + s)^{2} (n_{y} + s + 1)} .$

	$MSE (θ_{x \| y}^{ML}) = E_{θ} [E_{n} [{(\frac{n_{x y}}{n_{y}} - θ_{x \| y})}_{x \| y}^{2}]]$
	$= E_{θ} [\frac{1}{n_{y}^{2}} {V % ar}_{n} (n_{x y})] = \frac{n_{y}}{n_{y}^{2}} E_{θ} [θ_{x \| y} (1 - θ_{x \| y})]$
	$= \frac{1}{n_{y}} (E_{θ} [θ_{x \| y}] - E_{θ} {[θ_{x \| y}]}_{x \| y}^{2} - {V ar%}_{θ} (θ_{x \| y}))$
	$= \frac{1}{n_{y}} ({~ α}_{x} - {~ α}_{x}^{2} - \frac{{~ α}_{x} - {~ α}_{x}^{2}}{s + 1}) = \frac{s}{n_{y}} \frac{{~ α}_{x} - {~ α}_{x}^{2}}{s + 1} .$		(13)

Hierarchical Multinomial-Dirichlet model for the estimation of conditional probability tables

Authors

Learning Bayesian Networks with Low Rank Conditional Probability Tables

Accurate parameter estimation for Bayesian Network Classifiers using Hierarchical Dirichlet Processes

Categorical Aspects of Parameter Learning

User Interface Tools for Navigation in Conditional Probability Tables and Elicitation of Probabilities in Bayesian Networks

An Odds Ratio Based Inference Engine

MADE: Masked Autoencoder for Distribution Estimation

Learning from demonstration using products of experts: applications to manipulation and task prioritization

I Introduction

Ii Estimation Under Multinomial-Dirichlet Model

Iii Hierarchical Model

Iii-a Posterior moments for $θ_{X | y}$

Definition 1.

Theorem 1.

Iii-B Properties of the estimator ${^θ}_{X | y}$

Theorem 2.

Iv Experiments

Iv-a MSE analysis

Iv-B Joint distribution fitting

Iv-C Classification

V Conclusions

Lemma 1.

Proof of Lemma 1.

Proof of Theorem 1.

Proof of Theorem 2.

Hierarchical Multinomial-Dirichlet model for the estimation of conditional probability tables

Authors

Learning Bayesian Networks with Low Rank Conditional Probability Tables

Accurate parameter estimation for Bayesian Network Classifiers using Hierarchical Dirichlet Processes

Categorical Aspects of Parameter Learning

User Interface Tools for Navigation in Conditional Probability Tables and Elicitation of Probabilities in Bayesian Networks

An Odds Ratio Based Inference Engine

MADE: Masked Autoencoder for Distribution Estimation

Learning from demonstration using products of experts: applications to manipulation and task prioritization

This week in AI

I Introduction

Ii Estimation Under Multinomial-Dirichlet Model

Iii Hierarchical Model

Iii-a Posterior moments for θX|y

Definition 1.

Theorem 1.

Iii-B Properties of the estimator ^θX|y

Theorem 2.

Iv Experiments

Iv-a MSE analysis

Iv-B Joint distribution fitting

Iv-C Classification

V Conclusions

Lemma 1.

Proof of Lemma 1.

Proof of Theorem 1.

Proof of Theorem 2.

Iii-a Posterior moments for $θ_{X | y}$

Iii-B Properties of the estimator ${^θ}_{X | y}$