11institutetext: Rotman School of Management, University of Toronto
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
# Learning Consumer Preferences from Bundle Sales Data

Ningyuan Chen  Setareh Farajollahzadeh  Guan Wang

###### Abstract

Product bundling is a common selling mechanism used in online retailing. To
set profitable bundle prices, the seller needs to learn consumer preferences
from the transaction data. When customers purchase bundles or multiple
products, classical methods such as discrete choice models cannot be used to
estimate customers’ valuations. In this paper, we propose an approach to learn
the distribution of consumers’ valuations toward the products using bundle
sales data. The approach reduces it to an estimation problem where the samples
are censored by polyhedral regions. Using the EM algorithm and Monte Carlo
simulation, our approach can recover the distribution of consumers’
valuations. The framework allows for unobserved no-purchases and clustered
market segments. We provide theoretical results on the identifiability of the
probability model and the convergence of the EM algorithm. The performance of
the approach is also demonstrated numerically.

###### Keywords:

EM algorithm bundlecensored demand clustering

##  1 Introduction

Bundles are pervasive in online retailing, fast-food restaurants, and the
video game industry. Companies bundle products and services together at a
bargain price to boost the demand of customers. Many firms in various
markets—from combo deals in McDonald’s, discounts for home and auto policies
combined for insurance companies, to deals for flights as a package of a
price, the number of baggage and meals for airlines—use bundling as one of
their main promotion strategies. Due to the popularity of bundles, how to set
profitable bundle pricing strategies has been a focus of the literature since
Stigler (1963). A number of simple bundle pricing strategies have shown
promising theoretical and empirical performance, such as grand bundles
(bundling all products together at a discount) and bundle size pricing
(pricing the bundle based on its size only). See Chu, Leslie, and Sorensen
(2011) for a comprehensive empirical comparison.

A prerequisite for _almost all_ successful bundle pricing schemes is the
accurate estimation of consumers’ reservation prices or valuations. That is,
how consumers react to a price menu of products and bundles in the form of
aggregate market demand. A typical source of data to learn such information is
the past transaction data that the firm has collected. However, compared to
the number of studies on how to set prices for bundles assuming _known_ market
demand or valuation distributions, there is scant literature on how to learn
the valuations from the transaction data of bundle sales. This is a stark
contrast to the industry practice: according to Fisher et al. (2020),
companies make more than 90% of the effort estimating demand and 10% setting
optimal prices.

This study focuses on learning consumer valuations from the transaction data
of bundle sales, which allows the firm to set profitable bundle pricing
strategies. In a sense, our study can be viewed as a preprocessing step for
developing bundle pricing strategies. The firm observes the behavior of past
customers: the products and bundles that are offered at the time, the charged
prices, and the final purchase decisions. Based on the information, the firm
attempts to learn the valuations of the customers in the data to forecast the
shopping behavior of future customers.

The challenges of this learning problem arise from missing values, censored
data, and bundle promotion. We next summarize each of these challenges and the
reason for the insufficiency of existing approaches. First, customers’ value
for products is missing values in sales data. While customers’ purchase
decisions are observed in sales data, the realized valuations are not. That
is, the firm only observes that a past customer is willing to pay $5 for a
t-shirt, but it doesn’t know if the customer values the t-shirt at $7 or $10.
Second, demand censoring, which refers to the common issue that customers who
don’t make any purchases are not observed in the sales data, also plagues the
studied estimation problem. The first two challenges are common for estimation
problems based on sales data. The method used to uncensored data for discrete
choice models without bundles cannot be directly applied. For transaction data
without bundles, discrete choice models such as the multinomial logit (MNL)
model tackle the challenge by linking the purchase probabilities to the
parameters without observing the valuations directly. In the presence of
bundles, discrete choice models can be applied if one views a bundle as an
independent product. But such an arrangement neglects the connection between
bundles and the component products. This is another unique problem when the
transaction data include bundle sales. Therefore, the third challenge of
products commonly bundled together hurdle the use of existing methods. When
products are bundled together, the valuations for these products are
entangled. When a customer purchases the bundle, it is unclear which product
the customer has valued higher. Therefore, it is challenging for the firm to
learn products’ valuations among customers separately.

In this paper, we construct a simple model for consumer valuations and provide
an algorithm to estimate the parameters from the data. In particular, we
assume the valuations customers have toward the products have a multivariate
Gaussian distribution. For a bundle, the valuation is the sum of the
valuations of the component products. Such an additive utility model is
commonly used in the bundle pricing literature. We formulate the problem as a
maximum likelihood estimation with incomplete data where customers’ valuation
samples are not observed. We then design an algorithm combining the
expectation-maximization (EM) algorithm and Monte Carlo simulation to estimate
the mean and covariance matrix of the Gaussian distribution. The contribution
of the paper is threefold.

  * •

We show the equivalence between the problem and another estimation problem
that is of independent interest. We find the transaction data can be converted
to censored samples of customer valuations. More precisely, the transaction
data does not reveal the exact location of the valuations of each customer,
but record a polyhedral region with linear boundaries in which the sample
falls in. A special case of this problem, interval-censored data, is commonly
seen in survival analysis (Lindsey and Ryan 1998). Our formulation generalizes
the interval to polyhedra in high dimensions.

  * •

The framework and algorithm we propose are flexible enough to accommodate many
realistic considerations. We extend the framework to allow for Gaussian
mixture models (GMM). Therefore, customers may have clustered valuations due
to market segmentation and any distribution can be approximated. Our framework
can tackle censored demand. Computationally, we demonstrate the algorithm on a
synthetic dataset with 6 products, 42 bundles, and 4000 transactions. This is
because Monte Carlo simulation allows us to provide close-form updates in the
M-step. No known algorithms can handle the estimation problem for bundle sales
of this scale.

  * •

We answer two theoretical questions: when is the model identifiable, and
whether the proposed algorithm converges to the true parameters. Both
questions have not been answered before for the estimation problem with bundle
sales. For the first question, we identify a simple pricing policy used by
many firms that leads to identifiable transaction data. For the second
question, we leverage the recent advances (Balakrishnan, Wainwright, and Yu
2017) and show that the algorithm is guaranteed to converge linearly to a
neighborhood of the true parameters (global maximum of the likelihood
function) under a number of technical conditions. The convergence to the
global maximum is only established recently and we are the first to show such
a result for the problem we study.

##  2 Literature Review

Bundle pricing mechanisms have been studied extensively in the literature,
including their theoretical and computational properties. We list a number of
representative papers below. Bakos and Brynjolfsson (1999), Abdallah,
Asadpour, and Reed (2021) analyze the asymptotic performance of pure bundling
and bundle size pricing. Chu, Leslie, and Sorensen (2011) conduct a
comprehensive numerical experiment to compare the performance of four commonly
studied bundle pricing mechanisms. Hanson and Martin (1990), Wu, Hitt, Chen,
and Anandalingam (2008), Jiang, Shang, Kemerer, and Liu (2011), Li, Sun, and
Teo (2020) provide mixed-integer programming or convex optimization
formulations to compute the optimal bundle pricing policies. Hart and Nisan
(2017), Babaioff, Immorlica, Lucier, and Weinberg (2020) derive theoretical
performance bounds for simple mechanisms. A few recent studies (Ma and Simchi-
Levi 2021, Sun, Li, and Teo 2021, Chen, Li, Li, and Wang 2022) propose new
implementable pricing mechanisms and analyze their properties. These papers
usually rely on one crucial assumption: the distribution of consumers’
valuations is known or a number of realized samples are given. In practice,
however, such information needs to be learned from the sales data. Our work
provides a framework to achieve this goal and further justifies the practical
feasibility of the pricing mechanisms studied in these papers. As pointed out
by Fisher, Olivares, and Staats (2020), the estimation of demand (customers’
valuations) is $90\%$ of the work in practice.

Our work is related to the stream of literature that focuses on the estimation
of discrete choice models, especially the well-known multinomial logit (MNL)
model, from transaction data. Kök et al. (2008) provide a comprehensive
literature review of the earlier studies on such models in revenue management.
Building on Talluri and Van Ryzin (2004), Vulcano et al. (2012) combine the
MNL model with a non-homogeneous Poisson arrival process over multiple periods
and use the EM algorithm to estimate model parameters as well as the number of
no-purchases. Newman et al. (2014) incorporate price and product features into
their estimation approach. Abdallah and Vulcano (2021) study the
identifiability of this type of model and propose a computational approach to
optimize the likelihood function. Subramanian and Harsha (2021) use a mixed-
integer program to estimate the parameters under the loss-minimization
objective function. The estimation of other types of discrete choice models is
also studied, including the rank-based choice model (Farias, Jagabathula, and
Shah 2013, Van Ryzin and Vulcano 2015) and the Markov chain choice model
(Şimşek and Topaloglu 2018). The major contribution of this paper to the
literature is that we consider customers may purchase bundles that a firm
offer at a discount. Such a setting doesn’t fit into the classic discrete
choice models, and this estimation problem has not been studied in the
literature.

This study proposes a framework to estimate the customers’ valuations from
transaction data for multiple products when the firm may offer a price menu
for products and bundles. There are three closely related papers (Jedidi,
Jagpal, and Manchanda 2003, Letham, Sun, and Sheopuri 2014, Ma and Simchi-Levi
2020) in this domain. Jedidi et al. (2003) study a similar problem and propose
a hierarchical Bayesian approach to infer the parameters of a multivariate
Gaussian distribution that is used to model the valuations for the products.
However, they propose to use off-the-shelf Bayesian computational tools that
seem to suffer from computational issues when the number of products is
moderate or large, according to our experiments. Indeed, their numerical
experiments are conducted for only two products and one bundle. Instead, we
provide an analytical framework to study this problem and use the EM algorithm
and Monte Carlo simulation to estimate the model parameters. We show that the
problem can be translated to an estimation problem when the samples are
“censored by regions”, i.e., the exact locations of the samples are not
observed but a polyhedral region in which each sample falls is given. In
another related paper, Letham et al. (2014) approximate the polyhedral region
with rectangular regions, which allows them to simplify the likelihood
function. In contrast, we use the EM algorithm to solve the original
estimation problem without approximations. Ma and Simchi-Levi (2020) consider
a slightly different setting: the bundle discount is only offered when a
customer purchases all products. They do not assume a parametric form for the
valuation distribution and focus on a few quantities of interest, in order to
reconstruct linear demand curves of each product. They show their estimators
are consistent. Our setting assumes the Gaussian distribution and allows for
an arbitrary price menu. In addition to the different setups and
methodologies, we extend our approach to the censored demand (unobserved no-
purchases). This is not considered in Jedidi et al. (2003), Letham et al.
(2014), Ma and Simchi-Levi (2020).

The methodology used in this paper, the EM algorithm combined with Monte Carlo
simulation, has been used in previous studies, including Natarajan et al.
(2000), Lee and Scott (2012). In particular, Monte Carlo simulation is used to
approximate the E-step, when the expectation doesn’t have a closed form.
Because our study is focused on a specific practical problem, we provide the
concrete steps of the EM algorithm and the Monte Carlo simulation. Moreover,
the theoretical properties are derived for the bundle application. In general,
the EM algorithm is only guaranteed to converge to the local maximum of the
sample likelihood function. Balakrishnan et al. (2017) provide conditions and
statistical guarantees for the EM algorithm to converge to the global maximum
of the likelihood function i.e., the true parameters, when initialized in a
neighborhood. We prove that the conditions proposed by Balakrishnan et al.
(2017) hold for our problem.

##  3 The Data Model and Estimation

In this section, we first introduce the format of the dataset recording the
bundle transactions. Suppose there are $I$ products and we use
$[I]\triangleq\left\\{1,\dots,I\right\\}$ to denote the product space. There
are $J$ candidate bundles offered on the menu and we use
$j\in\left\\{0,1\right\\}^{I}$ to denote a bundle. That is, product $i\in[I]$
is included in bundle $j$ if and only if $j_{i}=1$. Note that not all the
bundles are offered to each customer and a bundle may include only one
product. With a little abuse of notation, we also use $j$ as the bundle index
in the menu $[J]\triangleq\left\\{1,\dots,J\right\\}$ and $i\in j$ when bundle
$j$ includes product $i$.

The dataset records the purchase of $N$ customers: customer
$n\in[N]\triangleq\left\\{1,\dots,N\right\\}$ chooses bundle
$c_{n}\in[J]\cup\left\\{0\right\\}$, where the choice $0$ stands for no
purchases. The case of censored demand is investigated in Section 4.1.
Moreover, the price of bundle $j$ faced by customer $n$ is denoted by
$p^{j}_{n}$, which is also recorded in the dataset. Therefore, the dataset
consists of the tuples $\\{(p^{1}_{n},\dots,p^{J}_{n},c_{n})\\}_{n=1}^{N}$.
Note that we do not require the offering of all the bundles in the menu to all
the customers. When a customer is only offered a subset of the available
bundles, we can artificially set the prices of the bundles that are not
offered to infinity.

To be able to utilize the dataset and estimate customer preferences, we
specify a utility model. Suppose customer $t$ is endowed with a vector of
valuations, $\bm{v}_{n}=(v_{1},\dots,v_{I})$, for the products. Under additive
utility, the valuation for bundle $j$ is $\sum_{i\in j}v_{i}$. Note that
additive utility is a standard assumption in the literature, although other
types of valuation functions have been proposed and studied recently. For
additive utility, the complementary or substitutability of products can be
captured by the correlation of the valuations.

Under this utility model, the choice $c_{n}=j\neq 0$ corresponds to the
following region of the valuation vector $\bm{v}$, given by the incentive-
compatible constraints:

$c_{n}=j\iff\left\\{\sum_{i\in
j}v_{i}-p^{j}_{n}\geq\max_{j^{\prime}\in[J]}\sum_{i\in
j^{\prime}}v_{i}-p^{j^{\prime}}_{n},\;\sum_{i\in j}v_{i}-p^{j}_{n}\geq
0\right\\}.$ (1)

Similarly, we have

$c_{n}=0\iff\left\\{\max_{j^{\prime}\in[J]}\sum_{i\in
j^{\prime}}v_{i}-p^{j^{\prime}}_{n}\leq 0\right\\}.$ (2)

Suppose the valuation $\bm{v}$ of each customer is independently drawn from an
$I$-dimensional multivariate normal distribution:

$\bm{v}\sim\mathcal{N}(\bm{\mu},\bm{\Sigma}).$ (3)

Note that we focus on the somewhat restrictive normal distribution for the
purpose of exposition. In Section 4.2, we relax it to a Gaussian mixture
model, which is flexible enough to approximate any probability distributions
(for example, see (Goodfellow et al. 2016, p.65)).

With the presented data model, now we can state the goal of the firm. Given
the dataset $\\{(p^{1}_{n},\dots,p^{J}_{n},c_{n})\\}_{n=1}^{N}$, the firm
wants to estimate the distribution of the customer valuations, in particular
$\bm{\mu}$ and $\bm{\Sigma}$, with the knowledge that the incentive-compatible
constraints (1) and (2) relate the valuations to the observations in the
dataset. Note that precisely because of this indirect relationship, this
problem is not standard in the statistics literature. For example, if
$\bm{v}_{n}$ were observed, then the estimation of $\bm{\mu}$ and
$\bm{\Sigma}$ would be trivial using the empirical average and covariance
matrix of $\\{\bm{v}_{n}\\}_{n=1}^{N}$. Next we explain how the problem can be
viewed as an estimation problem with missing data and present an algorithm for
the problem.

###  3.1 An Estimation Problem with Missing Data

The key observation to simplify the estimation problem is that given
$(p^{1}_{n},\dots,p^{J}_{n},c_{n})$, the incentive-compatible constraints (1)
and (2) specify a polyhedron in which the valuation vector $\bm{v}_{n}$ falls
in. More precisely, define

$\displaystyle R^{j}_{n}$ $\displaystyle\triangleq\left\\{\bm{v}:\sum_{i\in
j}v_{i}-p^{j}_{n}\geq\max_{j^{\prime}\in[J]}\sum_{i\in
j^{\prime}}v_{i}-p^{j^{\prime}}_{n},\;\sum_{i\in j}v_{i}-p^{j}_{n}\geq
0\right\\},\quad j\neq 0,$ $\displaystyle R^{0}_{n}$
$\displaystyle\triangleq\left\\{\bm{v}:\max_{j^{\prime}\in[J]}\sum_{i\in
j^{\prime}}v_{i}-p^{j^{\prime}}_{n}\leq 0\right\\},$ (4)

