Bayesian Variational Inference for Mixed Data Mixture Models

Junyang Wang1 James Bennett2 Victor Lhoste2 Sarah Filippi1

July 23, 2025

Abstract.
Heterogeneous, mixed type datasets including both continuous and categorical variables are ubiq-
arXiv:2507.16545v1 [[Link]] 22 Jul 2025

uitous, and enriches data analysis by allowing for more complex relationships and interactions to be
modelled. Mixture models offer a flexible framework for capturing the underlying heterogeneity and
relationships in mixed type datasets. Most current approaches for modelling mixed data either forgo
uncertainty quantification and only conduct point estimation, and some use MCMC which incurs a very
high computational cost that is not scalable to large datasets. This paper develops a coordinate ascent
variational inference algorithm (CAVI) for mixture models on mixed (continuous and categorical) data,
which circumvents the high computational cost of MCMC while retaining uncertainty quantification. We
demonstrate our approach through simulation studies as well as an applied case study of the NHANES
risk factor dataset. In addition, we show that the posterior means from CAVI for this model converge
to the true parameter value as the sample size n tends to infinity, providing theoretical justification for
our method.
Keywords. Mixture model • Mixed type data • Variational inference

1 Introduction
Mixture models provide a probabilistic framework for modelling latent subpopulations within a popu-
lation, without needing to know which subpopulation each individual data point belongs to. Mixture
models are ubiquitous in statistical modelling and are commonly used for clustering and density esti-
mation. Popular applications include image analysis [Niknejad et al., 2015, Beecks et al., 2015], finance
[Alexander and Lazar, 2006, Durham, 2007] and genomics [Pritchard et al., 2000, Rau et al., 2015] .
Commonly in data analysis, datasets include both continuous and categorical variables, and such
datasets are known as mixed type data. For instance, healthcare datasets might include continuous
variables such as a patient’s blood pressure, weight and height, but also categorical variables such as
smoking status. In such mixed datasets, the ability to model multimodal continuous and categorical
data simultaneously is critical in understanding the relationship between the variables of interest, and
the underlying subpopulations in the data.
Existing mixture models for mixed type data have been developed in various forms. Hunt and Jor-
gensen [1999] proposed a model where each mixture component is a product of multivariate Gaussian
distributions for the continuous data, and a categorical distribution for the categorical data. The authors
use an Expectation-Maximisation (EM) algorithm to calculate the maximum likelihood estimate of the
parameters of interest, which was later extended in Hunt and Jorgensen [2003] to include missing data.
Lawrence and Krzanowski [1996] uses a similar parametrisation, though allowing the mean of the Gaus-
sian components to depend on the value of the categorical data, but the covariance of the components
are assumed to be independent of the categorical data as well as the mixture components in order to
reduce the number of parameters estimated with EM.
There are also a number of works that assume categorical variables are generated from latent con-
tinuous variables, and the values of the categorical variables correspond to certain cutoffs in the latent
continuous variable. Examples of these include Ranalli and Rocci [2017] and McParland and Gormley
[2016], both of which use EM for parameter estimation.
1 Department of Mathematics, Imperial College London, London, United Kingdom.
2 School of Public Health, Imperial College London, London, United Kingdom.

1
 Bayesian Variational Inference for Mixed Data Mixture Models

Approaches using Bayesian inference to obtain a posterior distribution over the model parameters
rather than a point estimate have also been used widely in mixture models, especially for clustering
(Celeux et al. [2020], Celeux et al. [2000], Frühwirth-Schnatter [2001], Miller and Harrison [2018], Richard-
son and Green [1997], Zhang et al. [2004], Teh et al. [2006]). However, mixed type data mixture models
using Bayesian inference appear to be much more limited (Wang et al. [2019], Zhang and Berhane [2014]).
MCMC for mixture models are however not without significant drawbacks. MCMC is computationally
intensive especially in higher dimensions, and diagnosing convergence can be challenging. Due to the
non-identifiability of mixture components under symmetric priors, MCMC samples can exhibit label
switching, where component labels permute during sampling Jasra et al. [2005]. This makes direct inter-
pretation of component-specific parameters difficult without additional post-processing or constraints.
Variational inference (VI) is a popular method used in Bayesian statistics and machine learning to
approximate the posterior distribution when it is difficult or intractable to compute (Blei et al. [2017]).
VI approximates the posterior distribution using optimisation, and is typically faster and more scalable
than MCMC.
Variational inference has been utilised for mixture models previously. For finite mixtures and con-
tinuous data, it is well known that a Normal-inverse Wishart prior on the mixture component mean
and covariance and a Dirichlet prior on the mixture component probabilities, coupled with Gaussian
data, leads to closed form updates using the Coordinate Ascent Variational Inference (CAVI) algorithm
[Bishop, 2006]. For finite mixtures and categorical data, Ahlmann-Eltze and Yau [2018] proposed a
CAVI algorithm for multinomially distributed discrete data. Infinite mixture models using Dirichlet
Process and variational inference include the seminal work of Blei and Jordan [2006] as well as further
developments such as Wang et al. [2011], Hughes and Sudderth [2013] and Huynh et al. [2016].
In this paper, we propose a VI algorithm for finite mixture models on mixed type data, which allows
both the continuous and categorical variables to influence the make up of each mixture component. We
employ a mean-field variational approximation of the posterior with conjugate priors in order to obtain
analytical updates of the variational posterior using the CAVI algorithm, with the derivation equations
provided. We also perform a convergence analysis of the proposed CAVI algorithm. In particular, we
demonstrate that the posterior mean produced by the algorithm converges to the true parameter value
as the sample size approaches infinity, provided the data are generated according to the true model. We
investigate the empirical performance of our algorithm extensively on a number of simulated datasets,
where both point estimation and uncertainty quantification of the variational posterior are assessed and
compared with a Gibbs sampler and existing EM methods. Finally, we analyse data from the National
Health and Nutrition Examination Survey (NHANES), a program conducted by the National Center for
Health Statistics (NCHS) that assesses the health and nutritional status of the U.S. population through
interviews and physical examinations, from which we use continuous and categorical risk factor variables
such as body mass index (BMI) and smoking status. These risk factors frequently co-occur and are
often etiologically linked; for example, obesity is a known risk factor for dyslipidemia, elevated blood
pressure, and hyperglycemia, making it difficult to accurately assess health risks when considering each
factor in isolation. We show our method can produce medically interpretable clusters of participants,
and we visualise said clusters in a way that takes into account uncertainty using the posterior predictive
distribution, something that is missed in most frequentist and non model based clustering methods.

2 Mixed data mixture model

Consider a set of n independent observations {(xi , ci )}1≤i≤n where xi ∈ Rq contains continuous observa-
tions while ci ∈ Zp is a p-dimensional vector of categorical data such that cij ∈ {1, 2, . . . , dj }. Note that
dj denotes the number of categorical values for the j-th categorical variable. In this paper we focus on
a finite mixture model of K components, where K is a constant specified by the user. In particular, we
consider a model where the likelihood for a single data point (xi , ci ) can be written as follows:
K
X p
Y
p(xi , ci |µ, Λ, ψ, π) = πk N (xi |µk , Λ−1
k ) ψk,j,cij . (1)
k=1 j=1

Here, πk are the mixture weights, denoting the probability of each component 1 ≤ k ≤ K; and ψk,j,cij
denotes the probability of observing cij for the jth categorical variable under component k. In addition,

2
 Bayesian Variational Inference for Mixed Data Mixture Models

µk and Λk denotes the component centre and component precision matrix respectively for component k,
and we will also use Σk = Λ−1
k to denote the covariance matrix for component k. Note that, for a given
mixture component k, conditional independence is assumed between the p categorical variables in ci as
well as between ci and xi .
Often for computational reasons, it is beneficial to introduce a latent variable z of size n denoting
the realised mixture components of each observations. Let x ∈ Rn×q and c ∈ Zn×p represent the matrix
containing respectively the continuous and categorical data. Incorporating z, the model can be written
as:
n
Y Y p
n Y n
Y
p(z|π) = πzi , p(c|ψ, z) = ψzi ,j,cij , p(x|c, z, µ, Λ) = N (xi |µzi , Λ−1
zi ) (2)
i=1 i=1 j=1 i=1

As we are presenting a Bayesian approach to infer the parameters of this mixture model, we now define
the prior distributions over the model parameters, for which we consider the following:

µk ∼ N (mk , β −1 Λ−1
k ), π = (π1 , . . . , πK ) ∼ Dir(α, . . . , α)
−1
Λk ∼ W(ν, Φ ), ψk,j = (ψk,j,1 , . . . , ψk,j,dj ) ∼ Dir(ηj , . . . , ηj ) (3)

The prior of µk and Σk = Λ−1 k is a joint Normal-Inverse-Wishart distribution while π and ψk,j take
Dirichlet priors. The prior hyperparameters mk , β, ν, Φ, α, ηj are specified by the user. The posterior
distribution
p(x|c, z, µ, Λ)p(c|ψ, z)p(z|π)p(µ, Λ, ψ, π)
p(µ, Λ, ψ, π, z|x, c) = (4)
p(x, c)
induced by this prior is not available in closed form. In this paper, we derive a variational inference
procedure to approximate the posterior.

3 Variational Inference
3.1 Background on Variational Inference and Coordinate Ascent
In this section, we introduce the key concepts of Variational Inference and the Coordinate Ascent Varia-
tional Inference (CAVI) procedure that are essential for understanding the remainder of the paper. For
more detailed explanations, we refer the reader to Bishop [2006]. Suppose that we have some data x (not
necessarily the data described in the previous section) generated by a likelihood function p(x|θ) where
θ denotes some unknown parameters of interest, which have a prior distribution p(θ). In variational
inference, the posterior distribution p(θ|x) is approximated using a variational distribution q, which is
chosen from a family of distributions Q that are computationally more tractable than the true posterior.
The variational distribution q is chosen by selecting the distribution from the variational family which
minimises the Kullback-Leibler divergence (KL-divergence) between the variational distribution and the
true posterior. Specifically:
q = arg min KL(q(θ) ∥ p(θ|x)) (5)
q∈Q

In practice, the variational family Q is often assumed to be parametrised by variational parameters λ, and
the optimisation problem becomes equivalent to finding the parameter λ which minimises KL(qλ (θ) ∥ p(θ|x)).
One key advantage in this optimisation approach to Bayesian inference is that one no longer needs to eval-
uate the marginal distribution p(x), which is typically the main computational challenge in conducting
Bayesian inference through sampling. To see this, note that:

p(x|θ)p(θ)
λ∗ = arg min KL(qλ (θ) ∥ p(θ|x)) = arg min Eqλ [ln(qλ (θ)) − ln( )]
λ λ p(x)
= arg min Eqλ [ln(qλ (θ)) − ln(p(x|θ)) − ln(p(θ))] (6)
λ

where ln(p(x)) can be dropped from the optimisation objective as it is independent of qλ . The negative
of the remaining optimisation objective Eqλ [−ln(qλ (θ)) + ln(p(x|θ)) + ln(p(θ))] is commonly referred to
as the Evidence Lower Bound (ELBO). The nomenclature here refers to the fact that the ELBO is upper

3
 Bayesian Variational Inference for Mixed Data Mixture Models

bound by ln(p(x)), which follows from the non-negativity of the KL divergence KL(qλ (θ) ∥ p(θ|x)).
Note minimising KL(qλ (θ) ∥ p(θ|x)) is equivalent to maximising the ELBO.
A common assumption for the variational family is that the variational posterior factorises into a
product of independent distributions:
I
Y
q(θ) = qi (θi ) (7)
i=1

where qi is the posterior for the ith block of variables θi . This is known as the mean field assumption.
One popular and convenient algorithm for solving the optimisation problem in Equation (6) under the
mean field assumption is Coordinate Ascent Variational Inference (CAVI). CAVI maximises the ELBO
by iteratively updating each factor in the variational distribution while treating all other factors as fixed,
until the ELBO has converged to a local maximum. Given the variational posterior in Equation (7), the
ELBO becomes:
Z Y I I
X
ELBO = qi (θi )[ln(p(x, θ)) − ln(qi (θi ))]dθ
i=1 i=1
Z Y I
Z Y I
X
= qj (θj )[ln(p(x, θ)) qi (θi )dθi ]dθj − qi (θi ) ln(qi (θi ))dθ
i̸=j i=1 i=1
Z Y I Z
X
= qj (θj )[ln(p(x, θ)) qi (θi )dθi ]dθj − qi (θi )ln(qi (θi ))dθi
i̸=j i=1
Z Z
= qj (θj )[Eθ−j [ln(p(x, θ))] + C1 ]dθj − qj (θj )ln(qj (θj ))dθj + C2 (8)
R Q
where in the last line, Eθ−j [ln(p(x, θ)] = ln(p(x, θ)) i̸=j qi (θi )dθi . C1 and C2 are constants with
respect to qj . So to maximise the ELBO with respect to each qj while keeping the other qi , i ̸= j fixed,
one recognises the last line is a negative KL divergence between qj (θj ) and exp(C1 )exp(Eθ−j [ln(p(x, θ))]),
it follows that the ELBO is maximised when qj (θj ) ∝ exp(Eθ−j [ln(p(x, θ))]). Therefore, at each iteration,
CAVI updates the qj in turn via
qj (θj ) ∝ exp(Eθ−j [ln(p(x, θ))]) (9)
while keeping all qi , i ̸= j fixed.

3.2 CAVI for mixed data mixture model

In this section, we derive a variational inference procedure for the mixture model defined in Equation (1),
building on the general concepts and techniques introduced in the previous section. To approximate the
posterior distribution p(µ, Λ, ψ, π, z|x, c) from Equation (4), we consider a variational distribution of the
following form:
q(µ, Λ, ψ, π, z) = q(µ, Λ, ψ, π)q(z) . (10)
This variational posterior assumes a general joint distribution for the ‘global’ variables µ, Λ, ψ, π, and a
separate distribution for the ‘local’ variables z that is independent from the ‘global’ variables. The exact
distributional families of the variational posterior q(µ, Λ, ψ, π) and q(z) need not to be specified in advance
and will fall out during the derivation. Independence between q(µ, Λ, ψ, π) and q(z) is a standard mean
field assumption that allows us to split the variables into two blocks for the Coordinate Ascent Variational
Inference (CAVI) algorithm. In other words,the algorithm alternates between updating q(µ, Λ, ψ, π) and
q(z). By applying the CAVI update from equation (9) to our variational posterior in (10), we obtain
(details of the derivations can be found in the supporting information):
q(µ, Λ, ψ, π) ∝ exp(E−µ,−Λ,−ψ,−π lnp(x, c, z, µ, Λ, ψ, π) ) (11)
∝ exp(Ez [ln(p(x|c, z, µ, Λ)p(µ|Λ)p(Λ)p(c|ψ, z)p(ψ)p(z|π)p(π))])
∝ exp(Ez [ln(p(x|c, z, µ, Λ)p(µ|Λ)p(Λ))]) exp(Ez [ln(p(z|π)p(π))]) exp(Ez [ln(p(c|ψ, z)p(ψ))])
| {z }| {z }| {z }
∝q(µ,Λ) ∝q(π) ∝q(ψ)

4
 Bayesian Variational Inference for Mixed Data Mixture Models

The variational distribution can therefore be written as the product of three terms which can be derived
separately: q(µ, Λ, ψ, π) = q(µ, Λ)q(ψ)q(π). The variational distribution over π, q(π), can be shown to
take the form:
n
X
q(π) = Dir(π|α̂1 , . . . , α̂K ) with α̂k = α + rik (12)
i=1

where rik = Ezi [zik ] = q(zi = k) is the variational posterior probability of the i-th data being in
component k under the variational posterior probability. The term q(µ, Λ) can be shown to take the
following form:
Y
q(µ, Λ) = N (µk |m̂k , βk−1 Λ−1 −1
k )W(Λk |ν̂k , Φ̂k ))
k

where
Pn n
βmk + rjk xj
Pnj=1
X
m̂k = , β̂k−1 = (β + rjk )−1
β+ j=1 rjk j=1
n
X n
X n
X
ν̂k = ν + rik , Φ̂k = Φ − (β + rik )m̂k m̂Tk + βmk mTk + rik xi xTi (13)
i=1 i=1 i=1

This gives for each k, independent Normal-Inverse-Wishart distributions µk |Λk ∼ N (µk |m̂k , β̂k−1 Λ−1
k )
and Λk ∼ W(Λk |ν̂k , Φ̂−1
k )). Lastly, q(ψ) can be shown to take the form:

YY n
X
q(ψ) = Dir(ψk,j |η̂k,j,1 , . . . , η̂k,j,dj ) with η̂k,j,g = ηj + I{cij =g} rik (14)
k j i=1

After updating q(µ, Λ, ψ, π), we use its newly updated distribution to in turn update q(z), which can be
shown to take the form:
n o
q(z) ∝ exp E−z lnp(x, c, z, µ, Λ, ψ, π)
n X
nX K n X
X p
K X n X
X K o
∝ exp zik Eµ,Λ [ln(N (xi |µk , Λk−1 ))] + zik Eψ [ln(ψk,j,cij )] + zik Eπ [ln(πk )]
i=1 k=1 i=1 k=1 j=1 i=1 k=1
n Y
Y K
∝ ρzikik (15)
i=1 k=1

where
p
n X o ρik
ρik = exp Eµ,Λ [ln(N (xi |µk , Λ−1
k ))] + Eψ [ln(ψk,j,cij )] + Eπ [ln(πk )] , rik = PK (16)
j=1 j=1 ρij

The exact forms of each of the expectations inside ρik can be found in the supplement. Algorithm 1
summarises the CAVI updates for all the parameters. Details including the full derivation for each update
and the exact forms of the ELBO is reserved for the supplement.

3.3 Posterior Predictive Distribution

Now that we have the variational posterior, we can use the variational posterior predictive distribution to
predict possible future (unobserved) data points given the observed data, while accounting for uncertainty
in the model parameters. This can be used for purposes such as density estimation. We derive the
variational posterior predictive distribution of a new observation x̃, c̃, which conveniently can be obtained

5
 Bayesian Variational Inference for Mixed Data Mixture Models

Algorithm 1 CAVI algorithm for mixed data mixture model

Initialise α̂k , m̂k , β̂k , Φ̂k , ν̂k , η̂k,j,g , rik
while ELBO(α̂k , m̂k , β̂k , Φ̂k , ν̂k , η̂k,j,g , rik ) not converged do
Update α̂k for k = 1, . . . K using Equation (12)
Update m̂k , β̂k , Φ̂k , ν̂k for k = 1, . . . K using Equation (13)
Update η̂k,j,g for k = 1, . . . K using Equation (14)
Update rik for k = 1, . . . K using Equation (16)
compute ELBO
return α̂k , m̂k , β̂k , Φ̂k , ν̂k , η̂k,j,g , rik

in closed form in this case:

Z
q(x̃, c̃|x, c) = p(x̃, c̃|µ, Λ, ψ, π)q(µ, Λ)q(ψ)q(π)dµdΛdψdπ
K
X p
Y
= Eπ [πk ] Eµ,Λ [N (x̃|µk , Λ−1
k )] Eψ [ ψk,j,c̃j ]
k=1 j=1
K p
X α̂k  Φ̂k (β̂k + 1)  Y η̂k,j,c̃j
= PK tν̂k −q+1 x̃ m̂k , Pdj (17)
k=1 k′ =1 αk′ β̂k (ν̂k − q + 1) j=1 gj =1 η̂k,j,gj

where expectations such as Eπ are with respect to the variational posterior distributions, tν (x̃|m, Φ)
denotes a multivariate student t density function with ν degrees of freedom, location m and scale matrix
Φ.
From the form of (17), we can also obtain the marginals of the variational posterior predictive
densities for each jth component of newly observed continuous data x̃j or categorical data c̃j given a
particular mixture component k, i.e. q(x̃j |z̃ = k, c̃, x, c) or q(c̃j |z̃ = k, x, c). In particular, a component
of a multivariate t distribution is an univariate t distribution with mean and variance equal to the
corresponding components of the mean and covariance of the multivariate t distribution, so q(x̃j |z̃ =
k, c̃, x, c) and q(c̃j |z̃ = k, x, c) take the following forms:

Φ̂k,jj (β̂k + 1) η̂k,j,c̃j

q(x̃j |z̃ = k, c̃, x, c) = tν̂k −q+1 (x̃j |m̂k,j , ), q(c̃j |z̃ = k, x, c) = Pdj (18)
β̂k (ν̂k − q + 1) gj =1 η̂k,j,gj

These marginal distributions of the variational posterior predictive densities enable useful visualisations
of the variational posterior predictive density in high dimensions that will be presented in Section 5.

