1 Introduction↩︎

Bundles are pervasive in online retailing, fast-food restaurants, insurance, airlines, and the video game industry. Firms often combine products or services into discounted packages to stimulate demand and increase transaction value. Examples range from McDonald’s combo meals and bundled home-and-auto insurance policies to airline packages that combine fares, baggage allowances, and meals. Due to the popularity of bundles, how to set profitable bundle pricing strategies has been a focus of the literature since [1]. 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 [2] for a comprehensive empirical comparison.

A prerequisite for successful bundle pricing is an accurate estimate of consumers’ reservation prices, or valuations. These valuations determine how consumers respond to menus of individual products and bundles, and hence determine aggregate demand. A natural source for learning such information is the firm’s historical transaction data. 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 [3], companies make more than 90% of the effort estimating demand and 10% setting optimal prices.

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 valuations, censored demand, and the bundled products. We summarize each of these challenges and explain why existing approaches are insufficient. First, customers’ valuations for products are only partially observed in transaction data. While purchase decisions are recorded, the underlying willingness-to-pay (WTP) is not directly observed. Instead, each observed purchase provides only a lower bound on the customer’s valuation: if a customer buys a product at price $5, the firm can infer that the customer’s valuation is at least $5, but the exact valuation remains unknown. This leads to a form of interval-censored or partially identified data, which complicates the estimation of the valuation distribution. Second, the data suffer from demand censoring, as customers who do not make any purchase are typically not observed. As a result, the dataset over-represents purchasing customers and omits those with low valuations, introducing selection bias into the estimation problem. The first two challenges are common in estimation problems based on transaction data. In settings without bundles, discrete choice models such as the multinomial logit (MNL) model address these issues by linking observed purchase probabilities to model parameters, thereby bypassing the need to observe valuations directly. However, the presence of bundles introduces an additional layer of difficulty. A common simplification in the literature is to treat each bundle as an independent product and apply standard discrete choice models. While this approach enables tractable estimation, it ignores the structural relationship between bundles and their component products. This simplification is often adopted precisely because bundles introduce additional complexity: when products are sold jointly, their valuations become entangled, and it is no longer clear how to attribute the observed purchase decision to the valuation of individual products. As a result, the firm cannot separately identify product-level valuations using standard methods. These challenges highlight the limitations of existing approaches and motivate the need for a framework that can jointly handle partial observability, censoring, and the structural dependence induced by bundle offerings.

In this paper, we develop a parsimonious yet expressive model of consumer valuations and propose an estimation procedure based on transaction data. In particular, we assume that consumers’ valuations for products follow a multivariate Gaussian distribution, and that the valuation of a bundle is given by the sum of the valuations of its component products. This additive utility specification is standard in the bundle pricing literature and provides a tractable foundation for analysis. Building on this formulation, we cast the estimation problem as a maximum likelihood problem with incomplete data, where individual-level valuation realizations are unobserved. To address this challenge, we design an estimation procedure that integrates the expectation–maximization (EM) algorithm with Monte Carlo simulation, enabling the recovery of the mean vector and covariance matrix of the underlying valuation distribution. The contribution of the paper is threefold.

• We reformulate bundle transaction data as a missing-data estimation problem over IC polyhedra. Each purchase observation does not reveal the customer’s valuation vector; instead, it identifies a polyhedral region, defined by linear IC inequalities, in which that vector must lie. A special case of this problem, interval-censored data, is commonly seen in survival analysis [4]. Our formulation generalizes the interval to polyhedra in high dimensions.

• The framework and algorithm we propose are flexible enough to accommodate many practical considerations. We extend the framework to non-additive product valuations in bundles, in order to capture complementarity and substitutability. The model also allows for Gaussian mixture models (GMM) and thus clustered valuations due to market segmentation. Our framework can tackle censored demand. This is because Monte Carlo simulation allows us to provide closed-form updates in the M-step. No known algorithms can handle the estimation problem for bundle sales of this scale.

• We study two theoretical questions: when the model is identifiable, and under what conditions the EM 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.

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. [5], [6] analyze the asymptotic performance of pure bundling and bundle size pricing. [2] conduct a comprehensive numerical experiment to compare the performance of four commonly studied bundle pricing mechanisms. [7]–[10] provide mixed-integer programming or convex optimization formulations to compute the optimal bundle pricing policies. [11], [12] derive theoretical performance bounds for simple mechanisms. A few recent studies [13]–[16] 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 [3], 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. [17] provide a comprehensive literature review of the earlier studies on such models in revenue management. Building on [18], [19] 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. [20] incorporate price and product features into their estimation approach. [21] study the identifiability of this type of model and propose a computational approach to optimize the likelihood function. [22] 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 [23], [24] and the Markov chain choice model [25]. Our contribution relative to this stream is to model transaction data in which customers may choose discounted bundles whose utilities are linked through shared component products. This setting is not directly captured by standard discrete-choice specifications that treat alternatives as unrelated products. We therefore develop an estimation framework for product-level valuations under bundle menus.

This study proposes a framework to estimate customer valuations from transaction data for multiple products when the firm may offer a price menu for products and bundles.

In another related paper, [26] 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. [27] 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 [26]–[28].

The methodology used in this paper, the EM algorithm combined with Monte Carlo simulation, has been used in previous studies, including [29], [30]. In particular, Monte Carlo simulation is used to approximate the E-step when the expectation does not 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 specific application.

3 Bundle Sales Data 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 slight 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 choices of \(N\) customers: customer \(n\in [N]\triangleq \left\{1,\dots,N\right\}\) chooses alternative \(c_n\in [J]\cup\left\{0\right\}\), where \(c_n=0\) denotes no purchase. The case of censored demand is studied in Section 5.2. 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\). We do not require every bundle in the menu to be offered to every customer. When a customer is only offered a subset of the available bundles, we can artificially set the prices of 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 \(n\) is endowed with a vector of product valuations, \(\boldsymbol{v}_n=(v_{n1},\dots,v_{nI})\). Under additive utility, the valuation of bundle \(j\) for customer \(n\) is \(\sum_{i\in j} v_{ni}\). Note that additive utility is a standard assumption in the literature (e.g., see [2] and references therein), although other types of valuation functions have been proposed and studied recently [31]. Under additive utility, statistical dependence across products can be captured through correlations in product valuations. In Section 5.1, we introduce a more general non-additive model with product synergy, which directly captures complementarity and substitutability within bundles.

Under this utility model, the choice \(c_n=j\) for \(j\neq 0\) corresponds to the following region of the valuation vector \(\boldsymbol{v}\), given by the incentive-compatible constraints: \[\label{eq:cn61j} c_n=j\iff \left\{\sum_{i\in j}v_i-p^{j}_n\ge \max_{j'\in [J]}\sum_{i\in j'}v_i-p^{j'}_n ,\;\sum_{i\in j}v_i-p^{j}_n\ge 0\right\}.\tag{1}\] Similarly, we have \[\label{eq:cn610} c_n=0 \iff\left\{\max_{j'\in [J]}\sum_{i\in j'}v_i-p^{j'}_n\le 0 \right\}.\tag{2}\] Suppose the random valuation \(\boldsymbol{v}\sim \boldsymbol{V}\) of each customer is independently drawn from an \(I\)-dimensional multivariate normal distribution: \[\label{eq:pdf-normal} \boldsymbol{V}\sim \mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma}).\tag{3}\] Note that we focus on the somewhat restrictive normal distribution for the purpose of exposition. In Appendix 9, we relax it to a Gaussian mixture model, which is flexible enough to approximate any probability distributions (for example, see [32, p. 65]).