which we refer to as _IC polyhedra_. The observation $c_{n}=j$ is equivalent
to $\bm{v}_{n}\in R^{j}_{n}$. As a result, for customer $n$, the valuation
space $\mathbb{R}^{I}$ can be partitioned into the $J+1$ regions:
$\mathbb{R}^{I}=\cup_{j=1}^{J}R^{j}_{n}\cup R^{0}_{n}$. The observation
$(p^{1}_{n},\dots,p^{J}_{n},c_{n})$ characterizes the partition as well as the
IC polyhedron in which $\bm{v}_{n}$ falls. Therefore, instead of observing
$\bm{v}_{n}$ exactly, the firm observes a polyhedron $R_{n}^{c_{n}}$ for all
potential $\bm{v}_{n}$ that is consistent with the observation. This is
sometimes referred to as “censored” observations in statistics. To avoid
ambiguity, we do not use the term and reserve it for demand censoring, which
is the subject of Section 4.1. With this interpretation, we view the
estimation problem through the lens of inference with missing data.

Parameters:

we are interested in estimating $(\bm{\mu},\bm{\Sigma})$, the distribution of
the valuations for the products among consumers. By convention, we use
$\bm{\theta}\triangleq(\bm{\mu},\bm{\Sigma})$ to denote the parameters.

Observed data:

the firm observes the polyhedra $\left\\{R_{n}^{c_{n}}\right\\}_{n=1}^{N}$
given the data $\\{(p^{1}_{n},\dots,p^{J}_{n},c_{n})\\}_{n=1}^{N}$ in (3.1).
We use $\bm{X}\triangleq\left\\{R_{n}^{c_{n}}\right\\}_{n=1}^{N}$ for the
observed data. Note that by converting the data to polyhedra, we do not lose
any information regarding the possible valuations
$\left\\{\bm{v}_{n}\right\\}_{n=1}^{N}$.

Missing data:

we treat valuations $\bm{v}_{n}$, $n=1,\dots,N$, as the missing data. They are
related to the observed data by the simple fact that $\bm{v}_{n}\in
R_{n}^{c_{n}}$. We use $\bm{Z}\triangleq\\{\bm{v}_{n}\\}_{n=1}^{N}$ for the
missing data.

There are two common approaches to missing data: imputation and maximum
likelihood estimator (MLE). Since our estimation problem is model-based, we
adopt the latter. In particular, the likelihood function based on the observed
data is:

$L(\bm{\theta};\bm{X})=\int
p(\bm{X},\bm{Z}|\bm{\theta})\mathrm{d}\bm{Z}=\prod_{n=1}^{N}\int_{R^{c_{t}}_{t}}f(\bm{v}|\bm{\mu},\bm{\Sigma})\mathrm{d}\bm{v},$

where

$f(\bm{v}|\bm{\mu},\bm{\Sigma})=\frac{1}{\sqrt{(2\pi)^{I}|\bm{\Sigma}|}}\exp\left(-\frac{1}{2}(\bm{v}-\bm{\mu})^{\top}\bm{\Sigma}^{-1}(\bm{v}-\bm{\mu})\right)$

is the probability density function (PDF) of the multivariate normal
distribution (3).

To maximize the likelihood function over $\bm{\theta}$, a standard technique
is to apply logarithm to the likelihood function and maximize the log-
likelihood function instead. However, because of the integration, the
optimization is theoretically intractable and computationally challenging. As
a popular alternative, the EM algorithm iteratively maximizes the likelihood
function by focusing on the _complete-data_ likelihood. Next we elaborate on
the implementation of the EM algorithm in this problem.

###  3.2 The EM Algorithm

Instead of focusing on the (log-)likelihood function of the observed data with
missing values, the EM algorithm focuses on the complete-data log-likelihood
function, i.e.,

$l(\bm{\theta};\bm{X},\bm{Z})=\sum_{n=1}^{N}\log
f(\bm{v}_{n}|\bm{\mu},\bm{\Sigma}).$ (5)

Note that as an objective function, $l(\bm{\theta};\bm{X},\bm{Z})$ is a lot
easier to handle than $l(\bm{\theta};\bm{X})$: because of the form of $f$
(belonging to the exponential family), the complete-data log-likelihood is a
quadratic function, which is easy to optimize. However, $\bm{Z}$ is not
observed from the data.

The EM algorithm circumvents the issue by estimating $\bm{\theta}^{(t)}$
iteratively. In each iteration, suppose
$\bm{\theta}^{(t)}=\left(\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)$ is given. We
first evaluate the expected complete-data log-likelihood over $\bm{Z}$, i.e.,

$\displaystyle Q\left(\bm{\theta}|\bm{\theta}^{(t)}\right)$
$\displaystyle\triangleq\mathbb{E}\left[l(\bm{\theta};\bm{X},\bm{Z})\bigg{|}\bm{X},\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right]$
$\displaystyle=\sum_{n=1}^{N}\int_{\bm{v}\in
R_{n}^{c_{n}}}\frac{f\left(\bm{v}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)}{\int_{R_{n}^{c_{n}}}f\left(\bm{x}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)\mathrm{d}\bm{x}}\log
f(\bm{v}|\bm{\mu},\bm{\Sigma})\mathrm{d}\bm{v}.$ (6)

This is the so-called _E-step_ and the expectation is taken with respect to
the missing data $\bm{Z}$. Given
$\left(\bm{X},\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)$, the missing variable
$\bm{v}_{n}$ has a conditional Gaussian distribution inside the polyhedron
$R_{n}^{c_{n}}$ with the given mean and covariance. Its PDF is thus
$f\left(\bm{v}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)/\int_{R_{n}^{c_{n}}}f\left(\bm{x}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)\mathrm{d}\bm{x}$
for $\bm{v}\in R_{n}^{c_{n}}$, which is the first term in the integral. The
second term $\log f(\bm{v}|\bm{\mu},\bm{\Sigma})$ is the log-likelihood
function $l(\bm{\theta};\bm{X},\bm{Z})$ and is integrated with the conditional
PDF for its expected value.

To find a new estimator for $\bm{\theta}$ in the $t$-th iteration, the _M-
step_ in the EM algorithm maximizes
$Q\left(\bm{\theta}|\bm{\theta}^{(t)}\right)$ and let
$\bm{\theta}^{(t+1)}=\text{argmax}_{\bm{\theta}}Q\left(\bm{\theta}|\bm{\theta}^{(t)}\right)$.
However, the major difficulty is to deal with the integral and the denominator
which don’t have closed forms. To deal with this challenge, we use Monte Carlo
simulation. More precisely, for all $n=1,\dots,N$, we simulate $L$ samples
according to the same distribution as
$\bm{v}|R_{n}^{c_{n}},\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}$. This is the
conditional multivariate normal distribution inside $R_{n}^{c_{n}}$. The
standard acceptance-rejection method can be used and we provide the details in
Algorithm 1.

Algorithm 1 Acceptance-rejection method

1:Input: $R_{n}^{c_{n}}$, $\bm{\mu}^{(t)}$, $\bm{\Sigma}^{(t)}$, $L$

2:$l\leftarrow 1$

3:while $l\leq L$ do

4: Repeat $\bm{x}\sim\mathcal{N}(\bm{\mu}^{(t)},\bm{\Sigma}^{(t)})$ until
$\bm{x}\in R_{n}^{c_{n}}$

5: $\bm{x}_{n}^{(l)}\leftarrow\bm{x}$

6: $l\leftarrow l+1$

7:end while

8:return $\left\\{\bm{x}_{n}^{(1)},\dots,\bm{x}_{n}^{(L)}\right\\}$

Once the Monte Carlo samples $\\{\bm{x}_{n}^{(l)}\\}$ for $n=1,\dots,N$ and
$l=1,\dots,L$ have been generated, we may replace
$Q(\bm{\theta}|\bm{\theta}^{(t)})$ in (3.2) by

$\displaystyle\hat{Q}(\bm{\theta}|\bm{\theta}^{(t)})$
$\displaystyle=\frac{1}{L}\sum_{n=1}^{N}\sum_{l=1}^{L}\log
f(\bm{x}_{n}^{(l)}|\bm{\mu},\bm{\Sigma})$
$\displaystyle=C+\frac{1}{L}\sum_{n=1}^{N}\sum_{l=1}^{L}\left(-\frac{1}{2}\left(\bm{x}_{n}^{(l)}-\bm{\mu}\right)^{\top}\bm{\Sigma}^{-1}\left(\bm{x}_{n}^{(l)}-\bm{\mu}\right)-\frac{1}{2}\log|\bm{\Sigma}|\right),$
(7)

where $C$ is constant with respect to $\bm{\mu}$ and $\bm{\Sigma}$. Note that
this is the standard MLE for multivariate normal distribution with $NL$
samples. Therefore, we have

$\displaystyle\bm{\mu}^{(t+1)}=$
$\displaystyle\frac{1}{NL}\sum_{n=1}^{N}\sum_{l=1}^{L}\bm{x}_{n}^{(l)},$
$\displaystyle\bm{\Sigma}^{(t+1)}=$
$\displaystyle\frac{1}{NL}\sum_{n=1}^{N}\sum_{l=1}^{L}\left(\bm{x}_{n}^{(l)}-\bm{\mu}^{(t+1)}\right)\left(\bm{x}_{n}^{(l)}-\bm{\mu}^{(t+1)}\right)^{\top}.$

This completes the $t$-th iteration of the EM algorithm.

The EM algorithm is terminated when the estimation $\bm{\theta}^{(t)}$ has
converged. In practice, we may impose a small tolerance level $\epsilon>0$ and
terminate the algorithm once
$\|\bm{\theta}^{(t+1)}-\bm{\theta}^{(t)}\|\leq\epsilon$ for a chosen norm
$\|\cdot\|$. We summarize the steps in Algorithm 2.

Algorithm 2 The EM algorithm

1:Input: $\\{(p^{1}_{n},\dots,p^{J}_{n},c_{n})\\}_{n=1}^{N}$, $\epsilon$, $L$,
$\bm{\mu}^{(0)}$, $\bm{\Sigma}^{(0)}$

2:Calculate $\left\\{R_{n}^{c_{n}}\right\\}_{n=1}^{N}$ from the data according
to (3.1)

3:$error\leftarrow 100$, $t\leftarrow 0$

4:while $error>\epsilon$ do

5: Run Algorithm 1 for $\left\\{R_{n}^{c_{n}}\right\\}_{n=1}^{N}$

6: Update  $\displaystyle\bm{\mu}^{(t+1)}\leftarrow$
$\displaystyle\frac{1}{NL}\sum_{n=1}^{N}\sum_{l=1}^{L}\bm{x}_{n}^{(l)}$
$\displaystyle\bm{\Sigma}^{(t+1)}\leftarrow$
$\displaystyle\frac{1}{NL}\sum_{n=1}^{N}\sum_{l=1}^{L}\left(\bm{x}_{n}^{(l)}-\bm{\mu}^{(t+1)}\right)\left(\bm{x}_{n}^{(l)}-\bm{\mu}^{(t+1)}\right)^{\top}$

7:
$error\leftarrow\|\bm{\mu}^{(t+1)}-\bm{\mu}^{(t)}\|_{1}+\|\bm{\Sigma}^{(t+1)}-\bm{\Sigma}^{(t)}\|_{1}$,
$t\leftarrow t+1$

8:end while

9:return $\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}$

Sampling efficiency of Monte Carlo simulation. In practice, the major
computational complexity of Algorithm 2 is caused by the Monte Carlo
simulation. In particular, when the region $R_{n}^{c_{n}}$ is small or distant
from $\bm{\mu}^{(t)}$, Algorithm 1 may be inefficient as it requires a large
number of samples in Step 4 to draw one accepted in $R_{n}^{c_{n}}$. We
provide three remedies to improve the sampling efficiency.

First, in each iteration of Algorithm 2 (Step 4), we can generate a large
number of samples from
$\mathcal{N}\left(\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)$ and store them.
When running Algorithm 1 for each $R_{n}^{c_{n}}$, the algorithm can use the
same set of samples for acceptance/rejection. This saves the computation time
to generate a large number of samples for each $n=1,\dots,N$.

Second, we may not require the same number of samples $L$ in all regions
$R_{n}^{c_{n}}$. For example, we may simulate $L_{n}$ Monte Carlo samples for
region $R_{n}^{c_{n}}$ so that $L_{n}$ could be smaller for low-probability
$R_{n}^{c_{n}}$. As long as we use $\sum_{l=1}^{L_{n}}\log
f(\bm{x}_{n}^{(l)}|\bm{\mu},\bm{\Sigma})/L_{n}$ in
$\hat{Q}(\bm{\theta}|\bm{\theta}^{(t)})$, it is still unbiased. However, one
needs to caution that the Monte Carlo error is bottlenecked by the smallest
$L_{n}$.

Third, for low-probability $R_{n}^{c_{n}}$, we may use _importance sampling_
to increase the sampling efficiency. In particular, consider (3.2) and some
$n$ that $f\left(\bm{v}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)\approx 0$ for
$\bm{v}\in R_{n}^{c_{n}}$. In this case, it is computationally challenging as
Step 4 is repeated many times before acceptance. Instead, we can consider a
different normal distribution
$\mathcal{N}(\bm{\mu}^{\prime},\bm{\Sigma}^{\prime})$ such that $\bm{v}\in
R_{n}^{c_{n}}$ with high probability. This can be achieved, for example, by
choosing $\bm{\mu}^{\prime}$ close to the center of $R_{n}^{c_{n}}$ when
$R_{n}^{c_{n}}$ is bounded. As a result, we have

$\displaystyle\int_{\bm{v}\in
R_{n}^{c_{n}}}\frac{f\left(\bm{v}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)}{\int_{R_{n}^{c_{n}}}f\left(\bm{x}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)\mathrm{d}\bm{x}}\log
f(\bm{v}|\bm{\mu},\bm{\Sigma})\mathrm{d}\bm{v}$ $\displaystyle=$
$\displaystyle\int_{\bm{v}\in
R_{n}^{c_{n}}}\frac{f(\bm{v}|\bm{\mu}^{\prime},\bm{\Sigma}^{\prime})}{\int_{R_{n}^{c_{n}}}f(\bm{x}|\bm{\mu}^{\prime},\bm{\Sigma}^{\prime})\mathrm{d}\bm{x}}\frac{f\left(\bm{v}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)\int_{R_{n}^{c_{n}}}f(\bm{x}|\bm{\mu}^{\prime},\bm{\Sigma}^{\prime})\mathrm{d}\bm{x}}{f(\bm{v}|\bm{\mu}^{\prime},\bm{\Sigma}^{\prime})\int_{R_{n}^{c_{n}}}f\left(\bm{x}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)\mathrm{d}\bm{x}}\log
f(\bm{v}|\bm{\mu},\bm{\Sigma})\mathrm{d}\bm{v}.$

By simulating $L$ samples from a truncated normal distribution
$\mathcal{N}(\bm{\mu}^{\prime},\bm{\Sigma}^{\prime})$ on region
$R_{n}^{c_{n}}$ (run Algorithm 1 with the new distribution), denoted as
$\bm{x}_{n}^{(1)},\dots,\bm{x}_{n}^{(L)}$, we can approximate the above term
by

$\displaystyle\frac{1}{\frac{\int_{R_{n}^{c_{n}}}f\left(\bm{x}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)\mathrm{d}\bm{x}}{\int_{R_{n}^{c_{n}}}f\left(\bm{x}|\bm{\mu}^{\prime},\bm{\Sigma}^{\prime}\right)\mathrm{d}\bm{x}}}$
$\displaystyle\frac{1}{L}\sum_{l=1}^{L}\frac{f\left(\bm{x}_{n}^{(l)}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)}{f\left(\bm{x}_{n}^{(l)}|\bm{\mu}^{\prime},\bm{\Sigma}^{\prime}\right)}\log
f\left(\bm{x}_{n}^{(l)}|\bm{\mu},\bm{\Sigma}\right)$
$\displaystyle=\frac{1}{\sum_{l=1}^{L}\frac{f\left(\bm{x}_{n}^{(l)}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)}{f\left(\bm{x}_{n}^{(l)}|\bm{\mu}^{\prime},\bm{\Sigma}^{\prime}\right)}}\sum_{l=1}^{L}\frac{f\left(\bm{x}_{n}^{(l)}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)}{f\left(\bm{x}_{n}^{(l)}|\bm{\mu}^{\prime},\bm{\Sigma}^{\prime}\right)}\log
f\left(\bm{x}_{n}^{(l)}|\bm{\mu},\bm{\Sigma}\right),$ (8)

where the equality follows from approximating the denominator outside the sum:

$\frac{\int_{R_{n}^{c_{n}}}f\left(\bm{x}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)\mathrm{d}\bm{x}}{\int_{R_{n}^{c_{n}}}f(\bm{x}|\bm{\mu}^{\prime},\bm{\Sigma}^{\prime})\mathrm{d}\bm{x}}=\int_{R_{n}^{c_{n}}}\frac{f\left(\bm{x}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)}{f(\bm{x}|\bm{\mu}^{\prime},\bm{\Sigma}^{\prime})}\frac{f(\bm{x}|\bm{\mu}^{\prime},\bm{\Sigma}^{\prime})}{\int_{R_{n}^{c_{n}}}f(\bm{x}|\bm{\mu}^{\prime},\bm{\Sigma}^{\prime})\mathrm{d}\bm{x}}\mathrm{d}\bm{x}=\frac{1}{L}\sum_{l=1}^{L}\frac{f\left(\bm{x}_{n}^{(l)}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)}{f\left(\bm{x}_{n}^{(l)}|\bm{\mu}^{\prime},\bm{\Sigma}^{\prime}\right)}.$

We Define
$w_{nl}\triangleq\frac{f\left(\bm{x}_{n}^{(l)}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)}{f\left(\bm{x}_{n}^{(l)}|\bm{\mu}^{\prime},\bm{\Sigma}^{\prime}\right)}/\sum_{l=1}^{L}\frac{f\left(\bm{x}_{n}^{(l)}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)}{f\left(\bm{x}_{n}^{(l)}|\bm{\mu}^{\prime},\bm{\Sigma}^{\prime}\right)}$.
Using the Monte Carlo samples, and (8), the approximate $Q$ function (3.2) can
be cast in the following form

$\hat{Q}\left(\bm{\theta}|\bm{\theta}^{(t)}\right)=\sum_{n=1}^{N}\sum_{l=1}^{L}w_{nl}\log
f\left(\bm{x}_{n}^{(l)}|\bm{\mu},\bm{\Sigma}\right).$

The $M$-step thus leads to the following update

$\displaystyle\bm{\mu}^{(t+1)}=$
$\displaystyle\frac{\sum_{n=1}^{N}\sum_{l=1}^{L}w_{nl}\bm{x}_{n}^{(l)}}{\sum_{n=1}^{N}\sum_{l=1}^{L}w_{nl}},$
$\displaystyle\bm{\Sigma}^{(t+1)}=$
$\displaystyle\frac{\sum_{n=1}^{N}\sum_{l=1}^{L}w_{nl}\left(\bm{x}_{n}^{(l)}-\bm{\mu}^{(t+1)}\right)\left(\bm{x}_{n}^{(l)}-\bm{\mu}^{(t+1)}\right)^{\top}}{\sum_{n=1}^{N}\sum_{l=1}^{L}w_{nl}}.$

Therefore, using importance sampling results in a similar procedure to
Algorithm 2 but reduces the computation in Algorithm 1.

##  4 Model Extensions

In this section, we consider two extensions of the base model in Section 3.

###  4.1 Censored Demand

In practice, the transactions observed by the firm are typically censored: the
firm only records the transactions in which customers purchase at least one
item and the customers who have not purchased anything are censored. Next we
tackle this problem using our framework.

For ease of exposition, in this subsection, we consider that the firm propose
only one price menu $\bm{p}=(p^{1},\dots,p^{J})$ seen by all customers.
Therefore, we assume that the firm observes $(p^{1},\dots,p^{J},c_{n})$ for
customer $n$, where the choice $c_{n}>0$ for $n=1,\dots,N$. This data can be
converted to the IC polyhedron defined in (3.1),
$\\{R^{0},R^{1},\dots,R^{J}\\}$. Note that they are constant for all customers
because of the fixed price menu. Similarly we treat $\bm{v}_{n}$,
$n=1,\dots,N$, the valuation of customer $n$, as the missing data.

Because of the demand censoring, the total number of customers that have
visited the store, $N^{\prime}\geq N$, is also unobserved, as well as the
valuations of the censored customers $\bm{v}_{n}$, $n=N+1,\dots,N^{\prime}$.
Therefore, the missing data is
$\bm{Z}\triangleq\left\\{N^{\prime},\bm{v}_{1},\dots,\bm{v}_{N^{\prime}}\right\\}$.
It is known that $\bm{v}_{n}\in R^{0}$ for the censored customers
$n=N+1,\dots,N^{\prime}$. Moreover, we record the total number of customers
that buy bundle $j$ as $N_{j}$, where $\sum_{j=1}^{J}N_{j}=N$. The observed
data can thus be encoded as
$\bm{X}\triangleq\left\\{N_{1},\dots,N_{J}\right\\}$. Below we describe the
implementation of the EM algorithm to estimate
$\bm{\theta}=(\bm{\mu},\bm{\Sigma})$.

#### E-step.

We first express the complete-data log-likelihood as

$l(\bm{\theta};\bm{X},\bm{Z})=\sum_{n=1}^{N^{\prime}}\log
f(\bm{v}_{n}|\bm{\mu},\bm{\Sigma}).$ (9)

We then take the expectation of $l(\bm{\theta};\bm{X},\bm{Z})$ conditional on
$\bm{X}$ under the current estimate $\bm{\theta}^{(t)}$:

$\displaystyle Q(\bm{\theta}|\bm{\theta}^{(t)})$
$\displaystyle\triangleq\mathbb{E}\left[l(\bm{\theta};\bm{X},\bm{Z})\bigg{|}\bm{X},\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right]$
$\displaystyle=\mathbb{E}\left[\mathbb{E}\left[\sum_{n=1}^{N^{\prime}}\log
f(\bm{v}_{n}|\bm{\mu},\bm{\Sigma})\bigg{|}\bm{X},\bm{\theta}^{(t)},N^{\prime}\right]\bigg{|}\bm{X},\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right]$
$\displaystyle=\mathbb{E}\left[\sum_{n=1}^{N^{\prime}}\int_{\bm{v}\in
R_{n}^{c_{n}}}\frac{f\left(\bm{v}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)}{\int_{\bm{x}\in
R_{n}^{c_{n}}}f\left(\bm{x}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)\mathrm{d}\bm{x}}\log
f(\bm{v}|\bm{\mu},\bm{\Sigma})\mathrm{d}\bm{x}\bigg{|}\bm{X},\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right].$
(10)

In the second equation, we first take expectation of the valuations
$\bm{v}_{n}$ conditional on the total number of customers $N^{\prime}$ using
the tower property. This follows the same derivation as the base model in
Section 3.

Next we analyze the distribution $N^{\prime}\mid\bm{X},\bm{\theta}^{(t)}$.
Note that under $\bm{\theta}^{(t)}$, a customer is censored with probability
$\mathbb{P}(R^{0}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)})$, which represents the
probability of the IC polyhedron $R^{0}$. Therefore, the number of total
customers $N^{\prime}\mid\bm{X},\bm{\theta}^{(t)}$ can be described by the
following probability model: after $N^{\prime}$ i.i.d. Bernoulli trials with
success probability $1-\mathbb{P}(R^{0}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)})$ (an
uncensored customer is counted as a success), there are $N$ successful trials
being observed. Therefore, for $n=N,N+1,\dots$, we have

$\displaystyle\mathbb{P}\left(N^{\prime}=n|\bm{X},\bm{\theta}^{(t)}\right)=\mathbb{P}\left(N^{\prime}=n|N,\bm{\theta}^{(t)}\right)$
$\displaystyle=\frac{\mathbb{P}\left(N=m|N^{\prime}=n,\bm{\theta}^{(t)}\right)\mathbb{P}\left(N^{\prime}=n|\bm{\theta}^{(t)}\right)}{\mathbb{P}\left(N=m|\bm{\theta}^{(t)}\right)}$
$\displaystyle=\frac{\mathbb{P}\left(N|N^{\prime}=n,\bm{\theta}^{(t)}\right)\mathbb{P}(N^{\prime}=n)}{\sum_{\tilde{n}=N}^{\infty}\mathbb{P}(N|N^{\prime}=\tilde{n},\bm{\theta}^{(t)})\mathbb{P}(N^{\prime}=\tilde{n})}.$

This probability is unspecified unless we have a prior distribution for
$N^{\prime}$. For this purpose, we use the improper prior following marquis de
Laplace (1840) suggesting that one should apply a uniform prior to unknown
events from “the principle of insufficient reason”, i.e.,
$\mathbb{P}(N^{\prime}=n_{1})=\mathbb{P}(N^{\prime}=n_{2})$ for $n_{1}\neq
n_{2}$. Therefore, the above term can be expressed as

$\displaystyle\frac{\binom{n}{N}\mathbb{P}\left(R^{0}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)^{n-N}\left(1-\mathbb{P}\left(R^{0}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)\right)^{N}}{\sum_{\tilde{n}=N}^{+\infty}\binom{\tilde{n}}{N}\mathbb{P}\left(R^{0}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)^{\tilde{n}-N}(1-\mathbb{P}(R^{0}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}))^{N}}$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad=\binom{n}{N}\mathbb{P}\left(R^{0}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)^{n-N}\left(1-\mathbb{P}\left(R^{0}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)\right)^{N+1}.$
(11)

Note that it has the same distribution as the negative binomial distribution
with success rate
$1-\mathbb{P}\left(R^{0}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)$ and $N+1$
successes. Recall that the negative binomial distribution with $N+1$ successes
describes the probability mass function of the number of trials in independent
Bernoulli trials until the $(N+1)$th success is observed. The use of the
improper prior equates the total number of customers to that right before the
$(N+1)$th customer who makes a purchase.

With this formula, we next proceed to (4.1). Similarly we resort to Monte
Carlo simulation in the following steps:

  1. 1.

Consider $L$ instances. In instance $l$, simulate $N^{\prime(l)}$ according to
the negative binomial distribution (11).

  2. 2.

Then simulate $N_{j}$ samples from
$\mathcal{N}(\bm{\mu}^{(t)},\bm{\Sigma}^{(t)})$ conditional on $\bm{x}\in
R^{j}$, denoted $\\{\bm{x}^{(l)}_{j,s}\\}_{s=1}^{N_{j}}$, and
$N^{\prime(l)}-N$ samples in $R_{0}$, denoted
$\\{\bm{x}^{(l)}_{0,s}\\}_{s=1}^{N^{\prime(l)}-N}$.

  3. 3.

Approximate the expectation in (4.1) by

$\hat{Q}\left(\bm{\theta}|\bm{\theta}^{(t)}\right)=\frac{1}{L}\sum_{l=1}^{L}\left(\sum_{j=1}^{J}\sum_{s=1}^{N_{j}}\log
f\left(\bm{x}_{j,s}^{(l)}|\bm{\mu},\bm{\Sigma}\right)+\sum_{s=1}^{N^{\prime(l)}-N}\log
f\left(\bm{x}_{0,s}^{(l)}\big{|}\bm{\mu},\bm{\Sigma}\right)\right).$ (12)

Note that the information of $\bm{\theta}^{(t)}$ has been fully absorbed by
the simulated samples $N^{{}^{\prime}(l)}$ and $\\{\bm{x}_{j,s}^{(l)}\\}$. The
formula in (12) is a function of $(\bm{\mu},\bm{\Sigma})$, which will be
optimized in the M-step.

Given $\hat{Q}(\bm{\theta}|\bm{\theta}^{(t)})$ in (12), the M-step is similar
to (3.2). In particular, $\bm{\mu}^{(t+1)}$ and $\bm{\Sigma}^{(t+1)}$ are the
sample mean and variance in the standard MLE for multivariate normal
distributions:

$\displaystyle\bm{\mu}^{(t+1)}=$
$\displaystyle\frac{1}{\sum_{l=1}^{L}N^{\prime(l)}}\sum_{l=1}^{L}\left(\sum_{j=1}^{J}\sum_{s=1}^{N_{j}}\bm{x}_{j,s}^{(l)}+\sum_{s=1}^{N^{\prime(l)}-N}\bm{x}_{0,s}^{(l)}\right),$
$\displaystyle\bm{\Sigma}^{(t+1)}=$
$\displaystyle\frac{1}{\sum_{l=1}^{L}N^{\prime(l)}}\sum_{l=1}^{L}\Bigg{(}\sum_{j=1}^{J}\sum_{s=1}^{N_{j}}\left(\bm{x}_{j,s}^{(l)}-\bm{\mu}^{(t+1)}\right)\left(\bm{x}_{j,s}^{(l)}-\bm{\mu}^{(t+1)}\right)^{\top}$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad+\sum_{s=1}^{N^{\prime(l)}-N}\left(\bm{x}_{0,s}^{(l)}-\bm{\mu}^{(t+1)}\right)\left(\bm{x}_{0,s}^{(l)}-\bm{\mu}^{(t+1)}\right)^{\top}\Bigg{)}.$

Then Algorithm 2 can be adapted to the censored-demand data.

To conclude this section, we discuss how the EM algorithm can be generalized
when customers may observe different price menus provided by the firm, for
example, based on the arrival time in a day. Suppose that the firm offers $m$
price menus to customers. The transaction data consists of $N^{m}$ observed
customer purchases from each price menu and their choices. Note that the IC
polyhedra differ across the price menus. Using the same procedure stated
above, we can use Monte Carlo to simulate the number of total customers
$N^{\prime m}$ _for each price menu_ and simulate the valuations of the
customers accordingly. They can be used in the M-step to update the
parameters.

###  4.2 Gaussian Mixture Models

In this section, we generalize the base model in Section 3 by considering
valuation $\bm{v}$ of each customer independently drawn from an
$I$-dimensional Gaussian mixture model with $K$ components:

$\bm{v}\sim\sum_{k=1}^{K}\phi_{k}\mathcal{N}(\bm{\mu}_{k},\bm{\Sigma}_{k}),$

where $\phi_{k}\geq 0$ and $\sum_{k=1}^{K}\phi_{k}=1$ are the mixture
proportions, $\mathcal{N}(\bm{\mu}_{k},\bm{\Sigma}_{k})$ is the density of a
multivariate Gaussian distribution with mean $\bm{\mu}_{k}$ and covariance
matrix $\bm{\Sigma}_{k}$, and
$\bm{\theta}=(\phi_{1},\cdots,\phi_{k},\bm{\mu}_{1},\cdots,\bm{\mu}_{k},\bm{\Sigma}_{1},\cdots,\bm{\Sigma}_{k})$
denotes the vector of parameters that we need to estimate. Such a Gaussian
mixture model allows for clustering of consumer valuations and can flexibly
approximates any distribution.

We use $\bm{\pi}^{n}\in\\{0,1\\}^{K}$ to denote the component association of
customer $n$: $\pi_{k}^{n}=1$ if $\bm{v}_{n}$ is generated from
$\mathcal{N}(\bm{\mu}_{k},\bm{\Sigma}_{k})$ and zero otherwise. Note that in
this model, for each transaction there are two missing values: (i) the
unobserved valuations, and (ii) the associated component to this valuation.
Therefore, the missing data is
$\bm{Z}\triangleq\\{\bm{\pi}^{1},\dots,\bm{\pi}^{N},\bm{v}_{1},\dots,\bm{v}_{N}\\}$.
We express the complete-data log-likelihood function as

$l(\bm{\theta};\bm{X},\bm{Z})=\sum_{n=1}^{N}\sum_{k=1}^{K}\bm{\pi}_{k}^{n}[\log(\bm{\phi}_{k})+\log
f(\bm{v}_{n}|\bm{\mu}_{k},\bm{\Sigma}_{k})],$

where recall that $f(\cdot)$ is the PDF of the multivariate normal
distribution.

#### E-step.

Since the $\bm{\pi}_{k}^{n}$ and $\bm{v}_{n}$ are not observed, we consider
the expected complete-data log-likelihood. We first compute
$Q\left(\bm{\theta}|\bm{\theta}^{(t)}\right)$ by computing the expectation of
$l(\bm{\theta};\bm{X},\bm{Z})$ conditional on $\bm{X}$ under the current
estimate $\bm{\theta}^{(t)}$, and thus, we have