4 Theory
This section is dedicated to examining the convergence of the proposed CAVI algorithm to the true,
data generating parameters {πk∗ , µ∗k , Λ∗k , {ψk,g ∗
}g }k as the sample size n tends to infinity, assuming the
data have been generated according to Equation (1). For ease of presentation we assume there is p = 1
categorical variable with d categorical PK values. We also assume Pd d and K, the number of components K is
correctly specified, so πk∗ > 0, k=1 πk∗ = 1 and ψk,g ∗
> 0, g=1 ψk,g ∗
= 1 for 1 ≤ k ≤ K and 1 ≤ g ≤ d.
Our approach follows that of Titterington and Wang [2006], who developed convergence results for a
CAVI algorithm for a Gaussian mixture model. In particular, we extend Theorem 1 of Titterington and
Wang [2006] to the case of mixed type data involving mixtures of Gaussian distributed continuous data
and categorical data that are conditionally independent within each mixture.
The CAVI algorithm derived in the previous section iteratively updates the variational hyperparame-
ters α̂k , m̂k , β̂k , Φ̂k , ν̂k , η̂k,g . Denoting with an upper-script t the value of these hyperparameters at the
t-th iteration of the CAVI algorithm, α̂kt , m̂tk , β̂kt , Φ̂tk , ν̂kt , η̂k,g
t
are deterministic functions of α̂kt−1 , m̂t−1
k ,
t−1 t−1 t−1 t−1
β̂k , Φ̂k , ν̂k , η̂k,g where the transformation is computed by evaluating the intermediate quantities
(t−1)
rik . In this section, instead of directly working with the variational hyper-parameters, we will focus

6
 Bayesian Variational Inference for Mixed Data Mixture Models

on slightly different quantities defined as follows:

n Pn (t−1) Pn (t−1)
1 X (t−1) i=1 rik xi i=1 I{ci =g} rik
π̂kt = r , µ̂tk = Pn (t−1) ,
t
ψ̂k,g =
n i=1 ik
Pn (t−1)
i=1 rik i=1 rik
n n
(t−1) (t−1)
X X
Λ̂tk = ( rik )( rik (xi − µ̂tk )(xi − µ̂tk )T )−1 (19)
i=1 i=1

Note that these quantities are related to the variational hyperparameter estimates from CAVI:

nµ̂tk π̂kt + βmk

α̂kt = nπ̂kt + α, t
η̂k,g = nπ̂kt ψ̂k,g
t
+ η, m̂tk = , β̂kt = nπ̂kt + β
nπ̂kt + β
nπ̂kt
Φ̂tk = nπ̂kt (Λ̂tk )−1 + β(µ̂tk − mk )(µ̂tk − mk )T + Φ, ν̂kt = nπ̂kt + ν (20)
nπ̂kt + β

The newly introduced quantities π̂kt , µ̂tk , Λ̂tk , ψ̂k,g

t
are closely related to the variational posterior means of
Eπ [πk ], Eµ [µk ], EΛ [Λk ], Eψ [ψk ] respectively. More precisely, π̂kt is equal to the variational posterior mean
of q(πk ) parametrised by α̂kt in the case where the prior hyperparameter α = 0. Similar connection can
be done for each new quantity and the associated variational posterior mean in the case where the prior
hyperparameters β, ν, Φ, η are set to 0. It turns out that analysing these somewhat simplified quantities
are sufficient in the sense that if they converge as n → ∞, then the corresponding variational posterior
means will asymptotically converge to the same limits.
Similarly to the variational hyperparameters, the newly introduced quantities π̂kt , µ̂tk , Λ̂tk , ψ̂k,g t
can
(t−1)
be iteratively computed at each iteration by evaluating functions of the intermediate terms rik . De-
(t−1)
noting by Θt the set of all quantities {π̂kt , µ̂tk , Λ̂tk , {ψ̂k,g
t
}g }k at iteration t, then as rik are functions
of the variational hyperparameters at t − 1 and therefore Θ , we can write Θ = T (Θ(t−1) ) for some
t−1 t

deterministic transformation T which encapsulates our proposed CAVI algorithm. In the following we
define this function T more formally. To do so we introduce some vectorized notation and define the
required vector space for the domain of our function. We drop the upper script t so that the quantities
such as π̂k represent a generic input to the function T .
Let π̂k ∈ R, µ̂k ∈ Rq , Λ̂k ∈ S n (R), the set of real, symmetric q by q matrices, and ψ̂k ∈ Rd . Next, let
Sπ , Sµ , SΛ , Sψ be the K-fold direct sums of R, Rq , S n (R), Rd respectively. Then we define the following
vectors:
   ∗
 
π̂1
 
µ̂1
 
Λ̂1
 
ψ̂1 π̂ π
 ..   ..   ..   ..   µ̂  ∗  µ∗ 
π̂ =  . , µ̂ =  . , Λ̂ =  . , ψ̂ =  . , Θ =  Λ̂, Θ = Λ∗ 
   (21)
π̂K µ̂K Λ̂K ψ̂K ψ̂ ψ ∗

L L L
so that π̂ ∈ Sπ , µ̂ ∈ Sµ , Λ̂ ∈ SΛ , ψ̂ ∈ Sψ . Θ belongs to direct sum space Sπ Sµ SΛ Sψ . We also
let Θ∗ be the true data generating parameter of the likelihood Equation (1). As the notation implies,
Sπ , Sµ , SL
Λ , SψLare vector
L spaces which contain possible values of the corresponding parameters π̂, µ̂, Λ̂, ψ̂
and Sπ Sµ SΛ Sψ a vector space which contains the possible values of the overall parameter Θ.
(t−1)
As the rik are functions of the variational hyperparameters at t − 1, they are also functions of
(t−1)
Θ , so from here on we will use notation such as rik (Θt−1 ) = rik
t−1
or rik (Θ) when highlighting the
dependence of rik on Θ1 . Let:
n Pn Pn
1X rik (Θ)xi I{c =g} rik (Θ)
Πk (Θ) = rik (Θ), Mk (Θ) = Pn i=1
, Ψk (Θ)g = i=1 Pn i
n i=1 i=1 rik (Θ) i=1 rik (Θ)
n
X n
X
Sk (Θ) = ( rik (Θ))( rik (Θ)(xi − Mk (Θ))(xi − Mk (Θ))T )−1 (22)
i=1 i=1
1 Note r (Θ) is not defined for all Θ as some Θ outside the support of the prior results in variational hyperparameters
ik
outside the hyperparameter space, but for those Θ we can trivially map rik (Θ) to a value like 1/K.

7
 Bayesian Variational Inference for Mixed Data Mixture Models

These functions denote the maps which takes the simplified estimates Θ(t) from one iteration to the next.
Vectorising these maps:
       
Π1 (Θ) M1 (Θ) S1 (Θ) Ψ1 (Θ)
Π(Θ) =  ... , M (Θ) =  ..  .   . 
, S(Θ) =  .. , Ψ(Θ) =  ..  (23)
   
.
ΠK (Θ) MK (Θ) SK (Θ) ΨK (Θ)

we can define the iterative procedure Θ(t) = T (Θ(t−1) ), where:

 
Π(Θ)
M (Θ)
T (Θ) =  S(Θ) 
 (24)
Ψ(Θ)

For theoretical purposes we consider a generalised iterative procedure, the map T ϵ defined for 0 < ϵ < 2
as:

Θ(t) = T ϵ (Θ(t−1) ) = (1 − ϵ)Θ(t−1) + ϵT (Θ(t−1) ) (25)

Note T 1 (Θ(t−1) ) = T (Θ(t−1) ), so that ϵ = 1 is our CAVI algorithm. It is suggested in Titterington and
Wang [2006] that different choices of ϵ could improve convergence rate relative to the base case of ϵ = 1,
but we do not investigate this in this paper. L L L
For our convergence analysis, we need to define a suitable norm on Sπ Sµ SΛ Sψ . To define
this norm, consider the following inner products:
1
⟨u1 , u2 ⟩πk = u1 ( )u2 u1 , u2 ∈ R
πk∗
⟨v1 , v2 ⟩µk = v1T (πk∗ Λ∗k )v2 v1 , v2 ∈ Rq
π∗
⟨W1 , W2 ⟩Λk = k T r(W1 Λ∗−1 k W2 Λ∗−1
k ) W1 , W2 ∈ S n (R) :
2
π∗
⟨a2 , a2 ⟩ψk = aT1 diag( k∗ )a2 a1 , a2 ∈ Rd (26)
ψk
π∗ ∗
πk
where diag( ψk∗ ) is a d by d diagonal matrix with entries ∗
ψk,g , g = 1, . . . , d. We can easily show that
k
2
⟨WL1 , W2 ⟩L
Λk is indeed
L an inner product. Using the above inner products, we define ⟨B1 , B2 ⟩, on B1 , B2 ∈
Sπ Sµ SΛ Sψ by simply adding the inner products of the constituent spaces:
K
X K
X K
X K
X
⟨., .⟩ = ⟨., .⟩πk + ⟨., .⟩µk + ⟨., .⟩Λk + ⟨., .⟩ψk (27)
k=1 k=1 k=1 k=1

Sψ via ∥B∥2 = ⟨B, B⟩.

L L L
This inner product induces a vector norm on Sπ Sµ SΛ

We can now state the main result:

Theorem 1. With probability 1 as the sample size n → ∞, the iterative procedure defined in Equation (25)
converges locally to the true data generating parameter Θ∗ whenever 0 < ϵ < 2. Locally here means
whenever the starting values are sufficiently near Θ∗ .
∗−1
2 ⟨W , W ⟩
1 2 Λk is an inner product. Indeed, T r(W2 Λk W1 Λ∗−1
k ) = T r(W1 Λ∗−1
k W2 Λ∗−1
k ) by the cyclic property of trace
so symmetry is satisfied. For α, β ∈ R, W1 , W2 , W3 ∈ S n (R),
T r((αW1 + βW3 )Λ∗−1
k W2 Λ∗−1
k ) = T r(αW1 Λ∗−1
k W2 Λ∗−1
k + βW3 Λ∗−1
k W2 Λ∗−1
k )
= αT r(W1 Λ∗−1
k W2 Λ∗−1
k ) + βT r(W3 Λ∗−1
k W2 Λ∗−1
k )
so linearity in the first argument is satisfied. Lastly,
1
∗− 2 1
∗− 2 ∗− 1 ∗− 1 1
∗− 2 1
∗− 2 ∗− 1 1
∗− 2
T r(W Λ∗−1
k W Λ∗−1
k ) = T r(Λk W Λk Λk 2
W Λk 2
) = T r((Λk W Λk )(Λk 2
W Λk T
) )≥0
so we have positive definiteness.

8
 Bayesian Variational Inference for Mixed Data Mixture Models

The proof of Theorem 1 is contained in the supplement. The proof strategy involves showing that
with probability 1 as n → ∞, the map T ϵ (Θ) is locally contractive near Θ∗ . Specifically, there exists
0 ≤ λ < 1 such that:

∥T ϵ (Θ) − T ϵ (Θ∗ )∥ ≤ λ∥Θ − Θ∗ ∥ (28)

for Θ sufficiently near Θ∗ . To show Equation (28), by Taylor’s theorem (see P315 of Bhatia [1997]), we
have:

∥T ϵ (Θ) − T ϵ (Θ∗ )∥ ≤ ∥∇T ϵ (Θ∗ )∥op ∥Θ − Θ∗ ∥ + O(∥Θ − Θ∗ ∥2 ) (29)

It is therefore sufficient to show ∇T ϵ (Θ∗ ) converges with probability 1 (as n → ∞) to an operator that
has operator norm less than 1.
The vector norm that we use in the context of Equation (28) and Equation (29) is induced by the inner
product defined in Equation (27). Similarly, the operator norm used in Equation (29) on ∥∇T ϵ (Θ∗ )∥op
is defined as the operator norm with respect to the vector norm, so:
M M M
∥A∥op = sup{∥AB∥ : ∥B∥ ≤ 1, B ∈ Sπ Sµ SΛ Sψ } (30)

Theorem 1 gives rise to the following corollary, whose proof can be found in the supplement:

Corollary 1. Given Theorem 1, the variational posterior mean of the CAVI algorithm also converges to
the true data generating parameter Θ∗ .

5 Experimental Results
In this section, we conduct simulation experiments to assess and validate our method. For all simu-
lated data in this section, the continuous part of the data was standardised before fitting the model.
Throughout this section we take the following values for the prior hyperparameters: mk = 0 for all k,
β = 1, Φ = 0.25I, ν = q + K + 1, α = 1/K, ηj = 1/dj , where dj is the number of possible values of
the jth categorical variable. We also initialise the CAVI algorithm using the Kprototypes algorithm, by
setting the initial rik to be equal to 0.9 at the k that Kprototypes predicts for each data point (xi , ci ),
and 0.1 otherwise.3 In this section, we will generally use Σk = Λ−1
k to refer to the covariance matrix of
component k to more easily compare the size of the component covariance compared to how much the
component means {µk }k are separated.

5.1 An illustrative one-dimensional example

For ease of visualisation, we begin with an example where the continuous data is of dimension q = 1. We
simulate n = 5000 data with the following model parameters: p = q = 1, K ∗ = 3, π1∗ = 0.4, π2∗ = 0.35,
π3∗ = 0.25, µ∗1 = 0, µ∗2 = 2, µ∗3 = 5, Σ∗1 = 0.1, Σ∗2 = 0.2 and Σ∗3 = 0.5. In addition, ψk,k ∗
= 0.6
∗
for all 1 ≤ k ≤ K and ψk,g = 0.2 for all 1 ≤ g ≤ K with g ̸= k. Figures 1, 2 and 3 compare the
marginal posterior distributions obtained using our variational inference approach to the ones using a
Gibbs sampler for the same model to examine whether the variational distribution obtained using our
proposed CAVI algorithm is a suitable approximation of the posterior distribution. We observe that,
for each parameter, the marginal variational distribution matches the Gibbs marginal posterior well,
particularly in terms of posterior mode, but has slightly less variance than the Gibbs posterior, which is
most noticeably for Σ. Overall the variational posterior provides a suitable estimate for the true, data
generating value and performs similarly to the Gibbs sampler. From Figure 4, the variational posterior
predictive distribution is shown to be close to the true likelihood functions for both continuous and
categorical data, indicating that the variational posterior is a good fit for the data.
3 Note that the r
ik does not sum to 1 over k with this initialisation, but the initial rik is not required to be probability
vectors.

9
 Bayesian Variational Inference for Mixed Data Mixture Models

posterior of µ, k=1 posterior of Σ, k=1

120
Gibbs 500 Gibbs
100 variational posterior variational posterior
true value 400 true value
80
Density

Density
300
60
200
40

20 100

0 0
−0.965 −0.960 −0.955 −0.950 −0.945 −0.940 0.022 0.023 0.024 0.025 0.026 0.027
Value Value

posterior of µ, k=2 posterior of Σ, k=2

80
Gibbs 250 Gibbs
variational posterior variational posterior
60 true value 200 true value

150
Density

Density
40
100
20
50

0 0
0.005 0.010 0.015 0.020 0.025 0.030 0.035 0.040 0.038 0.040 0.042 0.044 0.046 0.048 0.050 0.052
Value Value

posterior of µ, k=3 posterior of Σ, k=3

40 Gibbs 80 Gibbs
variational posterior variational posterior
30 true value 60 true value
Density

Density
20 40

10 20

0 0
1.45 1.46 1.47 1.48 1.49 1.50 1.51 0.110 0.115 0.120 0.125 0.130 0.135 0.140 0.145 0.150
Value Value

Figure 1: Plot of marginal posterior distribution of µ and Σ for both variational inference and Gibbs

posterior of Ψ, k=1 g=1 posterior of Ψ, k=2 g=1

Gibbs 40 Gibbs
variational posterior variational posterior
30
true value 30 true value
Density

Density

20
20

10 10

0 0
0.56 0.58 0.60 0.62 0.64 0.18 0.19 0.20 0.21 0.22 0.23 0.24
Value Value

posterior of Ψ, k=1 g=2 posterior of Ψ, k=2 g=2

Gibbs Gibbs
40 30
variational posterior variational posterior
true value true value
30
20
Density

Density

20

10
10

0 0
0.18 0.19 0.20 0.21 0.22 0.23 0.56 0.57 0.58 0.59 0.60 0.61 0.62 0.63
Value Value

posterior of Ψ, k=1 g=3 posterior of Ψ, k=2 g=3

Gibbs 40 Gibbs
40 variational posterior variational posterior
true value true value
30
30
Density

Density

20
20

10 10

0 0
0.17 0.18 0.19 0.20 0.21 0.22 0.16 0.17 0.18 0.19 0.20 0.21 0.22 0.23
Value Value

Figure 2: Plot of marginal posterior distribution of ψ1,g , ψ2,g for both variational inference and Gibbs

5.2 Simulation study

We now conduct an extensive simulation study for three, higher dimensional scenarios, for data of various
sample sizes n. In each scenario, we fix a set of true parameter values (µ∗ , Σ∗ , ψ ∗ , π ∗ ) and K ∗ components,
and simulate both the component allocation zi∗ as well as the observation (x∗i , c∗i ) for each data point
i = 1, . . . n. Full detail of the data generating parameters for the three scenarios can be found in the
supplement.

10
 Bayesian Variational Inference for Mixed Data Mixture Models

posterior of Ψ, k=3 g=1 posterior of π, k=1

60
Gibbs Gibbs
variational posterior 50 variational posterior
30
true value true value
40
Density

Density
20 30

20
10
10

0 0
0.16 0.17 0.18 0.19 0.20 0.21 0.22 0.23 0.39 0.40 0.41 0.42 0.43
Value Value

posterior of Ψ, k=3 g=2 posterior of π, k=2

60
Gibbs Gibbs
30 variational posterior 50 variational posterior
true value true value
40
Density

Density
20
30

20
10
10

0 0
0.18 0.19 0.20 0.21 0.22 0.23 0.24 0.25 0.31 0.32 0.33 0.34 0.35
Value Value

posterior of Ψ, k=3 g=3 posterior of π, k=3

30
Gibbs Gibbs
60
25 variational posterior variational posterior
true value true value
20
40
Density

Density
15

10
20
5

0 0
0.56 0.58 0.60 0.62 0.24 0.25 0.26 0.27 0.28
Value Value

Figure 3: Plot of marginal posterior distribution of ψ3,g , π for both variational inference and Gibbs

Distribution of x given c=0 Distribution of x given c=1

1.75 True p(x|c, µ∗ , Λ∗ , ψ ∗ , π ∗ ) True p(x|c, µ∗ , Λ∗ , ψ ∗ , π ∗ )
Variational posterior predictive q(x̃|c̃, x, c) 1.0 Variational posterior predictive q(x̃|c̃, x, c)
1.50 Empirical Data Empirical Data

1.25 0.8

1.00
Density

Density

0.6

0.75
0.4
0.50

0.2
0.25

0.00 0.0

−1.5 −1.0 −0.5 0.0 0.5 1.0 1.5 2.0 2.5 −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5 2.0 2.5
x x

posterior predictive of c
0.5
Distribution of x given c=2 variational posterior predictive q(c̃|x, c)
true p(c)
0.7 True p(x|c, µ∗ , Λ∗ , ψ ∗ , π ∗ ) empirical distribution of c
Variational posterior predictive q(x̃|c̃, x, c) 0.4
0.6 Empirical Data

0.5
0.3
Probability

0.4
Density

0.3 0.2

0.2
0.1
0.1

0.0
0.0
−1 0 1 2 Category 1 Category 2 Category 3
x

Figure 4: Comparison of the variational posterior predictive distribution of x given c, and c to the
empirical data and the likelihood function evaluated at the data generating parameters.

5.2.1 Simulation set-ups and error metrics

We consider the three following scenarios:
• In the first scenario, we simulate data with µ∗k that are pairwise well separated, and Σ∗k with
standard deviations that are small compared to the distance between the µ∗k . The categorical
variable has a relatively weak effect on the mixture [Link] idea of this scenario is to

11
 Bayesian Variational Inference for Mixed Data Mixture Models

test the model in the case where the mixtures are largely driven by separation across different
components in the continuous variables rather than categorical variables.
• For the second scenario, we instead use Σ∗k with standard deviations that are larger than the pairwise
distance between µ∗k , so data generated from the mixture components substantially overlap. In
contrast, the categorical variables are generated in a way to heavily correlate with the mixture
component assignments. The idea of this scenario is to test the model in the case where the
mixtures are largely driven by the categorical variables rather than the continuous variables.
• For the third scenario, we simulate data using parameters that is a middle ground between sce-
narios 1 and 2. Specifically, we use µ∗k with separation that are similar in magnitude compared to
the standard deviations of Σ∗k , causing some continuous data overlap. We define a single categor-
ical variable which is significantly correlated with the mixture component z, but not as strongly
correlated compared to scenario 2.
In order to assess the empirical performance of our algorithm, we employ several performance mea-
sures. First, to test parameter estimation, we calculate for π, µ, Σ, ψ the mean absolute error between
the posterior mean of the variational posterior and the true parameter values:
∗ ∗
K K
1 X ∗ 1 X ∗ Φ̂k
Errorµ = ∥µk − m̂k ∥1 ErrorΣ = 2 ∗ Σk − (31)
qK ∗ q K ν̂k − q − 1
k=1 k=1 1
∗ dj ∗
K K Xp X
1 X ∗ α̂k 1 X
∗ η̂k,j,gj
Errorπ = ∗ |πk − PK | Errorψ = Pp ∗
|ψk,j,g − Pdj |
K k=1 α̂k
( d
j=1 j )K j=1 g =1
j
′ η̂k,j,g′
k=1 k=1 j gj =1 j