Given the dataset \(\{(p^1_n,\dots,p^J_n,c_n)\}_{n=1}^N\), the firm’s objective is to estimate the distribution of customer valuations, in particular \(\boldsymbol{\mu}\) and \(\boldsymbol{\Sigma}\), with the knowledge that the incentive-compatible constraints 1 and 2 relate the valuations to the observations in the dataset. Note that because the parameters \(\boldsymbol{\mu}\) and \(\boldsymbol{\Sigma}\) are linked to the observed data only indirectly through the IC constraints, this estimation problem is not standard in the statistics literature. For example, if \(\boldsymbol{v}_n\) were observed, then the estimation of \(\boldsymbol{\mu}\) and \(\boldsymbol{\Sigma}\) would be trivial using the empirical average and covariance matrix of \(\{\boldsymbol{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 \(\boldsymbol{v}_n\) falls in. More precisely, define \[\begin{align} \label{eq:ic-polytope} R^j_n&\triangleq \left\{\boldsymbol{v}: \sum_{i\in j}v_i-p^{j}_n\ge \max_{j'\in [J]}\sum_{i\in j'}v_i-p^{j'}_n ,\;\sum_{i\in j}v_i-p^{j}_n\ge 0\right\}, \quad j\neq 0,\notag\\ R^0_n&\triangleq \left\{\boldsymbol{v}: \max_{j'\in [J]}\sum_{i\in j'}v_i-p^{j'}_n\le 0 \right\}, \end{align}\tag{4}\] which we refer to as IC polyhedra. The observation \(c_n=j\) is equivalent to \(\boldsymbol{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 \(\boldsymbol{v}_n\) falls. Therefore, instead of observing \(\boldsymbol{v}_n\) exactly, the firm observes a polyhedron \(R_n^{c_n}\) for all potential \(\boldsymbol{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 5.2. With this interpretation, we view the estimation problem through the lens of inference with missing data.

Parameters:

we are interested in estimating \((\boldsymbol{\mu}, \boldsymbol{\Sigma})\), the distribution of the valuations for the products among consumers. By convention, we use \(\boldsymbol{\theta}\triangleq (\boldsymbol{\mu}, \boldsymbol{\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 4 . We use \(\boldsymbol{D} \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\{\boldsymbol{v}_n\right\}_{n=1}^N\).

Missing data:

we treat valuations \(\boldsymbol{v}_n\), \(n=1,\dots,N\), as the missing data. They are related to the observed data by the simple fact that \(\boldsymbol{v}_n\in R_n^{c_n}\). We use \(\boldsymbol{Z}\triangleq\{\boldsymbol{v}_n\}_{n=1}^N\) to denote the missing data.

There are two common approaches to missing data: imputation and maximum likelihood estimation (MLE). Since our estimation problem is model-based, we adopt the latter. In particular, the likelihood function based on the observed data is: \(\mathcal{L}(\boldsymbol{\theta};\boldsymbol{D}) = \int p(\boldsymbol{D}, \boldsymbol{Z}|\boldsymbol{\theta})\mathrm{d}\boldsymbol{Z} = \prod_{n=1}^N\int_{R^{c_n}_n} f(\boldsymbol{v}|\boldsymbol{\mu},\boldsymbol{\Sigma}) \mathrm{d}\boldsymbol{v}\), where \(f(\boldsymbol{v}|\boldsymbol{\mu},\boldsymbol{\Sigma}) = \tfrac{1}{\sqrt{(2\pi)^I|\boldsymbol{\Sigma}|}}\exp\left(-\frac{1}{2}(\boldsymbol{v}-\boldsymbol{\mu})^\top \boldsymbol{\Sigma}^{-1}(\boldsymbol{v}-\boldsymbol{\mu})\right)\) is the probability density function (PDF) of the multivariate normal distribution 3 .

Taking logarithms improves numerical scaling but does not remove the main difficulty: the objective still contains Gaussian integrals over IC polyhedra, which do not have closed-form expressions in general. 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., \[\label{eq:complete-data-llhx} \ell (\boldsymbol{\theta};\boldsymbol{D}, \boldsymbol{Z}) = \sum_{n=1}^{N} \log f(\boldsymbol{v}_n|\boldsymbol{\mu},\boldsymbol{\Sigma}).\tag{5}\] Note that as an objective function, \(\ell(\boldsymbol{\theta};\boldsymbol{D}, \boldsymbol{Z})\) is easier to handle than \(\ell(\boldsymbol{\theta};\boldsymbol{D})\). Given complete valuation data, the Gaussian log-likelihood has standard closed-form maximizers for \(\boldsymbol{\mu}\) and \(\boldsymbol{\Sigma}\). However, \(\boldsymbol{Z}\) is not observed in the transaction data.

The EM algorithm circumvents the issue by estimating \(\boldsymbol{\theta}^{(t)}\) iteratively. In each iteration, suppose \(\boldsymbol{\theta}^{(t)}=\left(\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)\) is given. We first evaluate the expected complete-data log-likelihood over \(\boldsymbol{Z}\), i.e., \[\begin{align} \label{eq:Q-func} Q\left(\boldsymbol{\theta}|\boldsymbol{\theta}^{(t)}\right)&\triangleq \mathbb{E}\left[\ell(\boldsymbol{\theta};\boldsymbol{D},\boldsymbol{Z})\bigg|\boldsymbol{D},\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right]\notag\\ &= \sum_{n=1}^{N} \int_{\boldsymbol{v}\in R_n^{c_n}} \frac{f\left(\boldsymbol{v}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)}{\int_{R_n^{c_n}}f\left(\boldsymbol{\xi}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)\mathrm{d}\boldsymbol{\xi} }\log f(\boldsymbol{v}|\boldsymbol{\mu},\boldsymbol{\Sigma})\mathrm{d}\boldsymbol{v}. \end{align}\tag{6}\] This is the so-called E-step and the expectation is taken with respect to the missing data \(\boldsymbol{Z}\). Given \(\left(\boldsymbol{D},\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)\), the missing variable \(\boldsymbol{v}_n\) follows a multivariate normal distribution with mean \(\boldsymbol{\mu}^{(t)}\) and covariance \(\boldsymbol{\Sigma}^{(t)}\), truncated to the polyhedron \(R_n^{c_n}\). Its PDF is thus \(f\left(\boldsymbol{v}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)/\int_{R_n^{c_n}}f\left(\boldsymbol{\xi}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)\mathrm{d}\boldsymbol{\xi}\) for \(\boldsymbol{v}\in R_n^{c_n}\), which is the first term in the integral. The second term \(\log f(\boldsymbol{v}|\boldsymbol{\mu},\boldsymbol{\Sigma})\) is the log-likelihood function \(\ell(\boldsymbol{\theta};\boldsymbol{D},\boldsymbol{Z})\) and is integrated with the conditional PDF for its expected value.

To update \(\boldsymbol{\theta}\) in the \(t\)-th iteration, the M-step maximizes \(Q\left(\boldsymbol{\theta}|\boldsymbol{\theta}^{(t)}\right)\) and sets \(\boldsymbol{\theta}^{(t+1)}=\argmax_{\boldsymbol{\theta}}Q\left(\boldsymbol{\theta}|\boldsymbol{\theta}^{(t)}\right)\). The main difficulty is that the truncated-normal expectations in 6 do not have closed forms. To deal with this challenge, we use Monte Carlo simulation. More precisely, for all \(n=1,\dots,N\), we generate \(L\) samples according to the same distribution as \(\boldsymbol{v}|R_n^{c_n},\boldsymbol{\mu}^{(t)},\boldsymbol{\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.

Once the Monte Carlo samples \(\{\boldsymbol{v}_{n}^{(l)}\}\) for \(n=1,\dots,N\) and \(l=1,\dots,L\) have been generated, we may replace \(Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(t)})\) in 6 by \[\begin{align} \label{eq:approx-Q} \hat{Q}(\boldsymbol{\theta}|\boldsymbol{\theta}^{(t)})&= \frac{1}{L}\sum_{n=1}^{N} \sum_{l=1}^L \log f(\boldsymbol{v}_{n}^{(l)}|\boldsymbol{\mu},\boldsymbol{\Sigma})\notag\\ &= C+ \frac{1}{L}\sum_{n=1}^N\sum_{l=1}^L \left(- \frac{1}{2}\left(\boldsymbol{v}_n^{(l)}-\boldsymbol{\mu}\right)^\top \boldsymbol{\Sigma}^{-1}\left(\boldsymbol{v}_{n}^{(l)}-\boldsymbol{\mu}\right) - \frac{1}{2}\log |\boldsymbol{\Sigma}| \right), \end{align}\tag{7}\] where \(C\) is constant with respect to \(\boldsymbol{\mu}\) and \(\boldsymbol{\Sigma}\). Note that this is the standard MLE for multivariate normal distribution with \(NL\) samples. Therefore, we have \[\begin{align} \boldsymbol{\mu}^{(t+1)} &= \frac{1}{NL}\sum_{n=1}^N\sum_{l=1}^L \boldsymbol{v}_{n}^{(l)},\tag{8}\\ \quad\boldsymbol{\Sigma}^{(t+1)} &= \frac{1}{NL} \sum_{n=1}^N\sum_{l=1}^L \left(\boldsymbol{v}_n^{(l)}-\boldsymbol{\mu}^{(t+1)}\right)\left(\boldsymbol{v}_n^{(l)}-\boldsymbol{\mu}^{(t+1)}\right)^{\top}\tag{9}. \end{align}\] This completes the \(t\)-th iteration of the EM algorithm.

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

Sampling efficiency of Monte Carlo simulation. In practice, the major computational complexity of Algorithm 2 is caused by the Monte Carlo simulation. when the region \(R_n^{c_n}\) is small or distant from \(\boldsymbol{\mu}^{(t)}\), Algorithm 1 may be inefficient as it requires a large number of samples in Step [step:acceptance] to draw an accepted sample in \(R_n^{c_n}\). We provide three remedies to improve the sampling efficiency.

First, in each iteration of Algorithm 2 (Step [step:em-iteration]), we can generate a large number of samples from \(\mathcal{N}\left(\boldsymbol{\mu}^{(t)},\boldsymbol{\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 generate \(L_n\) Monte Carlo samples for region \(R_n^{c_n}\), choosing \(L_n\) adaptively based on the acceptance rate or the empirical variance of the sampled log-likelihood terms. As long as we use \(L_n^{-1}\sum_{l=1}^{L_n} \log f(\boldsymbol{v}_{n}^{(l)}|\boldsymbol{\mu},\boldsymbol{\Sigma})\) in \(\hat{Q}(\boldsymbol{\theta}|\boldsymbol{\theta}^{(t)})\), the Monte Carlo approximation remains unbiased for the corresponding conditional expectation. However, small values of \(L_n\) can increase Monte Carlo noise and may slow or destabilize convergence.

Third, for low-probability \(R_n^{c_n}\), we may use importance sampling to increase the sampling efficiency. In particular, consider 6 and some \(n\) that \(f\left(\boldsymbol{v}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)\approx 0\) for \(\boldsymbol{v}\in R_n^{c_n}\). In this case, it is computationally challenging as Step [step:acceptance] is repeated many times before acceptance. Instead, we can consider a different normal distribution \(\mathcal{N}(\boldsymbol{\mu}',\boldsymbol{\Sigma}')\) such that \(\boldsymbol{v}\in R_n^{c_n}\) with high probability. This can be achieved, for example, by choosing \(\boldsymbol{\mu}'\) close to the center of \(R_n^{c_n}\) when \(R_n^{c_n}\) is bounded. As a result, we have \[\begin{align} &\int_{\boldsymbol{v}\in R_n^{c_n}} \frac{f\left(\boldsymbol{v}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)}{\int_{R_n^{c_n}}f\left(\boldsymbol{\xi}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)\mathrm{d}\boldsymbol{\xi} }\log f(\boldsymbol{v}|\boldsymbol{\mu},\boldsymbol{\Sigma})\mathrm{d}\boldsymbol{v}\\ =& \int_{\boldsymbol{v}\in R_n^{c_n}} \frac{f(\boldsymbol{v}|\boldsymbol{\mu}',\boldsymbol{\Sigma}')}{\int_{R_n^{c_n}}f(\boldsymbol{\xi}|\boldsymbol{\mu}',\boldsymbol{\Sigma}')\mathrm{d}\boldsymbol{\xi} }\frac{f\left(\boldsymbol{v}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)\int_{R_n^{c_n}}f(\boldsymbol{\xi}|\boldsymbol{\mu}',\boldsymbol{\Sigma}')\mathrm{d}\boldsymbol{\xi} }{f(\boldsymbol{v}|\boldsymbol{\mu}',\boldsymbol{\Sigma}')\int_{R_n^{c_n}}f\left(\boldsymbol{\xi}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)\mathrm{d}\boldsymbol{\xi} }\log f(\boldsymbol{v}|\boldsymbol{\mu},\boldsymbol{\Sigma})\mathrm{d}\boldsymbol{v}. \end{align}\]

By simulating \(L\) samples from a truncated normal distribution \(\mathcal{N}(\boldsymbol{\mu}',\boldsymbol{\Sigma}')\) on region \(R_n^{c_n}\) (run Algorithm 1 with the new distribution), denoted as \(\boldsymbol{v}_{n}^{(1)},\dots,\boldsymbol{v}_n^{(L)}\), we can approximate the above term by \[\begin{align} \frac{1}{\frac{\int_{R_n^{c_n}}f\left(\boldsymbol{\xi}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)\mathrm{d}\boldsymbol{\xi} }{\int_{R_n^{c_n}}f\left(\boldsymbol{\xi}|\boldsymbol{\mu}',\boldsymbol{\Sigma}'\right)\mathrm{d}\boldsymbol{\xi} }} &\frac{1}{L} \sum_{l=1}^L\frac{f\left(\boldsymbol{v}_n^{(l)}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)}{f\left(\boldsymbol{v}_n^{(l)}|\boldsymbol{\mu}',\boldsymbol{\Sigma}'\right)} \log f\left(\boldsymbol{v}_n^{(l)}|\boldsymbol{\mu},\boldsymbol{\Sigma}\right) \notag \\ &=\frac{1}{\sum_{l=1}^L \frac{f\left(\boldsymbol{v}_n^{(l)}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)}{f\left(\boldsymbol{v}_n^{(l)}|\boldsymbol{\mu}',\boldsymbol{\Sigma}'\right)}}\sum_{l=1}^L\frac{f\left(\boldsymbol{v}_n^{(l)}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)}{f\left(\boldsymbol{v}_n^{(l)}|\boldsymbol{\mu}',\boldsymbol{\Sigma}'\right)} \log f\left(\boldsymbol{v}_n^{(l)}|\boldsymbol{\mu},\boldsymbol{\Sigma}\right), \label{eq:imp95samp95approx} \end{align}\tag{10}\] where the equality follows from approximating the denominator outside the sum: \[\frac{\int_{R_n^{c_n}}f\left(\boldsymbol{\xi}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)\mathrm{d} \boldsymbol{\xi} }{\int_{R_n^{c_n}}f(\boldsymbol{\xi}|\boldsymbol{\mu}',\boldsymbol{\Sigma}') \mathrm{d} \boldsymbol{\xi}} = \int_{R_n^{c_n}} \frac{f\left(\boldsymbol{\xi}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right) }{f(\boldsymbol{\xi}|\boldsymbol{\mu}',\boldsymbol{\Sigma}') } \frac{f(\boldsymbol{\xi}|\boldsymbol{\mu}',\boldsymbol{\Sigma}') }{\int_{R_n^{c_n}}f(\boldsymbol{\xi}|\boldsymbol{\mu}',\boldsymbol{\Sigma}') \mathrm{d} \boldsymbol{\xi} } \mathrm{d}\boldsymbol{\xi} = \frac{1}{L} \sum_{l=1}^L \frac{f\left(\boldsymbol{v}_n^{(l)}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)}{f\left(\boldsymbol{v}_n^{(l)}|\boldsymbol{\mu}',\boldsymbol{\Sigma}'\right)}.\] By defining the normalized importance weights \(w_{nl} \triangleq \frac{f\left(\boldsymbol{v}_n^{(l)}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)}{f\left(\boldsymbol{v}_n^{(l)}|\boldsymbol{\mu}',\boldsymbol{\Sigma}'\right)}/\sum_{l=1}^L \frac{f\left(\boldsymbol{v}_n^{(l)}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)}{f\left(\boldsymbol{v}_n^{(l)}|\boldsymbol{\mu}',\boldsymbol{\Sigma}'\right)}\), the approximate \(Q\) function 6 can be cast in the following form using the Monte Carlo samples, and 10 : \(\hat{Q}\left(\boldsymbol{\theta}|\boldsymbol{\theta}^{(t)}\right)= \sum_{n=1}^{N} \sum_{l=1}^L w_{nl}\log f\left(\boldsymbol{v}_{n}^{(l)}|\boldsymbol{\mu},\boldsymbol{\Sigma}\right)\). The \(M\)-step thus leads to the following update \[\begin{align} \boldsymbol{\mu}^{(t+1)} = & \frac{\sum_{n=1}^N\sum_{l=1}^L w_{nl}\boldsymbol{v}_{n}^{(l)}}{\sum_{n=1}^N\sum_{l=1}^L w_{nl}},\tag{11}\\ \boldsymbol{\Sigma}^{(t+1)} = & \frac{\sum_{n=1}^N\sum_{l=1}^L w_{nl}\left(\boldsymbol{v}_n^{(l)}-\boldsymbol{\mu}^{(t+1)}\right)\left(\boldsymbol{v}_n^{(l)}-\boldsymbol{\mu}^{(t+1)}\right)^{\top}}{\sum_{n=1}^N\sum_{l=1}^L w_{nl}}.\tag{12} \end{align}\] Therefore, using importance sampling results in a similar procedure to Algorithm 2 but reduces the computation in Algorithm 1. We rewrite the complete EM algorithm with importance sampling in detail in Algorithm 3.

4 Theoretical Properties↩︎

In this section, we study statistical properties of the base model in Section 3 and the associated EM algorithm. We first establish conditions under which the valuation distribution is identifiable from bundle transaction data. We then analyze local convergence of the EM algorithm.

4.1 Identifiability↩︎

Identifiability is a prerequisite for consistent parameter recovery: if there are multiple parameters generating the same distribution of the data, then the parameter cannot be estimated consistently. However, identifiability is not a technical requirement that is always satisfied in this model. Consider the following two examples.

Therefore, we need to establish conditions on \(\mathscr P\) that guarantee identifiability. Intuitively, the model has \(I+I(I+1)/2\) free parameters: \(I\) mean parameters and \(I(I+1)/2\) covariance parameters. Thus, identifiability requires sufficiently many informative region probabilities. However, because these probabilities depend nonlinearly on \((\boldsymbol{\mu},\boldsymbol{\Sigma})\), it is nontrivial to verify either sufficiency or necessity from a simple counting argument.

Note that it 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 sale, i.e., \((p_1,\dots,p_{i-1},p_i',p_{i+1},\dots,p_I)\) for product \(i\). If the above \(I+1\) price menus are offered infinitely often, then \(\mathscr P\) includes the following regions, whose probabilities turn out to be sufficient for identifiability: \[\begin{align} &\left\{\boldsymbol{v}\in \mathbb{R}^{I}, v_i\le p_i\right\}, \; \left\{\boldsymbol{v}\in \mathbb{R}^{I}, v_i\le p_i'\right\}, \; \forall i=1,\dots,I\\ &\left\{\boldsymbol{v}\in \mathbb{R}^{I}, v_{i_1}\le p_{i_1}, v_{i_2}\le p_{i_2}\right\}, \; \left\{\boldsymbol{v}\in \mathbb{R}^{I}, v_{i_1}\le p'_{i_1}, v_{i_2}\le p_{i_2}\right\},\; \left\{\boldsymbol{v}\in \mathbb{R}^{I}, v_{i_1}\le p_{i_1}, v_{i_2}\le p'_{i_2}\right\}, \;\forall i_1,i_2=1,\dots,I \end{align}\]

We provide a brief outline of the proof for Proposition 1. For simplicity, consider a firm offering two products, indexed by \(1\) and \(2\). Customer valuations are drawn from a multivariate Gaussian distribution with mean \(\boldsymbol{\mu}= (\mu_1,\mu_2)^\top\) and covariance matrix \(\boldsymbol{\Sigma}= \begin{bmatrix} \Sigma_{11} & \Sigma_{12}\\ \Sigma_{12} & \Sigma_{22} \end{bmatrix}\).

The necessity is immediate. If product \(i\) is observed at only one price level \(p_i\), then the observable probability \(\mathbb{P}(V_i\le p_i)=\Phi((p_i-\mu_i)/{\sqrt{\Sigma_{ii}}})\) provides only one equation in the two unknowns \((\mu_i,\Sigma_{ii})\). Hence \((\mu_i,\Sigma_{ii})\) cannot be uniquely determined, and the model is not identifiable.

To prove sufficiency, we first identify the parameters of each product. Under separate selling, each product \(i\in\{1,2\}\) is observed at two distinct price levels \(p_i\) and \(p_i'\). Therefore, from the two observable probabilities \(\mathbb{P}(V_i\le p_i)\) and \(\mathbb{P}(V_i\le p_i')\), we obtain two equations involving \((\mu_i,\Sigma_{ii})\). Since \(\Phi(\cdot)\) is strictly increasing, these two equations uniquely determine \(\mu_i\) and \(\Sigma_{ii}\). Hence the marginal distributions of both products are identified. Once \(\mu_1,\mu_2,\Sigma_{11},\Sigma_{22}\) are known, the joint probability under the regular menu can be written as \[\mathbb{P}(V_1\le p_1,\;V_2\le p_2) = \Phi_2\left( \frac{p_1-\mu_1}{\sqrt{\Sigma_{11}}}, \frac{p_2-\mu_2}{\sqrt{\Sigma_{22}}}; \frac{\Sigma_{12}}{\sqrt{\Sigma_{11}\Sigma_{22}}} \right),\] where \(\Phi_2(\cdot,\cdot;\rho)\) denotes the bivariate standard normal CDF with correlation \(\rho\). For fixed thresholds, \(\Phi_2(c_1,c_2;\rho)\) is strictly increasing in \(\rho\). It follows that the observed joint probability uniquely determines the correlation parameter, and hence uniquely determines \(\Sigma_{12}\).

Therefore, all entries of \(\boldsymbol{\mu}\) and \(\boldsymbol{\Sigma}\) are uniquely determined. This establishes identifiability. The detailed proof for Proposition 1 can be found in the online appendix 12.2.

4.2 Convergence of the EM Algorithm↩︎

Classical EM theory establishes convergence to stationary points of the likelihood function under suitable regularity conditions [33], [34]. Convergence to the statistically meaningful fixed point corresponding to the true parameters requires additional conditions, especially because the likelihood can have multiple stationary points. Recent advances such as [35] provide a framework for proving local contractivity of the population EM operator around the true parameter. In this section, we use the approach in [35] to analyze our problem. To simplify the analysis, we focus on the population-level \(Q\)-function, following the first step in [35]; equivalently, we study the idealized limit in which the sample size is infinite. We also assume that all consumers face the same price menu, so the IC polyhedra are common across observations. As a result, we can focus on a partition of \(\mathbb{R}^I\), denoted \(\{\mathcal{P}_\rho\}_{\rho=1}^{P}\). For tractability, we treat the true covariance matrix \(\boldsymbol{\Sigma}^*\) as known and analyze the EM update for the mean parameter \(\boldsymbol{\mu}\). The \(Q\)-function in this simplified setting can be expressed as \[\begin{align} \label{eq:q-func-theory} Q(\boldsymbol{\mu}'\mid \boldsymbol{\mu}) = \sum_{\rho=1}^{P} \int_{\boldsymbol{\xi} \in \mathcal{P}_\rho}f(\boldsymbol{\xi}|\boldsymbol{\mu}^*) d\boldsymbol{\xi} \int_{\boldsymbol{v}\in \mathcal{P}_\rho} \frac{f(\boldsymbol{v} \mid \boldsymbol{\mu})}{\left(\int_{\boldsymbol{\xi}\in \mathcal{P}_\rho}f(\boldsymbol{\xi}|\boldsymbol{\mu}) d\boldsymbol{\xi}\right)} \log f(\boldsymbol{v}|\boldsymbol{\mu}') \mathrm{d}\boldsymbol{v}. \end{align}\tag{13}\] Compared to 6 , we remove the dependence on \(\boldsymbol{\Sigma}\) and use the same partition. Moreover, because the function is defined at the population level, the probability \(\int_{\boldsymbol{\xi} \in \mathcal{P}_\rho}f(\boldsymbol{\xi}|\boldsymbol{\mu}^*)\mathrm d\boldsymbol{\xi}\) in each region replaces the samples in 6 . The population-level EM algorithm updates \(\boldsymbol{\mu}\) iteratively according to the operator: \[\begin{align} M(\boldsymbol{\mu}) \triangleq \argmax_{\boldsymbol{\mu}'} Q(\boldsymbol{\mu}' \mid \boldsymbol{\mu}). \end{align}\] More precisely, given initialization \(\boldsymbol{\mu}^{(0)}\), the EM algorithm uses \(\boldsymbol{\mu}^{(t+1)}=M\left(\boldsymbol{\mu}^{(t)}\right)\) to obtain a sequence \(\left\{\boldsymbol{\mu}^{(t)}\right\}_{t=0}^{\infty}\). The theoretical question is whether we have \(\boldsymbol{\mu}^{(t)}\to\boldsymbol{\mu}^*\). Note that the true value \(\boldsymbol{\mu}^*\) always maximizes the population-level likelihood function [36]. Moreover, it satisfies the self-consistency condition [37]: \(\boldsymbol{\mu}^* = \argmax_{\boldsymbol{\mu}'} Q(\boldsymbol{\mu}' \mid \boldsymbol{\mu}^*)\). Therefore, \(\boldsymbol{\mu}^*\) is a fixed point of the EM operator \(M\). [35] establish conditions under which \(M\) is a contraction mapping in a neighborhood of \(\boldsymbol{\mu}^*\), denoted by \(\mathbb{B}(r;\boldsymbol{\mu}^*)\) with Euclidean radius \(r\). Since iterates of a contraction mapping converge to its fixed point, our goal is to verify analogous conditions for the bundle estimation problem.

Following [35], convergence is guaranteed if \(M(\cdot)\) is a contraction mapping in a neighborhood of \(\boldsymbol{\mu}^*\). This is ensured by two conditions: (1) Concavity of \(Q(\boldsymbol{\mu}\mid\boldsymbol{\mu}^*)\). There exists \(\lambda > 0\) such that for all \(\boldsymbol{\mu}_1, \boldsymbol{\mu}_2\) in a neighborhood \(\mathbb{B}(r;\boldsymbol{\mu}^*)\), we have \(q(\boldsymbol{\mu}_1) - q(\boldsymbol{\mu}_2) - \langle \nabla q(\boldsymbol{\mu}_2), \boldsymbol{\mu}_1 - \boldsymbol{\mu}_2 \rangle \le -\tfrac{\lambda}{2} \|\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2\|_2^2\), where \(q(\boldsymbol{\mu}) \triangleq Q(\boldsymbol{\mu}\mid \boldsymbol{\mu}^*)= \sum_{\rho=1}^{P} \int_{\boldsymbol{v}\in \mathcal{P}_\rho} f(\boldsymbol{v} \mid \boldsymbol{\mu}^*)\log f(\boldsymbol{v}|\boldsymbol{\mu}) \mathrm{d}\boldsymbol{v}\). (2) First-order stability. There exists \(\gamma < \lambda\) such that \(\| \nabla Q(M(\boldsymbol{\mu}) \mid \boldsymbol{\mu}^*) - \nabla Q(M(\boldsymbol{\mu}) \mid \boldsymbol{\mu}) \|_2 \le \gamma \| \boldsymbol{\mu}- \boldsymbol{\mu}^* \|_2\).

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.

To provide some intuition for the assumption, consider the case \(J=0\), i.e., no bundle is offered and the only feasible action is no-purchase. In this case, \(\mathbb{E}[\boldsymbol{V}' \mid \mathcal{P}']\) takes a single value and \(\mathrm{Var}(\mathbb{E}[\boldsymbol{V}' \mid \mathcal{P}'])\) is a degenerate \(I \times I\) zero matrix, so Assumption 1 is clearly not satisfied. It is not surprising that the EM algorithm cannot recover the true parameters in this degenerate case, because the model is not identifiable, i.e., no information about \(\boldsymbol{\mu}^*\) can be learned from observing only no-purchase decisions.

Now consider the other extreme. Suppose the price menu induces a partition \(\{{\mathcal{P}_\rho}'\}_{\rho=1}^P\) that is so fine that each \({\mathcal{P}_\rho}'\) almost degenerates to a single point. In this case, \(\mathbb{E}[\boldsymbol{V}' \mid \boldsymbol{V}' \in {\mathcal{P}_\rho}']\) effectively traces out the entire support of \(\boldsymbol{V}'\), so \(\mathrm{Var}(\mathbb{E}[\boldsymbol{V}' \mid \mathcal{P}']) \approx \boldsymbol{I}\), satisfying Assumption 1. The two examples illustrate the intuition that the EM algorithm performs well when the price menu induces many small regions, which is consistent with the discussion after Example 2.

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

5 Model Extensions↩︎

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

5.1 Product Synergy in a Bundle↩︎

In Section 3, the bundle utility is modeled as the sum of the utilities of the included products, with the covariance structure of product valuations capturing implicit relationships across products. In practice, however, products within a bundle may exhibit explicit complementarity or substitutability that renders the overall utility of the bundle non-additive. These synergy effects imply that a bundle’s utility is not necessarily equal to the sum of its standalone components; instead, the combination of specific products may amplify or diminish the perceived value. For example, pairing a printer with compatible ink cartridges often yields additional value relative to purchasing the items separately (positive synergy), whereas bundling two devices with overlapping functionality may reduce incremental benefit (negative synergy). Prior studies typically incorporate such effects through pairwise or higher-order interaction terms [38]–[40]. In line with this literature, we adopt a pairwise interaction structure, which offers a flexible yet tractable representation of inter-product synergies. We next investigate the estimation problem when bundle choices may exhibit such synergy effects.

Vectorized representation of bundles. To express the utility model compactly, we first introduce an equivalent vectorized representation of bundles. Let \(\boldsymbol{x}_j\in \left\{0,1\right\}^I\) denote bundle \(j\), where the \(i\)-th entry equals 1 if product \(i\in[I]\) is included in the bundle. Customer \(n\) is faced with a menu of bundles \(\boldsymbol{X}_n\in\{0,1\}^{I\times J_n}\), where each column of \(\boldsymbol{X}_n\) is a binary vector denoting an offered bundle and \(J_n\) is the number of bundles shown to customer \(n\). The corresponding bundle prices are collected in \(\boldsymbol{p}_n=(p_n^1,\dots,p_n^{J_n})\), which is observed in the data. Customer \(n\) chooses \(c_n\in[J_n]\cup\{0\}\), where \(c_n = 0\) denotes the no-purchase option. The dataset therefore consists of tuples \(\{(\boldsymbol{X}_n, \boldsymbol{p}_n,c_n)\}_{n=1}^N\). Note that customers do not necessarily observe the same menu, nor do they see the complete set of possible bundles. This vectorized representation generalizes the formulation in Section 3 and provides a convenient way to incorporate synergy effects.

Given a vector of individual product valuation \(\boldsymbol{v}\), the utility of bundle \(\boldsymbol{x}_j\) to customer \(n\) is specified as \[\label{eq:bundle-utility-synergy} u_{nj}\triangleq \boldsymbol{x}_j^\top \boldsymbol{v}-p_n^j+\boldsymbol{x}_j^\top \boldsymbol{A}\boldsymbol{x}_j,\tag{14}\] where the upper triangular matrix \(\boldsymbol{A}\in\mathbb{R}^{I\times I}\) captures the pairwise synergy between products that appear in the same bundle. In particular, if products \(i_1\) and \(i_2\) are both present in bundle \(j\), then the bundle utility includes an additional term \(A_{i_1 i_2}\). The sign of \(A_{i_1 i_2}\) reflects whether the two products are complementary (\(A_{i_1 i_2}>0\)) or substitutable (\(A_{i_1 i_2}<0\)). We impose \(A_{ii}=0\) for all \(i\) so that \(\boldsymbol{A}\) only encodes interactions between distinct products. When \(\boldsymbol{A}\equiv\boldsymbol{0}\), the model reduces to the baseline additive specification in Section 3.

IC polyhedra with synergy. Under the synergy specification 14 , the IC polyhedra in 4 can be rewritten in a vectorized form, comparing the utility of each bundle to that of the chosen one and incorporating the synergy terms that arise whenever two products appear together in a bundle. In particular, the bundle-utility vector for customer \(n\) is given by \[\label{eq:bundle-utility-synergy-vec} \boldsymbol{u}_n(\boldsymbol{v}) = \boldsymbol{X}_n^\top \boldsymbol{v}-\boldsymbol{p}_n+ diag(\boldsymbol{X}_n ^\top \boldsymbol{A}\boldsymbol{X}_n)\in\mathbb{R}^{J_n},\tag{15}\] where \(diag(\cdot)\) extracts the diagonal of a square matrix as a column vector. Then for customer \(n\), we can define the corresponding IC region in the bundle utility space by \[U_n^{c_n} = \left\{ \begin{array}{ll} \{ \boldsymbol{u}\in \mathbb{R}^{J_n} | \; \boldsymbol{u}\le \boldsymbol{1}_{J_n}u_{nc_n},\; u_{nc_n}\ge 0 \}, & c_n\neq0;\\ \{ \boldsymbol{u}\in \mathbb{R}^{J_n} | \; \boldsymbol{u}\le \boldsymbol{0}_{J_n}\}, & c_n=0. \end{array} \right.\] where \(\boldsymbol{1}_n\) and \(\boldsymbol{0}_n\) denote \(n\)-dimensional column vectors of ones and zeros, respectively. A valuation vector \(\boldsymbol{v}\) satisfies the IC constraints if and only if its induced utility vector \(\boldsymbol{u}_n(\boldsymbol{v})\) lies in this set. Hence, the corresponding IC polyhedra in the product valuation space is \[\label{eq:ic-polytope-synergy} K_n^{c_n}\triangleq\left\{\boldsymbol{v} \in \mathbb{R}^{I} \left| \; \boldsymbol{u}_n(\boldsymbol{v})\in U_n^{c_n}\right. \right\}.\tag{16}\]

We next articulate how to adapt the computational framework introduced in Section 3 to this setting.

5.1.0.1 E-step.

Following the standard EM framework, we evaluate the expectation of the complete-data log-likelihood with respect to the conditional distribution of the latent variables \(\boldsymbol{v}\). The likelihood can be written as \(\mathcal{L}(\boldsymbol{\theta};\boldsymbol{D}) = \prod_{n=1}^N\int_{K_n^{c_n}(\mathbf{A})} f(\boldsymbol{v}|\boldsymbol{\mu},\boldsymbol{\Sigma}) \mathrm{d}\boldsymbol{v} =\prod_{n=1}^N\int \mathbb{I}\{\boldsymbol{v} \in K_n^{c_n}(\mathbf{A})\} f(\boldsymbol{v}|\boldsymbol{\mu},\boldsymbol{\Sigma}) \mathrm{d}\boldsymbol{v}\), where \(f(\cdot)\) denotes the PDF of the multivariate normal distribution. The complete-data log-likelihood can therefore be written as \[\label{eq:llh95complete95synergy} \begin{align} \ell(\boldsymbol{\theta};\boldsymbol{D},\boldsymbol{Z}) =\sum_{n=1}^N \log f(\boldsymbol{v}_n\mid \boldsymbol{\mu},\boldsymbol{\Sigma}) +\sum_{n=1}^N \log \mathbb{I}\{\boldsymbol{v}_n\in K_n^{c_n}(\boldsymbol{A})\}. \end{align}\tag{17}\] In the E-step, the conditional law of \(\boldsymbol{v}_n\) is no longer a truncated Gaussian over a fixed region. Thus, a direct M-step based on 17 is difficult because the term \(\mathbb{I}\{\boldsymbol{v}_n\in K_n^{c_n}(\boldsymbol{A})\}\) is non-differentiable in \(\boldsymbol{A}\), and the feasible region itself changes with \(\boldsymbol{A}\). We address the non-differentiability by a differentiable surrogate that approximates this objective.

In particular, we approximate the hard IC condition by a product of sigmoid terms: \[\label{eq:sigmoid95approx} \mathbb{I}\{\boldsymbol{v}\in K_n^{c_n}(\boldsymbol{A})\} \approx \prod_{j=1}^{J_n}\sigma\!\left(\frac{\Delta U_{nj}(\boldsymbol{v};\boldsymbol{A})}{\lambda}\right),\tag{18}\] where \(\sigma(\cdot) = (1 + e^{-\cdot})^{-1}\) is the sigmoid function and \(\lambda>0\) is a smoothing parameter, and \(\Delta U_{nj} = \mathbb{I}\{c_n\neq 0\}\,u_{n c_n} - \mathbb{I}\{c_n=0\;\text{or}\;j\neq c_n\}\,u_{nj}\). Essentially, we turn the IC indicator \(u_{nc_n}\ge u_{nj}\) to a soft margin \(\sigma((u_{nc_n}-u_{nj})/\lambda)\), which measures how much the observed choice \(c_n\) is preferred to alternative \(j\) under \((\boldsymbol{v},\boldsymbol{A})\). This approximation gives a differentiable surrogate likelihood, so \(\boldsymbol{A}\) can be updated in the M-step using gradient-based optimization. As \(\lambda\to 0\), one can show that the surrogate approaches the original indicator.

Consequently, under the proposed smoothing approximation with synergy matrix \(\boldsymbol{A}\), the indicator function in E-step is replaced by the sigmoid approximation in 18 . As a result, the hard truncation of the feasible region is replaced by a continuous weighting term. The smoothed likelihood is therefore approximated by \(\mathcal{L}(\boldsymbol{\theta};\boldsymbol{D})\approx \prod_{n=1}^N\int \prod_{j=1}^{J_n} \sigma\left(\frac{\Delta U_{nj}(\boldsymbol{v}; \mathbf{A})}{\lambda}\right) f(\boldsymbol{v}|\boldsymbol{\mu},\boldsymbol{\Sigma}) \mathrm{d}\boldsymbol{v}.\) The complete-data log-likelihood can therefore be written as \[\label{eq:likelihood95sigmoid95complete-data} \begin{align} \ell(\boldsymbol{\theta};\boldsymbol{D},\boldsymbol{Z})=\sum_{n=1}^N\left(\log f(\boldsymbol{v}_n | \boldsymbol{\mu}, \boldsymbol{\Sigma}) +\sum_{j=1}^{J_n} \log \sigma\left(\frac{\Delta U_{nj}(\boldsymbol{v}_n; \boldsymbol{A})}{\lambda}\right)\right). \end{align}\tag{19}\] Substituting this conditional density into the expectation yields the smoothed \(Q'\) function \[\label{eq:smoothed95Q} \begin{align} Q'(\boldsymbol{\theta} \mid \boldsymbol{\theta}^{(t)}) &=\sum_{n=1}^{N}\int q_n(\boldsymbol{v};\boldsymbol{\theta}^{(t)}) \left( \log f(\boldsymbol{v} \mid \boldsymbol{\mu},\boldsymbol{\Sigma})+\sum_{j=1}^{J_n}\log \sigma\left( \frac{\Delta U_{nj}(\boldsymbol{v};\boldsymbol{A})}{\lambda} \right) \right) \mathrm{d}\boldsymbol{v}\\ \end{align},\tag{20}\] where \(q_n(\boldsymbol{v}; \boldsymbol{\theta}^{(t)})= \left( f(\boldsymbol{v} \mid \boldsymbol{\mu}^{(t)}, \boldsymbol{\Sigma}^{(t)}) \prod_{j=1}^{J_n}\sigma(\frac{\Delta U_{nj}(\boldsymbol{v}; \boldsymbol{A}^{(t)})}{\lambda}) \right) \big/ \left( {\int f(\boldsymbol{\xi} \mid \boldsymbol{\mu}^{(t)}, \boldsymbol{\Sigma}^{(t)})\prod_{j=1}^{J_n}\sigma(\frac{\Delta U_{nj}(\boldsymbol{\xi}; \boldsymbol{A}^{(t)})}{\lambda}) \mathrm{d}\boldsymbol{\xi}} \right).\)

Compared to the EM formulation in Section 3.2, the resulting \(Q'\) function differs in two important aspects. First, the truncation of the utility distribution is replaced by a smooth sigmoid term. Second, the log-likelihood now contains an additional contribution from the smoothed IC constraints, which introduces explicit dependence on the synergy parameter \(\boldsymbol{A}\) inside the integrand. This modification enables gradient-based optimization with respect to \(\boldsymbol{A}\) while preserving the overall EM structure.

Similar to the base model, substituting the Monte Carlo approximation into Eq. (20 ) yields the \(\hat{Q}'\) function \[\label{eq:smoothed95Qhat} \begin{align} \hat{Q}'(\boldsymbol{\theta} \mid \boldsymbol{\theta}^{(t)})&=\sum_{n=1}^N\sum_{l=1}^L \overline{\omega}_{nl}\left( \log f(\boldsymbol{v}_n^{(l)}|\boldsymbol{\mu},\boldsymbol{\Sigma})+\log \omega(\boldsymbol{v}_n^{(l)};\boldsymbol{A}) \right)\\ &=C-\frac{N}{2}\log|\boldsymbol{\Sigma}|+\sum_{n=1}^N\sum_{l=1}^L \overline{\omega}_{nl}\left( -\frac{1}{2}(\boldsymbol{v}_n^{(l)}-\boldsymbol{\mu})^\top\boldsymbol{\Sigma}^{-1}(\boldsymbol{v}_n^{(l)}-\boldsymbol{\mu})+\log \omega(\boldsymbol{v}_n^{(l)};\boldsymbol{A}) \right), \end{align}\tag{21}\] where \(C\) is a constant with respect to \(\boldsymbol{\theta}\) and \[\label{eq:sigmoid95vals} \omega(\boldsymbol{v};\boldsymbol{A}) \triangleq\prod_{j=1}^{J_n}\sigma\left( \frac{ \Delta U_{nj}(\boldsymbol{v}; \boldsymbol{A}) }{\lambda} \right),\quad \overline{\omega}_{nl} \triangleq \frac{ \omega(\boldsymbol{v}_{n}^{(l)};\boldsymbol{A}^{(t)})}{\sum_{l=1}^L \omega(\boldsymbol{v}_{n}^{(l)};\boldsymbol{A}^{(t)})}.\tag{22}\]

5.1.0.2 M-step.

In the M-step, we maximize \(\hat{Q'}(\boldsymbol{\theta}\mid\boldsymbol{\theta}^{(t)})\) with respect to \((\boldsymbol{\mu},\boldsymbol{A},\boldsymbol{\Sigma})\). The objective function is concave in each parameter block. First, the Gaussian log-likelihood term is concave in \(\boldsymbol{\mu}\) and in \(\boldsymbol{\Sigma}^{-1}\), which follows from standard properties of the multivariate normal distribution. Second, Lemma 1 below establishes that the contribution involving \(\boldsymbol{A}\) is also concave. The detailed proof is deferred to the Appendix 12.4.

Unlike \(\boldsymbol{\mu}\) and \(\boldsymbol{\Sigma}\), which admit closed-form updates from the first-order conditions, the synergy matrix \(\boldsymbol{A}\) does not admit a closed-form solution in the M-step and must be estimated using gradient-based optimization, where the gradient is given in Theorem 2.

The proof of Theorem 2, the closed-form updates for \(\boldsymbol{\mu}\) and \(\boldsymbol{\Sigma}\), and the corresponding update formulas under the importance sampling framework are provided in Appendix 8.1.

5.2 Censored Demand↩︎

In practice, transaction data are often censored because the firm records purchases but does not observe customers who visit without buying. We next extend our framework to account for this demand censoring.

For ease of exposition, suppose the firm offers a single price menu \(\boldsymbol{p} = (p^{1},\dots,p^J)\) to all customers. The observed data consist of \((p^{1},\dots,p^J,c_n)\) for each purchasing customer \(n\), where \(c_n>0\) for \(n=1,\dots,N\). These data can be converted to the IC polyhedra \(\{R^0,R^1,\dots,R^J\}\) defined in 4 . Because the price menu is fixed, these regions are common across customers. As in the base model, we treat each customer valuation \(\boldsymbol{v}_n\), \(n=1,\dots,N\), as missing data.

Because of the demand censoring, the total number of customers that have visited the store, \(N'\ge N\), is also unobserved, as well as the valuations of the censored customers \(\boldsymbol{v}_n\), \(n=N+1,\dots,N'\). Therefore, the missing data is \(\boldsymbol{Z}\triangleq \left\{N', \boldsymbol{v}_{1},\dots,\boldsymbol{v}_{N'}\right\}\). It is known that \(\boldsymbol{v}_n\in R^0\) for the censored customers \(n=N+1,\dots,N'\). 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 \(\boldsymbol{D} \triangleq \left\{N_1,\dots,N_J\right\}\). Below we describe the implementation of the EM algorithm to estimate \(\boldsymbol{\theta}= (\boldsymbol{\mu}, \boldsymbol{\Sigma})\).

Relative to the base model in Section 3, the E-step now includes an additional outer expectation over the unobserved total number of customers, \(N'\mid \boldsymbol{D},\boldsymbol{\theta}^{(t)}\). Under \(\boldsymbol{\theta}^{(t)}\), a customer is censored with probability \(\mathbb{P}(R^0\mid \boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)})\), the probability of the no-purchase region. Thus, conditional on \(N'\), observing \(N\) purchases corresponds to \(N\) successes in \(N'\) Bernoulli trials with success probability \(1-\mathbb{P}(R^0\mid \boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)})\).

Once the distribution of \(N'\) is identified, the E-step is approximated by Monte Carlo simulation as in Section 3.2: in each instance, we first draw \(N'^{(l)}\) from a negative binomial distribution, then draw the observed customer valuations from truncated Gaussians on the corresponding IC polyhedra \(R^j\), and the censored customer valuations from a truncated Gaussian on \(R^0\). The M-step can then be adapted using a \(\hat{Q}\) function computed from the simulated samples above, and the parameters are updated by maximizing this \(\hat{Q}\).

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'^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.

The detailed derivation of the E-step, the negative binomial form for \(N'\mid \boldsymbol{D}, \boldsymbol{\theta}^{(t)}\), the Monte Carlo procedure, and the resulting M-step update formulas are provided in Appendix 8.2.

6 Numerical Experiments↩︎

To assess estimation accuracy, we use the \(\ell_1\)-error \((\|\boldsymbol{\mu}^{(t+1)}-\boldsymbol{\mu}^*\|_1+\|\boldsymbol{\Sigma}^{(t+1)}-\boldsymbol{\Sigma}^*\|_1)/(I^2 + I)\), which measures the average absolute deviation of the estimated parameters from their true values.

6.1 EM Algorithm in the Base Model↩︎

We first evaluate the convergence of Algorithm 3 by tracking the log-likelihood and estimation errors across iterations. We also investigate how the estimation accuracy scales with the number of products and customers.

6.1.0.1 EM convergence.

We investigate the convergence of Algorithm 3 against the number of iterations. In this experiment, we consider two products (\(I=2\)) and one bundle. We randomly generate a sample of \(N=1000\) transactions and initialize the algorithm from three different points with varying distances from the true parameters: within one standard deviation (close), two standard deviations (midway), and five standard deviations of the true mean (distant). In each EM iteration, we calculate the log-likelihood using the estimated parameters and define the error as the difference between two subsequent values. We terminate the algorithm when either the error falls below \(0.005\) or the number of EM iterations reaches \(500\).

The left panel of Figure [fig:byinitials] shows that the average log-likelihood of Algorithm 3 converges toward that under the true parameter across the three initializations. The right panel of Figure [fig:byinitials] shows the trajectory of the \(\ell_1\)-error. When the initial point is far from the true parameters, the \(\ell_1\)-error may initially increase before decreasing.

We provide a more comprehensive comparison for different numbers of products and observations in Table 1. For each setting, we initialize the estimator from five standard deviations from the true mean (distant point). We report the final \(\ell_1\)-error after convergence for \(I\in\{2,3,4,5,6\}\) and \(N\in \{1000, 1500,2000,2500,3000, 5000, 10000\}\). Taken together, these results show that Algorithm 3 exhibits stable and consistent improvement as the sample size increases, demonstrating its scalability in both dimensions.

We additionally report the computational running time of the algorithm. All experiments were conducted on a Linux server equipped with 64 physical CPU cores (128 logical cores) running Python 3.10.19. To provide a hardware-independent measure of computational cost, we report the total CPU time, defined as the sum of the CPU seconds consumed by all worker processes. Most experimental settings finish within several minutes to tens of minutes, although a small number of instances require longer runtime. Since the computational cost depends on the realized observations and the corresponding convergence behavior of the EM iterations, the reported runtimes exhibit variability across datasets.

6.1.0.2 Comparison to Bayesian methods.

We compare Algorithm 3 with the method introduced in [28]. The latter conducts posterior inference by combining data augmentation, Gibbs sampling, and Metropolis–Hastings steps within a general MCMC framework. The key idea is to replace a single difficult draw from the joint posterior by a sequence of easier draws from the full conditional distributions of parameter blocks, thereby generating a Markov chain whose stationary distribution is the target joint posterior. We adopt their methodology by factoring in the posterior probability \(\mathbb{P}\left(\boldsymbol{\mu}^{(t)}, \boldsymbol{\Sigma}^{(t)} \mid R_n^{c_n}\right) \propto \mathbb{P}\left(R_n^{c_n} \mid \boldsymbol{\mu}^{(t)}, \boldsymbol{\Sigma}^{(t)}\right) \mathbb{P}\left(\boldsymbol{\mu}^{(t)}\right) \mathbb{P}\left(\boldsymbol{\Sigma}^{(t)}\right)\). We refer to Appendix 10 for further details on this implementation.

6.1.0.3 Comparison to the MNL model.

We conduct numerical experiments comparing our Gaussian-based model with the multinomial logit (MNL) model, one of the most widely used benchmarks in discrete choice analysis. We refer to Appendix 11 for further details.

We generate the training data according to the base model. Thus, the MNL model is misspecified. We estimate the parameters using our method and the standard maximum-likelihood estimator for the MNL model. To evaluate the goodness-of-fit, we construct a test dataset of \(m\) menus, where each menu consists of a bundle-offering matrix \(\boldsymbol{X}_n\) and an associated bundle-price vector \(\boldsymbol{p}_n\), using the same procedure as the training data. We let \(m=100\). For each menu, we compute the ground-truth choice probability \(\mathbb{P}^{\text{true}}\) and measure the root mean squared error (RMSE) of an estimator \[\text{RMSE}=\frac{1}{m}\sum_{n=1}^m\sqrt{ \frac{1}{J_n+1} \sum_{j=1}^{J_n+1} \left( \mathbb{P}_{nj}^{\text{true}}-\mathbb{P}_{nj}^{\text{model}} \right)^2 }\] over the \(J_n\) products/bundles in a menu and \(m\) menus in the test set. We report the out-of-sample log-likelihood and the RMSE in Table 3.

Because of misspecification, particularly the failure to incorporate the correlation between the valuations of bundles and the included products, the MNL model has lower predictive performance. Indeed, the MNL model treats the (Gumbel) random utilities of bundles and the included products as independent. This result illustrates the limitation of adapting MNL directly to the bundle setting.

6.2 Product Synergy↩︎

In this subsection, we test the algorithm in Section 5.1. The experimental setup is similar to that of the base model (see Section 6), including the product set, offered bundles, price menus, and customers’ individual-product valuation distribution. We consider \(I \in\{2, 3, 4, 5, 6\}\) products and \(N \in\{1000, 1500, 2000, 2500, 3000, 5000, 10000\}\) transactions. The key difference from the base model is that each consumer’s utility now includes pairwise product interactions represented by the synergy matrix \(\boldsymbol{A}\). Specifically, the upper-triangular entries of \(\boldsymbol{A}\) are randomly sampled from a uniform distribution \(\mathcal{U}[-3, 3]\), while the diagonal and lower-triangular entries are fixed to zero.

We run Algorithm 9 on these datasets to jointly estimate the parameters \((\boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{A})\). We use the \(\ell_1\)-error \(\left( \|\boldsymbol{\mu}^{(t+1)}-\boldsymbol{\mu}^*\|_1+\|\boldsymbol{A}^{(t+1)}-\boldsymbol{A}^*\|_1+\|\boldsymbol{\Sigma}^{(t+1)}-\boldsymbol{\Sigma}^*\|_1 \right)/ \left(\left(3I^2 + I\right)/2\right)\) to measure estimation accuracy and track algorithmic convergence.

6.2.0.1 EM Convergence.

Figure [fig:l1error95synergy] illustrates the convergence of Algorithm 9 under different problem settings. Algorithm 9 exhibits stable convergence in the reported instances. Although instances with larger \(I\) require more iterations to reach convergence, the overall convergence pattern is similar across product dimensions. Larger samples lead to smaller estimation errors and faster convergence. Panel [fig:EMperformance95trajectory] illustrates the trajectory of \((\boldsymbol{\mu}^{(t)},\boldsymbol{A}^{(t)})\) for an instance of \(I=2\) products and \(N=10000\) transactions during Algorithm 9. Starting from a distant initialization, the algorithm rapidly moves into a neighborhood of the true parameter values and then gradually refines the estimates until convergence. This pattern indicates that, in this instance, the algorithm moves toward the region around the true parameters and then exhibits stable local refinement.

6.2.0.2 EM Performance.

Panel 7 illustrates recovery of \((\boldsymbol{\mu},\boldsymbol{\Sigma})\), which determine the distribution of the two individual-product valuations. The contours of the estimated and ground-truth distributions are largely aligned, suggesting that the algorithm recovers both the marginal distribution of the utilities and the cross-product correlations. We also study the approximation quality of the sigmoid smoothing in 18 relative to the original Monte Carlo EM algorithm, Algorithm 3. When the synergy matrix \(\boldsymbol{A}=\boldsymbol{0}\), the extended model reduces exactly to the base model introduced in Section 3. This setting therefore provides a natural benchmark to evaluate the approximation error introduced by the sigmoid method. To evaluate the approximation quality of the sigmoid method, we generate synthetic datasets with \(N=5000\) transactions consisting of two products and one bundle under the base model with \(\boldsymbol{A}=\boldsymbol{0}\). Both Algorithm 3 and Algorithm 9 are then applied to estimate \((\boldsymbol{\mu},\boldsymbol{\Sigma})\). The data-generating process follows the setting of the base model described in Section 6. Since the true data-generating process contains no bundle synergy effects, any performance difference between the two methods reflects the approximation error caused by replacing the indicator function with the sigmoid smoothing. Both algorithms are initialized from the same starting point to ensure a fair comparison. Figure [fig:2D95Gaussian95noA] illustrates the estimation results for the case \(I=2\) and \(N=5000\). The figure plots the ground-truth distribution together with the distributions implied by the parameters estimated by MCEM and the sigmoid-based EM algorithm. The two estimated distributions are close to the ground-truth distribution, suggesting that the sigmoid approximation introduces little additional estimation error in this setting.

To further investigate the statistical properties of Algorithm 9, we evaluate the \(\ell_1\)-error across \(N\in \{1000,2500,5000,10000\}\) and \(I\in \{2,3,4,5,6\}\). As reported in Table 4, the \(\ell_1\)-error generally decreases as the number of transactions \(N\) increases for a fixed product dimension \(I\). For larger \(I\), the estimation problem becomes more challenging because the number of valuation, covariance, and synergy parameters increases.

6.3 Censored Data↩︎

In this section, we evaluate how well the proposed method recovers model parameters when no-purchase observations are censored from the sales data.

Starting from the data generation process for the base model experiments, we remove all no-purchase observations, resulting in a censored dataset. The proportion of removed no-purchase transactions ranges from approximately 10.84% to 17.54% across different datasets. We fix the number of products at \(I=2\) and the number of uncensored transactions at \(N=2500\). We generate five datasets with different parameters \((\boldsymbol{\mu},\boldsymbol{\Sigma})\) and estimate the parameters using the EM method for the censored data and the complete data. We compute the log-likelihood for both methods on a test set with \(N=2000\).

The results in Table 5 show that the censored EM algorithm remains close to the complete-data benchmark in test log-likelihood, although its parameter error is higher, as expected from the information loss induced by censoring. Across all datasets, both methods achieve high test log-likelihood scores, while the censored-data estimates have larger \(\ell_1\)-errors than the complete-data estimates. We show the converging trajectories of log-likelihood for both methods in Figure [fig:censored95llh]. Taken together, these findings suggest that censoring has a more visible effect on parameter recovery than on out-of-sample likelihood in these experiments.

6.4 Real Data↩︎

This subsection evaluates the performance of our proposed EM framework by comparing several variants of our model, including the base EM, Gaussian mixture model (GMM), and the bundle synergy model, against two benchmark approaches: the multinomial logit (MNL) model and an MH-based Bayesian method.

We conduct the empirical analysis using the JD.com dataset [41], which provides detailed transaction and clickstream information. The dataset consists of 486,928 transactions involving 9,159 unique stock-keeping units (SKUs) shipped to 60 districts in March 2018. Each transaction record includes product prices and various promotional mechanisms, such as direct discounts, bundle discounts, and gift incentives. In addition, the dataset contains customers’ click histories across SKUs during the same period, which allows us to infer customers’ consideration sets. Notably, only approximately 1% of customers in the dataset purchase bundles under promotional schemes.

To obtain initial values, we randomly sample 1,000 observations from the training set and run each algorithm for 50 iterations. We then use the estimated parameters as the starting point for the main algorithms. In the main training phase, we stop when the change in average log-likelihood between consecutive iterations is below \(10^{-6}\). Top-\(X\) accuracy measures the fraction of test observations for which the actual purchase is among the model’s top-\(X\) predicted alternatives.

Table 6 shows that all three proposed models (Base, Synergy, and GMM) outperform the MH and MNL benchmarks across the reported evaluation metrics, indicating that the proposed framework captures consumer choice behavior more effectively in this dataset. Among the proposed models, the Base model already achieves strong predictive performance and provides a significant improvement over the benchmark methods. The Synergy and GMM models further extend the Base model by incorporating additional behavioral structures, leading to modest but consistent improvements in selected performance metrics. In particular, the GMM model achieves the best out-of-sample log-likelihood and Top-3 accuracy, while the Synergy model attains the highest Top-1 and Top-5 accuracy. These results suggest that the additional flexibility introduced by the Synergy and GMM formulations can capture additional patterns in consumer choice behavior that are not fully represented in the Base model. Nevertheless, the relatively small performance gap among the three proposed models indicates that the Base model alone is already capable of explaining much of the observed choice behavior.

7 Conclusion↩︎

We investigate the practical problem of estimating customers’ valuations in the presence of bundle sales promotions, which is essential for designing bundle sales promotions and pricing strategies. We formulate this problem as a statistical problem of fitting a distribution (multivariate Gaussian distribution in the base model) when observations are censored within various partitions, and we propose an Expectation-Maximization algorithm. We study the properties of the problem and show that under separate selling and offering each single product at two different prices, the estimation problem is identifiable. We also provide sufficient conditions under which the population EM operator converges locally to the true parameters when initialized within a basin of attraction. Building upon the base model, we develop a Product Synergy model to capture complementarity and substitutability among products within a bundle, allowing the utility of a bundle to depend not only on the valuations of its constituent products but also on their interactions. Furthermore, we extend the base model to censored-demand settings, where no-purchase observations may be unobserved, and to Gaussian mixture models that capture latent customer segments. We demonstrate the efficacy of our algorithm through synthetic numerical experiments and sales data from JD.com.

However, there are still unexplored problems in the literature, such as considering competitive prices when estimating the demand distribution given that the sales data of the competitor is censored for one’s sales data. Additionally, statistical questions related to calculating the error of Monte Carlo samples and the convergence error of the EM algorithm need further exploration. Moreover, identifying the number of components when performing the Gaussian Mixture model is an interesting direction for future research.

8 Detailed EM Algorithm for Model Extensions↩︎

This appendix provides the detailed derivations of the EM algorithm for the model extensions in Section 5.

8.1 Product Synergy in a Bundle↩︎

We use the same notation as in Section 5.1: \(\boldsymbol{D} = \{(p_n^1,\dots,p_n^{J_n},c_n)\}_{n=1}^N\) denotes the observed data, and \(\boldsymbol{\theta}=(\boldsymbol{\mu},\boldsymbol{\Sigma},\boldsymbol{A})\) denotes the parameters to be estimated, where \(\boldsymbol{A}\) is the bundle synergy matrix.

8.1.0.1 E-step.

Following the standard EM framework, we evaluate the expectation of the complete-data log-likelihood with respect to the conditional distribution of the latent variables \(\boldsymbol{v}\). Compared with the base model, the key difference arises from the replacement of the indicator function defining the IC polyhedron. The likelihood can be written as \(\mathcal{L}(\boldsymbol{\theta};\boldsymbol{D}) = \prod_{n=1}^N\int_{K_n^{c_n}(\mathbf{A})} f(\boldsymbol{v}|\boldsymbol{\mu},\boldsymbol{\Sigma}) \mathrm{d}\boldsymbol{v} =\prod_{n=1}^N\int \mathbb{I}\{\boldsymbol{v} \in K_n^{c_n}(\mathbf{A})\} f(\boldsymbol{v}|\boldsymbol{\mu},\boldsymbol{\Sigma}) \mathrm{d}\boldsymbol{v}\), where \(f(\cdot)\) denotes the PDF of the multivariate normal distribution. The complete-data log-likelihood can therefore be written as \[\label{eq:likelihood95synergy95complete-data} \begin{align} \ell(\boldsymbol{\theta};\boldsymbol{D},\boldsymbol{Z}) =\sum_{n=1}^N \log f(\boldsymbol{v}_n\mid \boldsymbol{\mu},\boldsymbol{\Sigma}) +\sum_{n=1}^N \log \mathbb{I}\{\boldsymbol{v}_n\in K_n^{c_n}(\boldsymbol{A})\}. \end{align}\tag{23}\]

Under the proposed sigmoid smoothing approximation with synergy matrix \(\boldsymbol{A}\), the indicator function in E-step is replaced by the sigmoid product approximation in Eq. (18 ). As a result, the hard truncation of the feasible region is replaced by a continuous weighting term. The conditional density used in the E-step therefore becomes \(\mathcal{L}(\boldsymbol{\theta};\boldsymbol{D})\approx \prod_{n=1}^N\int \prod_{j=1}^{J_n} \sigma\left(\frac{\Delta U_{nj}(\boldsymbol{v}; \mathbf{A})}{\lambda}\right) f(\boldsymbol{v}|\boldsymbol{\mu},\boldsymbol{\Sigma}) \mathrm{d}\boldsymbol{v}.\) The complete-data log-likelihood can therefore be written as \[\label{eq:app-likelihood95sigmoid95complete-data} \begin{align} \ell(\boldsymbol{\theta};\boldsymbol{D},\boldsymbol{Z})=\sum_{n=1}^N\left(\log f(\boldsymbol{v}_n \mid \boldsymbol{\mu}, \boldsymbol{\Sigma}) +\sum_{j=1}^{J_n} \log \sigma\left(\frac{\Delta U_{nj}(\boldsymbol{v}_n; \boldsymbol{A})}{\lambda}\right)\right). \end{align}\tag{24}\]

This replaces the truncated Gaussian posterior in the base model with a smoothly weighted Gaussian density, where the sigmoid terms softly enforce the incentive compatibility constraints. Substituting this conditional density into the expectation yields the smoothed \(Q'\) function \[\label{eq:app-smoothed95Q} \begin{align} Q'(\boldsymbol{\theta} \mid \boldsymbol{\theta}^{(t)}) &=\sum_{n=1}^{N}\int q_n(\boldsymbol{v};\boldsymbol{\theta}^{(t)}) \left( \log f(\boldsymbol{v} \mid \boldsymbol{\mu},\boldsymbol{\Sigma})+\sum_{j=1}^{J_n}\log \sigma\left( \frac{\Delta U_{nj}(\boldsymbol{v};\boldsymbol{A})}{\lambda} \right) \right) \mathrm{d}\boldsymbol{v}\\ \end{align},\tag{25}\] where \[q_n(\boldsymbol{v}; \boldsymbol{\theta}^{(t)})= \frac{ f(\boldsymbol{v} \mid \boldsymbol{\mu}^{(t)}, \boldsymbol{\Sigma}^{(t)}) \prod_{j=1}^{J_n}\sigma(\frac{\Delta U_{nj}(\boldsymbol{v}; \boldsymbol{A}^{(t)})}{\lambda}) }{ {\int f(\boldsymbol{\xi} \mid \boldsymbol{\mu}^{(t)}, \boldsymbol{\Sigma}^{(t)})\prod_{j=1}^{J_n}\sigma(\frac{\Delta U_{nj}(\boldsymbol{\xi}; \boldsymbol{A}^{(t)})}{\lambda}) \mathrm{d}\boldsymbol{\xi}} }.\]

Relative to the original EM formulation, the resulting \(Q'\) function differs in two important aspects. First, the truncation of the latent utility distribution is replaced by a smooth probabilistic weighting determined by the sigmoid terms. Second, the log-likelihood now contains an additional contribution from the smoothed IC constraints, which introduces explicit dependence on the synergy parameter \(\boldsymbol{A}\) inside the integrand. This modification enables gradient-based optimization with respect to \(\boldsymbol{A}\) while preserving the overall EM structure.

Similar to the base model, we substituting the Monte Carlo approximation into Eq. (25 ) yields the empirical smoothed \(\hat{Q'}\) function \[\label{eq:app-smoothed95Qhat} \begin{align} \hat{Q'}(\boldsymbol{\theta} \mid \boldsymbol{\theta}^{(t)})&=\sum_{n=1}^N\sum_{l=1}^L \overline{\omega}_{nl}\left( \log f(\boldsymbol{v}_n^{(l)}|\boldsymbol{\mu},\boldsymbol{\Sigma})+\log \omega(\boldsymbol{v}_n^{(l)};\boldsymbol{A}) \right)\\ &=C-\frac{N}{2}\log|\boldsymbol{\Sigma}|+\sum_{n=1}^N\sum_{l=1}^L \overline{\omega}_{nl}\left( -\frac{1}{2}(\boldsymbol{v}_n^{(l)}-\boldsymbol{\mu})^\top\boldsymbol{\Sigma}^{-1}(\boldsymbol{v}_n^{(l)}-\boldsymbol{\mu})+\log \omega(\boldsymbol{v}_n^{(l)};\boldsymbol{A}) \right), \end{align}\tag{26}\] where \(C\) is constant with respect to \(\boldsymbol{\theta}\) and \[\label{eq:app-sigmoid95vals} \begin{align} \omega(\boldsymbol{v};\boldsymbol{A}) &\triangleq\prod_{j=1}^{J_n}\sigma\left( \frac{ \Delta U_{nj}(\boldsymbol{v}; \boldsymbol{A}) }{\lambda} \right),\\ \overline{\omega}_{nl} &\triangleq \frac{ \omega(\boldsymbol{v}_{n}^{(l)};\boldsymbol{A}^{(t)})}{\sum_{l=1}^L \omega(\boldsymbol{v}_{n}^{(l)};\boldsymbol{A}^{(t)})}. \end{align}\tag{27}\]

8.1.0.2 M-step.

In the M-step, we maximize the smoothed Monte Carlo objective \(\hat{Q'}(\boldsymbol{\theta}\mid\boldsymbol{\theta}^{(t)})\) obtained in the E-step with respect to \((\boldsymbol{\mu},\boldsymbol{A},\boldsymbol{\Sigma})\). The objective function is concave in each parameter block.

Update for \(\boldsymbol{\mu}\) and \(\boldsymbol{\Sigma}\). Taking the gradient of \(\hat{Q'}\) with respect to \(\boldsymbol{\mu}\) gives \[\nabla_{\boldsymbol{\mu}}\hat{Q'} = \sum_{n=1}^N\sum_{l=1}^L \overline{\omega}_{nl}\boldsymbol{\Sigma}^{-1} (\boldsymbol{v}_n^{(l)}-\boldsymbol{\mu}).\] The first order condition yields the closed-form update of \(\boldsymbol{\mu}\) and \(\boldsymbol{\Sigma}\) \[\label{eq:mut195sigmoid} \boldsymbol{\mu}^{(t+1)}=\sum_{n=1}^N\sum_{l=1}^L\overline{\omega}_{nl}\boldsymbol{v}_n^{(l)}.\tag{28}\] \[\label{eq:Sigmat195sigmoid} \boldsymbol{\Sigma}^{(t+1)}=\sum_{n=1}^N\sum_{l=1}^L\overline{\omega}_{nl}(\boldsymbol{v}_n^{(l)}-\boldsymbol{\mu}^{(t+1)})(\boldsymbol{v}_n^{(l)}-\boldsymbol{\mu}^{(t+1)})^\top.\tag{29}\]

Gradient with respect to \(\boldsymbol{A}\). Unlike \(\boldsymbol{\mu}\) and \(\boldsymbol{\Sigma}\), the synergy matrix \(\boldsymbol{A}\) does not admit a closed-form solution in the M-step. This is because \(\boldsymbol{A}\) appears inside the nonlinear term \[\log \sigma\left(\frac{\Delta U_{nj}(\boldsymbol{v};\boldsymbol{A})}{\lambda}\right),\] which prevents the first-order optimality condition from yielding an explicit solution for \(\boldsymbol{A}\). Consequently, the maximization with respect to \(\boldsymbol{A}\) must be performed using gradient-based optimization, where the gradient is given in Theorem 2.

In summary, the M-step produces closed-form updates for \(\boldsymbol{\mu}\) and \(\boldsymbol{\Sigma}\), while \(\boldsymbol{A}\) is updated using the gradient of the smoothed objective. Algorithm 9 summarizes the full EM procedure. If the importance sampling framework used in the E-step is applied, the parameter updates can be written in terms of the normalized importance weights \(w_{nl}\).The normalized weights are defined as \[\label{eq:IS95sigmoid95weights} \tilde{w}_{nl}\triangleq\frac{w_{nl}\; \overline{\omega}_{nl} }{\sum_{l=1}^{L} w_{nl}\; \overline{\omega}_{nl}},\tag{30}\] where \(\omega_{nl}\) is the importance sampling density ratio and \(w_{nl}\) is the sigmoid feasibility weight.

\[\label{eq:muSigA95sigmoid95IS} \begin{align} \boldsymbol{\mu}^{(t+1)} &= \sum_{n=1}^{N}\sum_{l=1}^{L}\tilde{w}_{nl}\boldsymbol{v}_n^{(l)} \\ \boldsymbol{\Sigma}^{(t+1)} &= \sum_{n=1}^{N}\sum_{l=1}^{L}\tilde{w}_{nl} (\boldsymbol{v}_n^{(l)}-\boldsymbol{\mu}^{(t+1)})(\boldsymbol{v}_n^{(l)}-\boldsymbol{\mu}^{(t+1)})^\top \\ \nabla_{\boldsymbol{A}}\hat{Q}' &= \frac{1}{\lambda}\sum_{n=1}^{N}\sum_{j=1}^{J_n}\left( \frac{\partial \Delta U_{nj}(\boldsymbol{A})}{\partial \boldsymbol{A}} \sum_{l=1}^{L}\tilde{w}_{nl} \left( 1-\sigma\left(\frac{\Delta U_{nj}(\boldsymbol{v}_n^{(l)};\boldsymbol{A})}{\lambda}\right) \right) \right). \end{align}\tag{31}\]

8.2 Censored Demand↩︎

This appendix provides the detailed derivation of the EM algorithm for the censored demand extension in Section 5.2. We adopt the same notation: \(\boldsymbol{D} = \{N_1,\dots,N_J\}\) denotes the observed data, \(\boldsymbol{Z} = \{N',\boldsymbol{v}_1,\dots,\boldsymbol{v}_{N'}\}\) denotes the missing data, and \(\boldsymbol{\theta}=(\boldsymbol{\mu},\boldsymbol{\Sigma})\) denotes the parameters to be estimated.

8.2.0.1 E-step

We first express the complete-data log-likelihood as \[\ell(\boldsymbol{\theta};\boldsymbol{D}, \boldsymbol{Z}) = \sum_{n=1}^{N'} \log f(\boldsymbol{v}_n\mid \boldsymbol{\mu},\boldsymbol{\Sigma}).\] We then take the expectation of \(\ell(\boldsymbol{\theta};\boldsymbol{D}, \boldsymbol{Z})\) conditional on \(\boldsymbol{D}\) under the current estimate \(\boldsymbol{\theta}^{(t)}\): \[\begin{align} \label{eq:e-step-censored} Q(\boldsymbol{\theta}\mid \boldsymbol{\theta}^{(t)}) &\triangleq \mathbb{E}\left[\ell(\boldsymbol{\theta};\boldsymbol{D},\boldsymbol{Z})\,\Big|\,\boldsymbol{D},\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right] \notag\\ &= \mathbb{E}\left[\mathbb{E}\left[\sum_{n=1}^{N'} \log f(\boldsymbol{v}_n\mid \boldsymbol{\mu},\boldsymbol{\Sigma})\,\Big|\,\boldsymbol{D},\boldsymbol{\theta}^{(t)},N'\right]\,\Big|\,\boldsymbol{D},\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right] \notag\\ &= \mathbb{E}\left[\sum_{n=1}^{N'} \int_{\boldsymbol{v}\in R_n^{c_n}} \frac{f(\boldsymbol{v}\mid \boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)})}{\int_{\boldsymbol{\xi}\in R_n^{c_n}} f(\boldsymbol{\xi}\mid \boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)})\,\mathrm{d}\boldsymbol{\xi}} \log f(\boldsymbol{v}\mid \boldsymbol{\mu},\boldsymbol{\Sigma})\,\mathrm{d}\boldsymbol{v}\,\Big|\,\boldsymbol{D},\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right]. \end{align}\tag{32}\] In the second equality, we apply the tower property by first taking expectation of \(\boldsymbol{v}_n\) conditional on \(N'\), which mirrors the derivation in Section 3.2 for the base model. The outer expectation in 32 is taken with respect to the distribution of \(N'\mid \boldsymbol{D},\boldsymbol{\theta}^{(t)}\), which we derive next.

Under \(\boldsymbol{\theta}^{(t)}\), a customer is censored with probability \(\mathbb{P}(R^0\mid \boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)})\), the probability mass of the IC polyhedron \(R^0\). Therefore, \(N'\mid \boldsymbol{D},\boldsymbol{\theta}^{(t)}\) admits the following probability model: after \(N'\) i.i.d.Bernoulli trials with success probability \(1-\mathbb{P}(R^0\mid \boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)})\) (an uncensored customer counts as a success), exactly \(N\) successes are observed. By Bayes’ rule, for \(n\ge N\), \[\mathbb{P}\!\left(N'=n\,\big|\,\boldsymbol{D},\boldsymbol{\theta}^{(t)}\right) =\mathbb{P}\!\left(N'=n\,\big|\,N,\boldsymbol{\theta}^{(t)}\right) = \frac{\mathbb{P}(N\mid N'=n,\boldsymbol{\theta}^{(t)})\,\mathbb{P}(N'=n)}{\sum_{\tilde{n}=N}^{\infty}\mathbb{P}(N\mid N'=\tilde{n},\boldsymbol{\theta}^{(t)})\,\mathbb{P}(N'=\tilde{n})}.\]

Therefore, the above term can be expressed as \[\begin{align} &\frac{{n \choose N} \mathbb{P}\left(R^0| \boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)^{n-N}\left(1-\mathbb{P}\left(R^0| \boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)\right)^N}{\sum_{\tilde{n}=N}^{+\infty}{\tilde{n} \choose N} \mathbb{P}\left(R^0| \boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)^{\tilde{n}-N}(1-\mathbb{P}(R^0| \boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}))^N}\notag \\ &\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad= {n \choose N} \mathbb{P}\left(R^0| \boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)^{n-N}\left(1-\mathbb{P}\left(R^0| \boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)\right)^{N+1}.\label{eq:neg-binomial} \end{align}\tag{33}\] Note that it has the same distribution as the negative binomial distribution with success rate \(1-\mathbb{P}\left(R^0| \boldsymbol{\mu}^{(t)},\boldsymbol{\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 33 , the expectation in 32 is approximated by Monte Carlo simulation:

1. Consider \(L\) instances. In instance \(l\), generate \(N^{\prime(l)}\) from the negative binomial distribution 33 .

2. Generate \(N_j\) samples from \(\mathcal{N}(\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)})\) conditional on \(\boldsymbol{v}\in \mathcal{P}_\rho\) for \(j=1,\dots,J\), denoted \(\{\boldsymbol{v}^{(l)}_{j,s}\}_{s=1}^{N_j}\), and \(N^{\prime(l)}-N\) samples from \(\mathcal{N}(\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)})\) conditional on \(\boldsymbol{v}\in R^0\), denoted \(\{\boldsymbol{v}^{(l)}_{0,s}\}_{s=1}^{N^{\prime(l)}-N}\).

3. Approximate 32 by \[\label{eq:approx-Q-censor} \hat{Q}(\boldsymbol{\theta}\mid \boldsymbol{\theta}^{(t)}) = \frac{1}{L}\sum_{l=1}^L \left(\sum_{j=1}^J \sum_{s=1}^{N_j}\log f(\boldsymbol{v}_{j,s}^{(l)}\mid \boldsymbol{\mu},\boldsymbol{\Sigma}) + \sum_{s=1}^{N^{\prime(l)}-N}\log f(\boldsymbol{v}_{0,s}^{(l)}\mid \boldsymbol{\mu},\boldsymbol{\Sigma})\right).\tag{34}\]

The information of \(\boldsymbol{\theta}^{(t)}\) has been fully absorbed by the simulated \(\{N'^{(l)}\}\) and \(\{\boldsymbol{v}_{j,s}^{(l)}\}\), so 34 is a function of \((\boldsymbol{\mu},\boldsymbol{\Sigma})\) to be optimized in the M-step.

8.2.0.2 M-step

Given \(\hat{Q}(\boldsymbol{\theta}\mid \boldsymbol{\theta}^{(t)})\) in 34 , the M-step parallels the M-step in the base model. The updates \(\boldsymbol{\mu}^{(t+1)}\) and \(\boldsymbol{\Sigma}^{(t+1)}\) are the sample mean and sample covariance of the simulated valuations: \[\begin{align} \boldsymbol{\mu}^{(t+1)} &= \frac{1}{\sum_{l=1}^L N'^{(l)}}\sum_{l=1}^L \left(\sum_{j=1}^J \sum_{s=1}^{N_j}\boldsymbol{v}_{j,s}^{(l)} + \sum_{s=1}^{N^{\prime(l)}-N}\boldsymbol{v}_{0,s}^{(l)}\right), \\ \boldsymbol{\Sigma}^{(t+1)} &= \frac{1}{\sum_{l=1}^L N'^{(l)}}\sum_{l=1}^L \Bigg(\sum_{j=1}^J \sum_{s=1}^{N_j}(\boldsymbol{v}_{j,s}^{(l)}-\boldsymbol{\mu}^{(t+1)})(\boldsymbol{v}_{j,s}^{(l)}-\boldsymbol{\mu}^{(t+1)})^\top \\ &\qquad\qquad\qquad\qquad\qquad + \sum_{s=1}^{N^{\prime(l)}-N}(\boldsymbol{v}_{0,s}^{(l)}-\boldsymbol{\mu}^{(t+1)})(\boldsymbol{v}_{0,s}^{(l)}-\boldsymbol{\mu}^{(t+1)})^\top\Bigg). \end{align}\] With these E- and M-steps, Algorithm 2 can be adapted directly to censored-demand data.

9 Gaussian Mixture Models↩︎

9.1 Model↩︎

In this section, we generalize the base model in Section 3 by allowing each customer’s valuation vector \(\boldsymbol{v}\) to be drawn independently from an \(I\)-dimensional Gaussian mixture model with \(K\) components: \[\boldsymbol{v}\sim \sum_{k=1}^K \phi_k \mathcal{N}(\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k),\] where \(\phi_k \geq 0\) and \(\sum_{k=1}^{K} \phi_k =1\) are the mixture proportions, \(\mathcal{N}(\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)\) is the density of a multivariate Gaussian distribution with mean \(\boldsymbol{\mu}_k\) and covariance matrix \(\boldsymbol{\Sigma}_k\), and \(\boldsymbol{\theta} = (\phi_1, \cdots, \phi_K, \boldsymbol{\mu}_1, \cdots, \boldsymbol{\mu}_K, \boldsymbol{\Sigma}_1, \cdots, \boldsymbol{\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 approximate a wide range of valuation distributions.

We use \(\boldsymbol{\pi}^n \in \{0,1\}^K\) to denote the component association of customer \(n\): \(\pi_k^n = 1\) if \(\boldsymbol{v}_n\) is generated from \(\mathcal{N}(\boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)\) and zero otherwise. Note that in this model, for each transaction there are two types of latent quantities: (i) the unobserved valuations, and (ii) the mixture component associated with that valuation. Therefore, the missing data are \(\boldsymbol{Z} \triangleq \{ \boldsymbol{\pi}^1, \dots, \boldsymbol{\pi}^N, \boldsymbol{v}_1, \dots, \boldsymbol{v}_{N} \}\). We express the complete-data log-likelihood function as \[\ell(\boldsymbol{\theta}; \boldsymbol{D}, \boldsymbol{Z} ) = \sum_{n=1}^N \sum_{k=1}^K \boldsymbol{\pi}_k^n [\log (\boldsymbol{\phi}_k) +\log f(\boldsymbol{v}_n |\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)],\] where \(f(\cdot)\) denotes the PDF of the multivariate normal distribution.

Compared with the base model in Section 3, the conditional expectation of the complete-data log-likelihood now involves an outer expectation over the component memberships \(\boldsymbol{\pi}^n\mid \boldsymbol{D}, \boldsymbol{\theta}^{(t)}\). Given the choice \(c_n\) of customer \(n\), the probability that customer \(n\) is associated with component \(k\) is, by Bayes’ rule, proportional to \(\boldsymbol{\phi}_k^{(t)}\,\mathbb{P}(R_n^{c_n}\mid \boldsymbol{\mu}_k^{(t)},\boldsymbol{\Sigma}_k^{(t)},\pi_k^n=1)\), where \(\mathbb{P}(R_n^{c_n}\mid \boldsymbol{\mu}_k^{(t)},\boldsymbol{\Sigma}_k^{(t)},\pi_k^n=1)\) is the mass that component \(k\) assigns to the IC polyhedron \(R_n^{c_n}\). Intuitively, components that place more probability mass on the IC polyhedron consistent with the observed choice \(c_n\) are more likely to have generated customer \(n\).

Once the posterior of \(\boldsymbol{\pi}^n\) is identified, the E-step is approximated by Monte Carlo simulation as in Section 3.2: in each instance, we first estimate the membership posterior \(\hat{\pi}_k^n\) by drawing samples from each component \(\mathcal{N}(\boldsymbol{\mu}_k^{(t)},\boldsymbol{\Sigma}_k^{(t)})\) and compute the fraction falling in \(R_n^{c_n}\), then draw customer valuations from truncated Gaussians on the corresponding IC polyhedra under each component. The M-step can then be adapted using a \(\hat{Q}\) function computed from the simulated samples above, and the parameters \((\boldsymbol{\phi},\boldsymbol{\mu},\boldsymbol{\Sigma})\) are updated by maximizing this \(\hat{Q}\), which yields a weighted version of the standard MLE for Gaussian mixtures. In this way, Algorithm 2 can be adapted to GMM with only minor modifications. The detailed derivation of the E-step, the closed form for the membership posterior \(\hat{\pi}_k^n\), the Monte Carlo procedure, and the resulting M-step update formulas are provided in Appendix 9.

9.2 EM Algorithm↩︎

9.2.0.1 E-step.

We first express the complete-data log-likelihood as \[\ell(\boldsymbol{\theta};\boldsymbol{D},\boldsymbol{Z}) = \sum_{n=1}^N \sum_{k=1}^K \pi_k^n\big[\log\phi_k + \log f(\boldsymbol{v}_n\mid \boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)\big],\] where \(f(\cdot)\) is the PDF of the multivariate normal distribution. We then take the expectation of \(\ell(\boldsymbol{\theta};\boldsymbol{D}, \boldsymbol{Z})\) conditional on \(\boldsymbol{D}\) under the current estimate \(\boldsymbol{\theta}^{(t)}\): \[\begin{align} \label{eq:GMM-Q} Q(\boldsymbol{\theta}\mid \boldsymbol{\theta}^{(t)}) &\triangleq \mathbb{E}_{\boldsymbol{Z}\mid \boldsymbol{D},\boldsymbol{\theta}^{(t)}}\!\left[\ell(\boldsymbol{\theta};\boldsymbol{D},\boldsymbol{Z})\,\Big|\,\boldsymbol{D},\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)},\boldsymbol{\phi}^{(t)}\right] \notag\\ &= \mathbb{E}_{\boldsymbol{\pi}\mid \boldsymbol{D},\boldsymbol{\theta}^{(t)}}\!\left[\mathbb{E}_{\boldsymbol{v}\mid \boldsymbol{D},\boldsymbol{\theta}^{(t)},\boldsymbol{\pi}}\!\left[\sum_{n=1}^N\sum_{k=1}^K \pi_k^n\log\phi_k + \pi_k^n\log f(\boldsymbol{v}_n\mid \boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)\,\Big|\,\boldsymbol{D},\boldsymbol{\theta}^{(t)}\right]\right] \notag\\ &= \sum_{n=1}^N \sum_{k=1}^K \mathbb{P}(\pi_k^n=1\mid \boldsymbol{D},\boldsymbol{\theta}^{(t)})\,\mathbb{E}\!\left[\log\phi_k + \log f(\boldsymbol{v}_n\mid \boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)\,\Big|\,\boldsymbol{D},\boldsymbol{\theta}^{(t)},\pi_k^n=1\right], \end{align}\tag{35}\] where the inner conditional expectation is \[\mathbb{E}\!\left[\log f(\boldsymbol{v}_n\mid \boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)\,\Big|\,\boldsymbol{D},\boldsymbol{\theta}^{(t)},\pi_k^n=1\right] = \int_{\boldsymbol{v}\in R_n^{c_n}}\frac{f(\boldsymbol{v}\mid \boldsymbol{\mu}_k^{(t)},\boldsymbol{\Sigma}_k^{(t)})}{\int_{\boldsymbol{\xi}\in R_n^{c_n}} f(\boldsymbol{\xi}\mid \boldsymbol{\mu}_k^{(t)},\boldsymbol{\Sigma}_k^{(t)})\,\mathrm{d}\boldsymbol{\xi}}\log f(\boldsymbol{v}\mid \boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)\,\mathrm{d}\boldsymbol{v}.\]

To evaluate 35 , we first derive the posterior membership probability. Given the choice \(c_n\), by Bayes’ rule, \[\label{eq:posterior} \hat{\pi}_k^n \triangleq \mathbb{P}(\pi_k^n=1\mid \boldsymbol{D},\boldsymbol{\theta}^{(t)}) = \frac{\boldsymbol{\phi}_k^{(t)}\,\mathbb{P}(R_n^{c_n}\mid \boldsymbol{\mu}_k^{(t)},\boldsymbol{\Sigma}_k^{(t)},\pi_k^n=1)}{\sum_{k'=1}^K \boldsymbol{\phi}_{k'}^{(t)}\,\mathbb{P}(R_n^{c_n}\mid \boldsymbol{\mu}_{k'}^{(t)},\boldsymbol{\Sigma}_{k'}^{(t)},\pi_{k'}^n=1)}.\tag{36}\] The polyhedral mass \(\mathbb{P}(R_n^{c_n}\mid \boldsymbol{\mu}_k^{(t)},\boldsymbol{\Sigma}_k^{(t)},\pi_k^n=1)\) in 36 can be estimated by Monte Carlo simulation: generate \(L\) samples from \(\mathcal{N}(\boldsymbol{\mu}_k^{(t)},\boldsymbol{\Sigma}_k^{(t)})\) and compute the fraction falling in \(R_n^{c_n}\).

The expectation in 35 is then approximated by the following Monte Carlo procedure:

1. For each component \(k=1,\dots,K\), generate \(L\) samples \(\boldsymbol{v}_k^{(l)}\), \(l=1,\dots,L\), from \(\mathcal{N}(\boldsymbol{\mu}_k^{(t)},\boldsymbol{\Sigma}_k^{(t)})\) and use them to estimate the posterior membership probabilities \(\hat{\pi}_k^n\) via 36 : \[\hat{\mathbb{P}}(R_n^{c_n}\mid \boldsymbol{\mu}_k^{(t)},\boldsymbol{\Sigma}_k^{(t)}) = \frac{1}{L}\sum_{l=1}^L \mathbb{I}\{\boldsymbol{v}_k^{(l)}\in R_n^{c_n}\}.\]

2. For all \(n=1,\dots,N\) and \(k=1,\dots,K\), generate \(L\) samples \(\boldsymbol{v}_{nk}^{(l)}\), \(l=1,\dots,L\), from \(\mathcal{N}(\boldsymbol{\mu}_k^{(t)},\boldsymbol{\Sigma}_k^{(t)})\) conditional on \(\boldsymbol{v}_{nk}^{(l)}\in R_n^{c_n}\), using Algorithm 1.

3. Approximate 35 by \[\label{eq:Estep-gmm} \hat{Q}(\boldsymbol{\theta}\mid \boldsymbol{\theta}^{(t)}) = \sum_{n=1}^N \sum_{k=1}^K \hat{\pi}_k^n \log\phi_k + \frac{1}{L}\sum_{l=1}^L\sum_{n=1}^N\sum_{k=1}^K \hat{\pi}_k^n \log f(\boldsymbol{v}_{nk}^{(l)}\mid \boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k).\tag{37}\]

9.2.0.2 M-step.

Given \(\hat{Q}(\boldsymbol{\theta}\mid \boldsymbol{\theta}^{(t)})\) in 37 , we maximize over \(\boldsymbol{\theta}=(\boldsymbol{\phi},\boldsymbol{\mu},\boldsymbol{\Sigma})\). For the mixture proportions, applying the constraint \(\sum_k\phi_k=1\) to the first term in 37 yields \[\phi_k^{(t+1)} = \frac{1}{N}\sum_{n=1}^N \hat{\pi}_k^n.\] For the component means and covariances, the second term in 37 is a weighted Gaussian log-likelihood with weights \(\hat{\pi}_k^n\), so the updates are weighted versions of the standard MLE: \[\begin{align} \boldsymbol{\mu}_k^{(t+1)} &= \frac{\sum_{l=1}^L\sum_{n=1}^N \hat{\pi}_k^n\,\boldsymbol{v}_{nk}^{(l)}}{L\sum_{n=1}^N \hat{\pi}_k^n}, \\ \boldsymbol{\Sigma}_k^{(t+1)} &= \frac{\sum_{l=1}^L\sum_{n=1}^N \hat{\pi}_k^n\,(\boldsymbol{v}_{nk}^{(l)}-\boldsymbol{\mu}_k^{(t+1)})(\boldsymbol{v}_{nk}^{(l)}-\boldsymbol{\mu}_k^{(t+1)})^\top}{L\sum_{n=1}^N \hat{\pi}_k^n}. \end{align}\] With these E-steps and M-steps, Algorithm 2 can be adapted directly to the Gaussian mixture setting.

9.3 Numerical Experiments↩︎

In this experiment, we study the convergence of the EM algorithm for the Gaussian mixture model with two components and compare its performance with that of the base EM algorithm. We generate synthetic datasets using a procedure similar to that in the base-model experiments. Specifically, each consumer’s product valuation vector is drawn from a two-component Gaussian mixture distribution. The mixture weight \(\phi\) is randomly generated, while the remaining experimental settings follow the base-model experiments. To systematically evaluate the model, we fix the number of products at \(I=2\) and the transaction size at \(N=2500\). We generate five datasets with different parameters \((\boldsymbol{\mu},\boldsymbol{\Sigma})\) and estimate each dataset using both the GMM and the base model.

Notes:

1. Log-likelihood Score is defined as \(1 - (\mathcal L_{\text{model}} - \mathcal L_{\text{exact}})/\mathcal L_{\text{exact}}\), so a value closer to \(100\%\) indicates that the fitted model is closer to the oracle test log-likelihood.

As shown in Figure [fig:gmm95llh], the GMM model provides a substantially better fit to the observed data compared to the single-component baseline model, which fails to capture the underlying heterogeneity. Table 7 shows that the GMM consistently achieves lower RMSE and higher test log-likelihood scores than the single-component base EM model. This improvement is expected because the data are generated from a two-component mixture, whereas the base model imposes a single Gaussian distribution. The results highlight the value of the GMM extension when the underlying valuation distribution contains latent customer heterogeneity.

10 Hierarchical Bayesian Model↩︎

In this section, we introduce the Metropolis-Hastings algorithm for estimating the posterior distribution. We assume a normal proposal distribution and uniform prior distributions for both \(\boldsymbol{\mu}\) and a covariance factor \(\boldsymbol{S}\), where \(\boldsymbol{\Sigma}=\boldsymbol{S}\boldsymbol{S}^\top\).

11 Multinomial Logit Model for Bundle Choice↩︎

This appendix describes the MNL model used in the numerical experiments. The purpose of this benchmark is to evaluate how a classical discrete-choice model performs when the true data-generating process follows the proposed Gaussian bundle-choice model. The resulting estimator therefore serves as a misspecified benchmark against the proposed Gaussian-based estimator.

Suppose there are \(I\) products. The true latent product valuation vector is generated from a multivariate Gaussian distribution: \(\boldsymbol{v}_n\sim \mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})\), where \(\boldsymbol{\Sigma}\) is a diagonal matrix. Although the data are generated from a Gaussian model, the MNL estimator assumes the standard multinomial logit structure. Specifically, the MNL model assumes the utility specification \[\begin{align} U_{n0}^{\text{MNL}} &= \varepsilon_{n0}, \\ U_{nj}^{\text{MNL}} &= \boldsymbol{x}_{nj}^\top \boldsymbol{\mu}- p_{nj} + \varepsilon_{nj}, \qquad j=1,\dots,J_n, \end{align}\] where \(\boldsymbol{\mu}\in\mathbb{R}^I\) is the vector of product-level mean valuations to be estimated, and \(\varepsilon_{nj} \stackrel{\text{i.i.d.}}{\sim} \mathrm{Gumbel}(0,1).\) Therefore, under the MNL model, bundle utilities are simply the sums of constituent product mean valuations plus an independent Gumbel error term.

Under the MNL assumptions, the standard logit formula gives: \[\begin{align} \mathbb{P}(c_n = 0) &= \frac{1}{1+\sum_{j=1}^{J_n}\exp\!\left(\boldsymbol{x}_{nj}^\top \boldsymbol{\mu}- p_{nj}\right)},\\ \mathbb{P}(c_n = j) &= \frac{\exp\!\left(\boldsymbol{x}_{nj}^\top \boldsymbol{\mu} - p_{nj}\right)}{1+\sum_{j'=1}^{J_n}\exp\!\left(\boldsymbol{x}_{nj'}^\top \boldsymbol{\mu}- p_{nj'}\right)}. \end{align}\]

The MNL benchmark estimates only the product-level mean valuations \(\boldsymbol{\mu}\) by maximizing the log-likelihood implied by the logit probabilities: \[\widehat{\boldsymbol{\mu}}_{\mathrm{MNL}} = \arg\max_{\boldsymbol{\mu}} \sum_{n=1}^N \log \mathbb{P}(c_n \mid \boldsymbol{\mu}).\] Optimization is performed using the L-BFGS-B algorithm.

12 Proofs↩︎

12.1 The Derivation of the M-step↩︎

\[\begin{align} Q\left(\boldsymbol{\theta}|\boldsymbol{\theta}^{(t)}\right)&\triangleq \mathbb{E}\left[l(\boldsymbol{\theta};\boldsymbol{D},\boldsymbol{Z})\bigg|\boldsymbol{D},\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right]\notag\\ &= \sum_{n=1}^{N} \int_{\boldsymbol{v}\in R_n^{c_n}} \frac{f\left(\boldsymbol{v}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)}{\int_{R_n^{c_n}}f\left(\boldsymbol{\xi}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)\mathrm{d}\boldsymbol{\xi} } \log f(\boldsymbol{v}|\boldsymbol{\mu},\boldsymbol{\Sigma})\mathrm{d}\boldsymbol{v}. \end{align}\] To find a new estimation for \(\boldsymbol{\theta}\), the M-step in the EM algorithm maximizes \(Q\left(\boldsymbol{\theta}|\boldsymbol{\theta}^{(t)}\right)\) and let \(\boldsymbol{\theta}^{(t+1)}=\argmax_{\boldsymbol{\theta}}Q\left(\boldsymbol{\theta}|\boldsymbol{\theta}^{(t)}\right)\). Thus, we first solve \[\begin{align} \partial \mathbb{E}\left[l(\boldsymbol{\theta};\boldsymbol{D},\boldsymbol{Z})\mid\boldsymbol{D},\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right]/\partial \boldsymbol{\mu}= 0, \label{eq:deriv-mu} \end{align}\tag{38}\] and then we solve \[\begin{align} \partial \mathbb{E}\left[l(\boldsymbol{\theta};\boldsymbol{D},\boldsymbol{Z})\mid\boldsymbol{D},\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right]/\partial \boldsymbol{\Sigma} = 0. \label{eq:deriv-sigma} \end{align}\tag{39}\] Equality 38 equates to \[\begin{align} \frac{\partial \mathbb{E}\left[l(\boldsymbol{\theta};\boldsymbol{D},\boldsymbol{Z})\mid\boldsymbol{D},\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right]}{\partial \boldsymbol{\mu}} =& \frac{\partial }{\partial \boldsymbol{\mu}} \sum_{n=1}^{N} \int_{\boldsymbol{v}\in R_n^{c_n}} \frac{f\left(\boldsymbol{v}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)}{\int_{R_n^{c_n}}f(\boldsymbol{D}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)})\mathrm{d}\boldsymbol{\xi} } \log f(\boldsymbol{v}|\boldsymbol{\mu},\boldsymbol{\Sigma})\mathrm{d}\boldsymbol{v} \notag \\ =&\sum_{n=1}^{N} \int_{\boldsymbol{v}\in R_n^{c_n}} \frac{f\left(\boldsymbol{v}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)}{\int_{R_n^{c_n}}f\left(\boldsymbol{\xi}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)\mathrm{d}\boldsymbol{\xi} } \frac{\partial }{\partial \boldsymbol{\mu}} \log f(\boldsymbol{v}|\boldsymbol{\mu},\boldsymbol{\Sigma})\mathrm{d}\boldsymbol{v}\notag\\ =&\sum_{n=1}^{N} \int_{\boldsymbol{v}\in R_n^{c_n}} \frac{f\left(\boldsymbol{v}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)}{\int_{R_n^{c_n}}f\left(\boldsymbol{\xi}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)\mathrm{d}\boldsymbol{\xi} } \frac{1}{f(\boldsymbol{v}|\boldsymbol{\mu},\boldsymbol{\Sigma})} \boldsymbol{\Sigma}^{-1}(\boldsymbol{v} - \boldsymbol{\mu}) f(\boldsymbol{v}|\boldsymbol{\mu},\boldsymbol{\Sigma}) \mathrm{d}\boldsymbol{v} \notag\\ = &\sum_{n=1}^{N} \int_{\boldsymbol{v}\in R_n^{c_n}} \boldsymbol{\Sigma}^{-1}(\boldsymbol{v} - \boldsymbol{\mu}) \frac{f\left(\boldsymbol{v}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)}{\int_{R_n^{c_n}}f\left(\boldsymbol{\xi}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)\mathrm{d}\boldsymbol{\xi} } \mathrm{d}\boldsymbol{v} \notag\\ =&\sum_{n=1}^{N} \int_{\boldsymbol{v}\in R_n^{c_n}} \boldsymbol{v} \frac{f\left(\boldsymbol{v}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)}{\int_{R_n^{c_n}}f\left(\boldsymbol{\xi}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)\mathrm{d}\boldsymbol{\xi} } \mathrm{d}\boldsymbol{v} \notag \\ &- \sum_{n=1}^{N} \int_{\boldsymbol{v}\in R_n^{c_n}} \boldsymbol{\mu} \frac{f\left(\boldsymbol{v}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)}{\int_{R_n^{c_n}}f\left(\boldsymbol{\xi}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)\mathrm{d}\boldsymbol{\xi} } \mathrm{d}\boldsymbol{v} \notag\\ =&\sum_{n=1}^{N} \int_{\boldsymbol{v}\in R_n^{c_n}} \boldsymbol{v} \frac{f\left(\boldsymbol{v}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)}{\int_{R_n^{c_n}}f\left(\boldsymbol{\xi}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)\mathrm{d}\boldsymbol{\xi} } \mathrm{d}\boldsymbol{v} - N \boldsymbol{\mu} =0. \label{eq:Q-firstorder-mu} \end{align}\tag{40}\] Thus, we have \[\begin{align} \boldsymbol{\mu} = \frac{\sum_{n=1}^{N} \int_{\boldsymbol{v}\in R_n^{c_n}} \boldsymbol{v} \frac{f\left(\boldsymbol{v}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)}{\int_{R_n^{c_n}}f\left(\boldsymbol{\xi}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)\mathrm{d}\boldsymbol{\xi} } \mathrm{d}\boldsymbol{v}}{ N} = \frac{\sum_{n=1}^{N} \boldsymbol{\mathbb{E}}[\boldsymbol{v} |\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}, R_n^{c_n}]}{N}. \end{align}\] Equality 39 equates to \[\begin{align} \frac{\partial \mathbb{E}\left[l(\boldsymbol{\theta};\boldsymbol{D},\boldsymbol{Z})\mid\boldsymbol{D},\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right]}{\partial \boldsymbol{\Sigma}} &=\sum_{n=1}^{N}\frac{\partial }{\partial \boldsymbol{\Sigma}} \int_{\boldsymbol{v}\in R_n^{c_n}} \frac{f\left(\boldsymbol{v}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)}{\int_{R_n^{c_n}}f(\boldsymbol{\xi}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)})\mathrm{d}\boldsymbol{\xi} } \log f(\boldsymbol{v}|\boldsymbol{\mu},\boldsymbol{\Sigma}) \mathrm{d}\boldsymbol{v} \notag \\ &= \sum_{n=1}^{N} \int_{\boldsymbol{v}\in R_n^{c_n}} \frac{f\left(\boldsymbol{v}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)}{\int_{R_n^{c_n}}f(\boldsymbol{\xi}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)})\mathrm{d}\boldsymbol{\xi} }(-\frac{1}{2} \boldsymbol{\Sigma}^{-1} + \frac{1}{2} \boldsymbol{\Sigma}^{-1} (\boldsymbol{v}- \boldsymbol{\mu}) (\boldsymbol{v}- \boldsymbol{\mu}) ^\top \boldsymbol{\Sigma}^{-1}) \mathrm{d}\boldsymbol{v} \label{eq:Q-firstorder-sigma} \\ &=-\frac{1}{2} \boldsymbol{\Sigma}^{-1} \sum_{n=1}^{N} \int_{\boldsymbol{v}\in R_n^{c_n}} \frac{f\left(\boldsymbol{v}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)}{\int_{R_n^{c_n}}f(\boldsymbol{\xi}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)})\mathrm{d}\boldsymbol{\xi} } \mathrm{d}\boldsymbol{v} \notag \\ &\quad +\frac{1}{2} \boldsymbol{\Sigma}^{-1}\sum_{n=1}^{N} \int_{\boldsymbol{v}\in R_n^{c_n}} \frac{f\left(\boldsymbol{v}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)}{\int_{R_n^{c_n}}f(\boldsymbol{\xi}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)})\mathrm{d}\boldsymbol{\xi} } (\boldsymbol{v}- \boldsymbol{\mu}) (\boldsymbol{v}- \boldsymbol{\mu}) ^\top \boldsymbol{\Sigma}^{-1} \mathrm{d}\boldsymbol{v} = 0. \notag \end{align}\tag{41}\] Thus, we have \[\begin{align} \boldsymbol{\Sigma} &= \frac{ \sum_{n=1}^{N} \int_{\boldsymbol{v}\in R_n^{c_n}} (\boldsymbol{v}- \boldsymbol{\mu}) (\boldsymbol{v}- \boldsymbol{\mu}) ^\top \frac{f\left(\boldsymbol{v}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)}\right)}{\int_{R_n^{c_n}}f(\boldsymbol{\xi}|\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)})\mathrm{d}\boldsymbol{\xi} } \mathrm{d}\boldsymbol{v} }{N} \end{align}\]

12.2 Identifiability↩︎

We divide the proof into two parts.

12.2.0.1 Necessity.

We first show that having at least two distinct price points for each product is necessary. Suppose, to the contrary, that there exists a product \(i\) that is observed only at a single price \(p_i\) under all separate-selling menus. We construct two distinct parameter vectors that generate exactly the same region probabilities.

Fix the parameters of all products other than \(i\), and assume that product \(i\) is independent of the remaining products, that is, \[\Sigma_{ij}=0 \qquad \text{for all } j \neq i.\] Let \(\mathbb{P}_i\triangleq\mathbb{P}(_i \le p_i)\). Since only one threshold \(p_i\) is observed for product \(i\), the observable information about its marginal distribution is reduced to the single equation \[\mathbb{P}_i = \Phi\!\left(\frac{p_i-\mu_i}{\sqrt{\Sigma_{ii}}}\right).\] For any scalar \(\tau>0\), define \[\mu_i(\tau) \triangleq p_i - \tau \Phi^{-1}(\mathbb{P}_i), \qquad \Sigma_{ii}(\tau)\triangleq \tau^2.\] Then \[\Phi\!\left(\frac{p_i-\mu_i(\tau)}{\sqrt{\Sigma_{ii}(\tau)}}\right) = \mathbb{P}_i \qquad \text{for every } \tau>0.\] Choose two distinct values \(\tau_1 \neq \tau_2\), and let \(\boldsymbol{\theta}^{(1)}\) and \(\boldsymbol{\theta}^{(2)}\) be the corresponding full parameter vectors obtained by replacing only \((\mu_i,\Sigma_{ii})\) with \((\mu_i(\tau_1),\Sigma_{ii}(\tau_1))\) and \((\mu_i(\tau_2),\Sigma_{ii}(\tau_2))\), respectively, while keeping the remaining coordinates of \(\boldsymbol{\mu}\) and \(\boldsymbol{\Sigma}\) fixed.

Because product \(i\) is independent of the remaining products and appears only through the single threshold \(p_i\), the probability of every region induced by the separate-selling menus factors into a term depending on \(\mathbb{P}_i\) and a term depending only on the remaining products. Therefore \(\boldsymbol{\theta}^{(1)}\) and \(\boldsymbol{\theta}^{(2)}\) induce exactly the same collection of region probabilities, although \(\boldsymbol{\theta}^{(1)} \neq \boldsymbol{\theta}^{(2)}\). This contradicts identifiability.

Hence, at least two distinct price points for each product are necessary.

12.2.0.2 Sufficiency.

We now show that the menu design in Proposition 1 is sufficient for identification.

12.2.0.3 Step 1: identification of the diagonal terms.

Fix any product \(i\in\{1,\dots,I\}\). Under the regular menu \(\boldsymbol{p}^{(0)}\) and the corresponding single-sale menu \(\boldsymbol{p}^{(i)}\), product \(i\) is offered at two distinct price levels \(p_i < p_i'\). In the separate-selling context, the observable data allow us to estimate the probability mass of the random vector \(V_i \sim \mathcal{N}(\mu_i,\Sigma_{ii})\) over the intervals induced by these two prices. Let \[\mathcal{P}_{i,1} \triangleq (-\infty,p_i], \qquad \mathcal{P}_{i,2} \triangleq(p_i,p_i'], \qquad \mathcal{P}_{i,3} \triangleq (p_i',+\infty).\] By observing the fraction of customers who do not purchase product \(i\) at prices \(p_i\) and \(p_i'\), we can consistently estimate \[\mathbb{P}_i \triangleq \mathbb{P}\left(V_i\in \mathcal{P}_{i,1} \right), \qquad \mathbb{P}_i' \triangleq \mathbb{P}\left(V_i\in \left( \mathcal{P}_{i,1}\cup \mathcal{P}_{i,2} \right) \right).\] Since \(p_i<p_i'\), we have \(\mathcal{P}_{i,1}\subsetneq \left( \mathcal{P}_{i,1}\cup \mathcal{P}_{i,2} \right)\), and because the Gaussian density is strictly positive on \(\mathbb{R}\), it follows that \[\mathbb{P}_i < \mathbb{P}_i'.\] The identifiability of \((\mu_i,\Sigma_{ii})\) is equivalent to showing that the mapping from \((\mu_i,\Sigma_{ii})\) to the observed probabilities \((\mathbb{P}_i,\mathbb{P}_i')\) is injective. By the Gaussian CDF \(\Phi(\cdot)\), \[\mathbb{P}_i = \Phi\!\left(\frac{p_i-\mu_i}{\sqrt{\Sigma_{ii}}}\right), \qquad \mathbb{P}_i' = \Phi\!\left(\frac{p_i'-\mu_i}{\sqrt{\Sigma_{ii}}}\right).\] Since \(\Phi(\cdot)\) is a strictly increasing bijection from \(\mathbb{R}\) to \((0,1)\) and \(\mathbb{P}_i < \mathbb{P}_i'\), the values \(\Phi^{-1}(\mathbb{P}_i)\) and \(\Phi^{-1}(\mathbb{P}_i')\) are uniquely determined by the observed data. Therefore the system can be rearranged as \[\begin{cases} \mu_i + \Phi^{-1}(\mathbb{P}_i)\sqrt{\Sigma_{ii}} = p_i,\\ \mu_i + \Phi^{-1}(\mathbb{P}_i')\sqrt{\Sigma_{ii}} = p_i'. \end{cases}\] Equivalently, \[\begin{bmatrix} 1 & \Phi^{-1}(\mathbb{P}_i)\\ 1 & \Phi^{-1}(\mathbb{P}_i') \end{bmatrix} \begin{bmatrix} \mu_i\\ \sqrt{\Sigma_{ii}} \end{bmatrix} = \begin{bmatrix} p_i\\ p_i' \end{bmatrix}.\] To establish identifiability, we verify that this linear system has a unique solution. Consider the coefficient matrix \[\mathcal{A}_i \triangleq \begin{bmatrix} 1 & \Phi^{-1}(\mathbb{P}_i)\\ 1 & \Phi^{-1}(\mathbb{P}_i') \end{bmatrix}\] and the augmented matrix \[\left[\;\mathcal{A}_i\;|\;\boldsymbol{b}_i\;\right] \triangleq \left[ \begin{array}{cc|c} 1 & \Phi^{-1}(\mathbb{P}_i) & p_i\\ 1 & \Phi^{-1}(\mathbb{P}_i') & p_i' \end{array} \right],\] where \(\boldsymbol{b}_i=(p_i,p_i')^\top\). First, \[\det(\mathcal{A}_i)=\Phi^{-1}(\mathbb{P}_i')-\Phi^{-1}(\mathbb{P}_i)\neq 0,\] because \(\mathbb{P}_i<\mathbb{P}_i'\) and \(\Phi^{-1}\) is strictly increasing. Hence \[\operatorname{rank}(\mathcal{A}_i)=2.\] Since the augmented matrix has the same number of rows, its rank cannot exceed \(2\), and because it contains \(\mathcal{A}_i\) as a submatrix, its rank is also \(2\). Therefore \[\operatorname{rank}(\mathcal{A}_i)=\operatorname{rank}([\mathcal{A}_i\mid \boldsymbol{b}_i])=2.\] Thus the linear system has a unique solution. Hence \(\mu_i\) and \(\Sigma_{ii}\) are uniquely identified.

Since \(i\) was arbitrary, all entries of \(\boldsymbol{\mu}\) and all diagonal entries of \(\boldsymbol{\Sigma}\) are identified.

12.2.0.4 Step 2: identification of the off-diagonal terms.

Fix any pair \(1 \le i < j \le I\). Under the regular menu, the joint no-purchase probability \[\mathbb{P}_{ij} \triangleq \mathbb{P}(V_i \le p_i,\;V_j \le p_j)\] is observed from the data.

By Step 1, \(\mu_i,\mu_j,\Sigma_{ii},\Sigma_{jj}\) are already known. Define \[c_i \triangleq \frac{p_i-\mu_i}{\sqrt{\Sigma_{ii}}}, \qquad c_j \triangleq \frac{p_j-\mu_j}{\sqrt{\Sigma_{jj}}}, \qquad \rho_{ij} \triangleq \frac{\Sigma_{ij}}{\sqrt{\Sigma_{ii}\Sigma_{jj}}}.\] Then \[\mathbb{P}_{ij} = \Phi_2(c_i,c_j;\rho_{ij}),\] where \(\Phi_2(\cdot,\cdot;\rho)\) denotes the CDF of the standard bivariate normal distribution with correlation \(\rho\).

For fixed \((c_i,c_j)\), the map \[\rho \mapsto \Phi_2(c_i,c_j;\rho)\] is strictly increasing on \((-1,1)\). Indeed, by Plackett’s identity, \[\frac{\partial}{\partial \rho}\Phi_2(c_i,c_j;\rho)=\phi_2(c_i,c_j;\rho)>0,\] where \(\phi_2(\cdot,\cdot;\rho)\) is the corresponding bivariate normal density. Therefore \(\mathbb{P}_{ij}\) uniquely determines \(\rho_{ij}\), and hence \[\Sigma_{ij} = \rho_{ij}\sqrt{\Sigma_{ii}\Sigma_{jj}}\] is uniquely determined.

Since the pair \((i,j)\) was arbitrary, every off-diagonal entry \(\Sigma_{ij}\) is identified.

12.2.0.5 Conclusion.

Step 1 identifies all entries of \(\boldsymbol{\mu}\) and all diagonal entries of \(\boldsymbol{\Sigma}\). Step 2 identifies all off-diagonal entries of \(\boldsymbol{\Sigma}\). Therefore the full parameter pair \((\boldsymbol{\mu},\boldsymbol{\Sigma})\) is uniquely determined by the probabilities of the regions induced by the \(I+1\) separate-selling menus. This proves sufficiency.

Combining the necessity and sufficiency parts, we conclude that under separate selling, having at least two distinct price points for each product is necessary, and the menu design in Proposition 1 is sufficient, for the identification of \((\boldsymbol{\mu},\boldsymbol{\Sigma})\).

12.3 Proofs for Section 4.2: EM Convergence↩︎

12.3.1 Preliminaries↩︎

Recall the definition of the population-level \(Q\)-function: \[\label{eq:pop-Q} Q(\boldsymbol{\mu}'\mid \boldsymbol{\mu}) = \sum_{\rho=1}^{P} \left(\int_{\boldsymbol{\xi} \in {\mathcal{P}_\rho}} f(\boldsymbol{\xi} \mid \boldsymbol{\mu}^*) \, \mathrm d\boldsymbol{\xi} \right) \int_{\boldsymbol{v} \in {\mathcal{P}_\rho}} \frac{f(\boldsymbol{v} \mid \boldsymbol{\mu})}{\int_{\boldsymbol{\xi} \in {\mathcal{P}_\rho}} f(\boldsymbol{\xi} \mid \boldsymbol{\mu}) \, \mathrm d\boldsymbol{\xi}} \log f(\boldsymbol{v} \mid \boldsymbol{\mu}') \, \mathrm{d}\boldsymbol{v}.\tag{42}\] By definition, \(q(M(\boldsymbol{\mu})) = Q(M(\boldsymbol{\mu})\mid \boldsymbol{\mu}^*)\). Therefore, \[\begin{align} q(M(\boldsymbol{\mu})) &= Q(M(\boldsymbol{\mu}) \mid \boldsymbol{\mu}^*) \notag \\ &= \sum_{\rho=1}^{P} \int_{\boldsymbol{v} \in {\mathcal{P}_\rho}} f(\boldsymbol{v} \mid \boldsymbol{\mu}^*) \log f(\boldsymbol{v} \mid M(\boldsymbol{\mu})) \, \mathrm{d}\boldsymbol{v}\notag \\ &= \mathbb{E}_{\boldsymbol{\mu}^*}\left[\log f(\boldsymbol{V} \mid M(\boldsymbol{\mu}))\right]. \label{q-func-simp} \end{align}\tag{43}\] We rewrite \(Q(M(\boldsymbol{\mu})\mid \boldsymbol{\mu})\) as \[\begin{align} Q(M(\boldsymbol{\mu})\mid \boldsymbol{\mu}) &= \sum_{\rho=1}^{P} \left(\int_{\boldsymbol{\xi} \in {\mathcal{P}_\rho}} f(\boldsymbol{\xi} \mid \boldsymbol{\mu}^*) \, \mathrm d\boldsymbol{\xi} \right) \int_{\boldsymbol{v} \in {\mathcal{P}_\rho}} \frac{f(\boldsymbol{v} \mid \boldsymbol{\mu})}{\int_{\boldsymbol{\xi} \in {\mathcal{P}_\rho}} f(\boldsymbol{\xi} \mid \boldsymbol{\mu}) \, \mathrm d\boldsymbol{\xi}} \log f(\boldsymbol{v} \mid M(\boldsymbol{\mu})) \, \mathrm{d}\boldsymbol{v}\notag \\ &= \sum_{\rho=1}^{P} \frac{\int_{\boldsymbol{\xi} \in {\mathcal{P}_\rho}} f(\boldsymbol{\xi} \mid \boldsymbol{\mu}^*) \,\mathrm d\boldsymbol{\xi}}{\int_{\boldsymbol{\xi} \in {\mathcal{P}_\rho}} f(\boldsymbol{\xi} \mid \boldsymbol{\mu}) \, \mathrm d\boldsymbol{\xi}} \int_{\boldsymbol{v} \in {\mathcal{P}_\rho}} f(\boldsymbol{v} \mid \boldsymbol{\mu}) \log f(\boldsymbol{v} \mid M(\boldsymbol{\mu})) \, \mathrm{d}\boldsymbol{v}. \label{Q-func-simp} \end{align}\tag{44}\] For the Gaussian density, we have \[\nabla_{\boldsymbol{\mu}} \log f(\boldsymbol{v} \mid \boldsymbol{\mu}) = (\boldsymbol{\Sigma}^*)^{-1} (\boldsymbol{v} - \boldsymbol{\mu}).\] Using the above identity, we obtain \[\begin{align} \nabla Q(M(\boldsymbol{\mu})\mid \boldsymbol{\mu}) &=\sum_{\rho=1}^{P} \frac{\int_{\boldsymbol{\xi} \in {\mathcal{P}_\rho}} f(\boldsymbol{\xi} \mid \boldsymbol{\mu}^*) \, \mathrm d\boldsymbol{\xi}}{\int_{\boldsymbol{\xi} \in {\mathcal{P}_\rho}} f(\boldsymbol{\xi} \mid \boldsymbol{\mu}) \, \mathrm d\boldsymbol{\xi}} \int_{\boldsymbol{v} \in {\mathcal{P}_\rho}} f(\boldsymbol{v} \mid \boldsymbol{\mu}) (\boldsymbol{\Sigma}^*)^{-1} (\boldsymbol{v} - M(\boldsymbol{\mu})) \, \mathrm{d}\boldsymbol{v}\notag \\ &= \sum_{\rho=1}^{P} \frac{\int_{\boldsymbol{\xi} \in {\mathcal{P}_\rho}} f(\boldsymbol{\xi} \mid \boldsymbol{\mu}^*) \, \mathrm d\boldsymbol{\xi}}{\int_{\boldsymbol{\xi} \in {\mathcal{P}_\rho}} f(\boldsymbol{\xi} \mid \boldsymbol{\mu}) \, \mathrm d\boldsymbol{\xi}} \int_{\boldsymbol{v} \in {\mathcal{P}_\rho}} f(\boldsymbol{v} \mid \boldsymbol{\mu}) (\boldsymbol{\Sigma}^*)^{-1} \boldsymbol{v} \, \mathrm{d}\boldsymbol{v} - (\boldsymbol{\Sigma}^*)^{-1} M(\boldsymbol{\mu}). \label{deltatheta} \end{align}\tag{45}\] Similarly, for \(q(\boldsymbol{\mu}) = Q(\boldsymbol{\mu}\mid \boldsymbol{\mu}^*)\), we have \[\begin{align} \nabla q(\boldsymbol{\mu}) = \mathbb{E}_{\boldsymbol{\mu}^*}\left[(\boldsymbol{\Sigma}^*)^{-1}(\boldsymbol{V} - \boldsymbol{\mu})\right] = (\boldsymbol{\Sigma}^*)^{-1}(\boldsymbol{\mu}^* - \boldsymbol{\mu}). \label{deltathetastar} \end{align}\tag{46}\]

12.3.2 EM Convergence↩︎

Recall that \[q(\boldsymbol{\mu}) = \mathbb{E}_{\boldsymbol{\mu}^*}[\log f(\boldsymbol{V} \mid \boldsymbol{\mu})].\] For the Gaussian density, we have \[\log f(\boldsymbol{V} \mid \boldsymbol{\mu}) = -\frac{1}{2} (\boldsymbol{V} - \boldsymbol{\mu})^\top (\boldsymbol{\Sigma}^*)^{-1} (\boldsymbol{V} - \boldsymbol{\mu}) + C,\] where \(C\) is a constant independent of \(\boldsymbol{\mu}\).

Therefore, \[q(\boldsymbol{\mu}) = -\frac{1}{2} \mathbb{E}_{\boldsymbol{\mu}^*}\left[(\boldsymbol{V} - \boldsymbol{\mu})^\top (\boldsymbol{\Sigma}^*)^{-1} (\boldsymbol{V} - \boldsymbol{\mu})\right] + C.\]

Taking gradient and Hessian with respect to \(\boldsymbol{\mu}\), we obtain \[\nabla q(\boldsymbol{\mu}) = (\boldsymbol{\Sigma}^*)^{-1}(\boldsymbol{\mu}^* - \boldsymbol{\mu}), \qquad \nabla^2 q(\boldsymbol{\mu}) = -(\boldsymbol{\Sigma}^*)^{-1}.\]

Since \((\boldsymbol{\Sigma}^*)^{-1}\) is positive definite, we have \[\nabla^2 q(\boldsymbol{\mu}) \preceq -\frac{1}{\lambda_{\max}(\boldsymbol{\Sigma}^*)} I.\]

The result then follows from the standard characterization of strong concavity.

Using 45 , we can express the gradient difference as \[\begin{align} \nabla Q(M(\boldsymbol{\mu})\mid \boldsymbol{\mu}) - \nabla Q(M(\boldsymbol{\mu})\mid \boldsymbol{\mu}^*) = (\boldsymbol{\Sigma}^*)^{-1} \left( \sum_{\rho=1}^{P} \mathbb{P}_{\boldsymbol{\mu}^*}(\mathcal{P}_\rho) \, g_j(\boldsymbol{\mu}) - \boldsymbol{\mu}^* \right), \end{align}\] where \[g_j(\boldsymbol{\mu}) := \mathbb{E}_{\boldsymbol{\mu}}[\boldsymbol{V} \mid \boldsymbol{V} \in \mathcal{P}_\rho].\] A direct calculation shows that \[\nabla g_j(\boldsymbol{\mu}) = \operatorname{Var}_{\boldsymbol{\mu}}(\boldsymbol{V} \mid \mathcal{P}_\rho)(\boldsymbol{\Sigma}^*)^{-1}.\]

Applying a first-order Taylor expansion around \(\boldsymbol{\mu}^*\), \[g_j(\boldsymbol{\mu}) = g_j(\boldsymbol{\mu}^*) + \operatorname{Var}_{\boldsymbol{\mu}^*}(\boldsymbol{V} \mid \mathcal{P}_\rho)(\boldsymbol{\Sigma}^*)^{-1}(\boldsymbol{\mu}- \boldsymbol{\mu}^*) + r_j(\boldsymbol{\mu}),\] where \(\|r_j(\boldsymbol{\mu})\| = O(\|\boldsymbol{\mu}- \boldsymbol{\mu}^*\|^2)\).

Note that \(\sum_{\rho=1}^J \mathbb{P}_{\boldsymbol{\mu}^*}(\mathcal{P}_\rho) g_j(\boldsymbol{\mu}^*) = \mathbb{E}_{\boldsymbol{\mu}^*}[\boldsymbol{V}] = \boldsymbol{\mu}^*\). Summing the Taylor expansion over \(j\) weighted by \(\mathbb{P}_{\boldsymbol{\mu}^*}(\mathcal{P}_\rho)\), and using the law of total covariance \[\begin{align} \sum_{\rho=1}^{P} \mathbb{P}_{\boldsymbol{\mu}^*}(\mathcal{P}_\rho)\operatorname{Var}_{\boldsymbol{\mu}^*}(\boldsymbol{V} \mid \mathcal{P}_\rho) = \boldsymbol{\Sigma}^* - \operatorname{Var}_{\boldsymbol{\mu}^*}(\mathbb{E}_{\boldsymbol{\mu}^*}[\boldsymbol{V} \mid R']), \end{align}\] where \(R'\) is the categorical region label introduced in Section 4.2, we obtain \[\begin{align} \nabla Q(M(\boldsymbol{\mu})\mid \boldsymbol{\mu}) - \nabla Q(M(\boldsymbol{\mu})\mid \boldsymbol{\mu}^*) = (\boldsymbol{\Sigma}^*)^{-1} \left( \boldsymbol{\Sigma}^* - \operatorname{Var}_{\boldsymbol{\mu}^*}(\mathbb{E}_{\boldsymbol{\mu}^*}[\boldsymbol{V} \mid R']) \right)(\boldsymbol{\Sigma}^*)^{-1}(\boldsymbol{\mu}- \boldsymbol{\mu}^*) + O(\|\boldsymbol{\mu}- \boldsymbol{\mu}^*\|^2). \end{align}\]

Recall the whitening transformation \(\boldsymbol{V}' = (\boldsymbol{\Sigma}^*)^{-1/2}(\boldsymbol{V} - \boldsymbol{\mu}^*)\). The conditional means satisfy \(\mathbb{E}_{\boldsymbol{\mu}^*}[\boldsymbol{V}\mid R'] = (\boldsymbol{\Sigma}^*)^{1/2}\mathbb{E}[\boldsymbol{V}'\mid R']+\boldsymbol{\mu}^*\), so \[\operatorname{Var}_{\boldsymbol{\mu}^*}(\mathbb{E}_{\boldsymbol{\mu}^*}[\boldsymbol{V} \mid R']) = (\boldsymbol{\Sigma}^*)^{1/2} \operatorname{Var}(\mathbb{E}[\boldsymbol{V}' \mid R']) (\boldsymbol{\Sigma}^*)^{1/2}.\] Substituting back yields \[\begin{align} \nabla Q(M(\boldsymbol{\mu})\mid \boldsymbol{\mu}) - \nabla Q(M(\boldsymbol{\mu})\mid \boldsymbol{\mu}^*) = (\boldsymbol{\Sigma}^*)^{-1/2} \left( I - \operatorname{Var}(\mathbb{E}[\boldsymbol{V}' \mid R']) \right) (\boldsymbol{\Sigma}^*)^{-1/2}(\boldsymbol{\mu}- \boldsymbol{\mu}^*) + O(\|\boldsymbol{\mu}- \boldsymbol{\mu}^*\|^2). \end{align}\]

By the law of total variance, \(\operatorname{Var}(\boldsymbol{V}') = \operatorname{Var}(\mathbb{E}[\boldsymbol{V}'\mid R']) + \mathbb{E}[\operatorname{Var}(\boldsymbol{V}'\mid R')]\), and since \(\operatorname{Var}(\boldsymbol{V}')=I\) and \(\mathbb{E}[\operatorname{Var}(\boldsymbol{V}'\mid R')] \succeq 0\), we have \(0 \preceq \operatorname{Var}(\mathbb{E}[\boldsymbol{V}'\mid R']) \preceq I\). Combined with Assumption 1, \[0 \le \lambda_{\max}\!\left(I - \operatorname{Var}(\mathbb{E}[\boldsymbol{V}' \mid R'])\right) \le 1 - \epsilon.\]

By sub-multiplicativity of the spectral norm, \[\big\|(\boldsymbol{\Sigma}^*)^{-1/2}\big(I - \operatorname{Var}(\mathbb{E}[\boldsymbol{V}'\mid R'])\big)(\boldsymbol{\Sigma}^*)^{-1/2}\big\|_2 \le (1-\epsilon)\,\big\|(\boldsymbol{\Sigma}^*)^{-1/2}\big\|_2^{\,2} = \frac{1-\epsilon}{\lambda_{\min}(\boldsymbol{\Sigma}^*)}.\]

Choose \(r\) small enough so that for all \(\boldsymbol{\mu}\in \mathbb{B}_2(r;\boldsymbol{\mu}^*)\), the higher-order term satisfies \(\|O(\|\boldsymbol{\mu}-\boldsymbol{\mu}^*\|^2)\| \le \frac{\epsilon/2}{\lambda_{\min}(\boldsymbol{\Sigma}^*)}\|\boldsymbol{\mu}-\boldsymbol{\mu}^*\|\). Then \[\|\nabla Q(M(\boldsymbol{\mu})\mid \boldsymbol{\mu}) - \nabla Q(M(\boldsymbol{\mu})\mid \boldsymbol{\mu}^*)\|_2 \le \frac{1-\epsilon/2}{\lambda_{\min}(\boldsymbol{\Sigma}^*)} \|\boldsymbol{\mu}- \boldsymbol{\mu}^*\|_2.\]

The result follows from [35]. By Lemma 2, the function \(q(\boldsymbol{\mu})\) is strongly concave with parameter \[\lambda = \frac{1}{\lambda_{\max}(\boldsymbol{\Sigma}^*)}.\] By Lemma 3, the first-order stability condition holds with \[\gamma = \frac{1-\epsilon/2}{\lambda_{\min}(\boldsymbol{\Sigma}^*)}.\] Since \(\gamma < \lambda\) for any \(0 < \epsilon \le 1\), [35] implies that the EM operator \(M(\cdot)\) is a contraction mapping on \(\mathbb{B}_2(r;\boldsymbol{\mu}^*)\). In particular, \[\begin{align} \| M(\boldsymbol{\mu}) - \boldsymbol{\mu}^* \|_2 \le (1 - \epsilon/2) \|\boldsymbol{\mu}- \boldsymbol{\mu}^*\|_2, \quad \forall \boldsymbol{\mu}\in \mathbb{B}_2(r;\boldsymbol{\mu}^*). \end{align}\]

Applying this inequality recursively, we obtain for any initialization \(\boldsymbol{\mu}^{(0)} \in \mathbb{B}_2(r;\boldsymbol{\mu}^*)\), \[\|\boldsymbol{\mu}^{(t)} - \boldsymbol{\mu}^*\|_2 \le (1 - \epsilon/2)^t \|\boldsymbol{\mu}^{(0)} - \boldsymbol{\mu}^*\|_2,\] which establishes the linear convergence of the population EM sequence.

12.4 Concavity↩︎

Firstly, consider the scalar function \[g(z)=\log\sigma(z),\] where \(\sigma(z)=(1+e^{-z})^{-1}\) is the sigmoid function. Since \(\sigma'(z)=\sigma(z)(1-\sigma(z))\), we have \[g'(z)=\frac{\sigma'(z)}{\sigma(z)}=1-\sigma(z),\] and \[g''(z)=-\sigma(z)(1-\sigma(z))\le0.\] Therefore \(g(z)\) is concave in \(z\).

Then observe that the utility difference \(\Delta U_{nj}(\boldsymbol{v};\boldsymbol{A})\) is affine in \(\boldsymbol{A}\). In particular, \[\frac{\partial \Delta U_{nj}(\boldsymbol{v};\boldsymbol{A})}{\partial\boldsymbol{A}} = (1-0^{c_n})\,\boldsymbol{x}_{nc_n}\boldsymbol{x}_{nc_n}^\top - \left[0^{c_n}+(1-0^{c_n})(1-0^{j-c_n})\right]\boldsymbol{x}_{nj}\boldsymbol{x}_{nj}^\top ,\] which does not depend on \(\boldsymbol{A}\). Hence \(\Delta U_{nj}(\boldsymbol{v};\boldsymbol{A})\) is an affine function of \(\boldsymbol{A}\).

Since \(g(z)\) is concave and \(\Delta U_{nj}(\boldsymbol{v};\boldsymbol{A})\) is affine in \(\boldsymbol{A}\),the composition rule for concave functions implies that \[\log\sigma\left(\frac{\Delta U_{nj}(\boldsymbol{v};\boldsymbol{A})}{\lambda}\right)\] is concave in \(\boldsymbol{A}\). Since sums and expectations preserve concavity, the objective \(\hat{Q'}(\boldsymbol{\theta}|\boldsymbol{\theta}^{(t)})\) is concave in \(\boldsymbol{A}\).

12.5 Gradient↩︎

Consider the Monte Carlo approximation of the smoothed objective \[\hat{Q}'(\boldsymbol{\theta}|\boldsymbol{\theta}^{(t)}) =\sum_{n=1}^{N}\sum_{l=1}^{L}\overline{w}_{nl} \left( \log f(\boldsymbol{v}_n^{(l)}|\boldsymbol{\mu},\boldsymbol{\Sigma})+\sum_{j=1}^{J_n}\log \sigma\left( \frac{\Delta U_{nj}(\boldsymbol{v}_n^{(l)};\boldsymbol{A})}{\lambda} \right) \right).\] When differentiating with respect to \(\boldsymbol{A}\), the Gaussian log-likelihood term \(\log f(\boldsymbol{v}_n^{(l)}|\boldsymbol{\mu},\boldsymbol{\Sigma})\) does not depend on \(\boldsymbol{A}\) and therefore vanishes. Hence only the sigmoid contribution needs to be differentiated.

Using the identity \[\frac{d}{dz}\log\sigma(z)=1-\sigma(z),\] and applying the chain rule yields \[\frac{\partial}{\partial \boldsymbol{A}}\log\sigma\left( \frac{\Delta U_{nj}(\boldsymbol{v}_n^{(l)};\boldsymbol{A})}{\lambda} \right) =\frac{1}{\lambda}\left( 1-\sigma\left( \frac{\Delta U_{nj}(\boldsymbol{v}_n^{(l)};\boldsymbol{A})}{\lambda} \right) \right) \frac{\partial \Delta U_{nj}(\boldsymbol{v}_n^{(l)};\boldsymbol{A})}{\partial \boldsymbol{A}}.\]

Substituting this derivative into the objective and exchanging summation with differentiation gives \[\nabla_{\boldsymbol{A}}\hat{Q'} = \frac{1}{\lambda} \sum_{n=1}^{N}\sum_{j=1}^{J_n} \left( \frac{\partial \Delta U_{nj}(\boldsymbol{A})}{\partial \boldsymbol{A}} \sum_{l=1}^{L}\overline{\omega_{nl}} \left( 1-\sigma\left( \frac{\Delta U_{nj}(\boldsymbol{v}_n^{(l)};\boldsymbol{A})}{\lambda} \right) \right) \right).\]

Finally, since \(\Delta U_{nj}(\boldsymbol{v};\boldsymbol{A})\) is affine in \(\boldsymbol{A}\), its derivative admits the closed form stated in the theorem. This proves the theorem.

References↩︎