$\displaystyle Q\left(\bm{\theta}|\bm{\theta}^{(t)}\right)\triangleq$
$\displaystyle\mathbb{E}_{\bm{Z}|\bm{X},\bm{\theta}^{(t)}}\left[l(\bm{\theta};\bm{X},\bm{Z})\bigg{|}\bm{X},\bm{\mu}^{(t)},\bm{\Sigma}^{(t)},\bm{\phi}^{(t)}\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}_{\bm{\pi},\bm{v}|\bm{X},\bm{\theta}^{(t)}}\left[\sum_{n=1}^{N}\sum_{k=1}^{K}\pi_{k}^{n}\log(\phi_{k})+\pi_{k}^{n}\log
f(\bm{v}_{n}|\bm{\mu}_{k},\bm{\Sigma}_{k})\bigg{|}\bm{X},\bm{\mu}^{(t)},\bm{\Sigma}^{(t)},\bm{\phi}^{(t)}\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}_{\bm{\pi}|\bm{X},\bm{\theta}^{(t)}}\left[\mathbb{E}_{\bm{v}|\bm{X},\bm{\theta}^{(t)},\bm{\pi}}\left[\sum_{n=1}^{N}\sum_{k=1}^{K}\pi_{k}^{n}\log(\phi_{k})+\pi_{k}^{n}\log
f(\bm{v}_{n}|\bm{\mu}_{k},\bm{\Sigma}_{k})\bigg{|}\bm{X},\bm{\mu}^{(t)},\bm{\Sigma}^{(t)},\bm{\phi}^{(t)}\right]\right]$
$\displaystyle=$
$\displaystyle\sum_{n=1}^{N}\sum_{k=1}^{K}\mathbb{P}\left(\pi_{k}^{n}=1\big{|}\bm{X},\bm{\theta}^{(t)}\right)\mathbb{E}\left[\log(\phi_{k})+\log
f(\bm{v}_{n}|\bm{\mu}_{k},\bm{\Sigma}_{k})\bigg{|}\bm{X},\bm{\theta}^{(t)},\pi_{k}^{n}=1\right]$
(13)

In the last line of (13), we further have

$\mathbb{E}\left[\log
f(\bm{v}_{n}|\bm{\mu}_{k},\bm{\Sigma}_{k})\bigg{|}\bm{X},\bm{\theta}^{(t)},\pi_{k}^{n}=1\right]=\int_{\bm{v}\in
R_{n}^{c_{n}}}\frac{f\left(\bm{v}|\bm{\mu}^{(t)}_{k},\bm{\Sigma}^{(t)}_{k}\right)}{\int_{R_{n}^{c_{n}}}f\left(\bm{x}|\bm{\mu}^{(t)}_{k},\bm{\Sigma}^{(t)}_{k}\right)\mathrm{d}\bm{x}}\log
f(\bm{v}|\bm{\mu}_{k},\bm{\Sigma}_{k})\mathrm{d}\bm{v}.$

To calculate (13), we first calculate the probability
$\mathbb{P}\left(\pi_{k}^{n}=1|\bm{X},\bm{\mu}^{(t)}_{k},\bm{\Sigma}^{(t)}_{k},\bm{\phi}_{k}^{(t)}\right)$,
i.e., given the choice $c_{n}$ of customer $n$, what is the probability that
s/he is associated with component $k$:

$\hat{\pi}^{n}_{k}\triangleq\mathbb{P}\left(\pi_{k}^{n}=1\big{|}\bm{X},\bm{\theta}^{(t)}\right)=\frac{\mathbb{P}\left(R_{n}^{c_{n}}\big{|}\bm{\mu}^{(t)}_{k},\bm{\Sigma}^{(t)}_{k},\pi_{k}^{n}=1\right)\bm{\phi}_{k}^{(t)}}{\sum_{k^{\prime}=1}^{K}\mathbb{P}\left(R_{n}^{c_{n}}\big{|}\bm{\mu}^{(t)}_{k^{\prime}},\bm{\Sigma}^{(t)}_{k^{\prime}},\pi_{k^{\prime}}^{n}=1\right)\bm{\phi}_{k^{\prime}}^{(t)}}.$
(14)

Note that (14) allows us to directly simulate $\hat{\pi}_{k}^{n}$ for all $k$
and $n$ as an independent subroutine. In particular, given
$\bm{\theta}^{(t)}$, to estimate
$\mathbb{P}\left(R_{n}^{c_{n}}\big{|}\bm{\mu}^{(t)}_{k^{\prime}},\bm{\Sigma}^{(t)}_{k^{\prime}},\pi_{k^{\prime}}^{n}=1\right)$,
we can generate many Monte Carlo samples from the Gaussian distribution
$\mathcal{N}\left(\bm{\mu}^{(t)}_{k},\bm{\Sigma}^{(t)}_{k}\right)$ and record
the fraction in $R_{n}^{c_{n}}$. As a result, we can use (14) to estimate
$\hat{\pi}_{k}^{n}$.

Combining this additional step to estimate $\hat{\pi}^{n}_{k}$, we can proceed
with the EM algorithm for (13):

  1. 1.

_Estimating $\hat{\pi}^{n}_{k}$:_ for $k=1,...,K$, generate $L$ samples
$\bm{x}_{k}^{(l)}$, $l=1,...,K$, from
$\mathcal{N}\left(\bm{\mu}^{(t)}_{k},\bm{\Sigma}^{(t)}_{k}\right)$. Let

$\hat{\mathbb{P}}\left(R_{n}^{c_{n}}\big{|}\bm{\mu}^{(t)}_{k},\bm{\Sigma}^{(t)}_{k}\right)=\frac{1}{L}\sum_{l=1}^{L}\mathbb{I}\left\\{x_{k}^{(l)}\in
R_{n}^{c_{n}}\right\\}$

and plug it into (14) to estimate $\hat{\pi}^{n}_{k}$.

  2. 2.

For all observations $n=1,\dots,N$ and $K$ components $k=1,\dots,K$, generate
$L$ samples $\bm{x}_{nk}^{(l)}$, $l=1,\dots,L$ from
$\mathcal{N}\left(\bm{\mu}_{k}^{(t)},\bm{\Sigma}_{k}^{(t)}\right)$ conditional
on $\bm{x}_{nk}^{(l)}\in R_{n}^{c_{n}}$, using Algorithm 1.

  3. 3.

_E-step:_ We represent (13) as

$\sum_{n=1}^{N}\sum_{k=1}^{K}\hat{\pi}^{n}_{k}\log(\phi_{k})+\frac{1}{L}\sum_{l=1}^{L}\sum_{n=1}^{N}\sum_{k=1}^{K}\hat{\pi}^{n}_{k}\log
f\left(\bm{x}_{nk}^{(l)}|\bm{\mu}_{k},\bm{\Sigma}_{k}\right).$ (15)

  4. 4.

_M-step:_ from the first term in (15) we can maximize $\phi_{k}$ for
$k=1,\dots,K$:

$\phi^{(t+1)}_{k}=\sum_{n=1}^{N}\frac{1}{N}\hat{\pi}^{n}_{k}.$

For the mean and variance, we have

$\displaystyle\bm{\mu}^{(t+1)}_{k}=$
$\displaystyle\frac{\sum_{l=1}^{L}\sum_{n=1}^{N}\hat{\pi}_{k}^{n}x_{nk}^{(l)}}{L\sum_{n=1}^{N}\hat{\pi}^{n}_{k}},$
$\displaystyle\bm{\Sigma}^{(t+1)}_{k}=$
$\displaystyle\frac{\sum_{l=1}^{L}\sum_{n=1}^{N}\hat{\pi}^{n}_{k}\bm{\left
missing}(\bm{x}_{nk}^{(l)}-\bm{\mu}_{k}^{(t+1)}\right)\left(\bm{x}_{nk}^{(l)}-\bm{\mu}_{k}^{(t+1)}\right)^{\top}}{L\sum_{n=1}^{N}\hat{\pi}^{n}_{k}}.$

  5. 5.

We repeat the above steps until the parameters converge.

##  5 Theoretical Properties

In this section, we present the statistical properties of the base model in
Section 3 and the EM algorithm. We first establish conditions for the
identifiability of the problem. We then show some properties of the
convergence of the EM algorithm.

###  5.1 Identifiability

The identifiability is concerned with the notion that whether the parameter
$\bm{\theta}=(\bm{\mu},\bm{\Sigma})$ can ever be learned from the data. For
example, suppose the firm has collected a large data set after keeping the
price menu for a long period. Equivalently, the probabilities of the IC
polyhedra in (3.1) can be accurately estimated, from the fraction of customers
purchasing each bundle. Identifiability guarantees that knowing the
probability of the regions (IC polyhedra) uniquely determines the parameters.

To simplify the notation and formulate the problem in a way of independent
theoretical interest, consider a class of regions
$\mathcal{R}\triangleq\\{R_{1},\dots,R_{K}\\}$ where
$R_{k}\subset\mathbb{R}^{I}$. The identifiability of our model is defined as

######  Definition 1 (Identifiability)

For a family of multivariate normal distribution in $\mathbb{R}^{I}$

$\mathcal{P}=\left\\{f(\bm{v}|\bm{\theta})=\frac{1}{\sqrt{(2\pi)^{I}|\bm{\Sigma}|}}\exp\left(-\frac{1}{2}(\bm{v}-\bm{\mu})^{\top}\bm{\Sigma}^{-1}(\bm{v}-\bm{\mu})\right)\Big{|}\bm{\theta}=(\bm{\mu},\bm{\Sigma})\right\\},$

we say that the model is identifiable over $\mathcal{R}$ if we have

$\int_{\bm{v}\in R_{i}}f(\bm{v}|\bm{\theta})dv=\int_{\bm{v}\in
R_{i}}f(\bm{v}|\bm{\theta})d\bm{v}\quad\text{for all
}R_{i}\in\mathcal{R}\implies\bm{\theta}_{1}=\bm{\theta}_{2}.$

In other words, for the collection of IC polyhedra $\mathcal{R}$ induced by
the prices of the $J$ bundles, the model is identifiable when the
probabilities of the regions in $\mathcal{R}$ uniquely determines the
parameter $\bm{\theta}$. Note that identifiability is almost a prerequisite of
any statistical procedure: if there are multiple parameters generating the
same distribution of the data, then there is no hope to estimate the parameter
accurately. Identifiability is not a trivial requirement that is always
satisfied in this model. Consider the following two examples.

###### Example 1

The firm always bundles two products together in the price menu. That is, they
are never sold separately. In this case, the model is not identifable over the
resulting IC polyhedra because it is not possible to estimate the marginal
mean of either product.

###### Example 2

Consider a single product $I=J=1$ always offered at price $\$10$. In this
case, $\mathcal{R}$ can be represented by two regions (intervals):
$R_{1}=(-\infty,10)$ and $R_{2}=[10,+\infty)$. The model is not identifiable
because $\bm{\theta}_{1}$ and $\bm{\theta}_{2}$ satisfying
$(\mu_{1}-10)/\sigma_{1}=(\mu_{2}-10)/\sigma_{2}$ always leads to the same
probability
$\mathbb{P}(R_{1}|\bm{\theta}_{1})=\mathbb{P}(R_{1}|\bm{\theta}_{2})$ and
$\mathbb{P}(R_{2}|\bm{\theta}_{1})=\mathbb{P}(R_{2}|\bm{\theta}_{2})$.

Therefore, we need to establish conditions for $\mathcal{R}$ that guarantee
identifiability. Intuitively, since there are $K$ regions and $(I+1)^{2}/2$
effective parameters in the model, we require $K\geq(I+1)^{2}/2$. However,
because of the nonlinearity in the probabilities of the regions, it is hard to
verify the sufficiency or the necessity of this condition.

Based on the application, we provide the following sufficient condition that
can be easily satisfied in practice. Suppose the firms are selling products
separately. Technically, the firm determines the prices for the $I$ individual
prices and offering $J=2^{I}$ bundles; the price of a bundle $J$ is the sum of
the prices of the included individual products. Separate selling is prevalent
in retailing, when firms don’t bundle products together.

######  Proposition 1 (Identifiability under separate selling)

When the firm sells products separately, as long as at least two prices are
set for each individual products infinitely often, then the model is
identifiable.

Note that Proposition 1 does not require all combinations of the two prices of
the products to be offered infinitely often, which would result in $2^{I}$
price menus. Instead, a simple and realistic scenario with $I+1$ price menus
would suffice: a price menu of regular prices of the products
$(p_{1},\dots,p_{I})$ and $I$ price menus in which only one product is on
discount, i.e., $(p_{1},\dots,p_{i-1},p_{i}^{\prime},p_{i+1},\dots,p_{I})$ for
product $i$. If the above $I+1$ price menus are offered infinitely often, then
$\mathcal{R}$ includes the following regions, whose probabilities turn out to
be sufficient for identifiability:

$\displaystyle\left\\{\bm{v}\in\mathbb{R}^{I},\bm{v}_{i}\leq
p_{i}\right\\},\;\left\\{\bm{v}\in\mathbb{R}^{I},\bm{v}_{i}\leq
p_{i}\right\\},\;\forall i=1,\dots,I$
$\displaystyle\left\\{\bm{v}\in\mathbb{R}^{I},\bm{v}_{i_{1}}\leq
p_{i_{1}},\bm{v}_{i_{2}}\leq
p_{i_{2}}\right\\},\;\left\\{\bm{v}\in\mathbb{R}^{I},\bm{v}_{i_{1}}\leq
p^{\prime}_{i_{1}},\bm{v}_{i_{2}}\leq
p_{i_{2}}\right\\},\;\left\\{\bm{v}\in\mathbb{R}^{I},\bm{v}_{i_{1}}\leq
p_{i_{1}},\bm{v}_{i_{2}}\leq p^{\prime}_{i_{2}}\right\\}\;\forall
i_{1},i_{2}=1,\dots,I$

By Example 2, just setting one price for a product is not sufficient for
identifiability. Therefore, Proposition 1 provides a sufficient and necessary
condition for identifiability under separate selling.

###  5.2 Convergence of the EM Algorithm

The EM algorithm has been shown to converge to a local maximum of the
likelihood function in the seminal papers by Dempster et al. (1977), Wu
(1983). However, the convergence to the global maximum of the likelihood
function, which corresponds to the true parameters with large samples, remains
an active research area. The new theoretical advances Balakrishnan et al.
(2017) in the study of EM algorithms allow one to characterize the conditions
for the convergence to the global maximum of the likelihood function. In this
section, we use the approach in Balakrishnan et al. (2017) to analyze our
problem. To simplify the analysis, we focus on the population-level
$Q$-function, following the first step in Balakrishnan et al. (2017). That is,
the EM algorithm is run on a dataset with an infinite sample size. We also
assume the IC polyhedra are identical for all consumers, i.e., they face the
same price menu. As a result, we can focus on a partition of $\mathbb{R}^{n}$,
denoted $\\{R_{k}\\}_{k=1}^{K}$. Moreover, we assume the true
$\bm{\Sigma}^{*}$ is known and the EM algorithm finds the maximum likelihood
estimator for $\bm{\mu}$. The $Q$-function in this simplified setting can be
expressed as

$\displaystyle Q(\bm{\mu}^{\prime}\mid\bm{\mu})=\sum_{k=1}^{K}\int_{\bm{x}\in
R_{k}}f(\bm{x}|\bm{\mu}^{*})d\bm{x}\int_{\bm{v}\in
R_{k}}\frac{f(\bm{v}\mid\bm{\mu})}{\left(\int_{\bm{x}\in
R_{k}}f(\bm{x}|\bm{\mu})d\bm{x}\right)}\log
f(\bm{v}|\bm{\mu}^{\prime})d\bm{v}.$ (16)

Compared to (3.2), we remove the dependence on $\bm{\Sigma}$ and use the same
partition. Moreover, because the function is defined at the population level,
the probability $\int_{x\in R_{k}}f(x|\bm{\mu}^{*})$ in each region replaces
the samples in (3.2). The population-level EM algorithm updates $\bm{\mu}$
iteratively according to the following operator:

$\displaystyle
M(\bm{\mu})\triangleq\text{argmax}_{\bm{\mu}^{\prime}}Q(\bm{\mu}^{\prime}\mid\bm{\mu}).$

More precisely, given initialization $\bm{\mu}^{(0)}$, the EM algorithm uses
$\bm{\mu}^{(t+1)}=M\left(\bm{\mu}^{(t)}\right)$ to obtain a sequence
$\left\\{\bm{\mu}^{(t)}\right\\}_{t=0}^{\infty}$. The theoretical question is
whether we have $\bm{\mu}^{(t)}\to\bm{\mu}^{*}$. Note that the true value
$\bm{\mu}^{*}$ always maximizes the population-level likelihood function (Van
der Vaart 2000). Moreover, it satisfies the _self-consistency_ condition
(McLachlan and Krishnan 2007):
$\bm{\mu}^{*}=\text{argmax}_{\bm{\mu}^{\prime}}Q(\bm{\mu}^{\prime}\mid\bm{\mu}^{*})$.
Therefore, $\bm{\mu}^{*}$ is a fixed point of $M(\bm{\mu})=\bm{\mu}$.
Balakrishnan et al. (2017) establish a set of conditions that guarantee
$M(\bm{\mu})$ to be a contraction mapping when $\bm{\mu}^{(0)}$ is in a
neighborhood of $\bm{\mu}^{*}$, denoted by $\mathbb{B}(r;\bm{\mu}^{*})$ with
an Euclidean radius $r$. Since the contraction mapping converges to the fixed
point, our goal in this section is to verify the set of conditions in our
problem.