Secondly, to test clustering capability, we compare for each data point the true component label zi∗
with the highest variational posterior probability (‘hard clustering’), given by arg maxk rik . We use
P ropz to denote the proportion of data points for which the component label matches the true com-
ponent label. Third, to assess the density estimation, we sample additional sample size test data
ntest = min(0.4n, 2000) and calculate for each data point (x̃i , c̃i ) in the test dataset, the mean ab-
solute error between the loglikelihood for Equation (1) evaluated at the true data generating parameters
(µ∗ , Λ∗ , ψ ∗ , π ∗ ), and log of the posterior predictive distribution in Equation (17) using the final parameter
estimates obtained from running Algorithm 1 in the training set. We calculate the average of this error
over the test set Errorlogppd :
nX
test
1
Errorlogppd = |ln(p(x̃i , c̃i |µ∗ , Λ∗ , ψ ∗ , π ∗ )) − ln(q(x̃i , c̃i |x, c))| (32)
ntest i=1

Lastly, to assess the uncertainty quantification over the parameters (µ, Σ, ψ, π), we examine the fre-
quentist coverage of the posterior highest density intervals (HDI) of µ, Σ, ψ, π by counting the proportion
of 95% variational posterior HDI containing the data generating parameter values over 100 runs of the
variational inference procedure for independent datasets. Conveniently, the variational posterior density
q(µ, Λ, ψ, π) as well as the marginal variational posterior densities q(µ), q(Λ), q(ψ), and q(π) are all
available in closed form, which allows us to easily calculate the HDI and check the frequentist coverage
of µ, Σ, ψ, π jointly as well as individually.
Note that the model likelihood in Equation (1) is invariant under any permutation σ on 1, . . . , K
K
X p
Y K
X p
Y
πσ(k) N (xi |µσ(k) , Λ−1
σ(k) ) ψσ(k),j,cij = πk N (xi |µk , Λ−1
k ) ψk,j,cij (33)
k=1 j=1 k=1 j=1

When the model is correctly specified, this non identifiability of the likelihood usually causes the estimated
model parameters from Algorithm 1 to correspond to the true parameters up to a permutation. Because
of this, when reporting errors metrics listed or checking the frequentist coverage of posterior HDI on
simulated data, we match the permuted parameters with the true parameters via linear sum assignment
between the estimated m̂k and the true µ∗k , and we use the resulting permutation σ to relabel Λ, ψ, and
π.
We calculate the aforementioned error metrics and report the average over 100 runs for the three
scenarios, for sample sizes n = 2500, 5000, 10000, 20000 assuming that the number of component K = K ∗

12
 Bayesian Variational Inference for Mixed Data Mixture Models

is known a priori. For n = 5000, we also compare the results to the case where the inference procedure
is run with K = 2K ∗ to see how the algorithm performs when the number of clusters is misspecified. To
assess the competitiveness of our method, we compare the point estimation errors with the Expectation-
Maximization (EM) algorithm of Hunt and Jorgensen [1999], which unlike our method does not provide
uncertainty quantification. We also compare both the point estimation error and coverage with MCMC,
in this case a Gibbs sampler. MCMC is in a sense the gold standard in Bayesian inference as the empirical
distribution of its samples converges to the true posterior with infinite computational power. For the
Gibbs sampler, the posterior mean is used as the point estimate when calculating the errors, and for
P ropz , the last sample of z is used. Errorlogppd is not calculated for the Gibbs sampler as the posterior
predictive density is not available in closed form, and for coverage only the overall HDI is calculated,
using the joint density instead of the posterior as it is proportional to the posterior, while the true
posterior is unavailable in closed form.

5.2.2 Summary of simulation study results

Errorµ ErrorΣ Errorψ Errorπ P ropz Errorlogppd

n = 2500, VI 0.0138 0.00759 0.0158 0.00636 0.999 0.190
n = 5000, VI 0.00977 0.00515 0.0110 0.00447 0.999 0.127
n = 5000, K = 2K ∗ , VI 0.00988 0.00520 0.0110 0.00444 0.998 0.128
n = 5000, Gibbs 0.00978 0.00515 0.0110 0.00447 0.999
n = 5000, EM 0.00974 0.00496 0.0110 0.00447 0.999
n = 10000, VI 0.00722 0.00361 0.00792 0.00306 0.999 0.0873
n = 20000, VI 0.00506 0.00252 0.00566 0.00218 0.999 0.0596

Table 1: Various error metrics averaged over 100 runs for various sample sizes n, scenario 1

overall π Σ µ ψ
n = 2500, VI 0.00 0.92 0.00 0.97 0.92
n = 5000, VI 0.08 0.97 0.02 0.96 0.98
n = 5000, Gibbs 0.49
n = 10000, VI 0.37 0.95 0.2 0.91 0.96
n = 20000, VI 0.61 0.97 0.57 0.92 0.96

Table 2: Frequentist coverage of 95% HDI over 100 runs for various sample sizes n, scenario 1

For scenario 1, from Table 1, it can be seen that our approach accurately estimates the parameters
(µ, Σ, ψ, π), and the error decreases as the sample size n increases for all parameters. Our approach
manages to predict the true mixture components for almost all the data points, which is perhaps un-
surprising as the data is well separated in continuous space. The true log-likelihood is reasonably well
estimated by the log posterior predictive density, which, for comparison with the size of the error, has an
average absolute value of 4.6 in the test datasets. The errors for the case with the misspecified number
of clusters K are only marginally higher than the case with correctly specified K with the same sample
size, with P ropz still very close to 1, indicating that the variational posterior is still able to converge to
the true number of clusters. In general, our approach performs very similar to the EM algorithm of Hunt
and Jorgensen [1999] as well as the Gibbs sampler in terms of the point estimation error. From Table 2,
we observed that the frequentist coverage is close to 95% for all parameters except Σ, which nevertheless
achieves reasonable coverage as sample size is increased to 20000. This could be due to Σ having higher
dimension compared to the other parameters, making it more difficult to estimate. Similarly, the overall
coverage appears to be bottlenecked by the coverage of Σ. The Gibbs sampler achieves considerably
higher overall coverage for the case N = 5000, which is perhaps not unexpected, as variational inference
approaches are known to underestimate variance, although still considerably lower than 95
For scenario 2, from Table 3 we observe that our approach again achieves good point estimation of
the parameters (µ, Λ, ψ, π) for scenario 2, with the error decreasing as the sample size n increases for all
parameters. The proportion of data points with the correctly predicted hard component is somewhat
lower than scenario 1, which may be due to the highly overlapping nature of the data across different

13
 Bayesian Variational Inference for Mixed Data Mixture Models

Errorµ ErrorΣ Errorψ Errorπ P ropz Errorlogppd

n = 2500, VI 0.0398 0.0409 0.0136 0.00675 0.941 0.205
n = 5000, VI 0.0287 0.0289 0.00955 0.00504 0.942 0.142
n = 5000, K = 2K ∗ , VI 0.0286 0.0290 0.00959 0.00505 0.941 0.144
n = 5000, Gibbs 0.0287 0.0289 0.00953 0.005 0.920
n = 5000, EM 0.0287 0.0290 0.00955 0.00498 0.942
n = 10000, VI 0.0197 0.0208 0.00677 0.00338 0.943 0.0995
n = 20000, VI 0.0135 0.0147 0.00467 0.00265 0.944 0.0694

Table 3: Various error metrics averaged over 100 runs for various sample sizes n, scenario 2

overall π Σ µ ψ
n = 2500, VI 0.23 0.90 0.57 0.85 0.47
n = 5000, VI 0.38 0.84 0.73 0.83 0.54
n = 5000, Gibbs 0.63
n = 10000, VI 0.39 0.9 0.78 0.84 0.51
n = 20000, VI 0.47 0.86 0.75 0.91 0.55

Table 4: Frequentist coverage of 95% HDI over 100 runs for various sample sizes n, scenario 2

components in continuous space, but is still very high. The log posterior predictive density has absolute
value 11.6 on the test dataset on average, so the log posterior predictive distribution is able to estimate
the true log likelihood with lower error on average compared to scenario 1. Like in scenario 1, when
n = 5000, the errors for our VI approach is almost identical to the error for the EM algorithm and for
the Gibbs sampler showing our approach performs very competitively in terms of point estimation. In
addition the errors in the case when the number of clusters K is misspecified are very similar to the
correctly specified case.
In terms of frequentist coverage, from Table 4, it can be seen that the coverage of π and µ are close to
95%, albeit somewhat worse compared to scenario 1. Interestingly, Σ has better coverage than scenario
1, but ψ has considerably worse coverage, which could be due to ψ having much higher dimension in
scenario 2 compared to scenario 1. Again the Gibbs sampler achieves higher overall frequentist coverage
as expected.

Errorµ ErrorΣ Errorψ Errorπ P ropz Errorlogppd

n = 2500, VI 0.0266 0.01935 0.0111 0.00660 0.968 0.167
n = 5000, VI 0.0185 0.0137 0.00793 0.00487 0.969 0.116
n = 5000, K = 2K ∗ , VI 0.0192 0.0138 0.00818 0.00566 0.965 0.119
n = 5000, Gibbs 0.0186 0.0137 0.00795 0.00486 0.954
n = 5000, EM 0.0186 0.0138 0.00795 0.00486 0.969
n = 10000, VI 0.0132 0.00966 0.00553 0.00323 0.970 0.0808
n = 20000, VI 0.00950 0.00682 0.00396 0.00224 0.969 0.0574

Table 5: Various error metrics averaged over 100 runs for various sample sizes n, scenario 3

overall π Σ µ ψ
n = 2500, VI 0.48 0.90 0.49 0.87 0.85
n = 5000, VI 0.52 0.92 0.58 0.88 0.86
n = 5000, Gibbs 0.60
n = 10000, VI 0.68 0.90 0.74 0.82 0.90
n = 20000, VI 0.64 0.95 0.73 0.87 0.86

Table 6: Frequentist coverage of 95% HDI over 100 runs for various sample sizes n, scenario 3

For scenario 3, from Table 5, as for the previous scenarios we observe that good point estimation of
the parameters (µ, Λ, ψ, π) is achieved with the error decreasing as the sample size n increases for all
parameters. The proportion of data points with the correctly predicted hard mixture component is close

14
 Bayesian Variational Inference for Mixed Data Mixture Models

to 100%, and the error in log posterior predictive density is low compared to the average absolute log
posterior predictive density on the test set of 7.3. The errors for the misspecified K case is somewhat
higher but still close to the correctly specified K case, EM and the Gibbs sampler, which agree strongly.
In terms of frequentist coverage, from Table 6 the posterior HDI is able to achieve reasonably good
coverage for all parameters, though Σ like in scenario 1 is the most difficult parameter to achieve good
coverage. The difference in overall coverage between VI and Gibbs is also the smallest in this case.

5.3 Application: Identifying groups of individuals with comparable health

phenotypes based on survey responses
We apply our VI algorithm on the National Health and Nutrition Examination Survey (NHANES)
dataset, an ongoing program conducted by the National Center for Health Statistics (NCHS) to monitor
the health and nutritional status of the U.S. population. For this analysis we take a subset of the
NHANES dataset from 1988 to 2018, for middle aged males aged between 40 and 59, as both sex and
age group have strong epidemiological effects that can strongly drive cluster make up, resulting in less
epidemiologically interesting clusters. We remove rows containing missing data, and details of the data
cleaning procedure can be found in Lhoste et al. [2024]. Each row of the dataset represents a participant
of the survey and contains a range of continuous and categorical risk factors including anthropometric
measurements, blood pressure and heart rate, blood sugar level, cholesterol level and kidney function
measurements. Specifically, for anthropometric measurements, we have: BMI, defined as a participant’s
weight divided by height squared (kgm−2 ); height (cm) ; WHtR, defined as waist circumference divided
by height, which is used as a measure of abdominal obesity. For blood pressure and heart rate, we have
systolic blood pressure (SBP), diastolic blood pressures (DBP) and pulse rate, which is defined as the
resting heart rate over one minute. For cholesterol level, we have both HDL and non-HDL cholesterol
levels, where non-HDL is associated with higher risk of ischemic heart disease and stroke. Lastly, we have
HbA1c which is a measurement of glucose levels in the blood in recent weeks and glomerular filtration
rate (eGFR) as a measure of kidney function, where higher eGFR indicates better kidney function. For
categorical variables, we include the participant’s smoking status, classified as current, former, or never
smoker.
We investigate whether or not our proposed approach can produce epidemiologically heterogeneous
and meaningful clusters that could then be used for planning prevention and clinical care, and to monitor
population health over time. For the model likelihood, we use Equation (1), and we make the same prior
assumptions and hyperparameter choice as in Section 2. We run the proposed variational inference
procedure with K = 10 as this was the number of clusters identified as most suitable by Lhoste et al.
[2024] when analysing the continuous part of the same dataset with a K-means algorithm.
To examine the features of the clusters produced by our method, we use the variational posterior
predictive marginal distributions, as described earlier in Equation (18), because they are both suitable
for visualisation in high-dimensional settings and effectively summarize the composition of each cluster
while accounting for posterior uncertainty. For each cluster, the variational posterior predictive marginal
distributions are represented using radar plots for the continuous data in Figure 5, and using bar plots
for the categorical data in Figure 6.
In Figure 5, one can see that clusters 1 and 2 have very small posterior mixture probabilities, meaning
that the model has effectively settled on less clusters than originally specified, indicating some ability
to compensate for the over specification of K even on real data. Cluster 5 is the largest cluster (in the
sense of posterior mean of πk ), and is representative of a median adult, as the posterior predictive means
of most continuous variables are close to the empirical median albeit with somewhat lower pulse, HbA1c
and SBP, and its posterior predictive distribution for smoking status also does not differ greatly from
the median. Clusters 9 and 10 also have posterior predictive means close to the empirical median for
most continuous features, but distinguishes from cluster 5 and each other by cluster 9 having somewhat
elevated BMI and healthier cholesterol levels (HDL and non-HDL) with mostly non smokers, and cluster
10 having elevated cholesterol levels and predominantly current smokers.
Clusters 4, 6 and 8 have posterior predictive means at close to or exceeding 75% quantile for multiple
continuous features, and few if any below 50% quantile. Clusters 4 and 6 are both characterized by very
high posterior means for HbA1c, suggesting they are characterized by diabetic people. However, they
differ in anthropometric measures and, to some extent, lipid profiles: Cluster 6 shows markedly elevated

15
 Bayesian Variational Inference for Mixed Data Mixture Models

Radar Plots for marginal posterior predictive means and 80.0% credible interval in empirical distribution probability, per cluster
Cluster 1 ,0.000807 Cluster 2 ,0.000578 Cluster 3 ,0.0833 Cluster 4 ,0.0742 Cluster 5 ,0.228
hba1c WHtR hba1c WHtR hba1c WHtR hba1c WHtR hba1c WHtR

hdl height hdl height hdl height hdl height hdl height
100% 100% 100% 100% 100%
75% 75% 75% 75% 75%
50% 50% 50% 50% 50%
25% 25% 25% 25% 25%
non_hdl 0% bmi non_hdl 0% bmi non_hdl 0% bmi non_hdl 0% bmi non_hdl 0% bmi

sbp pulse sbp pulse sbp pulse sbp pulse sbp pulse

dbp eGFR dbp eGFR dbp eGFR dbp eGFR dbp eGFR
Cluster 6 ,0.0534 Cluster 7 ,0.127 Cluster 8 ,0.109 Cluster 9 ,0.136 Cluster 10 ,0.189
hba1c WHtR hba1c WHtR hba1c WHtR hba1c WHtR hba1c WHtR

hdl height hdl height hdl height hdl height hdl height
100% 100% 100% 100% 100%
75% 75% 75% 75% 75%
50% 50% 50% 50% 50%
25% 25% 25% 25% 25%
non_hdl 0% bmi non_hdl 0% bmi non_hdl 0% bmi non_hdl 0% bmi non_hdl 0% bmi

sbp pulse sbp pulse sbp pulse sbp pulse sbp pulse

dbp eGFR dbp eGFR dbp eGFR dbp eGFR dbp eGFR
Mean 10.0 to 90.0% Quantile

Figure 5: Radar plots of the marginal posterior predictive distribution of each continuous variable,
conditioned on cluster component. The posterior mean of each πk is also added atop each radar plot. The
concentric circles and quantiles in each radar plot represents the quantiles of the empirical distribution of
the whole data. For the sake of visual consistency, the distribution of height, eGFR and hdl are reversed
so that higher values in the radar plot for all continuous variables correspond to higher risk of chronic
diseases.

k=1 k=2 k=3 k=4 k=5

1.0
0.889
0.8
0.687
0.6
Probability

Probability

Probability

Probability

Probability 0.544
0.415 0.385
0.4 0.346 0.347
0.292 0.293 0.269
0.2 0.180
0.132 0.109
0.056 0.056
0.0
0 1 2 0 1 2 0 1 2 0 1 2 0 1 2
k=6 k=7 k=8 k=9 k=10
1.0

0.8

0.6 0.571
Probability

Probability

Probability

Probability

Probability

0.496
0.450 0.422 0.438
0.4 0.319 0.323 0.338
0.255 0.240 0.278
0.231 0.224 0.227
0.189
0.2

0.0
0 1 2 never0 smoked1 2former smoker0 1current 2smoker 0 Empirical
1 proportion
2 0 1 2

Figure 6: Barplot of the marginal posterior predictive distribution of smoking status, conditioned on
cluster component.

posterior means for WHtR and BMI, along with relatively high HDL, whereas cluster 4 exhibits posterior
means for WHtR and BMI that are closer to the population median but still somewhat elevated, and close
to median level of HDL. Cluster 8 is characterized by very high posterior means for systolic and diastolic
blood pressure, BMI, and WHtR, but near-median levels of HbA1c and cholesterol, distinguishing it from
clusters 4 and 6. Smoking distributions for all three clusters are similar.

16
 Bayesian Variational Inference for Mixed Data Mixture Models

Clusters 7 and 3 are characterized by below-median values for most continuous variables, reflecting
generally healthier profiles. Cluster 3 differs from cluster 7 by having higher posterior means for systolic
and diastolic blood pressure and pulse. Individuals in cluster 3 are also more likely to be current smokers,
whereas smoking status in cluster 7 is more evenly distributed, with no clear dominance of one category.

6 Discussion
We have developed a scalable, coordinate ascent variational inference algorithm for mixed data mixture
models that allows both the continuous and categorical variables to influence the make up of each
component or cluster. We have conducted simulation study under scenarios where the heterogenenity
in the mixture components are driven by the continuous data, the categorical data, or a combination
of both respectively. Through these simulation studies, we demonstrated that our proposed approach is
able to conduct good density estimation, and estimate the underlying parameters well and competitively
with Gibbs sampler and an existing EM algorithm, even when the number of mixtures is incorrectly
specified. While our method performs somewhat worse than the Gibbs sampler in terms of frequentist
coverage, it is nevertheless able to produce adequate levels of overall coverage, and most of the marginal
95% HDI are able to higher levels of coverage, often close to 95%. Our proposed VI approach is also
more computationally scalable than the Gibbs sampler.
We demonstrated the real world utility of our approach by applying it on a challenging dataset of risk
factor data of US participants in the NHANES study. We showed our method can produce epidemiolog-
ically interpretable clusters, which is crucial when clustering epidemiological datasets. We visualized the
heterogeneity between the different clusters using marginal posterior predictive distributions, which takes
into account posterior uncertainty, something that is often missed. However, while the medians of each
continuous marginal posterior predictive distribution is distinguishable across clusters, the uncertainty
quantification becomes very conservative for large probability credible intervals, causing the uncertainty
quantification across different clusters to be less distinguishable.
Furthermore, we established a theoretical convergence result for our method, showing that the pos-
terior mean from the CAVI algorithm converge to the true parameter value as the sample size n tends
to infinity. This provides formal justification for the validity of the algorithm.
One future research direction could include developing VI approaches that enforces more separa-
tion between the clusters in continuous space, which could produce more distinguishable uncertainty
quantification across different clusters. Another research direction could be to extend the approach to
infinite mixtures such as Dirichlet Processes, making it easier to select a suitable number of mixture
components.

Acknowledgements This work was supported by the UK Medical Research Council (MRC), under
the grant ‘Worldwide phenotypes and transitions in obesity-related multimorbidity’, reference number
MR/V034057/1. The authors are grateful to Majid Ezzati and Bin Zhou for useful discussions.

References
C. Ahlmann-Eltze and C. Yau. Mixdir: Scalable bayesian clustering for high-dimensional categorical
data. In 2018 IEEE 5th International Conference on Data Science and Advanced Analytics (DSAA),
pages 526–539, 2018. doi: 10.1109/DSAA.2018.00068.
C. Alexander and E. Lazar. Normal mixture garch(1,1): applications to exchange rate modelling. Journal
of Applied Econometrics, 21(3):307–336, 2006. doi: [Link]
C. Beecks, M. S. Uysal, and T. Seidl. Content-based image retrieval with gaussian mixture models. In
X. He, S. Luo, D. Tao, C. Xu, J. Yang, and M. A. Hasan, editors, MultiMedia Modeling, pages 294–305,
Cham, 2015. Springer International Publishing. ISBN 978-3-319-14445-0.
R. Bhatia. Matrix analysis. Springer New York, 1997.
C. M. Bishop. Pattern Recognition and Machine Learning (Information Science and Statistics). Springer-
Verlag, Berlin, Heidelberg, 2006. ISBN 0387310738.
D. M. Blei and M. I. Jordan. Variational inference for Dirichlet process mixtures. Bayesian Analysis, 1
(1):121 – 143, 2006. doi: 10.1214/06-BA104.

17
 Bayesian Variational Inference for Mixed Data Mixture Models

D. M. Blei, A. Kucukelbir, and J. D. McAuliffe. Variational inference: A review for statisticians. Journal
of the American Statistical Association, 112(518):859–877, 2017. doi: 10.1080/01621459.2017.1285773.
G. Celeux, M. Hurn, and C. P. Robert. Computational and inferential difficulties with mixture posterior
distributions. Journal of the American Statistical Association, 95(451):957–970, 2000. doi: 10.1080/
01621459.2000.10474285.
G. Celeux, C. P. Robert, and S. Fruhwirth-Schnatter, editors. Handbook of mixture analysis. CRC Press,
London, England, Dec. 2020.
G. B. Durham. Sv mixture models with application to s&p 500 index returns. Journal of Financial
Economics, 85(3):822–856, 2007. ISSN 0304-405X. doi: [Link]
S. Frühwirth-Schnatter. Markov chain monte carlo estimation of classical and dynamic switching and
mixture models. Journal of the American Statistical Association, 96(453):194–209, 2001. doi: 10.1198/
016214501750333063.
M. C. Hughes and E. Sudderth. Memoized online variational inference for dirichlet process mixture
models. In C. Burges, L. Bottou, M. Welling, Z. Ghahramani, and K. Weinberger, editors, Advances
in Neural Information Processing Systems, volume 26. Curran Associates, Inc., 2013.
L. Hunt and M. Jorgensen. Theory & methods: Mixture model clustering using the multimix program.
Australian & New Zealand Journal of Statistics, 41(2):154–171, 1999. doi: [Link]
1467-842X.00071.
L. Hunt and M. Jorgensen. Mixture model clustering for mixed data with missing information. Com-
putational Statistics & Data Analysis, 41(3):429–440, 2003. ISSN 0167-9473. doi: [Link]
1016/S0167-9473(02)00190-1. Recent Developments in Mixture Model.
V. Huynh, D. Phung, and S. Venkatesh. Streaming variational inference for dirichlet process mixtures. In
G. Holmes and T.-Y. Liu, editors, Asian Conference on Machine Learning, volume 45 of Proceedings
of Machine Learning Research, pages 237–252, Hong Kong, 20–22 Nov 2016. PMLR.
A. Jasra, C. C. Holmes, and D. A. Stephens. Markov Chain Monte Carlo Methods and the Label
Switching Problem in Bayesian Mixture Modeling. Statistical Science, 20(1):50 – 67, 2005. doi:
10.1214/088342305000000016.
C. J. Lawrence and W. J. Krzanowski. Mixture separation for mixed-mode data. Statistics and Com-
puting, 6(1):85–92, Mar 1996. ISSN 1573-1375. doi: 10.1007/BF00161577.
V. P. F. Lhoste, B. Zhou, A. Mishra, J. E. Bennett, S. Filippi, P. Asaria, E. W. Gregg, G. Danaei,
and M. Ezzati. Cardiometabolic and renal phenotypes and transitions in the united states pop-
ulation. Nature Cardiovascular Research, 3(1):46–59, Jan 2024. ISSN 2731-0590. doi: 10.1038/
s44161-023-00391-y.
D. McParland and I. C. Gormley. Model based clustering for mixed data: clustmd. Advances
in Data Analysis and Classification, 10(2):155–169, Jun 2016. ISSN 1862-5355. doi: 10.1007/
s11634-016-0238-x.
J. W. Miller and M. T. Harrison. Mixture models with a prior on the number of components. Journal
of the American Statistical Association, 113(521):340–356, 2018. doi: 10.1080/01621459.2016.1255636.
PMID: 29983475.
M. Niknejad, H. Rabbani, and M. Babaie-Zadeh. Image restoration using gaussian mixture models with
spatially constrained patch clustering. IEEE Transactions on Image Processing, 24(11):3624–3636,
2015. doi: 10.1109/TIP.2015.2447836.
J. K. Pritchard, M. Stephens, and P. Donnelly. Inference of population structure using multilocus
genotype data. Genetics, 155(2):945–959, 06 2000. ISSN 1943-2631. doi: 10.1093/genetics/155.2.945.
M. Ranalli and R. Rocci. Mixture models for mixed-type data through a composite likelihood approach.
Computational Statistics & Data Analysis, 110:87–102, 2017. ISSN 0167-9473. doi: [Link]
10.1016/[Link].2016.12.016.
A. Rau, C. Maugis-Rabusseau, M.-L. Martin-Magniette, and G. Celeux. Co-expression analysis of high-
throughput transcriptome sequencing data with poisson mixture models. Bioinformatics, 31(9):1420–
1427, 01 2015. ISSN 1367-4803. doi: 10.1093/bioinformatics/btu845.
S. Richardson and P. J. Green. On bayesian analysis of mixtures with an unknown number of components
(with discussion). Journal of the Royal Statistical Society: Series B (Statistical Methodology), 59(4):
731–792, 1997. doi: [Link]
Y. W. Teh, M. I. Jordan, M. J. Beal, and D. M. Blei. Hierarchical dirichlet processes. Journal of the
American Statistical Association, 101(476):1566–1581, 2006. doi: 10.1198/016214506000000302.

18
 Bayesian Variational Inference for Mixed Data Mixture Models

D. M. Titterington and B. Wang. Convergence properties of a general algorithm for calculating variational
Bayesian estimates for a normal mixture model. Bayesian Analysis, 1(3):625 – 650, 2006. doi: 10.
1214/06-BA121.
C. Wang, J. Paisley, and D. M. Blei. Online variational inference for the hierarchical dirichlet process. In
G. Gordon, D. Dunson, and M. Dudı́k, editors, Proceedings of the Fourteenth International Conference
on Artificial Intelligence and Statistics, volume 15 of Proceedings of Machine Learning Research, pages
752–760, Fort Lauderdale, FL, USA, 11–13 Apr 2011. PMLR.
S. Wang, J. G. Yabes, and C.-C. H. Chang. A bayesian finite mixture model with variable selection for
data with mixed-type variables, 2019.
Y. Zhang and K. Berhane. Bayesian mixed hidden markov models: a multi-level approach to modeling
categorical outcomes with differential misclassification. Statistics in Medicine, 33(8):1395–1408, 2014.
doi: [Link]
Z. Zhang, K. L. Chan, Y. Wu, and C. Chen. Learning a multivariate gaussian mixture model with the
reversible jump mcmc algorithm. Statistics and Computing, 14(4):343–355, Oct 2004. ISSN 1573-1375.
doi: 10.1023/B:STCO.0000039484.36470.41.

19
 Bayesian Variational Inference for Mixed Data Mixture Models

Supplement to ‘Bayesian Variational Inference for Mixed Data

Mixture Models’
This supplement consists of several parts. Appendix A contains the derivation for the variational posterior
update equations of the CAVI algorithm in the main text of the paper. Appendix B contains the ELBO
calculations for the CAVI algorithm. Appendix C contains the Gibbs sampler updates for the model
considered in the paper. The proof of Theorem 1 in the main text of paper is given in Appendix D,
and auxiliary calculations needed for said proof in Appendix D.2. Appendix E contains detail of the
parameters used for the three scenarios in the Simulation study under the Experimental Results section
of the main text of the paper.

A Coordinate Ascent Variational Inference update equations

In this section we provide the CAVI (Coordinate Ascent Variational Inference) update equations for our
mixed data mixture model.
Recall the prior distributions are as follows:
µk ∼ N (mk , β −1 Λ−1 k ).
Λk ∼ W(ν, Φ−1 ) or equivalently Σk ∼ IW(ν, Φ), where Σk = (Λk )−1
π = (π1 , π2 , . . . , πK ) ∼ Dir(α, . . . , α)
ψk,j = (ψk,j,1 , ψk,j,2 , . . . , ψk,j,dj ) ∼ Dir(ηj , . . . , ηj ).
The model (likelihood) terms can be written as follows:
Qn Qn QK
x|c, z, µ, Λ ∼ i=1 N (xi |µzi , Λ−1 ) = i=1 k=1 N (xi |µk , Λ−1
k )
zik
Qn Qn QK zziik
z|π = i=1 πzi = i=1 k=1 πk
Qn Qp Qn Qp QK zik
c|ψ, z = i=1 j=1 ψzi ,j,cij = i=1 j=1 k=1 ψk,j,c ij
where zik = 1 if zi = k, zik = 0 otherwise.
The likelihood with z marginalised out (for a single datapoint) is:

K
X p
Y
p(xi , ci |µ, Λ, ψ, π) = πk N (xi |µk , Λ−1
k ) ψk,j,cij (34)
k=1 j=1

This leads to a joint distribution of the form:

p(x, c, z, µ, Λ, ψ, π) = p(x|c, z, µ, Λ)p(c|ψ, z)p(z|π)p(µ|Λ)p(Λ)p(ψ)p(π)

= p(x|c, z, µ, Λ)p(µ|Λ)p(Λ)p(c|ψ, z)p(ψ)p(z|π)p(π) (35)

We assume the variational posterior is of the form q(µ, Λ, ψ, π, z) = q(µ, Λ, ψ, π)q(z). The exact
distributional families of the variational posterior need not to be specified and will come out during
the derivation. To derive its form we use CAVI with two blocks, so we alternate between updating
q(µ, Λ, ψ, π) and q(z).

q(µ, Λ, ψ, π) ∝ exp(E−µ,−Λ,−ψ,−π lnp(x, c, z, µ, Λ, ψ, π) ) (36)

∝ exp(Ez [ln(p(x|c, z, µ, Λ)p(µ|Λ)p(Λ)p(c|ψ, z)p(ψ)p(z|π)p(π))])
∝ exp(Ez [ln(p(x|c, z, µ, Λ)p(µ|Λ)p(Λ))]) exp(Ez [ln(p(z|π)p(π))]) exp(Ez [ln(p(c|ψ, z)p(ψ))])
| {z }| {z }| {z }
∝q(µ,Λ) ∝q(π) ∝q(ψ)

From this it’s clear q(µ, Λ, ψ, π) = q(µ, Λ)q(ψ)q(π), each term in this product can therefore be derived
separately as follows:
For q(π):

20
 Bayesian Variational Inference for Mixed Data Mixture Models

q(π) ∝ exp(Ez [lnp(z|π) + lnp(π)])

n Y
Y K
∝ exp(Ez [ln( πkzik ) + ln(Dir(π|α, . . . , α))])
i=1 k=1
n X
X K K
X
∝ exp(Ez [ zik lnπk + (α − 1)ln(πk )])
i=1 k=1 k=1
K
X n
X
∝ exp( (α − 1 + rik )lnπk )
k=1 i=1
K
(α−1+ n
P
i=1 rik )
Y
∝ πk (37)
k=1
Pn Pn
which by observation is a Dirichlet distribution Dir(α + i=1 ri1 , . . . , α + i=1 riK ), where rik =
Ezi [zik ] = q(zi = k), the posterior probability of the ith data being in cluster k under the variational
posterior probability, which will be obtained in the CAVI update for q(z) later on.
For q(µ, Λ):

q(µ, Λ) ∝ exp(Ez [lnp(x|c, z, µ, Λ) + lnp(µ|Λ) + lnp(Λ)])

n Y
Y K K
Y K
Y
∝ exp(Ez [ln( N (xi |µk , Λ−1
k )
zik
)] + ln( N (µk |mk , β −1 Λ−1
k )) + ln( W(Λk |ν, Φ−1 )))
i=1 k=1 k=1 k=1
n X
X K K
X K
X
∝ exp( rik ln(N (xi |µk , Λ−1
k )) + ln(N (µk |mk , β −1 Λ−1
k )) + ln(W(Λk,g |ν, Φ−1 )))
i=1 k=1 k=1 k=1
K
Y Xn
∝ exp( rik ln(N (xi |µk , Λ−1
k )) + ln(N (µk |mk , β
−1 −1
Λk )) + ln(W(Λk |ν, Φ−1 ))) (38)
k=1 i=1

From here we first derive a conjugate Gaussian posterior for each µk |Λk , and subsequently the leftover
terms can be ‘absorbed’ into the Wishart prior to form a conjugate Wishart posterior for Λk . For each
µk , the relevant terms are:

21
 Bayesian Variational Inference for Mixed Data Mixture Models