Among the set of technical conditions in Balakrishnan et al. (2017), there are
two key conditions relevant to our model:

  * •

_Concavity of $Q(\bm{\mu}\mid\bm{\mu}^{*})$_. Define $q(\bm{\mu})\triangleq
Q(\bm{\mu}\mid\bm{\mu}^{*})=\sum_{k=1}^{K}\int_{\bm{v}\in
R_{k}}f(\bm{v}\mid\bm{\mu}^{*})\log f(\bm{v}|\bm{\mu})\mathrm{d}\bm{v}$ as the
$Q$-function when the current $\bm{\mu}$ is the true value $\bm{\mu}^{*}$. The
required condition states that there exists $\lambda>0$, such that

$\displaystyle q(\bm{\mu}_{1})-q(\bm{\mu}_{2})-\left\langle\nabla
q(\bm{\mu}_{2}),\bm{\mu}_{1}-\bm{\mu}_{2}\right\rangle\leq-\frac{\lambda}{2}\|\bm{\mu}_{1}-\bm{\mu}_{2}\|^{2}_{2},$
(17)

for all pairs of $\bm{\mu}_{1},\bm{\mu}_{2}\in\mathbb{B}(r;\bm{\mu}^{*})$,
where $\nabla$ is the gradient operator of a function. Represented in the
Taylor expansion, (17) ensures that $q(\bm{\mu})$ is strictly concave when
$\bm{\mu}$ is in the neighborhood.

  * •

_First-order stability_. Note that $\nabla Q(\bm{\mu}^{\prime}\mid\bm{\mu})$
is the gradient of $Q$ at $\bm{\mu}^{\prime}$ given the current $\mu$. This
condition requires that for some $\gamma<\lambda$, we have

$\displaystyle\|\nabla Q(M(\bm{\mu})\mid\bm{\mu}^{*})-\nabla
Q(M(\bm{\mu})\mid\bm{\mu})\|_{2}\leq\gamma\|\bm{\mu}-\bm{\mu}^{*}\|_{2}.$ (18)

In other words, the $Q$-function needs to be sufficiently stable around
$\bm{\mu}^{*}$ when taking the gradient. The relationship $\gamma<\lambda$ is
essential, as $\gamma/\lambda$ turns to be the contraction constant.

Note that neither condition can be easily checked in our model. The main
contribution of this section is to identify reasonable assumptions for the
conditions to hold.

Let $\bm{X}\sim\mathcal{N}(\bm{\mu}^{*},\bm{\Sigma}^{*})$ be the distribution
of the valuation and $\bm{Z}=(\bm{\Sigma}^{*})^{-1/2}(\bm{X}-\bm{\mu}^{*})$ be
the linear transformation of $\bm{X}$ so that
$\bm{Z}\sim\mathcal{N}(\bm{0},\bm{I})$. Let
$R_{k}^{\prime}=(\bm{\Sigma}^{*})^{-1/2}(R_{k}-\bm{\mu}^{*})$ be the
transformed regions in the $\bm{Z}$ space. Therefore,
$\mathbb{E}[\bm{Z}|R^{\prime}_{k}]$ is the conditional expectation of $\bm{Z}$
in $R_{k}^{\prime}$. We let $R^{\prime}$ be a categorical random variable so
that $\mathbb{P}(R^{\prime}=R_{k}^{\prime})=\mathbb{P}(\bm{Z}\in
R^{\prime}_{k})$. The assumption below provides a sufficient condition for the
EM algorithm to converge. We use $\lambda_{min}(\cdot)$ to denote the minimal
eigenvalue of a positive semi-definite matrix.

###### Assumption 1

$\lambda_{min}(\mathrm{Var}(\mathbb{E}[\bm{Z}\mid R^{\prime}]))\geq\epsilon$
for some $\epsilon>0$.

Note that $\mathbb{P}(\mathbb{E}[\bm{Z}\mid R^{\prime}]=\mathbb{E}[\bm{Z}\mid
R^{\prime}_{k}])=\mathbb{P}(\bm{Z}\in R^{\prime}_{k})$ and thus
$\mathbb{E}(\bm{Z}\mid R^{\prime})$ is an $I$-dimensional random vector
supported on $K$ points. Assumption 1 states that the eigenvalues of the
covariance matrix of the random vector are bounded away from zero.

To provide some intuition for the assumption, consider the case $K=1$, i.e.,
there is one region. In this case, $\mathbb{E}(\bm{Z}\mid R^{\prime})$ has a
single value and $\mathrm{Var}(\mathbb{E}[\bm{Z}\mid R^{\prime}])$ is a
degenerate $I\times I$ zero matrix. Assumption 1 is clearly not satisfied. It
is not surprising that the EM algorithm doesn’t converge to the true
parameters for this degenerate case because the model is not identifiable,
i.e., no information about the parameters can be learned from the data.

Now consider the other extreme. Suppose $R_{k}^{\prime}$ partitions the space
extremely finely so that $R_{k}^{\prime}$ almost degenerates to a single
point. In this case, $\mathbb{E}[Z|R_{k}^{\prime}]$ is itself approximately
$\mathcal{N}(0,\bm{I})$. Therefore, $\mathrm{Var}(\mathbb{E}[\bm{Z}\mid
R^{\prime}])\approx\bm{I}$ and apparently satisfies Assumption 1. The two
examples point to the intuition that the EM algorithm performs well when there
are many small regions. This is consistent with the intuition mentioned after
Example 2.

Assumption 1 guarantees the convergence of the EM algorithm to the true
parameter locally.

###### Theorem 1

Suppose Assumption 1 holds. There exists a neighborhood
$\mathbb{B}(r;\bm{\mu}^{*})$ of $\bm{\mu}^{*}$ such that for
$\mu^{(0)}\in\mathbb{B}(r;\bm{\mu}^{*})$ the EM algorithm guarantees

$\displaystyle\|\bm{\mu}^{(t)}-\bm{\mu}^{*}\|_{2}\leq\left(1-\epsilon/2\right)^{t}\|\bm{\mu}^{(0)}-\bm{\mu}^{*}\|_{2}.$

The detailed proof for Theorem 1 can be found in the online appendix.

##  6 Numerical Experiments

In this section, we conduct the numerical experiments to study the performance
of the EM algorithm. Because there are few methods that can be applied to
bundle sales data, our main benchmark is an adapted MNL model, treating each
bundle as an individual product. We also implement the algorithm in Jedidi et
al. (2003).

We generate multiple synthetic datasets with $I=2,4,6$ products and
$N=1000,2000,4000$ transactions. Each transaction records a customer’s
consideration set, their prices and the customer’s choice (including no-
purchase). Customers’ consideration set is randomly drawn with probability 0.5
from the product set and also includes random bundles of sizes 2 and 3,
consisting of the products viewed by the customer. We consider that prices for
individual products following a uniform distribution $[\mu_{i}-3,\mu_{i}+3]$,
where $\mu_{i}$ denote the mean of the valuation distribution for product $i$.
For bundles, we get the summation of individual prices within the bundle and
generate a discount from a uniform distribution [0,5]. For example, consider a
bundle including products $\\{1,2,3\\}$, with respective prices 10, 20, 30. If
the random discount is 2, then the bundle price is $10+20+30-2=58$. For the
ground truth, we consider two underlying demand models: _Gaussian
distribution_ , where its mean is drawn from a uniform distribution
$\bm{\mu}\sim\mathcal{U}[10,50]^{I}$, and the covariance matrix is
$\bm{\Sigma}=AA^{\top}$, where the component of $A$ follows
$\mathcal{U}[1,2]$; _Gumbel distribution_ , where we consider an independent
Gumbel distribution for each product $i$ with parameters $\bm{\mu}$ and the
diagonals of $\bm{\Sigma}$. We use an additive function to evaluate customers’
valuation for bundles.

###  6.1 The Base Model

In this subsection, we compare the performance of the EM algorithm with the
MNL model treating each bundle accounts as a different/separate product. We
use 80% of our dataset for training and the remaining 20% for testing the
model prediction. We measure the estimation by the $\ell_{1}$-error
$(\|\bm{\mu}^{(t+1)}-\bm{\mu}^{*}\|_{1}+\|\bm{\Sigma}^{(t+1)}-\bm{\Sigma}^{*}\|_{1})/\left(K^{2}+K\right)$
which is the average absolute difference of the estimated component and the
true parameter.

EM convergence. Figure 1 depicts the error of the EM algorithm over the
iterations. Figure 1(a) shows that the EM algorithm almost recovers the true
parameters (i.e., $error<0.1$) in 15 iterations. Figure 1(b) demonstrates that
the more products or transactions, the lower the error, and the faster the
algorithm converges to the true parameters.

Figure 1: Errors in the iterations of the EM algorithm

(a) 6 products, 36 bundles (b) 2000 observations

EM vs. MNL. In Figure 2, we measure the performance of the EM algorithm and
the adapted MNL model in terms of root-mean square error (RMSE). This is
because we cannot evaluate the error against the true parameter as the models
are misspecified. Figure 2 demonstrates that when the ground truth follows
either the Gaussian or Gumbel distributions, the EM algorithm achieves
significantly lower RMSE in predicting the shopping behavior of the customers
than the adapted MNL model. EM algorithm performs better because it considers
that bundles are a basket of multiple products while the MNL model ignores
this crucial feature and considers all choices are independent.

Figure 2: RMSE comparison of the EM algorithm and the MNL model

(a) Gaussian distribution (b) Gumbel distribution

EM vs. Hierarchical Bayesian Model. In this experiment, we compare the
efficacy of the proposed EM algorithm with the approach in Jedidi et al.
(2003), in which they form a multinomial probit likelihood which involves
integration over high-dimensional multivariate normal distribution, and
propose a hierarchical Bayesian method to infer the parameters based on the
posterior distribution using Metropolis Hasting method. We implement their
method in our setting by considering the posterior probability
$\mathbb{P}\left(\bm{\mu}^{(t)},\bm{\sigma}^{(t)}\mid
R_{n}^{c_{n}}\right)\propto\mathbb{P}\left(R_{n}^{c_{n}}\mid\bm{\mu}^{(t)},\bm{\sigma}^{(t)}\right)\mathbb{P}\left(\bm{\mu}^{(t)}\right)\mathbb{P}\left(\bm{\sigma}^{(t)}\right)$
(see the Metropolis Hasting algorithm in Appendix B.3). We run this model for
two products with fixed covariance matrix $\begin{bmatrix}4\ 1\\\ 1\
5\end{bmatrix}$. We generate 1000 random observations with mean
$[\mu_{1},\mu_{2}]=[0.5,0.5]$. Our EM algorithm converges to $[0.562,0.588]$
after 20 iterations in 2.368 minutes, which implies an error of $0.075$. For
the Metropolis Hasting algorithm, we iterate the model 10000 times for about 7
days, and it converges to $[0.391,\ 0.685]$, with an error $0.147$. This shows
the advantage of our model over the Metropolis Hasting algorithm to infer the
true parameters.

###  6.2 Censored Data

In this subsection, we show the efficacy of our model to recover true
parameters when no-purchases are censored from the sales data and there is no
record of such transactions. In this experiment, we randomly generate a
complete dataset (including no-purchase transactions) with $I=2,4,6$ products
and $N=1000,2000,4000$ transactions as explained in Subsection 6.1 with an
exception that the number of observation sets, price menus and bundles are
limited to 1, 3, and 6. As explained in Section 4.1, we need the seller to
collect data for the same price menu in order to estimate the model from
censored demand. This is why we use a smaller number of price menus than the
experiments for the base model.

We perform the base EM algorithm on the complete dataset and then, we remove
the no-purchase transactions (approx. 5-15% of the transactions) and perform
the EM algorithm for censored data. Figure 3(a) shows that the EM algorithm
for censored data converges, and the error is less when the number of
transactions increases. Moreover, we have that the difference between the RMSE
of complete and censored datasets is less than 0.002 (0.6%).

Figure 3: EM censored performance EM censored for 4 products, 3 bundles, and 3
price menus

## References

  * Abdallah and Vulcano [2021] T. Abdallah and G. Vulcano.  Demand estimation under the multinomial logit model from sales transaction data.  _Manufacturing Service Oper. Management_ , 23(5):1196–1216, 2021.
  * Abdallah et al. [2021] T. Abdallah, A. Asadpour, and J. Reed.  Large-scale bundle-size pricing: A theoretical analysis.  _Oper. Res._ , 69(4):1158–1185, 2021.
  * Babaioff et al. [2020] M. Babaioff, N. Immorlica, B. Lucier, and S. M. Weinberg.  A simple and approximately optimal mechanism for an additive buyer.  _Journal of the ACM (JACM)_ , 67(4):1–40, 2020\.
  * Bakos and Brynjolfsson [1999] Y. Bakos and E. Brynjolfsson.  Bundling information goods: Pricing, profits, and efficiency.  _Management science_ , 45(12):1613–1630, 1999\.
  * Balakrishnan et al. [2017] S. Balakrishnan, M. J. Wainwright, and B. Yu.  Statistical guarantees for the EM algorithm: From population to sample-based analysis.  _The Annals of Statistics_ , 45(1):77–120, 2017\.
  * Chen et al. [2022] N. Chen, X. Li, Z. Li, and C. Wang.  Component pricing with bundle size discount.  _Working Paper_ , 2022.
  * Chu et al. [2011] C. S. Chu, P. Leslie, and A. Sorensen.  Bundle-size pricing as an approximation to mixed bundling.  _The American Economic Review_ , pages 263–303, 2011.
  * Dempster et al. [1977] A. P. Dempster, N. M. Laird, and D. B. Rubin.  Maximum likelihood from incomplete data via the EM algorithm.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 39(1):1–22, 1977.
  * Farias et al. [2013] V. F. Farias, S. Jagabathula, and D. Shah.  A nonparametric approach to modeling choice with limited data.  _Management science_ , 59(2):305–322, 2013.
  * Fisher et al. [2020] M. Fisher, M. Olivares, and B. R. Staats.  Why empirical research is good for operations management, and what is good empirical operations management?  _Manufacturing & Service Operations Management_, 22(1):170–178, 2020.
  * Goodfellow et al. [2016] I. Goodfellow, Y. Bengio, and A. Courville.  _Deep learning_.  MIT press, 2016.
  * Hanson and Martin [1990] W. Hanson and R. K. Martin.  Optimal bundle pricing.  _Management Science_ , 36(2):155–174, 1990.
  * Hart and Nisan [2017] S. Hart and N. Nisan.  Approximate revenue maximization with multiple items.  _Journal of Economic Theory_ , 172:313–347, 2017.
  * Jedidi et al. [2003] K. Jedidi, S. Jagpal, and P. Manchanda.  Measuring heterogeneous reservation prices for product bundles.  _Marketing Science_ , 22(1):107–130, 2003.
  * Jiang et al. [2011] Y. Jiang, J. Shang, C. F. Kemerer, and Y. Liu.  Optimizing e-tailer profits and customer savings: Pricing multistage customized online bundles.  _Marketing Science_ , 30(4):737–752, 2011.
  * Kök et al. [2008] A. G. Kök, M. L. Fisher, and R. Vaidyanathan.  Assortment planning: Review of literature and industry practice.  _Retail supply chain management_ , pages 99–153, 2008.
  * Lee and Scott [2012] G. Lee and C. Scott.  EM algorithms for multivariate gaussian mixture models with truncated and censored data.  _Computational Statistics & Data Analysis_, 56(9):2816–2829, 2012.
  * Letham et al. [2014] B. Letham, W. Sun, and A. Sheopuri.  Latent variable copula inference for bundle pricing from retail transaction data.  In _International conference on machine learning_ , pages 217–225. PMLR, 2014.
  * Li et al. [2020] X. Li, H. Sun, and C. P. Teo.  Convex optimization for bundle size pricing problem.  In _Proceedings of the 21st ACM Conference on Economics and Computation_ , pages 637–638, 2020.
  * Lindsey and Ryan [1998] J. C. Lindsey and L. M. Ryan.  Methods for interval-censored data.  _Statistics in medicine_ , 17(2):219–238, 1998\.
  * Ma and Simchi-Levi [2020] W. Ma and D. Simchi-Levi.  Constructing demand curves from a single observation of bundle sales.  _Working Paper_ , 2020.
  * Ma and Simchi-Levi [2021] W. Ma and D. Simchi-Levi.  Reaping the benefits of bundling under high production costs.  In _International Conference on Artificial Intelligence and Statistics_ , pages 1342–1350. PMLR, 2021.
  * marquis de Laplace [1840] P. S. marquis de Laplace.  _Essai philosophique sur les probabilités_.  Bachelier, 1840.
  * McLachlan and Krishnan [2007] G. J. McLachlan and T. Krishnan.  _The EM algorithm and extensions_ , volume 382.  John Wiley & Sons, 2007.
  * Natarajan et al. [2000] R. Natarajan, C. E. McCulloch, and N. M. Kiefer.  A monte carlo EM method for estimating multinomial probit models.  _Computational statistics & data analysis_, 34(1):33–50, 2000.
  * Newman et al. [2014] J. P. Newman, M. E. Ferguson, L. A. Garrow, and T. L. Jacobs.  Estimation of choice-based models using sales data from a single firm.  _Manufacturing & Service Operations Management_, 16(2):184–197, 2014.
  * Şimşek and Topaloglu [2018] A. S. Şimşek and H. Topaloglu.  An expectation-maximization algorithm to estimate the parameters of the Markov chain choice model.  _Operations Research_ , 66(3):748–760, 2018.
  * Stigler [1963] G. J. Stigler.  United states v. loew’s inc.: A note on block-booking.  _The Supreme Court Review_ , 1963:152–157, 1963.
  * Subramanian and Harsha [2021] S. Subramanian and P. Harsha.  Demand modeling in the presence of unobserved lost sales.  _Management Science_ , 67(6):3803–3833, 2021\.
  * Sun et al. [2021] H. Sun, X. Li, and C.-P. Teo.  Product bundle recommendation and pricing: How to make it work?  _Available at SSRN 3874843_ , 2021.
  * Talluri and Van Ryzin [2004] K. Talluri and G. Van Ryzin.  Revenue management under a general discrete choice model of consumer behavior.  _Management Science_ , 50(1):15–33, 2004.
  * Van der Vaart [2000] A. W. Van der Vaart.  _Asymptotic statistics_ , volume 3.  Cambridge university press, 2000.
  * Van Ryzin and Vulcano [2015] G. Van Ryzin and G. Vulcano.  A market discovery algorithm to estimate a general class of nonparametric choice models.  _Management Science_ , 61(2):281–300, 2015.
  * Vulcano et al. [2012] G. Vulcano, G. Van Ryzin, and R. Ratliff.  Estimating primary demand for substitutable products from sales transaction data.  _Operations Research_ , 60(2):313–334, 2012.
  * Wu [1983] C. J. Wu.  On the convergence properties of the EM algorithm.  _The Annals of statistics_ , pages 95–103, 1983.
  * Wu et al. [2008] S.-y. Wu, L. M. Hitt, P.-y. Chen, and G. Anandalingam.  Customized bundle pricing for information goods: A nonlinear mixed-integer programming approach.  _Management Science_ , 54(3):608–622, 2008.