n
nX o
q(µk |Λk ) ∝ exp rik ln(N (xi |µk , Λ−1
k )) + ln(N (µk |mk , β
−1 −1
Λk ))
i=1
n
nX 1 1 1
∝ exp rik [ ln(det(Λk )) − (xi − µk )T Λk (xi − µk )] + ln(det(βΛk ))
i=1
2 2 2
1 o
− (µk − mk )T βΛk (µk − mk )
2
n 1 Xn X n
∝ exp − [µTk (β + rik )Λk µk + (−βmTk − rik xTi )Λk µk )
2 i=1 i=1
n
X n
X
+ µTk Λk (−βmk − rik xi ) + mTk βΛk mk + xTi rik Λk xi ]
i=1 i=1
n
X 1 1 o
+ rik ( ln(det(Λk )) + ln(det(βΛk ))
i=1
2 2
n 1 n
1 X
∝ exp ( ln(det(βΛk )) − (µk − m̂k )T (β + rik )Λk (µk − m̂k )
2 2 i=1
n n n
1 T
X
T T
X X o
− [−m̂k (β + rik )Λk m̂k + mk βΛk mk + xi rik Λk xi − rik ln(det(Λk ))] (39)
2 i=1 i=1 i=1

βmk + n
P
r xj
Where m̂k = Pnj=1 jk .
The first two terms in the last expression of the derivation above forms
β+ j=1 rjk
Pn
a conjugate posterior µk |Λk ∼ N (µk |m̂k , β̂k−1 Λ−1 −1
k ), where β̂k = (β + j=1 rjk )
−1
(the determinant term
is correct up to a constant that does not depend on Λk ). The second line in the last expression contains
the leftover terms that need to be ‘absorbed’ by the Wishart prior, which we will see in the derivation
for q(Λk ) as follows:
For q(Λk ):

n
n 1 X
q(Λk ) ∝ exp ln(W(Λk |ν, Φ−1 )) − [−m̂Tk (β + rik )Λk m̂k + mTk βΛk mk
2 i=1
n
X n
X o
+ xTi rik Λk xi − rik ln(det(Λk ))]
i=1 i=1
nν − q − 1 n
1 1 X
∝ exp ln(det(Λk )) − T r(ΦΛk ) − [−m̂Tk (β + rik )Λk m̂k + mTk βΛk mk
2 2 2 i=1
n
X n
X o
+ xTi rik Λk xi − rik ln(det(Λk ))]
i=1 i=1
n ν − q − 1 + Pn r n
i=1 ik 1  X
∝ exp ln(det(Λk )) − T r ΦΛk − (β + rik )m̂k m̂Tk Λk + βmk mTk Λk
2 2 i=1
n
X o
+ rik xi xTi Λk (40)
i=1

where in the last line we used the cyclic propertyPof trace. The last line is also in the form of a
n
Wishart distribution with degree of freedom ν̂k = ν + i=1 rik and scale matrix Φ̂−1
k , where
n
X n
X
Φ̂k = Φ − (β + rik )m̂k m̂Tk + βmk mTk + rik xi xTi
i=1 i=1

Lastly we have the update for q(ψ):

22
 Bayesian Variational Inference for Mixed Data Mixture Models

n o
q(ψ) ∝ exp Ez [lnp(c|ψ, z) + lnp(ψ)]
nX p
K X
n X p
K X
X o
∝ exp rik ln(ψk,j,cij ) + ln(Dir(ψk,j |ηj , . . . , ηj )) (41)
i=1 k=1 j=1 k=1 j=1

so for each q(ψk,j ) we have:

n
nX o
q(ψk,j ) ∝ exp rik ln(ψk,j,cij ) + ln(Dir(ψk,j |ηj , . . . , ηj ))
i=1
n
nX dj o
X
q(ψk,j ) ∝ exp rik ln(ψk,j,cij ) + (ηj − 1)ln(ψk,j,g )
i=1 g=1
dj
n X
nX dj o
X
q(ψk,j ) ∝ exp I{cij =g} rik ln(ψk,j,g ) + (ηj − 1)ln(ψk,j,g )
i=1 g=1 g=1
dj
nX n o
X
q(ψk,j ) ∝ exp (ηj − 1 + I{cij =g} rik )ln(ψk,j,g ) (42)
g=1 i=1

where I{cij =g} = 1 if cij = g, I{cij =g}P= 0 otherwise.

n Pn
By observation, q(ψk,j ) ∼ Dir(ηj + i=1 I{cij =1} rik , . . . , ηj + i=1 I{cij =dj } rik )

Now that we have the update equations for q(µ, Λ, ψ, π), we can derive the update for q(z). This
will require the forms of the variational posterior q(µ, Λ, ψ, π) we just derived above in order to evaluate
expectations with respect to q(µ, Λ, ψ, π), in order to obtain an expression for rik :

n o
q(z) ∝ exp E−z lnp(x, c, z, µ, Λ, ψ, π)
n o
∝ exp Eµ,Λ,ψ,π [ln(p(x|c, z, µ, Λ)p(µ|Λ)p(Λ)p(c|ψ, z)p(ψ)p(z|π)p(π))]
n o
∝ exp Eµ,Λ [ln(p(x|c, z, µ, Λ)p(µ|Λ)p(Λ))] + Eψ [ln(p(c|ψ, z)p(ψ))] + Eπ [ln(p(z|π)p(π))]
n o
∝ exp Eµ,Λ [ln(p(x|c, z, µ, Λ))] + Eψ [ln(p(c|ψ, z))] + Eπ [ln(p(z|π))]
nX K
n X n X
X p
K X X K
n X o
∝ exp zik Eµ,Λ [ln(N (xi |µk , Λk−1 ))] + zik Eψ [ln(ψk,j,cij )] + zik Eπ [ln(πk )]
i=1 k=1 i=1 k=1 j=1 i=1 k=1
(43)

where in the penultimate expression we dropped terms which do not depend on z. At this point the
expression of q(z) essentially boils down to evaluating the three expectations Eµ,Λ [ln(N (xi |µk , Λ−1
k ))],
Eψ [ln(ψk,j,cij )] and Eπ [ln(πk )].

n 1 o
Eµ,Λ [ln(N (xi |µk , Λ−1
k ))] = Eµ,Λ [ln (2π)
−q/2
det(Λk )1/2 exp(− (xi − µk )T Λk (xi − µk )) ]
2
q 1 1 T
= − ln(2π) + EΛk [ln(det(Λk ))] − Eµ,Λ [ (xi − µk ) Λk (xi − µk )] (44)
2 2 2
Pq
It is well known that if Λ ∼ W(ν, Φ−1 ), then EΛ [ln(det(Λ))] = qln(2)+ln(det(Φ−1 ))+ i=1 𭟋( ν+1−i 2 ),
where 𭟋(x) is the digamma function. Since Λk ∼ W(ν̂k , Φ̂−1 k ), we have
Pq
EΛk [ln(det(Λk ))] = qln(2) + ln(det(Φ̂−1 k )) + i=1 𭟋( ν̂k +1−i
2 )
As for Eµ,Λ [ 12 (xi − µk )T Λk (xi − µk )], we have:

23
 Bayesian Variational Inference for Mixed Data Mixture Models

1
Eµ,Λ [ (xi − µk )T Λk (xi − µk )]
Z2Z
1
= (xi − µk )T Λk (xi − µk )q(µk , Λk )dµk dΛk
2
1
ZZ
= (xi − m̂k + m̂k − µk )T Λk (xi − m̂k
2
+ m̂k − µk )q(µk , Λk )dµk dΛk
1
ZZ
= [(xi − m̂k )T Λk (xi − m̂k ) + 2(xi − m̂k )T Λk (m̂k − µk )+ (45)
2
(m̂k − µk )T Λk (m̂k − µk )]q(µk |Λk )q(Λk )dµk dΛk
1
ZZ
= [(xi − m̂k )T Λk (xi − m̂k )
2
+ T r((µk − m̂k )(µk − m̂k )T Λk )]q(µk |Λk )q(Λk )dµk dΛk
1
ZZ
= [(xi − m̂k )T Λk (xi − m̂k ) + T r(β̂k−1 Λ−1
k Λk )]q(Λk )dΛk
2
1
= [(xi − m̂k )T ν̂k Φ̂−1 −1
k (xi − m̂k ) + q β̂k ] (46)
2
Next we evaluate Eπ [ln(πk )] and Eψ [ln(ψk,j,cij )]. The expectation of the natural logarithm of a
PK
Dirichlet is well known, specifically if X ∼ Dir(α1 , . . . , αK ), then E(ln(Xi )) = 𭟋(αi ) − 𭟋( j=1 αj ).
With this, we have:

n
X K X
X n
Eπ [ln(πk )] = 𭟋(α + rik ) − 𭟋(Kα + rij ) (47)
i=1 j=1 i=1

n dj n
X X X
Eψ [ln(ψk,j,cij )] = 𭟋(ηj + I{csj =cij } rsk ) − 𭟋(dj ηj + I{csj =gj } rsk ) (48)
s=1 gj =1 s=1

Pn Pn
Recalling that q(π) ∼ Dir(α
Pn + i=1 ri1 , . . . , α + Pi=1 riK )
n
and q(ψk,j ) ∼ Dir(ηj + i=1 I{cij =1} rik , . . . , ηj + i=1 I{cij =dj } rik )
We return to the expression for q(z) now that we have all the required expectations:

n X
nX K n X
X p
K X n X
X K o
q(z) ∝ exp zik Eµ,Λ [ln(N (xi |µk , Λk−1 ))] + zik Eψ [ln(ψk,j,cij )] + zik Eπ [ln(πk )]
i=1 k=1 i=1 k=1 j=1 i=1 k=1
n Y
Y K n p
X o
q(z) ∝ exp zik Eµ,Λ [ln(N (xi |µk , Λk−1 ))] + zik Eψ [ln(ψk,j,cij )] + zik Eπ [ln(πk )]
i=1 k=1 j=1
n Y
Y K n Y
Y K
q(z) ∝ exp(zik ln(ρik )) = ρzikik (49)
i=1 k=1 i=1 k=1

where
n p
X o
ρik = exp Eµ,Λ [ln(N (xi |µk , Λ−1
k ))] + Eψ [ln(ψ k,j,c ij
)] + Eπ [ln(π k )] (50)
j=1

ρ
By observation, we have rik = PK ik
ρij
. With this, we can alternate between updating q(µ, Λ, ψ, π),
j=1
which depend on rik and the rik which depend on q(µ, Λ, ψ, π).

24
 Bayesian Variational Inference for Mixed Data Mixture Models

B Evaluating the ELBO

Though the CAVI algorithm update equations do not require the ELBO, the termination condition of
CAVI is based on checking whether the ELBO has converged (in practice, whether the difference in
ELBO between subsequent iterations has reduced beneath a pre-specified threshold), which require us to
evaluate it. This involves another long calculation though some of the terms have been calculated from
deriving the CAVI update equations. Also, unlike the CAVI updates, one needs to be careful not to drop
any ‘constant’ terms here as we are no longer deriving probability distributions up to a constant.

ELBO(q) = Eq [ln(p(x, c, z, µ, Λ, ψ, π)) − ln(q(µ, Λ, ψ, π, z))] (51)

= Eq [ln(p(x|c, z, µ, Λ)p(µ|Λ)p(Λ)p(c|ψ, z)p(ψ)p(z|π)p(π)) − ln(q(µ, Λ)q(ψ)q(π)q(z))]
= Eµ,Λ,z [ln(p(x|c, z, µ, Λ))] + Eµ,Λ [ln(p(µ|Λ))] + EΛ [ln(p(Λ))]
+ Eψ,z [ln(p(c|ψ, z))] + Eψ [ln(p(ψ))] + Eπ,z [ln(p(z|π))] + Eπ [ln(p(π))]
− Eµ,Λ [ln(q(µ, Λ))] − Eψ [ln(q(ψ))] − Eπ [ln(q(π))] − Ez [ln(q(z))] (52)
There are 11 expectations here, which we evaluate in order:

n X
X K
Eµ,Λ,z [ln(p(x|c, z, µ, Λ))] = rik Eµ,Λ [ln(N (xi |µk , Λ−1
k ))]
i=1 k=1
n X
K q
X n q 1 −1
X ν̂k + 1 − i
= rik − ln(2π) + [qln(2) + ln(det(Φ̂k )) + 𭟋( )]
i=1 k=1
2 2 i=1
2
1 o
− [(xi − m̂k )T ν̂k Φ̂−1
k (x i − m̂ k ) + q β̂ −1
k ] (53)
2

K
Y K
X
Eµ,Λ [ln(p(µ|Λ))] = Eµ,Λ [ln( N (µk |mk , β −1 Λ−1
k ))] = Eµ,Λ [ln(N (µk |mk , β −1 Λ−1
k ))]
k=1 k=1
K
X n 1 o
= Eµ,Λ [ln (2π)−q/2 det(βΛk )1/2 exp(− (µk − mk )T βΛk (µk − mk )) ]
2
k=1
K n
X q 1 1 o
= − ln(2π) + EΛ [ln(det(βΛk ))] + Eµ,Λ [− (µk − mk )T βΛk (µk − mk )]
2 2 2
k=1
K n
X q q 1
= − ln(2π) + ln(β) + EΛ [ln(det(Λk ))]
2 2 2
k=1
1 o
+ Eµ,Λ [− (µk − m̂k + m̂k − mk )T βΛk (µk − m̂k + m̂k − mk )]
2
K n
X q q 1 1
= − ln(2π) + ln(β) + EΛ [ln(det(Λk ))] + Eµ,Λ [− (µk − m̂k )T βΛk (µk − m̂k )]
2 2 2 2
k=1
1 o
− Eµ,Λ [(µk − m̂k )T βΛk (m̂k − mk )] − EΛ [ (m̂k − mk )T βΛk (m̂k − mk )]
2
K
X n q q 1 1
= − ln(2π) + ln(β) + EΛ [ln(det(Λk ))] − (T r(Eµ,Λ [(µk − m̂k )(µk − m̂k )T βΛk ]))
2 2 2 2
k=1
1 o
− EΛ [ (m̂k − mk )T βΛk (m̂k − mk )]
2
K n q
X q q 1 X ν̂k + 1 − i q
= − ln(2π) + ln(β) + [qln(2) + ln(det(Φ̂−1 k )) + 𭟋( )] − β̂k−1 β
2 2 2 i=1
2 2
k=1
1 o
− (m̂k − mk )T β ν̂k Φ̂−1
k (m̂k − mk ) (54)
2

25
 Bayesian Variational Inference for Mixed Data Mixture Models

K
Y K
X
EΛ [ln(p(Λ))] = EΛ [ln( W(Λk |ν, Φ−1 ))] = EΛk [ln(W(Λk |ν, Φ−1 ))]
i=1 i=1
K (ν−q−1)/2
X det(Λk ) exp(−T r(ΦΛk )/2)
= EΛk [ln( )]
i=1
2νq/2 det(Φ−1 )ν/2 Γq (ν/2)
K
X ν−q−1 1 νq ν
= EΛk [ ln(det(Λk )) − T r(ΦΛk ) − ln(2) − ln(det(Φ−1 )) − ln(Γq (ν/2))]
i=1
2 2 2 2
K n q
X ν−q−1 X ν̂k + 1 − i 1
= [qln(2) + ln(det(Φ̂−1
k )) + 𭟋( )] − T r(Φν̂k Φ̂−1
k )
i=1
2 i=1
2 2
νq ν −1
o
− ln(2) − ln(det(Φ )) − ln(Γq (ν/2)) (55)
2 2
q(q−1)/4
Qq
where Γq (a) = π j=1 Γ(a + (1 − j)/2)

n Y
Y p
K Y n X
X p
K X
zik
Eψ,z [ln(p(c|ψ, z))] = Eψ,z [ln( ψk,j,c ij
)] = rik Eψ [ln(ψk,j,cij )]
i=1 k=1 j=1 i=1 k=1 j=1
p
K X
n X n n dj n o
X X X X
= rik 𭟋(ηj + I{csj =cij } rsk ) − 𭟋(dj ηj + I{csj =gj } rsk ) (56)
i=1 k=1 j=1 s=1 gj =1 s=1

Y p
K Y p
K X
X
Eψ [ln(p(ψ))] = Eψ [ln( Dir(ψk,j |ηj , . . . , ηj ))] = Eψ [ln(Dir(ψk,j |ηj , . . . , ηj ))]
k=1 j=1 k=1 j=1
p
K X dj K Xp dj
X 1 Y ηj −1 X 1 X
= Eψ [ln( ψk,j,g )] = Eψ [ln( ) + (ηj − 1) ln(ψk,j,g )]
B(ηj ) g=1 B(ηj ) g=1
k=1 j=1 j=1 k=1
p
K X dj n dj n
X 1 Xn X X X o
= [ln( ) + (ηj − 1) 𭟋(ηj + I{csj =g} rsk ) − 𭟋(dj ηj + I{csj =gj } rsk ) ] (57)
j=1
B(ηj ) g=1 s=1 g =1 s=1
k=1 j

QK Qdj
Γ(αk ) g=1 Γ(ηj )
where B(α) = k=1
Γ( K
P , where with a slight abuse of notation we define B(ηj ) = Pdj and
k=1 αk ) Γ( g=1 ηj )
QK
Γ(α)
B(α) = k=1
Γ( K
P .
k=1 α)

n Y
Y K n X
X K
Eπ,z [ln(p(z|π))] = Eπ,z [ln( πkzik )] = rik Eπ [ln(πk )]
i=1 k=1 i=1 k=1
n X
X K n
X K X
X n
= rik [𭟋(α + rik ) − 𭟋(Kα + rij )] (58)
i=1 k=1 i=1 j=1 i=1

K
1 Y α−1
Eπ [ln(p(π))] = Eπ [ln(Dir(π|α, . . . , α))] = Eπ [ln( πk ))]
B(α)
k=1
K
1 X
= ln( )+ (α − 1)Eπ [ln(πk )]
B(α)
k=1
K n K X n
1 X X X
= ln( )+ (α − 1)[𭟋(α + rik ) − 𭟋(Kα + rij )] (59)
B(α) i=1 j=1 i=1
k=1

26
 Bayesian Variational Inference for Mixed Data Mixture Models

Eµ,Λ [ln(q(µ, Λ))] = Eµ,Λ [ln(q(µ|Λ)q(Λ))] = Eµ,Λ [ln(q(µ|Λ))] + EΛ [ln(q(Λ))] (60)

We derive the two terms here separately for easier presentation.

K
Y K
X
Eµ,Λ [ln(q(µ|Λ))] = Eµ,Λ [ln( N (µk |m̂k , β̂k−1 Λ−1
k )] = Eµ,Λ [ln(N (µk |m̂k , β̂k−1 Λ−1
k )]
k=1 k=1
K
X n 1 o
= Eµ,Λ [ln (2π)−q/2 det(β̂k Λk )1/2 exp(− (µk − m̂k )T β̂k Λk (µk − m̂k )) ]
2
k=1
K n
X q 1 1 o
= − ln(2π) + EΛ [ln(det(β̂k Λk ))] + Eµ,Λ [− (µk − m̂k )T β̂k Λk (µk − m̂k )]
2 2 2
k=1
K n
X q q 1 1 o
= − ln(2π) + ln(β̂k ) + EΛ [ln(det(Λk ))] + Eµ,Λ [− T r((µk − m̂k )(µk − m̂k )T β̂k Λk )]
2 2 2 2
k=1
K n
X q q 1 qo
= − ln(2π) + ln(β̂k ) + EΛ [ln(det(Λk ))] −
2 2 2 2
k=1
K n q
X q q 1 X ν̂k + 1 − i qo
= − ln(2π) + ln(β̂k ) + [qln(2) + ln(det(Φ̂−1
k )) + 𭟋( )] − (61)
2 2 2 i=1
2 2
k=1

K
Y K
X
EΛ [ln(q(Λ))] = EΛ [ ln(W(Λk |ν̂k , Φ̂−1
k ))] = EΛk [ln(W(Λk |ν̂k , Φ̂−1
k ))]
k=1 k=1
K
X det(Λk )(ν̂k −q−1)/2 exp(−T r(Φ̂k Λk )/2)
= EΛk [ln( )]
i=1 2ν̂k q/2 det(Φ̂−1
k )
ν̂k /2 Γ (ν̂ /2)
q k
K
X ν̂k − q − 1 1 ν̂k q ν̂k
= EΛk [ ln(det(Λk )) − T r(Φ̂k Λk ) − ln(2) − ln(det(Φ̂−1
k )) − ln(Γq (ν̂k /2))]
i=1
2 2 2 2
K n q
X ν̂k − q − 1 X ν̂k + 1 − i ν̂k q
= [qln(2) + ln(det(Φ̂−1
k )) + 𭟋( )] −
i=1
2 i=1
2 2
ν̂k q ν̂k o
− ln(2) − ln(det(Φ̂−1
k )) − ln(Γq (ν̂k /2)) (62)
2 2

Y p
K Y p
K X
X
Eψ [ln(q(ψ))] = Eψ [ln( Dir(ψk,j |η̂k,j,1 , . . . , η̂k,j,dj ))] = Eψ [ln(Dir(ψk,j |η̂k,j,1 , . . . , η̂k,j,dj ))]
k=1 j=1 k=1 j=1
p
K X dj K Xp dj
X 1 Y η̂k,j,g −1
X 1 X
= Eψ [ln( ψk,j,g )] = Eψ [ln( ) + (η̂k,j,g − 1) ln(ψk,j,g )]
B(ηk,j ) g=1 B(ηk,j ) g=1
k=1 j=1 j=1
k=1
p
K X dj n dj n
X 1 Xn X X X o
= [ln( ) + (η̂k,j,g − 1) 𭟋(ηj + I{csj =g} rsk ) − 𭟋(dj ηj + I{csj =gj } rsk ) ]
B(ηk,j ) g=1 s=1 g =1 s=1
k=1 j=1 j

(63)
Pn
where η̂k,j,g = ηj + i=1 I{cij =g} rik

27
 Bayesian Variational Inference for Mixed Data Mixture Models

n n K
X X 1 Y αk −1
Eπ [ln(q(π))] = Eπ [ln(Dir(π|α + ri1 , . . . , α + riK ))] = Eπ [ln( πk ))]
i=1 i=1
B(α)
k=1
K
1 X
= ln( )+ (αk − 1)Eπ [ln(πk )]
B(α)
k=1
K n K X n
1 X X X
= ln( )+ (αk − 1)[𭟋(α + rik ) − 𭟋(Kα + rij )] (64)
B(α) i=1 j=1 i=1
k=1
Pn
where α̂k = α + i=1 rik

n Y
Y K n X
X K n X
X K
zik
Ez [ln(q(z))] = Ez [ln( rik )] = Ez [zik ln(rik )] = rik ln(rik ) (65)
i=1 k=1 i=1 k=1 i=1 k=1

One can alternatively calculate some of the ELBO terms in a somewhat simplified way by pairing
them up and computing the KL divergence instead. Specifically:

ELBO(q) = Eq [ln(p(x, c, z, µ, Λ, ψ, π)) − ln(q(µ, Λ, ψ, π, z))] (66)

= Eq [ln(p(x|c, z, µ, Λ)p(µ|Λ)p(Λ)p(c|ψ, z)p(ψ)p(z|π)p(π)) − ln(q(µ, Λ)q(ψ)q(π)q(z))]
= Eµ,Λ,z [ln(p(x|c, z, µ, Λ))] + Eµ,Λ [ln(p(µ|Λ))] + EΛ [ln(p(Λ))]
+ Eψ,z [ln(p(c|ψ, z))] + Eψ [ln(p(ψ))] + Eπ,z [ln(p(z|π))] + Eπ [ln(p(π))]
− Eµ,Λ [ln(q(µ, Λ))] − Eψ [ln(q(ψ))] − Eπ [ln(q(π))] − Ez [ln(q(z))]
= Eµ,Λ,z [ln(p(x|c, z, µ, Λ))] + Eψ,z [ln(p(c|ψ, z))] + Eπ,z [ln(p(z|π))] − Ez [ln(q(z))]
+ Eµ,Λ [ln(p(µ|Λ))] − Eµ,Λ [ln(q(µ|Λ))] + EΛ [ln(p(Λ))] − EΛ [ln(q(Λ))]
+ Eπ [ln(p(π))] − Eπ [ln(q(π))] + Eψ [ln(p(ψ))] − Eψ [ln(q(ψ))]
= Eµ,Λ,z [ln(p(x|c, z, µ, Λ))] + Eψ,z [ln(p(c|ψ, z))] + Eπ,z [ln(p(z|π))] − Ez [ln(q(z))]
− EΛ [KL(q(µ|Λ) ∥ p(µ|Λ))] − KL(q(Λ) ∥ p(Λ)) − KL(q(π) ∥ p(π)) − KL(q(ψ) ∥ p(ψ)) (67)

The KL terms are all KL divergences of well known distributions with analytical expressions, which
can be evaluated as follows:

K
Y K
Y
EΛ [KL(q(µ|Λ) ∥ p(µ|Λ))] = EΛ [KL( N (µk |m̂k , β̂k−1 Λ−1
k )∥ N (µk |mk , β −1 Λ−1
k ))]
k=1 k=1
K
X
= EΛ [KL(N (µk |m̂k , β̂k−1 Λ−1
k ) ∥ N (µk |mk , β
−1 −1
Λk ))]
k=1
K
X 1n det(β −1 Λ−1
k ) T −1 −1
o
= EΛ [ln( ) − q + ( m̂k − m k ) βΛk (m̂ k − mk ) + T r(βΛ k β̂ k Λk )]
2
k=1
det(β̂k−1 Λ−1
k )
K
X 1n β −q det(Λ−1
k )
o
= EΛ [ln( −q −1
) − q + (m̂k − mk )T βΛk (m̂k − mk ) + qβ β̂k−1 ]
2 β̂k det(Λk )
k=1
K
X 1n o
= −qln(β) + qln(β̂k ) − q + (m̂k − mk )T β ν̂k Φ̂−1
k (m̂k − mk ) + qβ β̂ −1
k (68)
2
k=1

28
 Bayesian Variational Inference for Mixed Data Mixture Models

K
Y K
Y
KL(q(Λ) ∥ p(Λ)) = KL( W(Λk |ν̂k , Φ̂−1
k )∥ W(Λk |ν, Φ−1 ))
k=1 i=1
K
X
= KL(W(Λk |ν̂k , Φ̂−1
k ) ∥ W(Λk |ν, Φ
−1
))
k=1
K
X 1n
= νln(det(Φ−1 )) − νln(det(Φ̂−1 −1
k )) + ν̂k T r(ΦΦ̂k ) + 2ln(Γq (ν/2)) − 2ln(Γq (ν̂k /2))
2
k=1
q
X ν̂k + 1 − i o
+ (ν̂k − ν) 𭟋( ) − ν̂k q (69)
i=1
2

n
X n
X
KL(q(π) ∥ p(π)) = KL(Dir(π|α + ri1 , . . . , α + riK ) ∥ Dir(π|α, . . . , α))
i=1 i=1
K
X n
X K
X K
X K
X n
X
= ln(Γ( (α + rik ))) − ln(Γ( α)) + ln(Γ(α)) − ln(Γ(α + rik ))
k=1 i=1 k=1 k=1 k=1 i=1
K
X n
X n n
X K
X n
X o
+ (α + rik − α) 𭟋(α + rik )) − 𭟋( (α + rik ))
k=1 i=1 i=1 k=1 i=1
K X
X n K
X K
X n
X
= ln(Γ(Kα + rik )) − ln(Γ(Kα)) + ln(Γ(α)) − ln(Γ(α + rik ))
k=1 i=1 k=1 k=1 i=1
K X
X n n n
X K X
X n o
+ rik 𭟋(α + rik ) − 𭟋(Kα + rik ) (70)
k=1 i=1 i=1 k=1 i=1

Y p
K Y p
K Y
Y
KL(q(ψ) ∥ p(ψ)) = KL( Dir(ψk,j |η̂k,j,1 , . . . , η̂k,j,dj ) ∥ Dir(ψk,j |ηj , . . . , ηj ))
k=1 j=1 k=1 j=1
p
K X
X
= KL(Dir(ψk,j |η̂k,j,1 , . . . , η̂k,j,dj ) ∥ Dir(ψk,j |ηj , . . . , ηj ))
k=1 j=1
p
K X dj dj dj dj
X X X X X
= [ln(Γ( η̂k,j,gj )) − ln(Γ( ηj )) + ln(Γ(ηj )) − ln(Γ(η̂k,j,gj ))
k=1 j=1 gj =1 gj =1 gj =1 gj =1
dj n dj o
X X
+ (η̂k,j,gj − ηj ) 𭟋(η̂k,j,gj ) − 𭟋( η̂k,j,gj ) ]
gj =1 gj =1
p
K X dj n dj dj n
X X X X X X
= [ln(Γ( (ηj + I{csj =gj } rsk ))) − ln(Γ(dj ηj )) + ln(Γ(ηj )) − ln(Γ(ηj + I{csj =gj } rsk ))
k=1 j=1 gj =1 s=1 gj =1 gj =1 s=1
dj n n n dj n o
X X X X X
+ (ηj + I{csj =gj } rsk − ηj ) 𭟋(ηj + I{csj =gj } rsk ) − 𭟋( (ηj + I{csj =gj } rsk )) ]
gj =1 s=1 s=1 gj =1 s=1
p
K X dj n dj dj n
X X X X X X
= [ln(Γ(dj ηj + I{csj =gj } rsk )) − ln(Γ(dj ηj )) + ln(Γ(ηj )) − ln(Γ(ηj + I{csj =gj } rsk ))
k=1 j=1 gj =1 s=1 gj =1 gj =1 s=1
dj n n n dj n o
X X X X X
+ I{csj =gj } rsk 𭟋(ηj + I{csj =gj } rsk ) − 𭟋(dj ηj + I{csj =gj } rsk ) ] (71)
gj =1 s=1 s=1 gj =1 s=1
Pn
where η̂k,j,gj = ηj + s=1 I{csj =gj } rsk

29
 Bayesian Variational Inference for Mixed Data Mixture Models

C Gibbs sampler for the model

With a Gibbs sampler we want to calculate the conditional distribution of every parameter conditioned
on every other parameter. Recall the parameters are µk , Σk , πk , ψk,j , zi .
First we update zi |x, c, µ, Σ, ψ, π, for each 1 ≤ i ≤ n, which has conditional distribution:

p(zi = k|x, c, µ, Σ, ψ, π) ∝ p(zi = k|π)p(xi |ci , zi , µ, Σ)p(ci |ψ, zi = k)

p
Y
∝ πk N (xi |µk , Σk ) ψk,j,cij (72)
j=1

This is just a discrete distribution so you normalise this probability by dividing by the sum of the
probabilities over k.
Next we update ψk,j for each 1 ≤ k ≤ K, 1 ≤ j ≤ p

Y
p(ψk,j |x, c, z, µ, Σ, π) ∝ p(ψk,j ) p(cij |zi = k, ψ)
i:zi =k
Y
∝ Dir(ψk,j |ηj , . . . , ηj ) ψk,j,cij (73)
i:zi =k

This is a Dirichlet prior with a categorical data likelihood, which gives a conjugate Dirichlet posterior.
In particular, the posterior is of the form:

n
X n
X
p(ψk,j |x, c, z, µ, Σ, π) ∼ Dir(ηj + I(cij = 1, zi = k), . . . , ηj + I(cij = ηj , zi = k)) (74)
i=1 i=1

Update µk , Σk in a block:

Y
p(µk , Σk |x, c, z, π, ψ) ∝ p(µk |Σk )p(Σk ) p(xi |zi = k, ci , µ, Σ)
i:zi =k
Y
∝ N (µk |mk , β0−1 Σk )IW(Σk |ν, Φ) N (xi |µk , Σk ) (75)
i:zi =k

This is a Normal-inverse Wishart prior with a Gaussian likelihood, which is conjugate and gives
another Normal-inverse Wishart posterior of the form:

P
β0 mk + i:zi =k xi Σk
p(µk |Σk , x, c, z, π, ψ) ∼ N (µk | , )
β0 + n k β0 + n k
X β0 n k
p(Σk |x, c, z, π, ψ) ∼ IW(Σk |ν + nk , Φ + (xi − x̄k )(xi − x̄k )T + (x̄k − mk )(x̄k − mk )T )
β0 + n k
i:zi =k
(76)
1
P
where x̄k = nk i:zi =k xi
Lastly, we update π:

n
Y
p(π|x, c, z, µ, γ, Σ, ψ, s) ∝ p(π) p(zi |π)
i=1
n
Y
∝ Dir(α, . . . , α) πzi (77)
i=1

30
 Bayesian Variational Inference for Mixed Data Mixture Models

Again like the ψ update this is a Dirichlet prior with categorical likelihood so we have conjugate
posterior:

n
X n
X
p(π|x, c, z, µ, γ, Σ, ψ, s) ∼ Dir(α + I(zi = 1), . . . , α + I(zi = K)) (78)
i=1 i=1

D Proof of Theorem 1
We have
∇π̂ Π ∇µ̂ Π ∇Λ̂ Π ∇ψ̂ Π
 
∇ π̂ M ∇µ̂ M ∇Λ̂ M ∇ψ̂ M 
∇T ϵ (Θ) = (1 − ϵ)I + ϵ

 ∇π̂ S
 (79)
∇µ̂ S ∇Λ̂ S ∇ψ̂ S 
∇π̂ Ψ ∇µ̂ Ψ ∇Λ̂ Ψ ∇ψ̂ Ψ
(t−1) (t−1)
Because the gradient terms in ∇T ϵ (Θ) involves rik or the derivatives of rik , we need to calculate
(t−1) (t−1) (t−1)
the limit of rik and the derivatives of rik as n → ∞. We start with the limit of rik , recall that
(t−1) ρ (Θ(t−1) )
rik = rik (Θ(t−1) ) = PK ik (t−1) ) , where
j=1 ρij (Θ

n o
ρik (Θ) = exp Eµ,Λ [ln(N (xi |µk , Λ−1
k ))] + Eψ [ln(ψk,ci )] + Eπ [ln(πk )] (80)

Where the expectations Eµ,Λ , Eψ , Eπ are with respect to the variational posterior with hyperparameter
values computed at Θ, which recall is related to Θ through the following relationships:

α̂k = nπ̂k + α
m̂k = (nµ̂k π̂k + βmk )/(nπ̂k + β)
β̂k = nπ̂k + β
ν̂k = nπ̂k + ν
nπ̂k
Φ̂k = nπ̂k (Λ̂k )−1 + β(µ̂k − mk )(µ̂k − mk )T + Φ
nπ̂k + β
η̂k,g = nπ̂k ψ̂k,g + η (81)
We will make use of the following lemma:
Lemma 1 (Lemma 1 of Titterington and Wang [2006]). Suppose that pn (x) is the probability density
function of the Rm valued random vector Xn = (x1n , . . . , xm T
n ) , that E(Xn ) = µn → µ and Covij (Xn ) =
1
O( n ) as n → ∞. Then, for any function f (.) with continuous second order derivative near µ, it holds
that:
1
E(f (Xn )) = f (µn ) + O( )
n
Proof. See Titterington and Wang [2006] ■
∗
πk p(xi ,ci |zi =k,Θ∗ )
Lemma 2. As n → ∞, ρik (Θ∗ ) = πk∗ p(xi , ci |zi = k, Θ∗ )+O( n1 ). Consequently, rik (Θ∗ ) → PK
π̂ ∗ p(xi ,ci |zi =k′ ,Θ∗ )
k′ =1 k′
   ∗
π̂ π
 µ̂   µ∗ 
Taking Θ =   =  ∗ 
   = Θ∗
Λ̂ Λ 
ψ̂ ψ∗
Proof. Recall that
q 1 1
Eµ,Λ [ln(N (xi |µk , Λ−1 T −1 −1
k ))] = − ln(2π) + EΛk [ln(det(Λk ))] − [(xi − m̂k ) ν̂k Φ̂k (xi − m̂k ) + qβk ]
2 2 2
(82)

31
 Bayesian Variational Inference for Mixed Data Mixture Models

1
− (xi − m̂k )T ν̂k Φ̂−1
k (xi − m̂k )
2
1 1 1 1
= − (xi − µ∗k − O( ))T (Λ∗k + O( ))(xi − µ∗k − O( ))
2 n n n
1 ∗ T ∗ ∗ 1
= − (xi − µk ) Λk (xi − µk ) + O( ) (83)
2 n
where we used m̂k = nµ∗k πk∗ /(nπk∗ + β) + βmk /(nπk∗ + β) → µ∗k + O( n1 ) and

nπk∗
ν̂k Φ̂−1 ∗ ∗ ∗ −1
k = (nπk + ν)(nπk (Λk ) + β(µ∗k − mk )(µ∗k − mk )T + Φ)−1
nπk∗ + β
Λ∗k Λ∗
= (nπk∗ + ν)(nπk∗ (Λ∗k )−1 (I + ∗ β(µ∗k − mk )(µ∗k − mk )T + k∗ Φ))−1
nπk + β nπk
∗ ∗ ∗
Λk Λ nπ + ν
= (I + β(µ∗k − mk )(µ∗k − mk )T + k∗ Φ)−1 k ∗ Λ∗k
nπk∗ + β nπk nπk
1
= Λ∗k + O( ) (84)
n

1
q β̂k−1 = O( )
n
 
For 12 EΛk [ln(det(Λk ))], EΛk (Λk ) = ν̂k Φ̂−1 ∗ −1 −1 −1 −1
k → Λk , Cov([Λk ]ij , [Λk ]ab ) = ν̂k [Φ̂k ]ia [Φ̂k ]jb + [Φ̂k ]ib [Φ̂k ]ja =
O n1 . So using Lemma 1 on q(Λk ), we have:

1 1 1
EΛ [ln(det(Λk ))] = ln(det(Λ∗k )) + O( ) (85)
2 k 2 n
Therefore
q 1 1
Eµ,Λ [ln(N (xi |µk , Λ−1
k ))] = − ln(2π) + EΛ [ln(det(Λk ))] − [(xi − m̂k )T ν̂k Φ̂−1 −1
k (xi − m̂k ) + qβk ]
2 2 k 2
q 1 1 1
= − ln(2π) + ln(det(Λ∗k )) − (xi − µ∗k )T Λ∗k (xi − µ∗k ) + O( ) (86)
2 2 2 n
∗
nπk +α
For q(π), we note that Eπ (πk ) = PKα̂k
α̂k′
= n+Kα → πk∗ and
k′ =1

δkl α̂k α̂k α̂l 1

Cov(πk , πl ) = PK PK − PK PK PK = O( )
( k′ =1 α̂k′ )( k′ =1 α̂k′ + 1) ( k′ =1 α̂k′ )( k′ =1 α̂k′ )( k′ =1 α̂k′ + 1) n

So by Lemma 1, we have
1
Eπ [ln(πk )] = ln(πk∗ ) + O( ) (87)
n
∗ ∗
η̂k,ci nπk ψk,c +η ∗
Similarly, for q(ψk ), we have Eψ [ψk,ci ] = Pd = nπk
i
∗ +dη → ψk,c .
g=1 η̂k,g
i

δgh η̂k,g η̂k,g η̂k,h 1

Cov(η̂g,h ) = Pd Pd − Pd Pd Pd = O( )
( g′ =1 η̂k,g′ )( g′ =1 η̂k,g′ + 1) ( g′ =1 η̂k,g′ )( g′ =1 η̂k,g′ )( g′ =1 η̂k,g′ + 1) n

So by Lemma 1, we have

∗ 1
Eψ [ln(ψk,ci )] = ln(ψk,c ) + O( ) (88)
i
n
Therefore:

32
 Bayesian Variational Inference for Mixed Data Mixture Models

−q 1 1 1
ρik (Θ∗ ) = πk∗ ψk,c
∗
(2π) det(Λ∗k ) 2 exp(− (xi − µ∗k )T Λ∗k (xi − µ∗k ) + O( ))
2
i
2 n
∗ ∗ ∗ −q
∗ 12 1 ∗ T ∗ ∗ 1
ρik (Θ ) = πk ψk,ci (2π) det(Λk ) exp(− (xi − µk ) Λk (xi − µk )) + O( )
2
2 n
∗ ∗ 1
= πk p(xi , ci |zi = k, Θ ) + O( ) (89)
n
∗
πk p(xi ,ci |zi =k,Θ∗ )
So as n → ∞, ρik (Θ∗ ) → πk∗ p(xi , ci |zi = k, Θ∗ ). Also, rik (Θ∗ ) → PK
πk∗ p(x ,c |z =k ′ ,Θ∗ ) . ■
k′ =1 ′ i i i

Next we need to calculate the derivatives of rik (Θ) with respect to Θ as n → ∞, at Θ = Θ∗ . To do

this, we first examine how the derivative of rik looks like in general. Writing rik = PKexp(ln(ρik ))
, and
i=1 exp(ln(ρik ))
applying the quotient rule, we have:
PK PK ∂ ln(ρ )
∂rik ρik ∂ ln(ρik )
∂θk′ ( j=1ρij ) − j=1 ρij ∂θ ′ij ρik
= PK
k

∂θk′ ( j=1 ρij )2

K
∂ ln(ρik ) X ∂ ln(ρij )
= rik − rij rik (90)
∂θk′ j=1
∂θk′

Here θk′ is any scalar parameter within Θ belonging to component k ′ , so we could have for example µ̂k′ l ,
the lth coordinate of µ̂k′ . Note if ρik is independent of θk′ for k ′ ̸= k, then

∂rik ∂ ln(ρik′ ) n o
= δkk′ rik − rik rik′ (91)
∂θk′ ∂θk′

and we only need to compute ∂ ln(ρik′ )

∂θk′ . ln(ρik (Θ)) as a function of Θ is the same as the continuous data
case (see Titterington and Wang [2006], above Eq 8) except with the addition of the Eψ [ln(ψk,ci )] term,
and because
n
X d X
X n
Eψ [ln(ψk,ci )] = 𭟋(η + I{cs =ci } rsk ) − 𭟋(dη + I{cs =g} rsk )
s=1 g=1 s=1
d
X
= 𭟋(η + nπ̂k ψ̂k,ci ) − 𭟋(dη + nπ̂k ψ̂k,g ) (92)
g=1

∂ ln(ρik ) ∂ ln(ρik )
is independent of µ̂, Λ̂, we therefore have ∂ µ̂kl , and ∂ Λ̂kjl
are unchanged from the continuous data
case, which as n → ∞ are of the form (see Titterington and Wang [2006]):

q
∂ ln(ρik ) X
→ Λ̂klj (xij − µ̂kj ) (93)
∂ µ̂kl j=1

∂ ln(ρik ) 1 −1
→ [Λ̂ − (xij − µ̂kj )(xil − µ̂kl )] (94)
∂ Λ̂kjl 2 kjl

∂ ln(ρik )
for ∂ ψ̂k,g
, we use the following lemma:

Lemma 3. Let 𭟋 be the digamma function on the positive reals, and b1 > 0, b2 ≥ 0, b3 ≥ 0 be real
constants and n positive integer. Let λ = b1 nx + b2 n + b3 , then:

∂𭟋(λ) ∂λ
∼ (λ − 1)−1 = b1 n(b1 nx + b2 n + b3 − 1)−1 → (x + b2 /b1 )−1
∂x ∂x
as n → ∞.

33
 Bayesian Variational Inference for Mixed Data Mixture Models

Proof. See Titterington and Wang [2006], Appendix D, Pages 643-644, which proves the case for λ =
nx + b3 but the proof is the same except substituting λ instead for λ = b1 nx + b2 n + b3 and changing
the chain rule derivative ∂λ
∂x = b1 n accordingly. ■
∂ ln(ρik )
Using Lemma 3, for ∂ ψ̂k,g
, we have:

Pd
∂ ln(ρik ) ∂𭟋(η + nπ̂k ψ̂k,ci ) ∂𭟋(dη + nπ̂k g=1 ψ̂k,g )
= −
∂ ψ̂k,g ∂ ψ̂k,g ∂ ψ̂k,g
Pd d
∂(η + nπ̂k ψ̂k,ci ) −1
∂(dη + nπ̂k g=1 ψ̂k,g ) X
∼ (η + nπ̂k ψ̂k,ci − 1) − (dη + nπ̂k ψ̂k,g − 1)−1
∂ ψ̂k,g ∂ ψ̂k,g g=1
d
δgci X δgci 1
→ − (nπ̂k )(nπ̂k ψ̂k,g )−1 = − Pd (95)
ψ̂k,ci g=1 ψ̂k,g g=1 ψ̂k,g

where δgci = I{ci =g} . Lastly, ∂ ln(ρ

∂ π̂k
ik )
is a little more involved, as the additional Eψ [ln(ψk,ci )] term
Equation (92) does depend on π̂k . Again using Lemma 3, we have:
Pd
∂Eψ [ln(ψk,ci )] ∂𭟋(η + nπ̂k ψ̂k,ci ) ∂𭟋(dη + nπ̂k g=1 ψ̂k,g )
= −
∂ π̂k ∂ π̂k ∂ π̂k
Pd d
∂(η + nπ̂k ψ̂k,ci ) −1
∂(dη + nπ̂k g=1 ψ̂k,g ) X
∼ (η + nπ̂k ψ̂k,ci − 1) − (dη + nπ̂k ψ̂k,g − 1)−1
∂ π̂k ∂ π̂k g=1
d
X d
X
→(nψ̂k,ci )(nπ̂k ψ̂k,ci )−1 − n ψ̂k,g (nπ̂k ψ̂k,g )−1 = π̂k−1 − π̂k−1 = 0 (96)
g=1 g=1

∂E [ln(N (x |µ ,Λ−1 ))] ∂ ln(ρik )

It was shown in Titterington and Wang [2006] that µ,Λ i
∂ π̂k
k k
→ 0 as n → ∞. So ∂ π̂k as
in the continuous data case only depends on the Eπ [ln(πk )] term as n → ∞. Recall that

K
X K
X
Eπ [ln(πk )] = 𭟋(α̂k ) − 𭟋( α̂k ) = 𭟋(α + nπ̂k ) − 𭟋(Kα + n π̂k ) (97)
k=1 k=1

As in Titterington and Wang [2006], rather than Equation (97) we instead use 𭟋(α+nπ̂k )−𭟋(Kα+n),
PK
which coincides with Equation (97) when k=1 π̂k = 1. This is always satisfied by the CAVI algorithm
PK
after the first update (and satisfied by the initial values too if the rik are initialised to satisfy k=1 rik =
1), so T (Θ) still matches our proposed algorithm for Θ that lies in the support of the prior. This will
simplify some of the following calculations by ensuring that ∂ ln(ρ ik )
∂ π̂k′ = 0 if k ′ ̸= k 4 . By Lemma 3, we
have:
∂ ln(ρik ) ∂𭟋(α + nπ̂k ) ∂𭟋(Kα + n)
→ −
∂ π̂k ∂ π̂k ∂ π̂k
∂(α + nπ̂k )
∼ (α + nπ̂k − 1)−1
∂ π̂k
1
→n(nπ̂k )−1 = (98)
π̂k

Now that we have computed the limits of the derivatives of ln(ρik ) as n → ∞, we can return to
our original objective of deriving the matrix of derivatives in Equation (79) as n → ∞ at Θ = Θ∗ .
PK
Let ϕ∗k = p(x, c|z = k, Θ∗ ) and ϕ∗ = k′ =1 πk∗′ p(x, c|zi = k ′ , Θ∗ ) denote the likelihood evaluated at an
4 In
Pd
contrast, we did not take g=1 ψ̂k,g =1 in Equation (92), as this results in Equation (95), and its expectation
δ 1 ϕ∗
E[( ψ∗gc − Pd ∗ ) k
ϕ∗
] = 0, which will simplify some of the subsequent calculations.
k,g g=1 ψk,g

34
 Bayesian Variational Inference for Mixed Data Mixture Models

arbitrary point x, c with the true parameter values Θ∗ for a specific mixture component and with the
mixture component integrated out, respectively. Also let:
ϕ∗k
Rkπ =
ϕ∗
ϕ∗k πk∗
Rkµ = Λ∗k (x − µ∗k )
ϕ∗
1 ∗−1 ϕ∗ π ∗
RkΛ = [Λk − (x − µ∗k )(x − µ∗k )T ] k ∗ k
2 ϕ
∗ ∗
δg c
′ ϕ π
Rkψ : Rkg
ψ
′ = ( ∗ − 1) k ∗ k
ψk,g ′ ϕ
πk∗ ϕ∗k
Lπk =
ϕ∗
ϕ∗k
Lµk = (x − µ∗k )
ϕ∗
∗ ∗ ∗ ∗ T ∗ ϕ∗k
LΛ
k = [Λk − Λk (x − µk )(x − µk ) Λk ]
ϕ∗
∗ ∗
(δg′ c − ψk,g ′ )ϕk
Lψ ψ
k : Lkg ′ = (99)
ϕ∗
 θ
L1
θ θ T θ T θ  .. 
and R = [(R1 ) , . . . , (RK ) ] and L =  .  for θ ∈ {π, µ, Λ, ψ}.
LθK
One can then show the gradients inside Equation (79) take the following forms as n → ∞:
∇π̂ Π(Θ∗ ) → I − E[Lπ Rπ ]
∇µ̂ Π(Θ∗ ) → −E[Lπ Rµ ]
∇Λ̂ Π(Θ∗ ) → −E[Lπ RΛ ]
∇ψ̂ Π(Θ∗ ) → −E[Lπ Rψ ]
∇π̂ M (Θ∗ ) → −E[Lµ Rπ ]
∇µ̂ M (Θ∗ ) → I − E[Lµ Rµ ]
∇Λ̂ M (Θ∗ ) → −E[Lµ RΛ ]
∇ψ̂ M (Θ∗ ) → −E[Lµ Rψ ]
∇π̂ S(Θ∗ ) → −E[LΛ Rπ ]
∇µ̂ S(Θ∗ ) → −E[LΛ Rµ ]
∇Λ̂ S(Θ∗ ) → I − E[LΛ RΛ ]
∇ψ̂ S(Θ∗ ) → −E[LΛ Rψ ]
∇π̂ Ψ(Θ∗ ) → −E[Lψ Rπ ]
∇µ̂ Ψ(Θ∗ ) → −E[Lψ Rµ ]
∇Λ̂ Ψ(Θ∗ ) → −E[Lψ RΛ ]
∇ψ̂ Ψ(Θ∗ ) → I − C ψ − E[Lψ Rψ ] (100)
ψ ∗
where C ψ : RKd → RKd , Ckgk ′ g ′ = δkk ′ ψk,g . So if a ∈ R
Kd
, u ∈ RK :
K X
X d K X
X d d
X
ψ ∗ ∗
Ckgk ′ g ′ ak ′ g ′ = δkk′ ψk,g ak′ g′ = ψk,g akg′
k′ =1 g ′ =1 k′ =1 g ′ =1 g ′ =1

The expectation in Equation (100) is taken with respect to x, c, which has density of the likelihood
evaluated at the true data generating parameter Θ∗ . Calculating the gradient limits in Equation (100)
require very long and tedious calculations, so we reserve the details to Appendix D.2.

35
 Bayesian Variational Inference for Mixed Data Mixture Models

Returning to the main objective of showing the operator ∇T ϵ (Θ∗ ) converges as n → ∞ (with prob-
ability 1) to an operator Aϵ (Θ∗ ) with norm less than 1, where the norm is induced by ⟨., .⟩. To do
this, we show Aϵ (Θ∗ ) < I and Aϵ (Θ∗ ) > −I, where for example Aϵ (Θ∗ ) < I means I − Aϵ (Θ∗ ) is
positive definite. To see why this is sufficient to imply ∥Aϵ (Θ∗ )∥op < 1, note that Aϵ (Θ∗ ) < I im-
plies ⟨B, Aϵ (Θ∗ )B⟩ ≤ ⟨B, B⟩ for all B with equality iff B = 0, and similarly Aϵ (Θ∗ ) > −I implies
⟨B, Aϵ (Θ∗ )B⟩ ≥ −⟨B, B⟩ with equality iff B = 0. ϵ ∗
L Lthis implies |⟨B, A (Θ )B⟩| < |⟨B,ϵ B⟩|∗ = 1
L Combined
for all ∥B∥ = 1. Because the unit circle {B ∈ Sπ Sµ SΛ Sψ : ∥B∥ = 1} is compact, |⟨B, A (Θ )B⟩|
attains its supremum on the unit circle which is smaller than 1. Lastly note that Aϵ (Θ∗ ) is self adjoint
5
and recall from Hilbert space theory
L L that L
the operator norm for self-adjoint operators A is equal to
sup{|⟨B, AB⟩| : ∥B∥ = 1, B ∈ Sπ Sµ SΛ  Sψ}.
u
v L L L
Let u ∈ Sπ , v ∈ Sµ , W ∈ SΛ , a ∈ Sψ , B = W  ∈ Sπ
 Sµ SΛ Sψ .

L L L a
Let R : Sπ Sµ SΛ Sψ → R such that:
K
X K
X K
X K
X
R(B) = Rkπ uk + (Rkµ )T vk + T r(RkΛ Wk ) + (Rkψ )T ak (101)
k=1 k=1 k=1 k=1

Lπ
 Lµ 
LΛ . Note we deliberately constructed the inner product ⟨., .⟩ so that:
Similarly let L =  

Lψ
K K X π∗ K K
X 1 X X π∗
⟨B, L⟩ = uk ( ∗ )Lπk + (vk )T (πk∗ Λ∗k )Lµk + k
T r(Wk Λ∗−1
k LΛ ∗−1
k Λk )+ aTk diag( k∗ )Lψ
πk 2 ψk k
k=1 k=1 k=1 k=1
K
X K
X K
X K
X
= uk Rkπ + vkT Rkµ + T r(RkΛ Wk ) + aTk Rkψ
k=1 k=1 k=1 k=1
= R(B) (102)

For the case Aϵ (Θ∗ ) > −I. From Equation (100), we have that 1 ∗
= I − C ψ −LE[LR].
L AL(Θ )L L Here
L we
ψ
use a slight abuse of notation where we extend the map C : Sπ Sµ SΛ Sψ → Sπ Sµ SΛ Sψ
by simply mapping the extended dimensions to 0 so:

2
Aϵ (Θ∗ ) > −I ⇐⇒ (1 − ϵ)I + ϵA1 (Θ∗ ) > −I ⇐⇒ C ψ + E[LR] < I
ϵ
2
⇐⇒ ⟨B, C ψ B⟩ + E[⟨B, LR(B)⟩] ≤ ⟨B, B⟩ with equality iff B = 0 (103)
ϵ
where
K X
d d K d d K
X πk∗ X ∗ X
∗
X X X
ψ
⟨B, C B⟩ = akg ∗ ψk,g akg =
′ πk akg akg =
′ πk∗ (aTk 1)2 (104)
g=1
ψk,g ′ g=1 ′
k=1 g =1 k=1 g =1 k=1
K K X π∗K K
X 1 X
∗ ∗ ∗−1 ∗−1
X πk∗
⟨B, B⟩ = u2k ( ) + v T
k (π k Λ k )v k + k
T r(W k Λ W k Λ ) + a T
k diag( )ak
πk∗ 2 k k
ψk∗
k=1 k=1 k=1 k=1

where 1 is a vector of dimension d where all the entries are 1.

E[⟨B, LR(B)⟩] is a much more involved term, for which we will derive an upper bound for below:
E[⟨B, LR(B)⟩] = E[(R(B))R(B)] (105)
2
= E[R(B) ] (106)
5 One can easily show similar to Equation (102) that ⟨L, B⟩ = R(B), and so E[⟨LR(B ), B ⟩] = E[R(B )⟨L, B ⟩] =
1 2 1 2
E[R(B1 )R(B2 )] = E[⟨B1 , L⟩R(B2 )] = E[⟨B1 , LR(B2 )⟩]. Likewise, one can easily derive ⟨C ψ B1 , B2 ⟩ and ⟨B1 , C ψ B2 ⟩
ψ ψ
similar to Equation (104) and show ⟨C B1 , B2 ⟩ = ⟨B1 , C B2 ⟩

36
 Bayesian Variational Inference for Mixed Data Mixture Models

E[R(B)2 ]
K
nX πk∗ ϕ∗k 1 ∗ ∗ T 1 ∗−1 ∗ ∗ T δc T o2
= E[ ∗
(u k ( ∗ ) + (Λ k (x − µk )) v k + T r( [Λ k − (x − µk )(x − µk ) ]W k ) + ∗
ak − aTk 1) ]
ϕ πk 2 ψk
k=1
K
X π ∗ ϕ∗ 1 1 δc T
≤ E[ k k
(uk ( ∗ ) + (Λ∗k (x − µ∗k ))T vk + T r( [Λ∗−1 − (x − µ∗
k )(x − µ∗ T
k ) ]W k ) + ak − aTk 1)2 ]
ϕ ∗ πk 2 k ψk∗
k=1
K
X π ∗ ϕ∗ 1 1
= E[ k k
(u2k ( ∗ )2 + vkT Λ∗k (x − µ∗k )(x − µ∗k )T Λ∗k vk + T r( [Λ∗−1 − (x − µ∗k )(x − µ∗k )T ]Wk )2
ϕ ∗ πk 2 k
k=1
T
δc 1 1 1
+( ak − aTk 1)2 + 2uk ( )(x − µ∗k )T Λ∗k vk + 2uk ( ∗ )T r( [Λ∗−1 − (x − µ∗k )(x − µ∗k )T ]Wk )
ψk∗ πk∗ πk 2 k
T
1 δc 1
+ 2uk ( ∗ )( ak − aTk 1) + 2(x − µ∗k )T Λ∗k vk T r( [Λ∗−1 − (x − µ∗k )(x − µ∗k )T ]Wk )
πk ψk ∗ 2 k
T
δc 1 δc T
+ 2(x − µ∗k )T Λ∗k vk ( ak − aTk 1) + 2T r( [Λ∗−1 − (x − µ∗
k )(x − µ∗ T
k ) ]W k )( ak − aTk 1)] (107)
ψk∗ 2 k ψk∗
δ
Here ψδc∗ is a vector of dimension d with components ψ∗cg . The second inequality above is due to
k,g
PK k PK PK
| k=1 ξk ηk |2 ≤ k=1 ξk2 ηk for ηk ≥ 0 and k=1 ηk = 1, and any real ξk . One can show this inequality
√ √
by applying Cauchy-Schwarz on ξk ηk and ηk . There are a lot of terms here but it will turn out the
‘diagonal’ terms will simplify nicely, and the non diagonal terms will have expected value of 0.
K K K
X π ∗ ϕ∗
k k 2 1 2
X 1 2 X 2 1
E[ u ( ) ] = πk∗ u2k ( ) = uk ( ∗ )
ϕ∗ k πk∗ πk∗ πk
k=1 k=1 k=1
K K
X π ∗ ϕ∗ X
E[ k k
(vkT Λ∗k (x − µ∗k )(x − µ∗k )T Λ∗k vk )] = vkT πk∗ Λ∗k vk
ϕ∗
k=1 k=1
K K
X π ∗ ϕ∗ 1 X πk
E[ k k
(T r( [Λ∗−1 − (x − µ∗k )(x − µ∗k )T ]Wk )2 )] = T r(Λ∗−1 Wk Λ∗−1 Wk )
ϕ∗ 2 k 2 k k
k=1 k=1
K T K T T
X πk∗ ϕ∗k δc X πk∗ ϕ∗k T δc δc T δc
E[ (( ∗ ak − aTk 1)2 )] = E[ (ak ak − 2a k 1 ak + (aTk 1)2 )]
ϕ∗ ψk ϕ∗ ψk∗ ψk∗ ψk∗
k=1 k=1
K d d d
X πk∗ ϕ∗k X X δcg′ δcg T
X δcg T 2
= E[ ( akg′ ∗ ∗ akg − 2ak 1 ∗ akg + (ak 1) )]
ϕ∗ ψk,g′ ψk,g ψ
g=1 k,g
k=1 g ′ =1 g=1
K K
X πk∗ X
= aTk diag( ∗ )ak − πk∗ (aTk 1)2 (108)
ψk
k=1 k=1

The third expectation above can be derived by writing the Trace in summation convention and using
Isserlis’s theorem on
PKthe resulting fourth moment. Note the four ‘diagonal’ terms above adds up to
precisely ⟨B, B⟩ − k=1 πk∗ (aTk 1)2 , and the ‘off-diagonal’ terms add up to 0:
K
X π ∗ ϕ∗ 1
E[ k k
∗
2uk ( ∗ )(x − µ∗k )T Λ∗k vk ] = 0
ϕ πk
k=1
K
X π ∗ ϕ∗ 1 1
E[ k k
2uk ( )T r( [Λ∗−1 − (x − µ∗k )(x − µ∗k )T ]Wk )] = 0
ϕ∗ πk∗ 2 k
k=1
K T
X π ∗ ϕ∗ 1 δc
E[ k k
2uk ( ∗ )( ∗ ak − aTk 1)] = 0 (109)
ϕ∗ πk ψk
k=1

37
 Bayesian Variational Inference for Mixed Data Mixture Models

K
X π ∗ ϕ∗ 1
E[ k k
2(x − µ∗k )T Λ∗k vk T r( [Λ∗−1 − (x − µ∗k )(x − µ∗k )T ]Wk )] = 0 (110)
ϕ∗ 2 k
k=1

By Isserlis’s Theorem.
K T
X π ∗ ϕ∗ δc 1
E[ k k
2(x − µ∗k )T Λ∗k vk ( ak − aTk 1)] = 0
ϕ∗ ψk∗ 2
k=1
K T
X π ∗ ϕ∗ 1 δc 1
E[ k k
2T r( [Λ∗−1 − (x − µ∗k )(x − µ∗k )T ]Wk )( ak − aTk 1)] = 0 (111)
ϕ∗ 2 k ψk∗ 2
k=1

By independence between x and c. Putting together Equation (108) and Equation (109) we have:
K
X
E[⟨B, LR(B)⟩] ≤ ⟨B, B⟩ − πk∗ (aTk 1)2 (112)
k=1

With an upper bound on E[⟨B, LR(B)⟩] we can now return to Equation (103):

⟨B, C ψ B⟩ + E[⟨B, LR(B)⟩]

K
X K
X
≤ πk∗ (aTk 1)2 + ⟨B, B⟩ − πk∗ (aTk 1)2
k=1 k=1
= ⟨B, B⟩ (113)
So the desired inequality is satisfied by picking 0 < ϵ < 2, which is enough for our CAVI algorithm
(ϵ = 1).
For the case Aϵ (Θ∗ ) < I. Again from Equation (100), we have that A1 (Θ∗ ) = I − C ψ − E[LR], so:

Aϵ (Θ∗ ) < I ⇐⇒ (1 − ϵ)I + ϵA1 (Θ∗ ) − I < 0 ⇐⇒ C ψ + E[LR] > 0 (114)

which is obviously true as E[LR] is positive definite and C ψ positive semidefinite under ⟨, ⟩
We have established ∇T ϵ (Θ∗ ) converges with probability 1 as n → ∞ to an operator with operator
norm less than 1, therefore there exist 0 ≤ λ < 1 such that:

∥Θ(t+1) − Θ∗ )∥ ≤ ∥T ϵ (Θ(t) ) − T ϵ (Θ∗ )∥ + ∥T ϵ (Θ∗ ) − Θ∗ ∥

≤ λ∥Θ(t) − Θ∗ ∥ + ∥T ϵ (Θ∗ ) − Θ∗ ∥ (115)
Calculations in Appendix D.2 shows that T ϵ (Θ∗ ) → Θ∗ as n → ∞, therefore, the iterative procedure
Θ = T ϵ (Θ(t) ) converges locally to Θ∗ as n → ∞.
(t+1)

D.1 Proof of Corollary 1

Using Equation (81), at iteration t, the variational posterior means converge take the following forms,
and converge to the limits shown below as n → ∞:

α̂t nπ̂kt + α
Eπ [πk ] = PK k = PK → πk∗
t n k=1 π̂kt + Kα
k=1 α̂k
Eµ,Λ [µk ] = m̂tk = (nµ̂tk π̂kt + βmk )/(nπ̂kt + β) → µ∗k
1 Φ −1
Eµ,Λ [Λk ] = ν̂kt (Φ̂tk )−1 = (nπ̂kt + ν)(nπ̂kt )−1 ((Λ̂tk )−1 + β(µ̂tk − mk )(µ̂tk − mk )T + )
nπ̂kt+β nπ̂k
→ Λ∗k
t
η̂k,g nπ̂kt ψ̂k,g
t
+η ∗
Eψ [ψk,g ] = Pd = d
→ ψk,g (116)
t t t
P
η̂
g=1 k,g nπ̂ k ψ̂
g=1 k,g + dη

38
 Bayesian Variational Inference for Mixed Data Mixture Models

D.2 Details of gradient limit calculations

We start with the gradients of Ψ. As before, let θk′ be any scaler variable in Θ belonging to component
k ′ . We have:
Pn
∂Ψk (Θ)g ∂ I{c =g} rik
= ( i=1
Pn i )
∂θk′ ∂θk′ i=1 rik
Pn ∂rik Pn Pn ∂rik Pn
( i=1 I{ci =g} ∂θ )( i=1 rik ) − ( i=1 ∂θ )( i=1 I{ci =g} rik )
k′ k′
= Pn 2
(117)
( i=1 rik )

Starting with the case θk′ = ψ̂k′ ,g′ . We consider the two terms in the last line as n → ∞ at Θ = Θ∗ .
Dividing each term in the top and bottom by n:
n n
1 X ∂rik (Θ∗ ) 1 X ∂ ln(ρik′ (Θ∗ )) n o
= δkk′ rik (Θ∗ ) − rik (Θ∗ )rik′ (Θ∗ )
n i=1 ∂ ψ̂k′ ,g′ n i=1 dψ̂k′ ,g′
δg ′ c π ∗ ϕ∗ π ∗ ϕ∗ π ∗′ ϕ∗ ′
→ E[( ∗ − 1)(δkk′ k ∗ k − k k ∗2k k )]
ψk′ ,g′ ϕ ϕ
∗
ψk,g ′ δg ′ c π ∗ ϕ∗ π ∗′ ϕ∗ ′ δg ′ c π ∗ ϕ∗ π ∗′ ϕ∗ ′
=( ∗ − 1)δkk′ πk∗ − E[( ∗ − 1) k k ∗2k k ] = −E[( ∗ − 1) k k ∗2k k ] (118)
ψk′ ,g′ ψk′ ,g′ ϕ ψk′ ,g′ ϕ

where in the second line, we used the limit of rik (Θ∗ ) previously derived in Lemma 2, and the strong
law of large numbers, with the expectation here is taken with respect to x, c, which has density of the
δ π∗ δ ′ ϕ∗
likelihood evaluated at the true data generating parameter Θ∗ . In the last line we used ψkk∗ ′ k E[ gϕc∗ k ] =
k′ ,g ′
∗
δkk′ πk ∗
ψk∗ ψk,g ′ = δkk′ πk∗ . Similarly:
′ ,g ′

n n
1X ∂rik (Θ∗ ) 1X ∂ ln(ρik′ (Θ∗ )) n o
I{ci =g} = I{ci =g} δkk′ rik (Θ∗ ) − rik (Θ∗ )rik′ (Θ∗ )
n i=1 ∂ψk′ ,g′ n i=1 ∂ψk′ ,g′
δg ′ c π ∗ ϕ∗ π ∗ ϕ∗ π ∗′ ϕ∗ ′
→ E[δgc ( ∗ − 1)(δkk′ k ∗ k − k k ∗2k k )]
ψk′ ,g′ ϕ ϕ
δg ′ c π ∗ ϕ∗ π ∗′ ϕ∗ ′
→ δgg′ δkk′ πk∗ − δkk′ ψk,g
∗
πk∗ − E[δgc ( ∗ − 1) k k ∗2k k ] (119)
ψk′ ,g′ ϕ

where in the last line we used δgc δg′ c = δgg′ δg′ c . Additionally:
n
1X π ∗ ϕ∗
Πk (Θ∗ ) = rik (Θ∗ ) → E[ k ∗ k ] = πk∗
n i=1 ϕ
n
1X π ∗ ϕ∗
I{ci =g} rik (Θ∗ ) → E[δgc k ∗ k ] = πk∗ ψk,g
∗
(120)
n i=1 ϕ

Consequently, Ψk (Θ∗ )g → ψk,g

∗
. Putting everything together and plugging into Equation (117) and we
get:
∗
∂Ψk (Θ∗ )g ∗ δg ′ c ϕ∗ π ∗′ ϕ∗ ′ δg′ c ψk,g ϕ∗k πk∗′ ϕ∗k′
→ δgg′ δkk′ − δkk′ ψk,g − E[δgc ( ∗ − 1) k k∗2 k ] + E[( ∗ − 1) ]
∂ψk′ ,g′ ψk′ ,g′ ϕ ψk′ ,g′ ϕ∗2
∗
∗
(δgc − ψk,g )ϕ∗k δg′ c πk∗′ ϕ∗k′
= δgg′ δkk′ − δkk′ ψk,g − E[ ( − 1) ] (121)
ϕ∗ ψk∗′ ,g′ ϕ∗

∂Ψk (Θ)g ∂Ψk (Θ)g ∂Ψk (Θ)g

∂ π̂k′ , ∂ µ̂k′ l , ∂ Λ̂k′ jl
can be calculated similarly to Equation (117) and replacing θk′ with the
corresponding variable to be differentiated:

39
 Bayesian Variational Inference for Mixed Data Mixture Models

∂Ψk (Θ∗ )g 1 π ∗ ϕ∗ π ∗ ϕ∗ π ∗′ ϕ∗ ′ πk∗

→ E[δgc ( ∗ )(δkk′ k ∗ k − k k ∗2k k )] ∗2
∂ π̂k′ πk ′ ϕ ϕ πk
∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗
1 π ϕ π ϕ π ′ϕ ′ π ψ
k k,g
− E[( ∗ )(δkk′ k ∗ k − k k ∗2k k )] ∗2
πk ′ ϕ ϕ πk
∗ ∗ ∗
(δgc − ψk,g )ϕk ϕk′
= −E[ ] (122)
ϕ∗ ϕ∗
Pq
Using Einstein’s summation convention j=1 Λ∗k′ lj (xj − µ∗k′ j ) = Λ∗k′ lj (xj − µ∗k′ j ), we have:

∂Ψk (Θ∗ )g π ∗ ϕ∗ π ∗ ϕ∗ π ∗′ ϕ∗ ′ πk∗

→ E[δgc Λ∗k′ lj (xj − µ∗k′ j )(δkk′ k ∗ k − k k ∗2k k )] ∗2
∂ µ̂k′ l ϕ ϕ πk
∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗
π ϕ π ϕ π ′ϕ ′ π ψ
k k,g
− E[Λ∗k′ lj (xj − µ∗k′ j )(δkk′ k ∗ k − k k ∗2k k )] ∗2
ϕ ϕ πk
∗ ∗ ∗ ∗
(δgc − ψk,g )ϕk ∗ ϕ ′π ′
= −E[ Λk′ lj (xj − µ∗k′ j ) k ∗ k ] (123)
ϕ∗ ϕ

∂Ψk (Θ∗ )g 1 πk∗ ϕ∗k πk∗ ϕ∗k πk∗′ ϕ∗k′ πk∗

→ E[δgc (Λ∗−1 ∗ ∗
k′ jl − (xj − µk′ j )(xl − µk′ l ))(δkk
′ − )] ∗2
∂ Λ̂k′ jl 2 ϕ∗ ϕ∗2 πk
1 πk∗ ϕ∗k πk∗ ϕ∗k πk∗′ ϕ∗k′ πk∗ ψk,g
∗
− E[ (Λ∗−1
′ − (x j − µ ∗
′
k j )(x l − µ∗
′
k l ))(δkk ′ − )]
2 k jl ϕ∗ ϕ∗2 πk∗2
∗
(δgc − ψk,g )ϕ∗k 1 ∗−1 ϕ∗ ′ π ∗′
= −E[ ∗
[Λk′ jl − (xj − µ∗k′ j )(xl − µ∗k′ l )] k ∗ k ] (124)
ϕ 2 ϕ
Next, for the gradients of Π:

n
∂Πk (Θ∗ ) ∂ 1X
= ( rik (Θ∗ ))
∂ ψ̂k′ ,g′ ∂ ψ̂k′ ,g′ n i=1
δg ′ c π ∗ ϕ∗ π ∗′ ϕ∗ ′
→ −E[( ∗ − 1) k k ∗2k k ] (125)
ψk′ ,g′ ϕ

Recall this is the same as Equation (118).

n
∂Πk (Θ∗ ) ∂ 1X
= ( rik (Θ∗ ))
∂ π̂k′ ∂ π̂k′ n i=1
1 π ∗ ϕ∗ π ∗ ϕ∗ π ∗′ ϕ∗ ′
→ E[ ∗ (δkk′ k ∗ k − k k ∗2k k )]
πk ′ ϕ ϕ
∗ ∗ ∗
π ϕk ϕk ′
= δkk′ − E[ k ∗2 ] (126)
ϕ

n
∂Πk (Θ∗ ) ∂ 1X
= ( rik (Θ∗ ))
∂ µ̂k′ l ∂ µ̂k′ l n i=1
πk∗ ϕ∗k π ∗ ϕ∗ π ∗′ ϕ∗ ′
→ E[Λ∗k′ lj (xj − µ∗k′ j )(δkk′ ∗
− k k ∗2k k )]
ϕ ϕ
∗ ∗ ∗ ∗
π ϕ π ′ ϕ ′
= −E[Λ∗k′ lj (xj − µ∗k′ j ) k k ∗2k k )] (127)
ϕ