Online Appendix

##  A Proofs

###  A.1 The Derivation of the M-step

$\displaystyle Q\left(\bm{\theta}|\bm{\theta}^{(t)}\right)$
$\displaystyle\triangleq\mathbb{E}\left[l(\bm{\theta};\bm{X},\bm{Z})\bigg{|}\bm{X},\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right]$
$\displaystyle=\sum_{n=1}^{N}\int_{\bm{v}\in
R_{n}^{c_{n}}}\frac{f\left(\bm{v}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)}{\int_{R_{n}^{c_{n}}}f\left(\bm{x}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)\mathrm{d}\bm{x}}\log
f(\bm{v}|\bm{\mu},\bm{\Sigma})\mathrm{d}\bm{v}.$ (A.1)

To find a new estimation for $\bm{\theta}$, the _M-step_ in the EM algorithm
maximizes $Q\left(\bm{\theta}|\bm{\theta}^{(t)}\right)$ and let
$\bm{\theta}^{(t+1)}=\text{argmax}_{\bm{\theta}}Q\left(\bm{\theta}|\bm{\theta}^{(t)}\right)$.
Thus, we first solve

$\displaystyle\partial\mathbb{E}\left[l(\bm{\theta};\bm{X},\bm{Z})\mid\bm{X},\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right]/\partial\bm{\mu}=0,$
(A.2)

and then we solve

$\displaystyle\partial\mathbb{E}\left[l(\bm{\theta};\bm{X},\bm{Z})\mid\bm{X},\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right]/\partial\bm{\Sigma}=0.$
(A.3)

Equality (A.2) equates to

$\displaystyle\frac{\partial\mathbb{E}\left[l(\bm{\theta};\bm{X},\bm{Z})\mid\bm{X},\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right]}{\partial\bm{\mu}}=$
$\displaystyle\frac{\partial}{\partial\bm{\mu}}\sum_{n=1}^{N}\int_{\bm{v}\in
R_{n}^{c_{n}}}\frac{f\left(\bm{v}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)}{\int_{R_{n}^{c_{n}}}f(\bm{x}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)})\mathrm{d}\bm{x}}\log
f(\bm{v}|\bm{\mu},\bm{\Sigma})\mathrm{d}\bm{v}$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{N}\int_{\bm{v}\in
R_{n}^{c_{n}}}\frac{f\left(\bm{v}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)}{\int_{R_{n}^{c_{n}}}f\left(\bm{x}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)\mathrm{d}\bm{x}}\frac{\partial}{\partial\bm{\mu}}\log
f(\bm{v}|\bm{\mu},\bm{\Sigma})\mathrm{d}\bm{v}$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{N}\int_{\bm{v}\in
R_{n}^{c_{n}}}\frac{f\left(\bm{v}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)}{\int_{R_{n}^{c_{n}}}f\left(\bm{x}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)\mathrm{d}\bm{x}}\frac{1}{f(\bm{v}|\bm{\mu},\bm{\Sigma})}\bm{\Sigma}^{-1}(\bm{v}-\bm{\mu})f(\bm{v}|\bm{\mu},\bm{\Sigma})\mathrm{d}\bm{v}$
$\displaystyle=$ $\displaystyle\sum_{n=1}^{N}\int_{\bm{v}\in
R_{n}^{c_{n}}}\bm{\Sigma}^{-1}(\bm{v}-\bm{\mu})\frac{f\left(\bm{v}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)}{\int_{R_{n}^{c_{n}}}f\left(\bm{x}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)\mathrm{d}\bm{x}}\mathrm{d}\bm{v}$
$\displaystyle=$ $\displaystyle\sum_{n=1}^{N}\int_{\bm{v}\in
R_{n}^{c_{n}}}\bm{v}\frac{f\left(\bm{v}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)}{\int_{R_{n}^{c_{n}}}f\left(\bm{x}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)\mathrm{d}\bm{x}}\mathrm{d}\bm{v}$
$\displaystyle-\sum_{n=1}^{N}\int_{\bm{v}\in
R_{n}^{c_{n}}}\bm{\mu}\frac{f\left(\bm{v}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)}{\int_{R_{n}^{c_{n}}}f\left(\bm{x}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)\mathrm{d}\bm{x}}\mathrm{d}\bm{v}$
$\displaystyle=$ $\displaystyle\sum_{n=1}^{N}\int_{\bm{v}\in
R_{n}^{c_{n}}}\bm{v}\frac{f\left(\bm{v}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)}{\int_{R_{n}^{c_{n}}}f\left(\bm{x}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)\mathrm{d}\bm{x}}\mathrm{d}\bm{v}-N\bm{\mu}=0.$
(A.4)

Thus, we have

$\displaystyle\bm{\mu}=\frac{\sum_{n=1}^{N}\int_{\bm{v}\in
R_{n}^{c_{n}}}\bm{v}\frac{f\left(\bm{v}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)}{\int_{R_{n}^{c_{n}}}f\left(\bm{x}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)\mathrm{d}\bm{x}}\mathrm{d}\bm{v}}{N}=\frac{\sum_{n=1}^{N}\bm{\mathbb{E}}[\bm{v}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)},R_{n}^{c_{n}}]}{N}.$
(A.5)

Equality (A.3) equates to

$\displaystyle\frac{\partial\mathbb{E}\left[l(\bm{\theta};\bm{X},\bm{Z})\mid\bm{X},\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right]}{\partial\bm{\Sigma}}$
$\displaystyle=\sum_{n=1}^{N}\frac{\partial}{\partial\bm{\Sigma}}\int_{\bm{v}\in
R_{n}^{c_{n}}}\frac{f\left(\bm{v}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)}{\int_{R_{n}^{c_{n}}}f(\bm{x}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)})\mathrm{d}\bm{x}}\log
f(\bm{v}|\bm{\mu},\bm{\Sigma})\mathrm{d}\bm{v}$
$\displaystyle=\sum_{n=1}^{N}\int_{\bm{v}\in
R_{n}^{c_{n}}}\frac{f\left(\bm{v}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)}{\int_{R_{n}^{c_{n}}}f(\bm{x}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)})\mathrm{d}\bm{x}}(-\frac{1}{2}\bm{\Sigma}^{-1}+\frac{1}{2}\bm{\Sigma}^{-1}(\bm{v}-\bm{\mu})(\bm{v}-\bm{\mu})^{T}\bm{\Sigma}^{-1})\mathrm{d}\bm{v}$
(A.6) $\displaystyle=-\frac{1}{2}\bm{\Sigma}^{-1}\sum_{n=1}^{N}\int_{\bm{v}\in
R_{n}^{c_{n}}}\frac{f\left(\bm{v}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)}{\int_{R_{n}^{c_{n}}}f(\bm{x}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)})\mathrm{d}\bm{x}}\mathrm{d}\bm{v}$
$\displaystyle\quad+\frac{1}{2}\bm{\Sigma}^{-1}\sum_{n=1}^{N}\int_{\bm{v}\in
R_{n}^{c_{n}}}\frac{f\left(\bm{v}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)}{\int_{R_{n}^{c_{n}}}f(\bm{x}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)})\mathrm{d}\bm{x}}(\bm{v}-\bm{\mu})(\bm{v}-\bm{\mu})^{T}\bm{\Sigma}^{-1}\mathrm{d}\bm{v}=0.$

Thus, we have

$\displaystyle N\bm{\Sigma}$ $\displaystyle=\sum_{n=1}^{N}\int_{\bm{v}\in
R_{n}^{c_{n}}}(\bm{v}-\bm{\mu})(\bm{v}-\bm{\mu})^{T}\frac{f\left(\bm{v}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)}{\int_{R_{n}^{c_{n}}}f(\bm{x}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)})\mathrm{d}\bm{x}}\mathrm{d}\bm{v}$
$\displaystyle=\sum_{n=1}^{N}\int_{\bm{v}\in
R_{n}^{c_{n}}}(\bm{v}\bm{v}^{T}-2\bm{\mu}^{T}v+\bm{\mu}\bm{\mu}^{T})\frac{f\left(\bm{v}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)}{\int_{R_{n}^{c_{n}}}f(\bm{x}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)})\mathrm{d}\bm{x}}\mathrm{d}\bm{v}$
$\displaystyle=\sum_{n=1}^{N}\mathbb{E}[\bm{v}\bm{v}^{T}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)},R_{n}^{c_{n}}]-2\bm{\mu}\mathbb{E}[\bm{v}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)},R_{n}^{c_{n}}]+\bm{\mu}\bm{\mu}^{T}$
$\displaystyle=\sum_{n=1}^{N}(\mathbb{E}[\bm{v}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)},R_{n}^{c_{n}}]-\bm{\mu})(\mathbb{E}[\bm{v}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)},R_{n}^{c_{n}}]-\bm{\mu})^{T}$
$\displaystyle\quad+\mathbb{E}[\bm{v}v^{T}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)},R_{n}^{c_{n}}]-\mathbb{E}[\bm{v}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)},R_{n}^{c_{n}}]\mathbb{E}[\bm{v}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)},R_{n}^{c_{n}}]^{T},$

which is equivalent to

$\displaystyle\bm{\Sigma}=\frac{\sum_{n=1}^{N}(\mathbb{E}[\bm{v}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)},R_{n}^{c_{n}}]-\bm{\mu})(\mathbb{E}[\bm{v}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)},R_{n}^{c_{n}}]-\bm{\mu})^{T}+\bm{U}^{n}}{N},$
(A.7)

where
$\bm{U}^{n}\triangleq\mathbb{E}[\bm{v}\bm{v}^{T}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)},R_{n}^{c_{n}}]-\mathbb{E}[\bm{v}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)},R_{n}^{c_{n}}]E[\bm{v}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)},R_{n}^{c_{n}}]^{T}.$

###  A.2 Identifiability

###### Proof

Proof of Proposition 1: We prove the result using induction on $I$.

Base case: We first show the statement holds for $I=1$. If the single product
is offered at two prices $p_{1}<p_{2}$, then the probabilities of the
following regions are determined: $R_{1}=(-\infty,p_{1}]$,
$R_{2}=(p_{1},p_{2}]$ and $R_{3}=(p_{2},+\infty)$. Since
$\mathbb{P}(R_{1}|\bm{\theta})=\Phi((p_{1}-\mu)/\sigma)$, where $\Phi(\cdot)$
is the CDF of the standard normal distribution, the probability of $R_{1}$
determines $(p_{1}-\mu)/\sigma$. Similarly, the probability of $R_{1}\cup
R_{2}$ uniquely determines $(p_{2}-\mu)/\sigma$. Because $p_{1}<p_{2}$, $\mu$
and $\sigma$ are uniquely determined as a result.

Inductive step: Suppose we have proved the claim for $I$ and now we prove it
for $I+1$ dimensions. Note that if we focus on products two to $I+1$, they are
offered separately at two prices. By the inductive hypothesis, we can uniquely
determine $\mu_{2},\mu_{3},\dots,\mu_{I+1}$ and

$\begin{bmatrix}\sigma_{22}&...&\sigma_{2,I+1}\\\ \vdots&\ddots&\vdots\\\
\sigma_{I+1,1}&...&\sigma_{I+1,I+1}\end{bmatrix}.$

Suppose product one is offered at at least two prices $p_{1}<p_{1}^{\prime}$.
It implies the probabilities of the regions
$R^{\prime}=\\{\bm{v}\in\mathbb{R}^{I+1},\bm{v}_{1}\leq p_{1}\\}$ (the
fraction of customers who have bought product one at price $p_{1}$) and
$R^{\prime\prime}=\\{\bm{v}\in\mathbb{R}^{I+1},\bm{v}_{1}\leq
p_{1}^{\prime}\\}$ are determined. By the argument in the base case, we can
uniquely determine $\mu_{1}$ and $\sigma_{11}$.

Now it remains to show that $\sigma_{1k}=\sigma_{k1}$ for $k=2,\dots,I+1$ can
be uniquely determined. Without loss of generality, we focus on $\sigma_{12}$.
Suppose product two are offered at two prices $p_{2}$ and $p_{2}^{\prime}$.
Since products one and two are offered together, one of the combinations of
their individual prices, $(p_{1},p_{2})$, $(p_{1},p_{2}^{\prime})$,
$(p_{1}^{\prime},p_{2})$ and $(p_{1}^{\prime},p_{2}^{\prime})$, must have been
offered infinitely in the data. Without loss of generality, suppose
$(p_{1},p_{2})$ is offered. It implies the probability of the following region
is determined:

$R_{1}=\left\\{\bm{v}\in\mathbb{R}^{I+1},\bm{v}_{1}\geq p_{1},\bm{v}_{2}\geq
p_{2}\right\\},\;R_{2}=\left\\{\bm{v}\in\mathbb{R}^{I+1},\bm{v}_{1}\leq
p_{1},\bm{v}_{2}\leq p_{2}\right\\}.$

Indeed, the probability of $R_{1}$ ($R_{2}$) is the fraction of customers who
bought both (neither) products when they are offered at $p_{1}$ and $p_{2}$,
respectively. Note that the joint distribution of $(\bm{V}_{1},\bm{V}_{2})$ is
Gaussian with mean $(\mu_{1},\mu_{2})$ and covariance matrix
$[\sigma_{11},\sigma_{12};\sigma_{12},\sigma_{22}]$, all of which are uniquely
determined but $\sigma_{12}$. It remains to show that $\sigma_{12}$ is also
uniquely determined by $\mathbb{P}(R_{1})$ and $\mathbb{P}(R_{2})$.

For this purpose, without loss of generality, for the simplicity of notations,
we assume $\mu_{1}=\mu_{2}=0$ and $\sigma_{11}=\sigma_{22}=1$. In this case,
$\sigma_{12}\in[-1,1]$ is the correlation of $V_{1}$ and $V_{2}$. We next
prove the result by contradiction. That is, suppose there exist
$\sigma_{12}\neq\sigma_{12}^{\prime}$ such that
$\mathbb{P}(R_{1}|\sigma_{12})=\mathbb{P}(R_{1}|\sigma_{12}^{\prime})$ and
$\mathbb{P}(R_{2}|\sigma_{12})=\mathbb{P}(R_{2}|\sigma_{12}^{\prime})$.
Consider the following two linear transformations of $(V_{1},V_{2})$:

$\displaystyle\mathcal{T}:\;$
$\displaystyle(X_{1},X_{2})=\left(V_{1},-\frac{\sigma_{12}}{\sqrt{1-\sigma_{12}^{2}}}V_{1}+\frac{1}{\sqrt{1-\sigma_{12}^{2}}}V_{2}\right)$
$\displaystyle\mathcal{T}^{\prime}:\;$
$\displaystyle(X^{\prime}_{1},X^{\prime}_{2})=\left(V_{1},-\frac{\sigma^{\prime}_{12}}{\sqrt{1-\sigma_{12}^{\prime
2}}}V_{1}+\frac{1}{\sqrt{1-\sigma_{12}^{\prime 2}}}V_{2}\right).$

Note that both transformations keep $V_{3},\dots,V_{I+1}$ unchanged, which
don’t play a role in the subsequent proof. It is easy to verify that $X_{1}$
and $X_{2}$ ($X^{\prime}_{1}$ and $X^{\prime}_{2}$) are indepndent standard
normal random variables if the correlation between $V_{1}$ and $V_{2}$ is
$\sigma_{12}$ ($\sigma_{12}^{\prime}$). We next consider the transformed
region of $R_{1}$ under $\mathcal{T}$:

$\displaystyle\mathcal{T}(R_{1})=\left\\{\bm{x}\in\mathbb{R}^{I+1},\bm{x}_{1}\geq
p_{1},\sigma_{12}\bm{x}_{1}+\sqrt{1-\sigma_{12}^{2}}\bm{x}_{2}\geq
p_{2}\right\\}.$

The transformed region of $R_{2}$ and both transformed regions under
$\mathcal{T}^{\prime}$ can be defined analogously, by replacing $\geq$ by
$\leq$ and $\sigma_{12}$ by $\sigma_{12}^{\prime}$. We consider the two lines
of $\mathbb{R}^{2}$:
$\sigma_{12}\bm{x}_{1}+\sqrt{1-\sigma_{12}^{2}}\bm{x}_{2}=p_{2}$ and
$\sigma_{12}^{\prime}\bm{x}_{1}+\sqrt{1-\sigma_{12}^{\prime
2}}\bm{x}_{2}=p_{2}$. Because $\sigma_{12}\neq\sigma_{12}^{\prime}$, it is
easy to show that
$\sigma_{12}/\sqrt{1-\sigma_{12}^{2}}\neq\sigma_{12}^{\prime}/\sqrt{1-\sigma_{12}^{\prime
2}}$. As a result, the two lines are not parallel and intersect at, say,
$(a_{1},a_{2})$. Without loss of generality, suppose the first line is always
on top when $\bm{x}_{1}<a_{1}$, i.e.,

$\frac{1}{\sqrt{1-\sigma_{12}^{2}}}(p_{2}-\sigma_{12}\bm{x}_{1})>\frac{1}{\sqrt{1-\sigma_{12}^{\prime
2}}}(p_{2}-\sigma_{12}^{\prime}\bm{x}_{1})$ (A.8)

for $\bm{x}_{1}<a_{1}$ and vice versa for $\bm{x}_{1}>a_{1}$. If $a_{1}\geq
p_{1}$, then for
$(\bm{x}_{1},\bm{x}_{2},\dots,\bm{x}_{I+1})\in\mathcal{T}^{\prime}(R_{2})$, we
have

$\displaystyle(\bm{x}_{1},\bm{x}_{2},\dots,\bm{x}_{I+1})\in\mathcal{T}^{\prime}(R_{2})\Rightarrow$
$\displaystyle\bm{x}_{1}\leq
p_{1},\bm{x}_{2}\leq\frac{1}{\sqrt{1-\sigma_{12}^{\prime
2}}}(p_{2}-\sigma_{12}^{\prime}\bm{x}_{1})$ $\displaystyle\Rightarrow$
$\displaystyle\bm{x}_{1}\leq
p_{1},\bm{x}_{2}\leq\frac{1}{\sqrt{1-\sigma_{12}^{2}}}(p_{2}-\sigma_{12}\bm{x}_{1})$
$\displaystyle\Rightarrow$
$\displaystyle(\bm{x}_{1},\bm{x}_{2},\dots,\bm{x}_{I+1})\in\mathcal{T}(R_{2}).$

The second line is due to (A.8) and the fact that $\bm{x}_{1}\leq p_{1}\leq
a_{1}$. It implies that
$\mathcal{T}^{\prime}(R_{2})\subset\mathcal{T}(R_{2})$. Moreover, because
(A.8) holds strictly, we can find point such as that is in
$\mathcal{T}(R_{2})$ but not in $\mathcal{T}^{\prime}(R_{2})$. Therefore,
$\mathcal{T}^{\prime}(R_{2})\subsetneq\mathcal{T}(R_{2})$ This leads to a
contradiction because

$\displaystyle\mathbb{P}((V_{1},V_{2},V_{3},\dots,V_{I+1})\in
R_{2}|\sigma_{12}^{\prime})$
$\displaystyle=\mathbb{P}((X^{\prime}_{1},X^{\prime}_{2},V_{3},\dots,V_{I+1})\in\mathcal{T}^{\prime}(R_{2}))$
$\displaystyle<\mathbb{P}((X_{1},X_{2},V_{3},\dots,V_{I+1})\in\mathcal{T}(R_{2}))=\mathbb{P}((V_{1},V_{2},V_{3},\dots,V_{I+1})\in
R_{2}|\sigma_{12}),$

which contradicts the fact that $R_{2}$ has the same probability under
$\sigma_{12}$ and $\sigma_{12}^{\prime}$. The inequality is due to the fact
that $(X_{2},X_{2})$ and $(X_{1}^{\prime},X_{2}^{\prime})$ have the same
normal distribution and the fact that
$\mathcal{T}^{\prime}(R_{2})\subsetneq\mathcal{T}(R_{2})$. If $a_{1}<p_{1}$,
we can obtain a similar contradiction by focusing on $\mathcal{T}(R_{1})$ and
$\mathcal{T}^{\prime}(R_{1})$. This shows that the probability of $R_{1}$ and
$R_{2}$ uniquely determines $\sigma_{12}$, and we have completed the inductive
step.

##  B Proofs for Section 5.2: EM Convergence

###  B.1 Preliminaries

Recall the definition of the population-level $Q$-function:

$Q(\bm{\mu}^{\prime}\mid\bm{\mu})=\sum_{k=1}^{K}\left(\int_{\bm{x}\in
R_{k}}f(\bm{x}|\bm{\mu}^{*})\mathrm{d}\bm{x}\right)\int_{\bm{v}\in
R_{k}}\frac{f(\bm{v}\mid\bm{\mu})}{\left(\int_{\bm{x}\in
R_{k}}f(x|\bm{\mu})\mathrm{d}\bm{x}\right)}\log
f(\bm{v}|\bm{\mu}^{\prime})\mathrm{d}\bm{v}.$ (B.1)

By definition $q(M(\bm{\mu}))=Q(M(\bm{\mu})\mid\bm{\mu}^{*})$, we have

$\displaystyle
q(M(\bm{\mu}))=Q(M(\bm{\mu})\mid\bm{\mu}^{*})=\sum_{k=1}^{K}\int_{\bm{v}\in
R_{k}}f(\bm{v}\mid\bm{\mu}^{*})\log
f(\bm{v}|M(\bm{\mu}))\mathrm{d}\bm{v}=\mathbb{E}_{\bm{\mu}^{*}}[\log f(v\mid
M(\bm{\mu}))]$ (B.2)

Moreover, $Q(M(\bm{\mu})\mid\bm{\mu}))$ is equal to

$\displaystyle Q(M(\bm{\mu})\mid\bm{\mu})=\sum_{k=1}^{K}\left(\int_{x\in
R_{k}}f(x|\bm{\mu}^{*})\mathrm{d}\bm{x}\right)\int_{v\in
R_{k}}\frac{f(\bm{v}\mid\bm{\mu})}{\left(\int_{x\in
R_{k}}f(x|\bm{\mu})\mathrm{d}\bm{x}\right)}\log
f(\bm{v}|M(\bm{\mu}))\mathrm{d}\bm{v}$
$\displaystyle=\sum_{k=1}^{K}\frac{\left(\int_{x\in
R_{k}}f(x|\bm{\mu}^{*})\mathrm{d}\bm{x}\right)}{\left(\int_{x\in
R_{k}}f(x|\bm{\mu})\mathrm{d}\bm{x}\right)}\int_{v\in
R_{k}}f(\bm{v}\mid\bm{\mu})\log f(\bm{v}|M(\bm{\mu}))\mathrm{d}\bm{v}$ (B.3)

To take the gradients of $q(\bm{\mu})=Q(\bm{\mu}\mid\bm{\mu}^{*})$ and
$Q(\bm{\mu}\mid\bm{\mu})$, note that

$\nabla\log f(\bm{v}|\bm{\mu})=(\bm{\Sigma}^{*})^{-1}(\bm{v}-\bm{\mu})$

by the PDF of the Gaussian distribution. Therefore,

$\displaystyle\nabla Q(M(\bm{\mu})\mid\bm{\mu})$
$\displaystyle=\sum_{k=1}^{K}\frac{\left(\int_{\bm{x}\in
R_{k}}f(\bm{x}|\bm{\mu}^{*})\mathrm{d}\bm{x}\right)}{\left(\int_{\bm{x}\in
R_{k}}f(\bm{x}|\bm{\mu})\mathrm{d}\bm{x}\right)}\int_{v\in
R_{k}}f(\bm{x}\mid\bm{\mu})(\bm{\Sigma}^{*})^{-1}(\bm{x}-M(\bm{\mu}))\mathrm{d}\bm{x}$
$\displaystyle=\sum_{k=1}^{K}\frac{\left(\int_{\bm{x}\in
R_{k}}f(\bm{x}|\bm{\mu}^{*})\mathrm{d}\bm{x}\right)}{\left(\int_{\bm{x}\in
R_{k}}f(\bm{x}|\bm{\mu})\mathrm{d}\bm{x}\right)}\left(\int_{v\in
R_{k}}f(\bm{x}\mid\bm{\mu})(\bm{\Sigma}^{*})^{-1}\bm{x}\mathrm{d}\bm{x}-\int_{\bm{v}\in
R_{k}}f(\bm{v}|\bm{\mu})\mathrm{d}\bm{v}(\bm{\Sigma}^{*})^{-1}M(\bm{\mu})\right)$
$\displaystyle=\sum_{k=1}^{K}\frac{\left(\int_{\bm{x}\in
R_{k}}f(x|\bm{\mu}^{*})\mathrm{d}\bm{x}\right)}{\left(\int_{\bm{x}\in
R_{k}}f(\bm{x}|\bm{\mu})\mathrm{d}\bm{x}\right)}\int_{v\in
R_{k}}f(\bm{x}\mid\bm{\mu})(\bm{\Sigma}^{*})^{-1}\bm{x}\mathrm{d}\bm{x}-(\bm{\Sigma}^{*})^{-1}M(\bm{\mu})$
(B.4) $\displaystyle\nabla q(\bm{\mu})$
$\displaystyle=\mathbb{E}_{\bm{\mu}^{*}}\left[(\bm{\Sigma}^{*})^{-1}\left(\bm{x}-\bm{\mu}\right)\right]=(\bm{\Sigma}^{*})^{-1}(\bm{\mu}^{*}-\bm{\mu})$
(B.5)

###  B.2 EM Convergence

###### Lemma B.1

For all $\bm{\mu}_{1},\bm{\mu}_{2}$, we have

$\displaystyle q(\bm{\mu}_{1})-q(\bm{\mu}_{2})-\left\langle\nabla
q(\bm{\mu}_{2}),\bm{\mu}_{1}-\bm{\mu}_{2}\right\rangle\leq-\frac{1}{2\lambda_{min}(\bm{\Sigma}^{*})}||\bm{\mu}_{1}-\bm{\mu}_{2}||^{2}_{2}.$

###### Proof

Proof of Lemma B.1: Expanding the left-hand side, we have

$\displaystyle q(\bm{\mu}_{1})-q(\bm{\mu}_{2})-\left\langle\nabla
q(\bm{\mu}_{2}),\bm{\mu}_{1}-\bm{\mu}_{2}\right\rangle$ $\displaystyle=$
$\displaystyle\mathbb{E}_{\bm{\mu}^{*}}\left[\log
f(\bm{x}\mid\bm{\mu}_{1})\right]-\mathbb{E}_{\bm{\mu}^{*}}\left[\log
f(\bm{x}\mid\bm{\mu}_{2})\right]-\left\langle\left(\bm{\Sigma}^{*}\right)^{-1}(\bm{\mu}^{*}-\bm{\mu}_{2}),\bm{\mu}_{1}-\bm{\mu}_{2}\right\rangle$
(B.6) $\displaystyle=$
$\displaystyle\mathbb{E}\left[-\frac{1}{2}(\bm{x}-\bm{\mu}_{1})^{\top}\left(\bm{\Sigma}^{*}\right)^{-1}(\bm{x}-\bm{\mu}_{1})\right]-\mathbb{E}\left[-\frac{1}{2}(\bm{x}-\bm{\mu}_{2})^{\top}\left(\bm{\Sigma}^{*}\right)^{-1}(\bm{x}-\bm{\mu}_{2})\right]-\left\langle\left(\bm{\Sigma}^{*}\right)^{-1}(\bm{\mu}^{*}-\bm{\mu}_{2}),\bm{\mu}_{1}-\bm{\mu}_{2}\right\rangle$
$\displaystyle=$
$\displaystyle-\frac{1}{2}(\bm{\mu}_{1}-\bm{\mu}_{2})^{\top}\left(\bm{\Sigma}^{*}\right)^{-1}(\bm{\mu}_{1}-\bm{\mu}_{2}),$

where (B.6) follows from (B.2), and (B.5). Thus, we have

$\displaystyle q(\bm{\mu}_{1})-q(\bm{\mu}_{2})-\left\langle\nabla
q(\bm{\mu}_{2}),\bm{\mu}_{1}-\bm{\mu}_{2}\right\rangle\leq\frac{\|(\bm{\Sigma}^{*})^{-1}\|_{2}}{2}\|\bm{\mu}_{1}-\bm{\mu}_{2}\|_{2},$
(B.7)

by the definition of the 2-norm of a matrix. Note that
$\|(\bm{\Sigma}^{*})^{-1}\|_{2}=\lambda_{max}((\bm{\Sigma}^{*})^{-1})=\lambda_{min}(\bm{\Sigma}^{*})^{-1}$,
where $\lambda_{max}(\cdot)$ and $\lambda_{min}(\cdot)$ are the maximum and
minimum eigenvalues of a matrix. This proves the result of Lemma B.1.

###### Lemma B.2

Suppose Assumption 1 holds. We can find a neighborhood
$\mathbb{B}_{2}(r;\bm{\mu}^{*})$ such that for
$\bm{\mu}\in\mathbb{B}_{2}(r;\bm{\mu}^{*})$ we have

$\displaystyle\|\nabla Q(M(\bm{\mu})\mid\bm{\mu})-\nabla
Q(M(\bm{\mu})\mid\bm{\mu}^{*})\|_{2}\leq\frac{(1-\epsilon/2)}{\lambda_{min}(\bm{\Sigma}^{*})}\|\bm{\mu}-\bm{\mu}^{*}\|_{2}.$
(B.8)

###### Proof

Proof of Lemma B.2. Using (B.4), we have

$\displaystyle\nabla Q(M(\bm{\mu})\mid\bm{\mu}^{*})-\nabla
Q(M(\bm{\mu})\mid\bm{\mu})$
$\displaystyle=-(\bm{\Sigma}^{*})^{-1}(\bm{\mu}^{*}-M(\bm{\mu}))+\sum_{k=1}^{K}\frac{\int_{x\in
R_{k}}f(x|\bm{\mu}^{*})\mathrm{d}\bm{x}}{\int_{\bm{x}\in
R_{k}}f(\bm{x}|\bm{\mu})\mathrm{d}\bm{x}}\int_{v\in
R_{k}}f(\bm{v}\mid\bm{\mu})(\bm{\Sigma}^{*})^{-1}\bm{v}\mathrm{d}\bm{v}-(\bm{\Sigma}^{*})^{-1}M(\bm{\mu})$
$\displaystyle=(\bm{\Sigma}^{*})^{-1}\left(\sum_{k=1}^{K}\int_{\bm{x}\in
R_{k}}f(\bm{x}|\bm{\mu}^{*})\mathrm{d}\bm{x}\frac{\int_{\bm{x}\in
R_{k}}\bm{x}f(\bm{x}|\bm{\mu})\mathrm{d}\bm{x}}{\int_{\bm{x}\in
R_{k}}f(\bm{x}|\bm{\mu})\mathrm{d}\bm{x}}-\bm{\mu}^{*}\right)$ (B.9)