40
 Bayesian Variational Inference for Mixed Data Mixture Models

n
∂Πk (Θ∗ ) ∂ 1X
= ( rik (Θ∗ ))
∂ Λ̂k′ jl ∂ Λ̂k′ jl n i=1
1 πk∗ ϕ∗k π ∗ ϕ∗ π ∗′ ϕ∗ ′
→ E[ (Λ∗−1 ∗ ∗
k′ jl − (xj − µk′ j )(xl − µk′ l ))(δkk
′
∗
− k k ∗2k k )]
2 ϕ ϕ
∗ ∗ ∗ ∗
1 ∗−1 π ϕ π ′ ϕ ′
= −E[ (Λk′ jl − (xj − µ∗k′ j )(xl − µ∗k′ l )) k k ∗2k k )] (128)
2 ϕ
Next, for the gradients of M :

Pn
∂Mk (Θ) ∂ rik xi
= ( Pi=1
n )
∂θk′ ∂θk′ i=1 rik
Pn ∂rik Pn Pn ∂rik Pn
( i=1 ∂θ xi )( i=1 rik ) − ( i=1 ∂θk′ )( i=1 rik xi )
k′
= Pn (129)
( i=1 rik )2

Note
∗ ∗
Pn πk ϕk x
ϕ∗ ]
∗
∗rik (Θ )x i E[
Mk (Θ ) = Pi=1
n ∗)
→ πk∗ ϕ∗ = µ∗k (130)
r
i=1 ik (Θ E[ ϕ∗ ] k

For θk′ = ψ̂k′ ,g′ , we have:

∂Mk (Θ∗ ) 1 n δg ′ c πk∗ ϕ∗k πk∗ ϕ∗k πk∗′ ϕ∗k′ ∗ δg ′ c πk∗ ϕ∗k πk∗ ϕ∗k πk∗′ ϕ∗k′ ∗ ∗ o
→ E[( − 1)(δ kk ′ − )x]π k − E[( − 1)(δ kk ′ − )]πk µk
∂ ψ̂k′ ,g′ πk∗2 ψk∗′ ,g′ ϕ∗ ϕ∗2 ψk∗′ ,g′ ϕ∗ ϕ∗2
δg ′ c ϕ∗ (x − µ∗k )ϕ∗k δg′ c πk∗′ ϕ∗k′
= δkk′ E[( ∗ − 1) k∗ (x − µ∗k )] − E[ ( ∗ − 1) ]
ψk′ ,g′ ϕ ϕ∗ ψk′ ,g′ ϕ∗
(x − µ∗k )ϕ∗k δg′ c πk∗′ ϕ∗k′
= −E[ ( − 1) ] (131)
ϕ∗ ψk∗′ ,g′ ϕ∗

For θk′ = π̂k′ , we have:

∂Mk (Θ∗ ) 1 n 1 π ∗ ϕ∗ π ∗ ϕ∗ π ∗′ ϕ∗ ′ 1 π ∗ ϕ∗ π ∗ ϕ∗ π ∗′ ϕ∗ ′ o
→ ∗2 E[ ∗ (δkk′ k ∗ k − k k ∗2k k )x]πk∗ − E[ ∗ (δkk′ k ∗ k − k k ∗2k k )]πk∗ µ∗k
∂ π̂k′ πk πk ′ ϕ ϕ πk ′ ϕ ϕ
1 ϕ∗k (x − µ ∗
)ϕ ∗
1 π ∗ ∗
′ ϕ ′
= δkk′ E[ ∗ ∗ (x − µ∗k )] − E[ k k
( ∗ ) k ∗k ]
πk ′ ϕ ϕ∗ πk ′ ϕ
(x − µ∗k )ϕ∗k ϕ∗k′
= −E[ ] (132)
ϕ∗ ϕ∗
For θk′ = µ̂k′ l , we have:

∂Mk (Θ∗ )r 1 n π ∗ ϕ∗ π ∗ ϕ∗ π ∗′ ϕ∗ ′
→ ∗2 E[Λ∗k′ lj (xj − µ∗k′ j )(δkk′ k ∗ k − k k ∗2k k )xr ]πk∗
∂ µ̂k′ l πk ϕ ϕ
∗ ∗ ∗ ∗ ∗ ∗
π ϕ π ϕ π ′ϕ ′ o
− E[Λ∗k′ lj (xj − µ∗k′ j )(δkk′ k ∗ k − k k ∗2k k )]πk∗ µ∗kr
ϕ ϕ
∗
ϕ (xr − µ∗kr )ϕ∗k ∗ πk∗′ ϕ∗k′
= δkk′ E[Λ∗k′ lj (xj − µ∗k′ j ) k∗ (xr − µ∗kr )] − E[ (Λ ′
k lj (x j − µ∗
′
k j )) ]
ϕ ϕ∗ ϕ∗
(xr − µ∗kr )ϕ∗k ∗ ∗ πk∗′ ϕ∗k′
= δkk′ δlr − E[ (Λ k ′ lj (xj − µk ′ j )) ] (133)
ϕ∗ ϕ∗

41
 Bayesian Variational Inference for Mixed Data Mixture Models

For θk′ = Λ̂k′ jl , we have:

∂Mk (Θ∗ )r 1 n 1 ∗−1 ∗ ∗ πk∗ ϕ∗k πk∗ ϕ∗k πk∗′ ϕ∗k′

→ E[ (Λ ′ − (x j − µ ′
k j )(x l − µ ′
k l ))(δ kk ′ − )xr ]πk∗
∂ Λ̂k′ jl πk∗2 2 k jl ϕ∗ ϕ∗2
1 πk∗ ϕ∗k πk∗ ϕ∗k πk∗′ ϕ∗k′ ∗ ∗ o
− E[ (Λ∗−1 ′ − (x j − µ ∗
′
k j )(x l − µ ∗
′
k l ))(δ kk ′ − )]πk µkr
2 k jl ϕ∗ ϕ∗2
1 ϕ∗k
= δkk′ E[ (Λ∗−1 ∗ ∗
k′ jl − (xj − µk′ j )(xl − µk′ l )) ∗ (xr − µkr )]
∗
2 ϕ
(xr − µ∗kr )ϕ∗k 1 ∗−1 π ∗′ ϕ∗ ′
− E[ ∗
( [Λk′ jl − (xj − µ∗k′ j )(xl − µ∗k′ l )]) k ∗ k ]
ϕ 2 ϕ
∗ ∗
(xr − µkr )ϕk 1 ∗−1 π ∗′ ϕ∗ ′
= −E[ ∗
( [Λk′ jl − (xj − µ∗k′ j )(xl − µ∗k′ l )]) k ∗ k ] (134)
ϕ 2 ϕ
Finally, for the gradients of S:

n n
∂Sk (Θ) ∂ nX X o
= ( rik )( rik (xi − Mk )(xi − Mk )T )−1
∂θk′ ∂θk′ i=1 i=1
n n
X ∂rik X
=( )( rik (xi − Mk )(xi − Mk )T )−1
i=1
∂θk ′
i=1
n n n n
X X X ∂ n o X
−( rik )( rik (xi − Mk )(xi − Mk )T )−1 ( rik (xi − Mk )(xi − Mk )T )( rik (xi − Mk )(xi − Mk )T )−1
i=1 i=1 i=1
∂θk ′
i=1
(135)

where we used Mk = Mk (Θ) as a shorthand. Before continuing we need to evaluate some intermediate
limits:
n
1X
rik (Θ∗ )(xi − Mk (Θ∗ ))(xi − Mk (Θ∗ ))T
n i=1
n
1X
= rik (Θ∗ )(xi − µ∗k + µ∗k − Mk (Θ∗ ))(xi − µ∗k + µ∗k − Mk (Θ∗ ))T
n i=1
n n
1X 1X
= rik (Θ∗ )(xi − µ∗k )(xi − µ∗k )T + rik (Θ∗ )(xi − µ∗k )(µ∗k − Mk (Θ∗ ))T
n i=1 n i=1
n n
1X 1X
+ (µ∗k − Mk (Θ∗ )) rik (Θ∗ )(xi − µ∗k )T + (µ∗k − Mk (Θ∗ ))( rik (Θ∗ ))(µ∗k − Mk (Θ∗ ))T
n i=1 n i=1
πk∗ ϕ∗k
→ E[ (x − µ∗k )(x − µ∗k )T ] = πk∗ (Λ∗k )−1 (136)
ϕ∗
Pn
Consequently, Sk (Θ∗ ) → Λ∗k . Using a similar argument, 1
n
∗
i=1 rik (Θ )(xi − Mk (Θ∗ )) → 0

42
 Bayesian Variational Inference for Mixed Data Mixture Models

n
1X ∂ n o
rik (xi − Mk )(xi − Mk )T )
n i=1 ∂θk′
n n n
1 X ∂rik −∂Mk 1 X 1X −∂Mk T
= (xi − Mk )(xi − Mk )T + ( ) rik (xi − Mk )T + rik (xi − Mk )( )
n i=1 ∂θk′ ∂θk′ n i=1 n i=1 ∂θk′
n n
1 X ∂rik 1 X ∂rik
= (xi − µ∗k )(xi − µ∗k )T + (µ∗k − Mk ) (xi − µ∗k )T
n i=1 ∂θk′ n i=1 ∂θk′
n n
1 X ∂rik 1 X ∂rik ∗
+ (xi − µ∗k )(µ∗k − Mk )T + (µ∗k − Mk ) (µ − Mk )T
n i=1 ∂θk′ n i=1 ∂θk′ k
n n
−∂Mk 1 X 1X −∂Mk T
+( ) rik (xi − Mk )T + rik (xi − Mk )( )
∂θk′ n i=1 n i=1 ∂θk′
n
1 X ∂rik
∼ (xi − µ∗k )(xi − µ∗k )T at Θ = Θ∗ , n → ∞ (137)
n i=1 ∂θk′

Putting everything together, for the case θk′ = ψ̂k′ ,g′ , we therefore have:

∂Sk (Θ∗ ) δg ′ c π ∗ ϕ∗ π ∗ ϕ∗ π ∗′ ϕ∗ ′
→ E[( ∗ − 1)(δkk′ k ∗ k − k k ∗2k k )](πk∗−1 Λ∗k )
∂ ψ̂k′ ,g′ ψk′ ,g′ ϕ ϕ
δg ′ c πk∗ ϕ∗k πk∗ ϕ∗k πk∗′ ϕ∗k′
− πk∗ (πk∗−1 Λ∗k )E[( − 1)(δ kk ′ − )(x − µ∗k )(x − µ∗k )T ](πk∗−1 Λ∗k )
ψk∗′ ,g′ ϕ∗ ϕ∗2
δg ′ c ϕ∗k ∗
= δkk′ E[( − 1) (Λ − Λ∗k (x − µ∗k )(x − µ∗k )T Λ∗k )]
ψk∗′ ,g′ ϕ∗ k
ϕ∗k ∗ δg ′ c π ∗′ ϕ∗ ′
− E[ ∗
(Λk − Λ∗k (x − µ∗k )(x − µ∗k )T Λ∗k )( ∗ − 1)( k ∗ k )]
ϕ ψk′ ,g′ ϕ
ϕ∗k ∗ δg ′ c π ∗′ ϕ∗ ′
= −E[ ∗
(Λk − Λ∗k (x − µ∗k )(x − µ∗k )T Λ∗k )( ∗ − 1)( k ∗ k )] (138)
ϕ ψk′ ,g′ ϕ