We next use the Taylor expansion on the above quantity. To do so, we first
derive the gradient of
$g_{k}(\bm{\mu})\triangleq\frac{\int_{R_{k}}\bm{x}f(\bm{x}\mid\bm{\mu})\mathrm{d}\bm{x}}{\int_{R_{k}}\bm{f}(\bm{x}\mid\bm{\mu})\mathrm{d}\bm{x}}$.

$\displaystyle\nabla g_{k}(\bm{\mu})=$ $\displaystyle\frac{\int_{\bm{x}\in
R_{k}}\bm{x}(\nabla_{\bm{\mu}}f(\bm{x}|\bm{\mu}))^{\top}\mathrm{d}\bm{x}\int_{\bm{x}\in
R_{k}}f(\bm{x}|\bm{\mu})\mathrm{d}\bm{x}-\int_{\bm{x}\in
R_{k}}\bm{x}f(\bm{x}|\bm{\mu})\mathrm{d}\bm{x}\int_{\bm{x}\in
R_{k}}(\nabla_{\bm{\mu}}f(\bm{x}|\bm{\mu}))^{\top}\mathrm{d}\bm{x}}{(\int_{\bm{x}\in
R_{k}}f(\bm{x}|\bm{\mu})\mathrm{d}\bm{x})^{2}}$ $\displaystyle=$
$\displaystyle\frac{\left(\int_{R_{k}}\bm{x}(\bm{x}-\bm{\mu})^{\top}(\bm{\Sigma}^{*})^{-1}f(\bm{x}|\bm{\mu})\mathrm{d}\bm{x}\right)\left(\int_{R_{k}}f(\bm{x}|\bm{\mu})\mathrm{d}\bm{x}\right)}{(\int_{R_{k}}f(\bm{x}|\bm{\mu})\mathrm{d}\bm{x})^{2}}$
$\displaystyle\quad-\frac{\left(\int_{R_{k}}\bm{x}f(\bm{x}|\bm{\mu})\mathrm{d}\bm{x}\right)\left(\int_{R_{k}}(\bm{x}-\bm{\mu})^{\top}(\bm{\Sigma}^{*})^{-1}f(\bm{x}|\bm{\mu})\mathrm{d}\bm{x}\right)}{(\int_{R_{k}}f(\bm{x}|\bm{\mu})\mathrm{d}\bm{x})^{2}}$
$\displaystyle=$
$\displaystyle\frac{\int_{R_{k}}(\bm{x}-\bm{\mu})(\bm{x}-\bm{\mu})^{\top}(\bm{\Sigma}^{*})^{-1}f(\bm{x}|\bm{\mu})\mathrm{d}\bm{x}+\int_{R_{k}}\bm{\mu}(\bm{x}-\bm{\mu})^{\top}\bm{(}\Sigma^{*})^{-1}f(\bm{x}|\bm{\mu})\mathrm{d}\bm{x}}{\int_{R_{k}}f(\bm{x}|\bm{\mu})\mathrm{d}\bm{x}}$
$\displaystyle\quad-\frac{\left(\int_{R_{k}}\bm{x}f(\bm{x}|\bm{\mu})\mathrm{d}\bm{x}\right)\left(\int_{R_{k}}(\bm{x}-\bm{\mu})^{\top}(\bm{\Sigma}^{*})^{-1}f(\bm{x}|\bm{\mu})\mathrm{d}\bm{x}\right)}{(\int_{R_{k}}f(\bm{x}|\bm{\mu})\mathrm{d}\bm{x})^{2}}$
$\displaystyle=$
$\displaystyle\left(\mathbb{E}_{\bm{\mu}}\left[(\bm{X}-\bm{\mu})(\bm{X}-\bm{\mu})^{\top}\mid
R_{k}\right]-\frac{\left(\int_{R_{k}}(\bm{x}-\bm{\mu})f(\bm{x}|\bm{\mu})\mathrm{d}\bm{x}\right)\left(\int_{R_{k}}(\bm{x}-\bm{\mu})^{\top}f(\bm{x}|\bm{\mu})\mathrm{d}\bm{x}\right)}{(\int_{R_{k}}f(\bm{x}|\bm{\mu})\mathrm{d}\bm{x})^{2}}\right)(\bm{\Sigma}^{*})^{-1}$
$\displaystyle=$
$\displaystyle\left(\mathbb{E}_{\bm{\mu}}[(\bm{X}-\bm{\mu})(\bm{X}-\bm{\mu})^{\top}\mid
R_{k}]-\mathbb{E}_{\bm{\mu}}[\bm{X}-\bm{\mu}\mid
R_{k}]\mathbb{E}_{\bm{\mu}}[\bm{X}-\bm{\mu}\mid
R_{k}]^{\top}\right)(\bm{\Sigma}^{*})^{-1}$ $\displaystyle=$
$\displaystyle\left(\mathbb{E}_{\bm{\mu}}\left[\bm{X}\bm{X}^{\top}\mid
R_{k}\right]-\mathbb{E}_{\bm{\mu}}\left[\bm{X}\mid
R_{k}\right]\mathbb{E}_{\bm{\mu}}\left[\bm{X}\mid
R_{k}\right]^{\top}\right)(\bm{\Sigma}^{*})^{-1}$ $\displaystyle=$
$\displaystyle\mathrm{Var}_{\bm{\mu}}(R_{k})(\bm{\Sigma}^{*})^{-1}.$

In the last line, we use $\mathrm{Var}(R_{k}|\bm{\mu})$ to denote the
covariance matrix of the conditional distribution of
$\mathcal{N}(\bm{\mu},\bm{\Sigma}^{*})$ in the region $R_{k}$. Therefore, by
the Taylor expansion, we have

$\displaystyle\nabla Q(M(\bm{\mu})\mid\bm{\mu})-\nabla
Q(M(\bm{\mu})\mid\bm{\mu}^{*})$
$\displaystyle=(\bm{\Sigma}^{*})^{-1}\sum_{k=1}^{K}\int_{\bm{x}\in
R_{k}}f(\bm{x}|\bm{\mu}^{*})\mathrm{d}\bm{x}\nabla
g_{k}(\bm{\mu}^{*})(\bm{\mu}-\bm{\mu}^{*})+O(\|\bm{\mu}^{*}-\bm{\mu}\|_{2}^{2})$
$\displaystyle=(\bm{\Sigma}^{*})^{-1}\sum_{k=1}^{K}\int_{\bm{x}\in
R_{k}}f(\bm{x}|\bm{\mu}^{*})\mathrm{d}\bm{x}\mathrm{Var}_{\bm{\mu}^{*}}(R_{k})(\bm{\Sigma}^{*})^{-1}(\bm{\mu}-\bm{\mu}^{*})+O(\|\bm{\mu}^{*}-\bm{\mu}\|_{2}^{2}).$
(B.10)

Note that if we interpret $R$ as a categorical random variable, with
$\mathbb{P}(R=R_{k})=\sum_{k=1}^{K}\int_{\bm{x}\in
R_{k}}f(\bm{x}|\bm{\mu}^{*})\mathrm{d}\bm{x}$, then by the law of total
covariance we have

$\sum_{k=1}^{K}\int_{\bm{x}\in
R_{k}}f(\bm{x}|\bm{\mu}^{*})\mathrm{d}\bm{x}\mathrm{Var}(R_{k})=\mathbb{E}[\mathrm{Var}(R)]=\mathrm{Var}(\bm{X})-\mathrm{Var}(\mathbb{E}[\bm{X}|R])=\bm{\Sigma}^{*}-\mathrm{Var}(\mathbb{E}[\bm{X}|R]),$

where $\bm{X}\sim\mathcal{N}(\bm{\mu}^{*},\bm{\Sigma}^{*})$. Consider the
following transformation
$\bm{Z}=(\bm{\Sigma}^{*})^{-1/2}(\bm{X}-\bm{\mu}^{*})$. It is well known that
$\bm{Z}$ has a standard normal distribution whose covariance matrix is equal
to the identity matrix $\bm{I}$. Therefore, (B.10) can be further expanded to

$\displaystyle\eqref{eq:taylor}$
$\displaystyle=(\bm{\Sigma}^{*})^{-1}(\bm{\Sigma}^{*}-\mathrm{Var}(\mathbb{E}[\bm{X}|R]))(\bm{\Sigma}^{*})^{-1}(\bm{\mu}-\bm{\mu}^{*})+O(\|\bm{\mu}^{*}-\bm{\mu}\|_{2}^{2})$
$\displaystyle=(\bm{\Sigma}^{*})^{-1}(\bm{\Sigma}^{*}-(\bm{\Sigma}^{*})^{1/2}\mathrm{Var}(\mathbb{E}[\bm{Z}|R^{\prime}])(\bm{\Sigma}^{*})^{1/2})(\bm{\Sigma}^{*})^{-1}(\bm{\mu}-\bm{\mu}^{*})+O(\|\bm{\mu}^{*}-\bm{\mu}\|_{2}^{2})$
$\displaystyle=(\bm{\Sigma}^{*})^{-1/2}(\bm{I}-\mathrm{Var}(\mathbb{E}[\bm{Z}|R^{\prime}]))(\bm{\Sigma}^{*})^{-1/2}(\bm{\mu}-\bm{\mu}^{*})+O(\|\bm{\mu}^{*}-\bm{\mu}\|_{2}^{2}),$

where $R^{\prime}$ represents the transformed regions in the $\bm{Z}$-space.

Suppose Assumption 1 holds. In this case, we know that

$\lambda_{max}(\bm{I}-\mathrm{Var}(\mathbb{E}[\bm{Z}|R^{\prime}]))\leq
1-\epsilon.$

As a result,

$\|(\bm{\Sigma}^{*})^{-1/2}(I-\mathrm{Var}(\mathbb{E}[\bm{Z}|R^{\prime}]))(\bm{\Sigma}^{*})^{-1/2}\|_{2}\leq(1-\epsilon)\lambda_{max}((\bm{\Sigma}^{*})^{-1})=(1-\epsilon)\lambda_{min}(\bm{\Sigma}^{*})^{-1}$

This implies that there exists a small neighborhood
$\mathbb{B}_{2}(r;\bm{\mu}^{*})$, such that for
$\bm{\mu}\in\mathbb{B}_{2}(r;\bm{\mu}^{*})$, we have

$\displaystyle\|\nabla Q(M(\bm{\mu})\mid\bm{\mu})-\nabla
Q(M(\bm{\mu})\mid\bm{\mu}^{*})\|_{2}$
$\displaystyle\leq\|(\bm{\Sigma}^{*})^{-1/2}(I-\mathrm{Var}(\mathbb{E}[\bm{Z}|R^{\prime}]))(\bm{\Sigma}^{*})^{-1/2}\|_{2}\|\bm{\mu}-\bm{\mu}^{*}\|_{2}+o(\|\bm{\mu}-\bm{\mu}^{*}\|_{2})$
$\displaystyle\leq(1-\epsilon)\lambda_{min}((\bm{\Sigma}^{*}))^{-1}\|\bm{\mu}-\bm{\mu}^{*}\|_{2}+o(\|\bm{\mu}-\bm{\mu}^{*}\|_{2}).$

For a sufficiently small $r$, we can guarantee that (B.8) holds.

###### Proof

Proof of Theorem 1. The desired result is a direct implication of
[Balakrishnan et al., 2017, Theorem 4], which requires Lemmas B.1 and B.2,
respectively. We set
$\lambda=\frac{1}{\min(\operatorname{eig}(\bm{\Sigma}^{*}))}$ and
$\gamma=\frac{1-\epsilon/2}{\min(\operatorname{eig}(\bm{\Sigma}^{*}))}$, where
$\gamma<\lambda$ for any $0<\epsilon\leq 1$. Therefore, using [Balakrishnan et
al., 2017, Theorem 4], we have

$\displaystyle\|M(\bm{\mu})-\bm{\mu}^{*}\|_{2}\leq(1-\epsilon/2)\|\bm{\mu}-\bm{\mu}^{*}\|_{2},\quad\text{for
all }\bm{\mu}\in\mathbb{B}_{2}(r;\bm{\mu}^{*}).$ (B.11)

Then, we immediately derive that for any initial point
$\bm{\mu}^{(0)}\in\mathbb{B}_{2}(r;\bm{\mu}^{*})$, the population EM sequence
$\left\\{\bm{\mu}^{(t)}\right\\}_{t=0}^{\infty}$ exhibit linear convergence,
and establish the desired result.

###  B.3 Hierarchal Bayesian model

In this section, we first present the Metropolis Hasting algorithm to estimate
the posterior distribution. Then we optimize this method and use a Bayesian
grid search algorithm which performs better in terms of speed and accuracy.

Algorithm 3 Metropolis Hasting

1:Input: $R_{n}^{c_{n}}$, $\bm{\mu}^{(t)}$, $\bm{\Sigma}^{(t)}$

2:while $error\leq\epsilon$ do

3: simulate $\bm{\mu}^{(t)}=\bm{\mu}^{(t)}+\epsilon$, $\bm{\Sigma}^{(t)}$

4: $\epsilon\sim unif(-0.5,0.5)$

5: Calculate
$\mathbb{P}\left(R_{n}^{c_{n}}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)\times\mathbb{P}\left(\bm{\mu}^{(t)}\right)\times\mathbb{P}\left(\bm{\Sigma}^{(t)}\right)$

6: Calculate
$\mathbb{P}\left(R_{n}^{c_{n}}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)\times\mathbb{P}\left(\bm{\mu}^{(t)}\right)\times\mathbb{P}\left(\bm{\sigma}^{(t)}\right)$

7:
$A=\frac{\mathbb{P}\left(R_{n}^{c_{n}}|\bm{\mu}^{(0)},\bm{\Sigma}^{(0)}\right)\times\mathbb{P}\left(\bm{\mu}^{(0)}\right)\times\mathbb{P}\left(\bm{\Sigma}^{(0)}\right)}{\mathbb{P}\left(R_{n}^{c_{n}}|\bm{\mu}^{(t)},\bm{\sigma}^{(t)}\right)\times\mathbb{P}\left(\bm{\mu}^{(t)}\right)\times\mathbb{P}\left(\bm{\sigma}^{(t)}\right)}$

8: $a\sim Unif(1)$

9: if $a\leq A$ then

10: $\bm{\mu}^{(t)}=\bm{\mu}^{(t+1)}$, $\bm{\Sigma}^{(t)}=\bm{\Sigma}^{(t+1)}$

11: else

12: $\bm{\mu}^{(t)}=\bm{\mu}^{(t)}$, $\bm{\Sigma}^{(t)}=\bm{\Sigma}^{(t)}$

13: end if

14:
$error\leftarrow\|\bm{\mu}^{(t+1)}-\bm{\mu}^{(t)}\|_{1}+\|\bm{\Sigma}^{(t+1)}-\bm{\Sigma}^{(t)}\|_{1}$,
$t\leftarrow t+1$

15:end while

16:return $\bm{\mu}^{(t)}$, $\bm{\Sigma}^{(t)}$

Algorithm 4 Bayesian Grid Search

1:Input: $R_{n}^{c_{n}}$, $prior(\bm{\mu}^{(t)})$,
$prior\left(\bm{\Sigma}^{(t)}\right)$

2:while $error\leq\epsilon$ do

3: Simulate $\bm{\mu}^{(t)}$, $\bm{\Sigma}^{(t)}$

4: Calculate
$\mathbb{P}\left(R_{n}^{c_{n}}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)$

5: Calculate
$\mathbb{P}\left(R_{n}^{c_{n}}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)\times\mathbb{P}\left(\bm{\mu}^{(t)}\right)\times\mathbb{P}\left(\bm{\Sigma}^{(t)}\right)$

6: Calculate
$\mathbb{P}\left(\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}|R_{n}^{c_{n}}\right)\propto\mathbb{P}\left(R_{n}^{c_{n}}|\bm{\mu}^{(t)},\bm{\Sigma}^{(t)}\right)\times\mathbb{P}\left(\bm{\mu}^{(t)}\right)\times\mathbb{P}(\bm{\Sigma}^{(t)})$

7:
$error\leftarrow\|\bm{\mu}^{(t+1)}-\bm{\mu}^{(t)}\|_{1}+\|\bm{\Sigma}^{(t+1)}-\bm{\Sigma}^{(t)}\|_{1}$,
$t\leftarrow t+1$

8:end while

9:search the maximum posterior return maximum posterior $\bm{\mu}$ and
$\bm{\Sigma}$