For θk′ = π̂k′ , we have:

∂Sk (Θ∗ ) 1 π ∗ ϕ∗ π ∗ ϕ∗ π ∗′ ϕ∗ ′
→ E[ ∗ (δkk′ k ∗ k − k k ∗2k k )](πk∗−1 Λ∗k )
∂ π̂k′ πk ′ ϕ ϕ
∗ ∗ ∗ ∗ ∗ ∗
1 π ϕ π ϕ π ′ϕ ′
− πk∗ (πk∗−1 Λ∗k )E[( ∗ )(δkk′ k ∗ k − k k ∗2k k )(x − µ∗k )(x − µ∗k )T ](πk∗−1 Λ∗k )
πk ′ ϕ ϕ
∗ ∗
ϕ ϕ ′
= −E[(Λ∗k − Λ∗k (x − µ∗k )(x − µ∗k )T Λ∗k ) k ∗2k ] (139)
ϕ
For θk′ = µ̂k′ l , we have:

∂Sk (Θ∗ ) π ∗ ϕ∗ π ∗ ϕ∗ π ∗′ ϕ∗ ′
→ E[Λ∗k′ lj (xj − µ∗k′ j )(δkk′ k ∗ k − k k ∗2k k )](πk∗−1 Λ∗k )
∂ µ̂k′ l ϕ ϕ
∗ ∗ ∗ ∗ ∗ ∗
π ϕ π ϕ π ′ϕ ′
− πk∗ (πk∗−1 Λ∗k )E[Λ∗k′ lj (xj − µ∗k′ j )(δkk′ k ∗ k − k k ∗2k k )(x − µ∗k )(x − µ∗k )T ](πk∗−1 Λ∗k )
ϕ ϕ
ϕ∗ π ∗′ ϕ∗ ′
= −E[(Λ∗k − Λ∗k (x − µ∗k )(x − µ∗k )T Λ∗k )Λ∗k′ lj (xj − µ∗k′ j ) k k∗2 k ] (140)
ϕ

For θk′ = Λ̂k′ jl , we have:

43
 Bayesian Variational Inference for Mixed Data Mixture Models

∂Sk (Θ∗ )rs 1 πk∗ ϕ∗k πk∗ ϕ∗k πk∗′ ϕ∗k′

→ E[ (Λ∗−1
′ − (xj − µ∗
k
∗
′ j )(xl − µk ′ l ))(δkk ′ − )](πk∗−1 Λ∗krs )
∂ Λ̂k′ jl 2 k jl ϕ∗ ϕ∗2
1 πk∗ ϕ∗k πk∗ ϕ∗k πk∗′ ϕ∗k′
− πk∗ (πk∗−1 Λ∗krc )E[ (Λ∗−1
′ − (xj − µ ∗
′
k j )(x l − µ∗
′
k l ))(δ kk ′ − )(xc − µ∗kc )(xd − µ∗kd )](πk∗−1 Λ∗kds )
2 k jl ϕ∗ ϕ∗2
1
= −δkk′ Λ∗krc [Λ∗−1 Λ∗−1 − (Λ∗−1 ∗−1 ∗−1 ∗−1 ∗−1 ∗−1
kjl Λkcd + Λkjc Λkld + Λkjd Λklc )](Λkds )
∗
2 kjl kcd
1 ϕ∗k πk∗′ ϕ∗k′
− E[(Λ∗krs − Λ∗krc (xc − µ∗kc )(xd − µ∗kd )Λ∗kds ) (Λ∗−1 ∗ ∗
k′ jl − (xj − µk′ j )(xl − µk′ l )) ]
2 ϕ∗2
1
= δkk′ (δrj δls + δjs δrl )
2
1 ϕ∗k πk∗′ ϕ∗k′
− E[(Λ∗krs − Λ∗krc (xc − µ∗kc )(xd − µ∗kd )Λ∗kds ) (Λ∗−1 ∗ ∗
k′ jl − (xj − µk′ j )(xl − µk′ l )) ] (141)
2 ϕ∗2

E Detail of data generating parameters for simulated examples

E.1 Scenario 1: a continuous data dominated example
For this scenario, we simulate data as follows: π ∗ are generated from normalising a random vector of
size K ∗ = 5, where each component is generated from U (0.5, 2.0). This ensures the ratio between the
most likely component and the least likely component is capped to at most 4 and not too skewed. We
use a single categorical variable with integer values from 1 to d1 = 4, where the true probability ψ ∗
for the categorical variable is generated from a Dirichlet distribution with all d1 hyperparameters equal
to 5. To ensure the continuous data are well separated between components, we generate K ∗ vectors
µ∗k that are a distance of at least 4K ∗ and at most 8K ∗ apart. To generate Σ∗k , for each k we sample
diagonal matrices Dk with eigenvalues between k 2 /2 and k 2 (so ‘standard deviation’ of at most k ∗ ),
then transforming it via conjugation with a random orthogonal matrix Uk to obtain a non diagonal
covariance matrix Σ∗k = Uk Dk UkT . This procedure generates components where the separation between
the component centres is relatively large compared to the standard deviations of the covariance matrices,
so the data sampled from different components are likely to be well separated.

E.2 A categorical data dominated example

For this scenario, again π ∗ are generated from normalising a random vector of size K ∗ = 5, where
each component is generated from U (0.5, 2). There are p = 2K ∗ binary categorical variables, which we
generate by having them correlate strongly with the component label as follows: if the component label is
equal to k, then the (2k − 1)th and 2kth categorical variables have 0.9 probability of equalling 1, while all
∗ ∗
other categorical variables have 0.9 probability of equalling 0. In other words, ψk,2k−1,1 = 0.9, ψk,2k,1 =
∗ ∗
0.9, ψk,j,1 = 0.1 for j ̸= 2k − 1, 2k. For the continuous data, we set q = K , with the component centres
µ∗k = ek , where ek is the unit vector in the kth direction. Each Σ∗k is chosen to be a diagonal matrix,
with the kth diagonal set to 9√while other diagonals are set to 4. In this example, the distance between
any two component centres is 2, but the standard deviations of the covariance matrices are 2 or 3. This
leads to the points generated from different components to overlap significantly with each other while
still being distinct, to ensure the components are largely driven by the value of the categorical variables.

E.3 A third example

For this scenario, the true component probabilities π ∗ are generated as for the two previous scenarios.
We define a single categorical variable (i.e. p = 1) with integer values between 1 and d1 = K ∗ that
∗ ∗
significantly correlates with the component label. In particular, for each k, ψk,1,k = 0.75 while ψk,1,l =
0.25/(d1 − 1) for all l ̸= k, so under the likelihood the categorical variable satisfies P (c = z) = [Link]
the continuous data, we set q = 5, µ∗k = 4K ∗ ek . Σ∗k is generated as in scenario 1, except the eigenvalue
of the diagonal matrices Dk is sampled uniformly to be between (1.6K ∗ )2 /2 and (1.6K ∗ )2 , significantly
higher than in example 1, then transformed via conjugation with an orthogonal matrix to obtain a non
diagonal matrix. The overall idea of this example is to represent a middle ground between scenario 1 and

44
 Bayesian Variational Inference for Mixed Data Mixture Models

2: have components that are less well separated in continuous space than scenario 1, in the sense that
the standard deviations of the covariance matrices is relatively larger in size compared to the separation
between the component centres, but compensated by a categorical variable that is significantly correlated
with the component label z, though not as much as in scenario 2.

45