Discussion Paper No. 18
A Large-Scale Marketing Model using Variational
Bayes Inference for Sparse Transaction Data

Tsukasa Ishigaki
Nobuhiko Terui
Tadahiko Sato
Greg M. Allenby
January , 2014

Data Science and Service Research
Discussion Paper

Center for Data Science and Service Research
Graduate School of Economic and Management
Tohoku University27-1 Kawauchi, Aobaku
Sendai 980-8576, JAPANA Large-Scale Marketing Model using Variational Bayes Inference
for Sparse Transaction Data

Tsukasa Ishigaki*
Nobuhiko Terui*Tadahiko Sato**
Greg M. Allenby***January, 2014

*Graduate School of Economics and Management, Tohoku University, Sendai 980-8576, Japan** Graduate School of Business Science, University of Tsukuba, Tokyo 112-0012, Japan*** Fisher College of Business, Ohio State University, Columbus, OH, U.S.A.
Ishigaki acknowledges the financial support of KAKENHI Grant-in-Aid for Young Scientists(A) 24683012. Terui acknowledges the financial support of the Japanese Ministry of Education
Scientific Research Grants (A) 25245054.

1

Abstract
Large-scale databases in marketing track multiple consumers across multiple productcategories. A challenge in modeling these data is the resulting size of the data matrix, which
often has thousands of consumers and thousands of choice alternatives with prices andmerchandising variables changing over time. We develop a heterogeneous topic model forthese data, and employ variational Bayes techniques for estimation that are shown to beaccurate in a Monte Carlo simulation study. We find the model to be highly scalable and usefulfor identifying effective marketing variables for different consumers, and for predicting the
choices of infrequent purchasers.

Key words and phrases:
database marketing, dimension reduction, Bayesian analysis, choice models, topic model,
targeting2

1. Introduction
Modern analytic techniques in marketing are continuously confronted with the necessity of
extracting relevant information from large volumes of data by identifying important drivers ofconsumer behavior. It is common for datasets to record household purchases of products thatare orders of magnitude larger than what current models of behavior are currently capable.
Existing models of choice and demand, for example, are typically limited to less than twenty
or so product alternatives that are tracked across possibly hundreds of consumers (see Rossi etal. 2005 and Chintagunta and Nair 2011).
Increasing the number of products analyzed is problematic because of potential
complexities in the structure of demand and the accompanying increase in the required numberof model parameters. Increasing the number of respondents is also problematic because ofcomputational constraints arising from respondent heterogeneity that is found to be importantin describing demand and deriving policy implications. While a variety of dimension-reducing
techniques have been studied in the fields of statistics and data-mining, the presence ofheterogeneous consumers and heterogeneous purchase environments with prices and othervariables change over occasions requires the use of model-based inference as opposed tomethods applied directly to the marginal data (Chintagunta and Nair 2011).
Naik et al. (2008) discusses three solutions to the challenges in massive data analysis:
increasing computer power, employing alternative approaches for data analysis, and using
scalable estimation methods. In this paper, we combine the second and third options to obtain
improved inferences about consumer behavior in large datasets. We extend the voting bloc
model of Spirling and Quinn (2010) and Grimmer (2011) that are a variation of topic models
used to conduct large-scale analysis of text data (Blei et al. 2003).
The topic model is a generalization of a finite mixture model (Kamakura and Russell 1989)
in which each data point is associated with a draw from a mixing distribution (Teh and Jordon
2010). Models of voting blocs (Spirling and Quinn 2010), for example, track the votes oflegislators (aye or nay) across multiple bills, with each bill associated with a potentiallydifferent concern or issue. Similarly, the latent Dirichlet allocation (LDA) model of Blei et al.
(2003) allocates words within documents to a small number of latent topics whose patterns aremeaningful and interpretable. Each vote and each word is associated with a potentially different
issue or topic, and hence the mixing distribution is applied to the individual datum. In our3analysis of household purchases, we allow every purchase (and every non-purchase) in every
product category to be related to a potentially different latent context (topic, or issue) for which
the good is purchased. This allows us to view a consumer's purchases as responding to different
needs or occasions (e.g., family dinner, snacks, etc.), and allows us to identify the ensemble ofgoods that collectively define latent product segments across a large number of products.
We obtain a scalable estimation method by employing variational Bayes (VB) inference asin Jordan et al. (1999) and Bishop (2006), instead of the standard Markov chain Monte Carlo(MCMC) inference. MCMC methods can incur large computational cost in large-scale
problems. VB inference approximates a posterior distribution of target by variational
optimization in a computationally efficient manner.
Our approach combines variational Bayes (VB) methods, as in Jordan et al. (1999) andBishop (2006), with a topic-like probit model to obtain a computationally feasible model ofconsumer purchases that is scalable to large databases. Individual-level inference is possible inour model, where we can identify the marketing variables that are effective for specificindividuals and the products for which they are effective. Our model is therefore similar toadaptive personalization systems proposed by Ansari and Mela (2003), Rust and Chung (2006),
Chung et al. (2009) and Braun and McAuliffe (2010). However, it is different in that our model
structure facilitates analysis of a much larger array of product categories.
In the next section, we propose a model for consumer purchases in multiple productcategories. Section 3 describes a variational Bayes inference scheme for the models and asimulation study that verifies the scalability. The prediction performance of the proposedmodels is presented in Section 4. Section 5 applies the model to actual customer purchases ina general merchandise store. Discussion and concluding remarks are presented in Section 6.

2. ModelEstimating parameters of choice model for a large number of consumers and products isoften computationally infeasible. In addition, the actual sample size of transactional data isoften much smaller than the data space reflected by a data cube with dimensions corresponding
to the number of consumers, number of products and time, making fixed-effect estimates ofmodel parameters with heterogeneity unsuitable. We address this challenge by relating
consumer purchases to latent segments (similar to topics and blocs) that greatly reduces the4dimensionality of the model. Response parameters are then introduced in the reduced
dimensional space by connecting each choice to marketing variables with a hierarchical probit
model.
2.1 Dimensional Reduction by Topic Models
Dimensional reduction is an important technique in massive data analysis. Here webriefly introduce the idea of introducing a latent variable to common in topic models in thecontext of consumer purchases. We seek the probability p  i | c  that consumer c purchases
item i . Dataset includes C consumers and I product items through T periods. However, theprobabilities cannot be directly calculated because of computational difficulty imposed by thelarge-scale setting and data sparseness. The topic model calculates p  i | c  by introducing alatent class z  1, , Z  whose dimension is significantly smaller than the number ofconsumers and items ( Z  C, I ).
The latent variable is used to represent the sparse data matrix as a finite mixture ofvectors commonly found in topic models: p(i  1 | c  1)  p(i  1 | c  C )  Z  p(1 | z )     p( z | 1)  p( z | C ).




 


z 1
 p(i  I | c  1)  p(i  I | c  C )  p( I | z ) 

(1)

More specifically, we decompose a large probability matrix of size C  I to two small
probability matrices of sizes I  Z and Z  C based on the property of conditional
independence.
The main difference between voting blocs model and LDA is assumed distributions
for probabilities p  i | z  in the I  Z size probability matrix. The voting blocs model
supposes a Bernoulli distribution for the probability p  i | z  . LDA assumes a categorical
distribution for the probability matrix.
In the analysis of purchase behavior using topic models for large consumer transaction
data, Iwata et al. (2009) extracted dynamic patterns between purchased product items andconsumer interests. Ishigaki et al. (2010) fused heterogeneous transaction data and consumer
lifestyle questionnaire data, while Iwata et al. (2012) identified consumer purchase patterns byusing a topic model with price information on the purchased products. These approaches
identify patterns among consumers and product items. The labeled LDA proposed by Ramage
et al. (2009), and the supervised LDA of Blei and McAuliffe (2007) extend the topic models
5

by incorporating additional data in the analysis. However, none of these approaches are suitable
for relating marketing variables to individual consumer choices as explanation variables. In thefollowing sections, we construct a model that links marketing variables with consumers andproducts.
2.2 A Reduced Dimensional Choice Model
Let ycit denote consumer c’s purchase record of product i at time t, assigning y cit  1
if consumer c purchased the item, and y cit  0 otherwise. Denote ucit as the utility of
consumer c’s purchase record of product i at time t. We assume a binary probit model with
ucit  0 if ycit  1 , and ucit  0 if ycit  0 . We couple the topic model in (1) with the binary
choice probability as in a voting bloc model to obtain the choice probability:p ( u cit  0)   z 1 p ( u cit  0 | z ) p ( z | c ) .
Z

(2)

(z)
We denote the utility associated with the latent class z as ucit
, and then the choice(z)probability can be represented as p  u cit  0 | z   p  u cit 0  . Assuming a linear

Gaussian structure on the segment utility

(z)
ucit

with marketing variables as

(z)
u cit
 N  x Tit β zi ,1  , the right hand side of (1) is represented as, F ( x1Tt β z1 ) 



 
  p ( z |1)  p ( z | C ) z 1T F (x β ) 
It zI Z

(3)

where β zi    zi1 , ,  ziM  is a response coefficient vector of latent class z with respect to item
Ti, x it   xit1 , , xitM T

is a vector of M marketing variable for item i at time t, and F   is the

cumulative distribution function (CDF) of the standard normal distribution. In our empirical
study, x it includes price and promotional variables.
We next set a categorical distribution C cz for the probability p( z | c) that consumer

c belongs to the latent class z. The categorical distribution is multinomial with parameters Cc .
The Cc is determined so that the selection probability of consumer c with respect to item i isconditionally independent if the latent class z is given. Then, the right hand side of (1) isrepresented by: F ( x1Tt β z1 ) 



 
 C1 z  CCz z 1
 F ( xT β ) 
It zI 

Z

6

(4)

Finally, segment-level heterogeneity is introduced through a hierarchical model witha random effect for response coefficient β zi

β zi  N M  μi ,Vi  ,(5)where the prior distributions for μ i and Vi follow an M-dimensional multivariable normal
distribution N M  μ ,  Vi  and an inverse-Wishart distribution IW W , w  , respectively.

 are hyperparameters. That is, we assume that the M-dimensionalμ ,   , W and w
coefficient vector  zi for each segment, z, is a draw from a distribution with mean andcovariance that is item-specific.
We specify a prior distribution for Cc , assuming the Dirichlet distribution as thenatural conjugate prior distribution of categorical distribution:

Cc ~ Dirichlet  γ  ,(6)where γ is a hyperparameter vector of the Dirichlet distribution.
The likelihood is given asCZl { ycit } |{Cc },{β zi },{x it }    Ccz p  ycit | x it , β zi , z  ,c 1 iI c tTc z 1(7)where p  ycit | x it , β zi , z  denotes the kernel of the binary probit model conditional on z, T cdenotes a subset of t in which consumer c purchased any item in a store, and I c is a subset
of items i purchased by consumer c during the period

 and I i |  y  0 .

Tc  t |  i 1 ycit  0
I

t  1,, T , that is,

T

c

t 1

cit

Equation (7) is difficult to use directly because the likelihood includes summations
over latent class z. Instead, we employ a data augmentation approach by Tanner (1987) withrespect to latent variable z. We introduce variables zcit  1, , Z  denoting the label of thelatent class for each consumer c, each purchased item i, and each purchasing event t.
Conditioning on the zcit for each purchasing transaction, as in the LDA of Blei et al. (2003),
the likelihood in (7) simplifies to:C

l { ycit } |{Cc },{zcit },{β zi },{x it }    p  zcit | Cc  p  ycit | x it , β zi , zcit  ,c 1 iI c tTc

(8)where p  zcit | Cc  denotes a categorical distribution when Cc is given.
Our model for massive data analysis is different to the LDA model in that it only deals withthe presence of products appearing in the purchase basket of the consumer, and does not deal
7

with non-purchase of product. This is different to model encountered in the analysis of text
data where it is the presence of words, and not their absence, that characterizes the latent topics.
The co-occurrence of the products selected during a shopping trip is what gives meaning tosegments, as modified by the marketing variables.
(z)

The posterior distribution of parameters including latent variables {zcit },{ucit } is thengiven by(z)
p {Cc },{zcit },{ucit},{β zi },{μi },{Vi } | {x it },{ ycit } p C c  | {zcit } p zcit  | Cc , β zi , x it , ycit 





(z)
 p ucit | β zi , zcit , xit , ycit  p μi , Vi  | {β zi }





(z)
, μi , Vi , x it 
 p β zi  | ucit

(9)

(z)
},{β zi },{μi },{Vi },{x it },{ ycit } p {C c },{zcit },{ucit
Z
 C
 I

=   p  Cc     p  μi , Vi  p  β zi | μi , Vi  z 1
 c 1  i 1


 C

(z)
    p  zcit | Cc  p  ycit | ucit , β zi , zcit , x it   c 1 iIc tTc

3. Variational Bayes Inference
We introduce VB inference in order to achieve computational feasibility for large-scaletransaction data. VB inference approximates the posterior, or target distribution in a Bayesianmodel. The advantage of this method over MCMC is low computational cost. VB also takes
advantage of parameters that can be decomposed into several mutually independent groups.
This is necessary for our analysis using a large database.
The target and approximate distributions are denoted as p and q, respectively. The latter iscalled the variational distribution. Distributions p and q share a parameter set θ . In general,
when the data D is given, the log marginal likelihood p ( D ) of the target distribution isdecomposed into two components as

8

log p ( D )  L ( q )  KL ( q || p ),L ( q )   q ( θ ) log  p ( D , θ ) q ( θ ) 1dZ ,KL( q || p )    q ( θ ) log  p ( θ | D ) q ( θ )

1

(10)

dZ .

L( q) is called the variational lower bound in VB inference, andKL( q || p ) is the

Kullback–Leibler divergence of the target and variational distributions. As is well known,

KL( q || p ) is zero if p and q are the same distribution. Therefore, a reasonable solution toestimating the posterior distribution p is the variational distribution q for which KL( q || p ) is
minimized. However, it is difficult to evaluate the value of KL( q || p ) because the expression
involves a posterior distribution of p(θ | D) .
In contrast, L( q) involves a joint distribution p ( D , θ ) that is easily evaluated in manycases because it is obtained as the product of the prior and the likelihood in Bayesian models.
We note that maximizing L( q) is equivalent to minimizing KL( q || p ) because the log
marginal likelihood of the target distribution is constant for a given dataset. In this situation,
assuming that the distribution q and parameter set θ are decomposable for some groups asq  θ    j q j  θ( j )*  , where the parameters θ( j )* are called variational parameters, L(q)J

can be maximized by the following updating algorithm (Jordan et al., 1999):θ( j )*{new}  arg max L  j q j θ( j )*
θ( j )*



J





 exp Ei  j log p  D, θ   .

(11)

The Ei  j   are the expectation value associated with q j distributions over all parametersθ( j )* , where i  j . The variational parameters are updated for each variational parameter set

θ( j )* until convergence of the algorithm. The initial variational parameters are proper random
values. The VB is guaranteed to converge after several iterations because L( q) is convex withrespect to each q j  θ( j )*  (Bishop, 2006). The variational lower bound monotonicallyincreases as the iteration proceeds; therefore, convergence can be confirmed by checking thevalue of L( q) at each iteration.
3.1 VB for the Proposed Model
We introduce the variational distributions and parameters for the modes of proposedmodel. The parameters and variational parameters are denoted as(z)θ  {Cc },{zcit },{ucit},{β zi },{μi },{Vi } and

9





θ*  C*cit  , γ*c  , β(citz )*  , μ*iz  , Viz*  , μi*  ,  i*  , wi*  , Wi*  respectively, while thevariational distributions are configured asq  θ | θ* ,{xit },{ ycit }
 C

 C
 C
(z)
  qc  Cc | γ*c    qz  zcit | C*cit    qu  ucit
| β(citz )* , xit , zcit  c1
  c1 iIc tTc  c1 iIc tTc


(12)

I
 I Z
*
*
*
* 
*
* 
q
|
,
V
β
μ
   zi zi zi    q ,V  μi ,Vi | μi ,  i , wi , Wi   i 1 z 1
  i 1
*
where qc is a Dirichlet distribution with variational parameter γ c , qc is a categorical*
distribution with variational parameter Ccit , qu is a truncated normal distribution withparameter zcit and variational parameter βcit , q is an M-dimensional multivariable( z )*

normal distribution with two variational parameters (mean vector μzi and covariance matrix*

Vzi* ), and q,V is a multivariable normal–inverse Wishart distribution with variationalparameters μi ,  i , wi ,Wi . Here, to realize effective variational inference, we assume that
*

*

*

*

all variational parameters are independent. The update equation and the derivations of thevariational parameters are detailed in Appendix A.

4. Simulation StudyIn this section we examine the performance of the proposed VB estimator relative to
MCMC using simulated data. We show that MCMC becomes too computationally demanding
as the size of the dataset increases, and that VB provides a computationally efficient andaccurate approximation to the posterior with good predictive properties.
The simulation datasets are generated as follows:a) We determine the number of consumers, items and time period for the dataset.
b) Consumers randomly are assigned to a segment. Consumers assigned to same
segment have a same product set of items, which are a subset of items that consumers
evaluate when making a purchase decision. The number of consumer segments andthe number of items in a product set is set to 0.02 C and 0.2  I , respectively, where
C and I are varied in our analysis.
c) Items in a product set are randomly assigned from all available items, with a lower
bounds on the number of items needed in the set.
10d) The marketing variables are comprised of x it  1, Pi  , where P i is a price discountTrate generated from a uniform distribution with interval [0.1 1.0].
e) Each consumer randomly visits the simulated shop on five occasions during data
period T, and purchases N items with probability proportional to the discount fromproduction set.
Computational time and predictive results are calculated using data simulated from theabove steps. The computational results reported below were calculated in same computational
environment (64-bit version of Python 2.7.5, implemented on a 3.5 GHz processor (Quad-Core
Xeon; Intel Corp.) with 64 GB memory).
4.1 Scalability
The scalability is investigated under the condition as follows: C = {1000, 5000, 10000}, I= {100, 500, 1000}, T = 30, Z = {5, 10, 20} and N = 10. Thus, 27 different scenarios were
explored in the simulation study. The simulation times in hours are shown in Table 1. Here, weTT
set hyperparameters as γ   0.1, ,0.1 , μ   0,  , 0  ,   1, W   I M , w   10 in

VB and MCMC. I M means identity matrix of size M. In VB method, iterations are

terminated when the variational lower bound improves by less than 105 % of current value in
two consecutive iterations. The MCMC method uses Gibbs sampler, and its simulation timesof 6,000 MCMC samples are estimated from those of 10 samples, as is to be computationallyinfeasible. We also note that the selection of 6,000 MCMC samples is consistent with thesimulation study of Braun and Mcauliffe (2010). Three kinds of settings for hyperparameters,
stopping rule of VB iteration and the number of MCMC samples, defined above, are adopted
in all empirical studies hereafter.
Table 1 shows the computational time for VB and MCMC methods. In both algorithms, thecost increases approximately linearly with the size of the dataset specified in term of thenumber of consumers, items, and latent classes. In all scenarios, the computational time of
MCMC exceeds that of VB. The VB algorithm is around 20 to 50 times more efficient than
MCMC, depending on the scenario. The time of estimation using large-scale data (C = 10000,
I = 1000) by MCMC is estimated over 450 h, and thus we recognize that MCMC is notapplicable for our problem.
11

Table 1 : Simulation time by VB and MCMC

4.2 Data Sparseness and RMSE of Prediction
The predictive performance of the models is investigated by simulation, focusing onwhether the model can adapt to sparse data or not. Sparse data are commonly encountered indatasets containing many items. In this paper, we define the data density rate (DDR) as

DDR   C  I  T 

1

C

I

T

 y
c1 i 1 t 1

(13)

cit

The DDR specifies the rate of ycit  1 events in the data space. The N controls values of DDR
in the simulation dataset. Here we generated datasets for C = 500, I = 100, T = 100, Z = 20, andN = {2, 4, 6, 8, 10}; this specification of N implies that DDR = {0.1%, 0.2%, 0.3%, 0.4%,
0.5%}. DDRs of actual scan-panel datasets are, to our knowledge, always below 1 %.
The prediction performance is measured by the root mean square error (RMSE), given byCRSME   I   Tc 
 c 1 

1 CI  y  p ( y  1)
c 1 i 1 tTc

2

cit

(14)

cit

where Tc denotes the number of elements Tc for each consumer; that is, the number of store
visits within the specified data period. The p( ycit  1) is calculated by  z 1 C cz F ( x Tit β zi ) .
Z

The hold-out and hold-in samples are generated by the same procedure. Using the Gibbs
sampler, we generate 6,000 samples of each parameter, where the first 5,000 are discarded asburn-in samples.
Table 2 shows the average RMSE obtained in three simulation trials of DDR for each offour methods, namely Random, Homogeneous, VB and MCMC. In this simulation, we set fivelevels of DDR from 0.1% to 0.5%. “Random” means the situation that consumers are permitted
completely random choice of purchased items and shop visits, and its RMSE is approximately
0.577. “Homogeneous” implies the model with Z = 1. We first observe that the proposed models
works well compared to “Random”. Second, heterogeneity significantly improves thepredictive performance. Third, the performances of VB and MCMC are comparable. These
properties hold throughout every setting of DDRs.

Table 2 RMSE as a function of DDR in simulations12

5 Empirical Application Using Customer DatabaseIn this section, we apply the proposed model to a real large-scale customer dataset.
The results of the simulation study indicate that the VB estimator provides a close
approximation to the MCMC with a large improvement in computational speed. We use VB in
the empirical application in this section and report on estimation results and their managerial
implications.
5.1 Data and Variables
A customer database from a general merchandise store, recorded from April 1 to June
30 in 2002, is used in our analysis. A customer identifier, price, display, and feature wererecorded for a given purchase occasion. The dataset contains 162,775 transactions involving
1,647 consumers and 1,004 items. The 1,004 selected items were displayed or featured at least
once in the data period. The DDR of the scan-panel data is 0.31 %.
The marketing variables are price ( Pit ), display ( Dit ), and feature ( F it ); that is,
x it  1 Pit

Dit

Fit  . Pit is the price relative to the maximum price of item i in theTobservational period. Display and feature are binary entries, equal to one if the item i isdisplayed or featured at time t, and zero otherwise.

5.2 Model Comparisons
The proposed models are compared in terms of RMSE. The parameters are estimatedwith the number of segment Z = {2, 3, 4, 5, 10, 20}. The hold-out sample comprises records
from July 1 to September 30, 2002. The RMSEs for comparable models are shown in Table 3.
We observe that our proposed models (Z greater than two) have smaller value of RMSEs than
“Random” and “Homogeneous” models. The models with Z greater than five have the same
RMSEs and thus we understand that the model with Z = 5 is appropriate for the empiricalanalysis below.
The comparison of RMSE of (i) all customers with that of (ii) infrequent customers
provides useful information of the performance of our models. The largest number of purchasesby one customer in our data is 390 items and 88 visits to the store, and we define infrequentcustomers as those with fewer than five purchases and three visits to the store. The predictive
performance of infrequent consumers slightly decreases compared to that of all consumers,
13however, the RMSE is almost equivalent to that all customers. Thus, our models present atolerable prediction performance for even infrequent consumers.

Table 3 : RMSEs for real customer database – all and infrequent customers

5.3 Segment-Level Parameter Estimates
Table 4 displays the parameter estimates for “Price”, “Display” and “Feature” and“Intercept.” The rows of the tables correspond to products, and the columns to the segments.
The rows are ordered in terms of the differences of estimates among the segments. The firstrow of each table is for the product estimates with maximum variation, and the last row of eachtable corresponds to estimates with minimum variation. The products positioned near the topof the table have larger heterogeneity in the response parameters, while those at bottom of thetable have relatively similar values among segments. The products are identified in terms of
their sales rank, with the product named “No. 1” the most purchased product in our database.
We observe that “Price” coefficients are estimated negative and “Display” and“Feature” are positive for most products. This means that our proposed model produces
economically reasonable estimates. The coefficient estimates also indicate the effectiveness ofthe marketing mix variables for each product. In the “Price” table, for example, product No.
205 has the highest rank and consumers in segment 3 do not respond to variation in price forthis product. For product No. 111, segment 1 is the least price sensitive and for product No.
153 the price insensitive segment is segment 5. Similar results are found with the display andfeature portions of the table, where we see that the most responsive segment is product-specific.
We also find that the lower ranked products in each table show nearly uniform response in eachsegment. The results imply that marketers can perform effective promotions to specificsegment for higher ranked products, however, homogeneous promotions to any segments areenough for lower ranked products.

Table 4: Characteristics of β zi for the five segments of consumers

Next, we extract relative preferences of product category for five segments. First, we
define the relative preference score of segment z for product i by using estimated intercepts for14

product i as RPzi   zi 0  Z 1  z 1  zi 0 . The consumers in segment z with the highest value ofZ

RPzi relatively prefer product i than consumers in other four segments. Thus, we observe
preference of five segments by ordering of the score with respect to each segment.
For visualization of segment’s preferences, we count the number of product category
in top 100 products in order of RPzi with respect to each segment. Table 5 shows name andnumber of product categories containing over three kinds of products in the top 100 products
for each segment. It discloses consumer’s preferences of purchased products for estimatedsegments.

Table 5: Relative preferences of purchased product category for five segments

5.4 Individual Level Parameter Estimates
Individual level estimates of market response is obtained by taking expectation ofsegment level estimates with respect to Ccz
Z

αci   Ccz βzi .

(15)

z 1

We characterize individual consumer in terms of her estimated response coefficient α ci . The
empirical marginal distribution of individual consumer parameter estimates taking averageI  α 
1

I

i 1

ci

i 1, , I

of 1,647 products for each marketing variable are displayed byhistograms in Figure 1(a). On the other hand, the empirical marginal distributions of individual
product, taking average over 1004 consumers, i.e., of

C 
1

C
c 1

α ci



c  1, , C, are depicted in

Figure 1(b). The products that never displayed and featured in the data period have been
omitted.
The marginal distributions provide reasonable individual consumer estimates sincealmost intercept’s and price’s coefficients are negative and almost display’s and feature’s
coefficient are positive. The distributions of individual product estimates show that thefeature’s distribution has a sharp peak around zero, implying that promotions by display andfeature are effective for many products, however, there are a lot of products that they are not15effective for the feature.
Table 6 shows the results of testing the significance of estimated response coefficientsfor individual consumers and products by using the 95% HPD (highest probability density)
region. Owing to the space constraints, the results are shown for a portion of the dataset (fiveconsumers and five products). “*” signifies that the 95% HPD region does not contain zero; if
the HPD region includes zero (i.e., if the coefficient is insignificant), this square is left blank.
We call this graph “Customer-Promotion Diagram”. For example, in the “Price” diagram, thepurchase of No.253 by consumers (a) and (b) is highly influenced by price; however, the priceof No. 318 affects the purchasing behavior of consumers (a) only. We also report that
discounting No.18 will not promote consumer (b) to purchase it. This analysis informs retailers
and marketers which promotion of specified product is effective to any specific individual
consumer. Thus, our proposed model enables marketers to develop effective pricing andpromotional strategies for targeted consumers and products.

Table 6: Personalized effective marketing variables for individual consumers and products

Figure 1 Marginal distribution of parameter estimates of individual consumers and products

5.5 Precision of Approximation to Posterior Density

We examine the precision of VB by comparing estimates to those obtained with
MCMC, as is assumed that MCMC is a more correct than VB as long as there are sufficient
iterations to fully characterize the posterior distribution. For comparison purposes, we reduced
the size of data so that MCMC computations terminate within one day. We extracted 500
customers randomly from our dataset, and choose the top 100 products in terms of sales volume.
Table 7 shows estimates of response parameters for VB methods. The vector of
estimates for 100 products are displayed in row according to the order of number of purchases,that is, first row specifies the estimates for the most purchased product. The numbers in thetable are sample mean of segment level estimates. “-” mean the product which has neverdisplayed or featured in the record. From this result, we observe that price coefficients arereasonably estimated in the sense most of product have negative values. The same thing holds
for the coefficients of display and features. That is, our proposed models perform well.
16

Table 7 : Estimated parameters by VB and MCMC

6. Discussion

This paper addresses two challenges in estimating models of demand in large
databases: i) the large number of available products and ii) the large number of consumers whopurchase these products. Existing models in marketing and methods of estimation tend to focus
on a narrow set of products and a subset of consumers to understand the richness of thecompetitive environment within a product category among a random sample of consumers.
This goal, however, is often at odds with the goals of practitioners who want to score existingdatasets to identify a wide set of customers and products to allocate promotional budgets andincrease sales.
We propose a descriptive model of demand based on the idea of topic models where
products purchased by consumers take the place of words used by authors in creating
documents. We allow for a product's purchase probability to be affected by price, display andfeature advertising variables, but do not treat purchases to arise from a process of constrained
utility maximization. The advantage of this approach is that it allows us to side-stepcomplications associated with competitive effects and model a much larger set of products thanthat possible with existing economic models. By retaining prices and other marketing variables
in our model we can still predict the effect of these variables on own-sales. This tradeoff isinevitable in the analysis of large-scale databases where purchases are tracked across thousandsof products. The proposed model links the characteristics of consumer segments to marketing
variables, and it is applicable to both segment-level and individual-level marketing across alarge set of products.
The scalability and predictive performance of the proposed models were confirmed
through a simulation study involving variational Bayes inference. In our analysis, we imposeda fairly conservative convergence criteria for VB of 105 %, but also found that coarser
thresholds (for instance, 103 %) produced similar results. We therefore believe that estimation
times can be further reduced in practice from those reported in this paper.
Finally, we employed the RMSE criterion for choosing the number of segments. In the
VB framework, the variational lower bound is used for this criterion, as is shown in, for17example, Corduneanu and Bishop (2001). The variational lower bound in our models issomewhat sensitive to changes in the number of segments, Z. We identify this as a future
research.

18

Appendix A: Derivation of VB Algorithm for Proposed Model

This appendix details the variational inference of proposed model. The update procedure
derives from the analytical calculation of equation (13). The update equation for eachvariational parameter is obtained from the following expectation valuesE q j  log p  D, θ    Ei  j  log p  D, θ    log p  D, θ   qi  θ ( i )* dθ ( i )* ,

(A1)

i j

where D  x it  ,  y cit  .
The update procedures of variational parameters Ccit , γ c , βcit , μiz , Viz , μ i ,  i , wi ,**( z )*

*

*

*

*

*

*

and Wi are presented below.
A.1 Optimization of γ c*

The Dirichlet and categorical distributions are of the following forms:    C 
Dirichlet  C | γ      
Z

z 1

c

z

Z

z 1

Z

z

Z

z 1

z 1
cz

(A2)

Categorical  zcit | Cc    Ccz cit z zz 1where    is the gamma function and   zcit  z  is the Dirac delta function defined as  zcit  z   1 if zcit  z , and   zcit  z   0 otherwise. The expectation value
E  qc  log p  D , θ   is then calculated for each c asEqc log p  D, θ    log p  Cc   Eqz log p zcit  | Cc    const.
Z
Z 

Z
*
 log   z 1 z   log z     z   Ccitz
 1  log Ccz  const.
z 1
z 1 
iI c tTc

 





(A3)

Here and hereafter, const. denotes any terms not included in the relevant parameters. The
second line of the above equations describes a log-Dirichlet function with parameter*z   Ccitz. Therefore,iI c tTc*
γ *(new) γ   Ccitzc
iI c tTc

(A4)

*

A2. Optimization of Ccitz
Here we denote a digamma function as     , which will be useful for later discussion,and summarize the property of truncated normal distribution in the probit model. u c( itz )19

follows a normal distribution with mean x it β zi and variance 1. Moreover, u c( itz ) must satisfyTy cit  1 if u cit  0

and y cit  0 if u cit  0 . Therefore, u c( itz ) is generated from atruncated normal distribution as

(z)
cit

u

 TN(0,)  xTit βiz ,1 if ycit  1~.
TTN( ,0)  xit βiz ,1 if ycit  0(A5)

where TN( n1 ,n2 ) ,  denotes a normal distribution truncated from n 1 to n 2 . The
distribution of u c( itz ) is therefore expressed as(z)p  ucit| β zi , zcit , xit , ycit  with cit  1  F xit βzi
( z)

T



  ycit 1

2
1
 1 (z) T
exp   ucit
 xit β zi  
2
 2


1
(citz )F  x β 
T
it zi

(A6)

  ycit 0. In addition, the expectation value andvariance are expressed as( z)  xTit β(citz )*  cit( z )E ucit(z)  1  xTit β(citz )*cit( z )  cit( z ) Var ucit

where (z)cit

 f   x Tit β(citz )*   1  F   x Tit β(citz )*    ycit 1 f   x Tit β(citz )*   F   x Tit β(citz )*  



(A7)

2

  ycit  0 .

Thus, the expected value E q  log p  D , θ   is given asz

Eqz log p  D, θ    Eqc log p  zcit | Cc  

(A8)

(z)
 Equ ,q log p  ucit
| β zi , zcit , x it , ycit    const.

The first term in the right hand side of Equation (A8) is obtained as   cz*     z 1  cz*
Z



(Blei et al. 2003), while the second term is evaluated as2
1 (z) T(z)
Equ ,q log p  ucit
| β zi , zcit , xit , ycit    Equ ,q   log 2 (citz )   ucit
 xit β zi  
2


2
2
1
1
(z)
E
 (z) T 
 Eq log (citz )   Equ  ucit
Eq  xTit β zi    const.
qu , q ucit xit β zi  





2
2

(A9)

*

To solve Equation (A8) for Ccitz , we must evaluate the four terms of Equation (A9). The firstterm includes a CDF from which the expectation value is difficult to obtain analytically. Thus,20

we expand the term as a first-order Taylor expansion in terms of the CDF of normal distribution
and the logarithm function. In addition, we assume that the expectation of the third term canbe approximated by a linear approximation. Such bold approximations are standard strategies
for adapting topic models with VB to practical computation (for examples, zeroth-order Taylor
approximation by Asuncion et al. (2009) and Sato and Nakagawa (2012), and zeroth and firstorder delta approximation by Braun and McAuliffe (2010)). The four expectation values inEquation (A9) are then written asEq  log  (citz )   T *1
  y  0  x it μ zi,  1 cit22(z) 2 (z)(z) 2   x Tit β(citz )*   cit
Equ   ucit
 var  ucit

,




Equ ,q  u x β zi    x β(A10)

 x μ  ,E   x β    x V x   x μ  .


q

(z)
cit

T
it

T
it

2

T
it

T
it

zi

( z )*
cit



*

zi

it

(z)
cit

T
it

T
it

*
zi

* 2
zi

In appendix A.3, we find that μzi is updated to the optimized value of β cit . Finally, Ccitz is*( z )**

updated as

exp   citz 
*
Ccitz
 Z
,
 exp   citz 

(A11)

z 1

where





 citz    cz*     z 1  cz*  E q  log  (citz )  Z

1 T * (z)x it μ zi cit  x TitV zi *x it .
2

(A12)

A.3 Optimization of βcit

( z )*

Similar to equations (A3) and (A9), the expected value that optimizes βcit

( z )*is(z)E qu  log p  D, θ    Eqz ,q  log p  ucit
| β zi , zcit , x it , ycit    const.
2

1
(z)
 exp   *  ucit
 x Tit μ*zi    const.
 2Ccitz


(A13)

(z)

Here we seek the mean vector of the truncated normal distribution of ucit . Therefore, theupdate equation becomesβ(citz )*  μ*zi .

21

(A14)

A.4 Optimization of μzi and Vzi**

First, we derive an inverse Wishart distribution function and adopt some well-known
properties of multivariable normal and inverse Wishart distributions (Anderson 2003, Bishop
2006).
IW W , w  

w /2
W

2


wM

( w / 2)Viw  M 12

 1
 1   ,exp   tr WV
i
 2


M
 w*  1  m *1
EqV  log Vi      i
  M log 2  log Wi ,
2
m 1



(A15)

EqV Vi 1   wi*Wi *1 ,Eq ,qV   β zi  μi  Vi 1  β zi  μi     β zi  μi*  wi*Wi *1  β zi  μi *    i* .


T

T

We obtain the optimization procedures of μiz and Viz by the following expected value:**E q log p  D, θ    Eq ,qV  log p  β zi | μi , Vi  





(z)
 Equ ,qz log p ucit
 | β zi , zcit , xit , ycit    const.
1T  Eq ,qV  β zi  μi  Vi 1  β zi  μi  


2
C
2
1
(z)
  Equ ,qz   ucit
 x Tit β zi    const.

2 c 1 tTc

(A16)

The first and second terms of the second line are given by the last and third lines of equation(A10), while the third and fourth terms are given by Equations (A2) and (A3), respectively,
derived in a manner similar to (A9). μiz and Vzi are then arithmetically updated as**

μ*zi  wi*Wi*1  X zi X iT  wi*Wi*1μi*  X zi u zi 1Vzi*  wi*Wi*1  X zi X iT 

1

(A17)

where





T

( z)
 , X  x 
 *



u zi   E uciti
it c 1,C ,tTc  , X zi  Ccitz xit c 1,C ,tT  .
c 1,C ,tTc cThe uzi is vector and Xi and X zi are matrices. The number of elements in u i , Xi and

X zi are decided by the size of the consumer base and by Tc .
A.5 Optimization of μ i ,  i , wi , and Wi
*

*

*

*

Here we consider a joint distribution of a multivariable normal distribution of μ i and an
inverse Wishart distribution of Vi , and derive the update equations for four types of22variational parameters from this joint distribution. To this end, we require the following
expectation value from the joint distribution function:Eq ,qV log p  D, θ    log p  μi ,Vi   Eq log p β zi  | μi ,Vi   const.
T11w  M  1
1
 1 
log Vi  tr WV  log Vi   1  μi  μ   Vi 1  μi  μ   
i
2
2
2
2
1
1 Z
T
 Z  Eq log Vi    Eq  μi  β zi  Vi 1  μi  β zi   const.
2
2 z 1(A18)

First, we extract from this expectation value all terms linked to multivariable variational
*parameters μ iand  i ; that is*T1
E  q  log p  D, θ      1  μ i  μ   Vi 1  μ i  μ  
2
.
1 Z
T
1


  E q  μ i  β zi  Vi  μ i  β zi   const.

2 z 1(A19)

The second term in the above equation is obtained in the same manner as (A15). The
multivariable normal distribution function is then constructed in a straightforward manner asfollows:Z1 
μi*   1  Z    1μ    μ*zi  ,
z 1

 .

(A20)

 i     Z  .
*

1

1

*

*

Next, we optimize wi and Wi using Equation (A15) and the relationshiplog q Vi   log q  μi ,Vi   log q  μi | Vi  .
E  qV  log p  D , θ    E  q , qV  log p  D , θ    E  q  log p  D , θ  

(A21)

The expectation value E  q  log p  D , θ   is calculated in a straightforward manner by usingV**(A16) and (A17). Finally, we obtain the update equations for wi and Wi asZZ
Wi*  W   z1Vzi*   1μ  μ  T   z1 μ*zi μ*ziT   1  Z  μi*μi*T ,Vi  w  Z .
*Notice that  i

*

A(22)

*

and wi are constant if the hyperparameters and the number latent class are

given.

Appendix B: Gibbs Sampler
The joint posterior distribution, assuming conditional independence between variables,23

provides the full conditional posterior distributions:

Cc |  ~ p  Cc | zcit zcit |  ~ p  zcit | Cc ,{β zi },{xit },{ ycit }

(z)
(z)
ucit
|  ~ p  ucit
| zcit , β zi , xit , ycit 

(C1)

(z)
β zi |  ~ p  β zi | {ucit}, μi ,Vi , xit 

μi |  ~ p  μi | {β zi },Vi 
Vi |  ~ p Vi | {β zi }, μi 
where TN denotes a truncated normal distribution.
B.1 Sampling of C c
The C c is generated by a Dirichlet categorical relation. The Dirichlet distribution is aconjugate prior of a categorical distribution. For each consumer c, n c   nc1 ,  , ncZ Tdenotes the number of generated latent classes z c by categorical distribution of parameterC c in each MCMC step. A Dirichlet categorical relation gives the posterior distribution withrespect to C c asp  Cc |    p  Cc  p  zc | Cc   Diriclet  nc   

(C2)

B.2 Sampling of zcit | The posterior probability of ( zcit  j ) is given asPr zcit  j | Cc , xit  , β zi  ,  ycit  where Acitz  F  xTit β zi ycit

1  F  x β T
it

Ccj Acitj

 z1 Ccz Acitz
Z

,

(C3)

1 ycitzi.

B.3 Sampling of ucit | (z)The distribution of u c( itz ) is described in Appendix A.2. u c( itz ) is sampled from atruncated normal distribution in Equation (A5). This well-known sampling approach is called
data augmentation (Tanner, 1987).
B.4 Sampling of β zi , μ i , and V i
The full conditional posterior distribution of β iz , μ i , and V i

is derived from ahierarchical linear regression model. In our case, β zi for each i and each z is sampled from

 





βiz ~ N M R 1  X ziT u(ziz )   Vi 1μi , R 1 ,

24

(C4)

T
1
(z)
where R  X zi X zi  Vi , u (ziz )  ucit
czc  z ,tTc 


T

T

and X zi  xit cz z,tT  .



c

c



μ i is sampled from
Z
1
1


μ i ~ N M   Z      β zi , V i   Z     I M  ,
z 1



(C5)

~  0 0 0 0 .
for each i. Here, the hyperparameters are set to μ
TFinally, Vi for each i is sampled fromVi ~ IW  w  Z , W  B T B  ,

Z


1
where B    β zi  Z  β zi  .
z 1 
z 1

Z

25

(C6)

References[1] Ansari, A., and Mela, C. F. (2003). “E-Customization”. Journal of Marketing Research.
40, 131-145.
[2] Asuncion, A., Welling, M., Smyth, P. and Teh, Y.W. (2009). “On Smoothing andInference for Topic Models” Proceedings of the Twenty-Fifth Conference on Uncertainty
in Artificial Intelligence, 27-34.
[3] Bishop, C.M. (2006). Pattern Recognition and Machine Learning, Springer, U.S.A.
[4] Blattberg, R.C., Kim, B.D., and Neslin, S.A. (2009). Database Marketing: Analyzing andManaging Customers, Springer: PA.
[5] Blei, D.M., Ng, A.Y., and Jordan, M.I. (2003). “Latent Dirichlet Allocation,” Journal of
Machine Learning Research. 3, 993-1022.
[6] Blei, D., and McAuliffe. J. (2007). “Supervised Topic Models,” Proceedings of Neural
Information Processing System.
[7] Braun, M., and McAuliffe. J. (2010). “Variational Inference for Large-Scale Models of
Discrete Choice,” Journal of the American Statistical Association, 105, 324-335.
[8] Chintagunta, P.K., and Nair. H.S. (2011). “Discrete-Choice Models of Consumer
Demand in Marketing,” Marketing Science, 30, 977-996.
[9] Chung, T.S., Rust, R., and Wedel. M. (2009). “My Mobile Music: An Adaptive
Personalization System for Digital Audio Players,” Marketing Science, 28, 52-68.
[10] Corduneanu, A., and Bishop, C.M. (2001). “Variational Bayesian Model Selection for
Mixture Distributions. In: Jaakkola, T., Richardson, T. (Eds.)”, Artificial Intelligence andStatistics, Morgan Kaufmann: Los Altos, CA, 27–34.
[11] Grimmer, J. (2011). “An Introduction to Bayesian Inference via Variational
Approximations,” Political Analysis, 19, 32-47.
[12] Ishigaki, T., Takenaka T., and Motomura. Y. (2010). “Category Mining by
Heterogeneous Data Fusion Using PdLSI Model in a Retail Service,” Proceeding ofIEEE International Conference on Data Mining, 857-862[13] Iwata, T., Watanabe, S., Yamada, and T., Ueda, N,. (2009). “Topic Tracking Model for
Analyzing Consumer Purchase Behavior,” Proceeding of International Joint Conferenceon Artificial Intelligence, 1427-1432.
26

[14] Iwata, T., and Sawada, H., (2012). “Topic Model for Analyzing Purchase Data with
Price Information,” Data Mining and Knowledge Discovery, 26, 559-573.
[15] Jordan, M.I., Ghahramani, Z., Jaakkola, T.S., and Saul, L.K. (1999). An Introduction to
Variational Methods for Graphical Models,” Machine Learning, 37, 183-233.
[16] Naik, P., Wedel, M., Bacon, L., Bodapati, A., Bradlow, E., Kamakura, W., Kreulen, J.,
Lenk, P., Madigan and D.M., Montgomery, A. (2008). “Challenges and opportunities inhigh-dimensional choice data analyses,” Marketing Letter, 19, 201-213.
[17] Ramage, D., Hall, D., Nallapati, R., and Manning, C.D. (2009). “Labeled LDA: a
Supervised Topic Model for Credit Attribution in Multi-labeled Corpora,” Proceedingsof the 2009 Conference on Empirical Methods in Natural Language Processing, 248256.
[18] Rossi, P.E., Allenby, G.M., and McCulloch, R. (2005). Bayesian Statistics andMarketing. John Wiley & Sons: Chichester, UK.
[19] Rust, R.T. and Chung, T.S. (2005). “Marketing Models of Service and Relationships,”
Marketing Science, 25, 560–580.
[20] Spirling, A. and Quinn, K. (2010). “Identifying Intraparty Voting Blocs in the U.K.
House of Commons”, Journal of the American Statistical Association, 105, 447-457.
[21] Sato, I. and Nakagawa, H. (2012). “Rethinking Collapsed Variational Bayes Inference for
LDA,” Proceedings of International Conference on Machine Learning, 999-1006.
[22] Tanner, M.A. and Wong, W.H. (1987). “The Calculation of Posterior Distributions by
Data Augmentation,” Journal of the American Statistics Association, 82, 528-540.
[23] Teh, Y.W. and M. I. Jordan. (2010). “Hierarchical Bayesian nonparametric models withapplications,” Bayesian Nonparametrics: Principles and Practice, eds N. Hjort, C.
Holmes, P. Mueller, and S. Walker, Cambridge University Press, Cambridge, UK:, 2010.
[24] Tsiptsis, K. and Chorianopoulos, A. (2010). Data Mining Techniques in CRM: Inside
Customer Segmentation, Wiley: UK.
[25] Wedel, M. and Kamakura, W.A. (1999). Market Segmentation: Conceptual andMethodological Foundations. Kluwer Academic Publishers: U.S.A.

27

Table 1 Simulation time by VB and MCMC

I

Z

5

VB
10

100
500
1000

0.2
0.4
0.7

0.3
0.7
0.9

100
500
1000

0.6
1.3
2.4

0.9
2.1
4.1

100
500
1000

1.5
3.3
5.9

2.1
5.0
8.6

20
5
C =1000
0.4
4.3
1.2
17.6
1.5
37.8
C = 5000
1.5
17.8
3.9
50.5
6.4
99.6
C = 10000
3.7
38.8
8.4 164.3
14.7 333.3

MCMC
10

20

6.0
22.8
41.6

12.2
32.1
48.9

23.3
62.6
111.0

35.8
80.9
128.2

52.3
201.0
372.7

78.5
264.0
450.2

The number means hour.

Table 2 RMSE as a function of DDR in simulations

RandomHomogeneous
VB
MCMC

0.1%

0.2%

0.245
0.227
0.226

0.236
0.213
0.214

DDR
0.3%
0.577
0.238
0.212
0.211

0.4%

0.5%

0.236
0.211
0.212

0.237
0.212
0.212

Simulated data (C = 500, I = 100).

Table 3 RMSEs for real customer database – all and infrequent customers# of Z
2
RandomHomogeneous
All customers
Infrequentcustomers

3

4

5

10

20

0.383

0.383

0.393 0.393

0.393

0.404

0.389

0.577
0.412
0.385 0.383

0.410

0.400

0.395

Real customer database (C = 1647, I = 1004)

28

Table 4 Characteristics of β zi for the five segments of consumers
Intercept
segment 1 segment 2 segment 3 segment 4 segment 5
111
-0.82
3.56
4.17
3.76
4.09
205
2.38
2.48
-1.60
2.68
2.57
153
1.44
1.25
1.48
1.52
-1.46
253
-1.42
-3.19
-1.43
-1.48
-1.52
120
-0.46
-0.50
-1.34
-0.41
-1.85
...
...
...
...
...
...
853
-0.89
-0.90
-0.90
-0.90
-0.90
1002
-0.73
-0.72
-0.73
-0.73
-0.73
822
-0.03
-0.03
-0.03
-0.02
-0.02
479
-0.13
-0.13
-0.13
-0.12
-0.13
166
0.05
0.06
0.06
0.06
0.06
No.

No.
195
18
182
225
47
...
936
483
871
181
794

Display
segment 1 segment 2 segment 3 segment 4 segment 5
1.26
1.32
1.29
1.31
2.74
1.35
0.14
0.11
0.14
0.12
1.28
1.18
1.21
2.27
1.22
1.57
0.62
0.61
0.63
1.33
0.66
1.59
0.66
0.62
0.58
...
...
...
...
...
-0.09
-0.08
-0.10
-0.10
-0.10
-0.07
-0.07
-0.06
-0.06
-0.07
0.32
0.32
0.33
0.33
0.33
0.30
0.31
0.31
0.29
0.31
0.70
0.71
0.70
0.71
0.70

No.
205
111
153
120
147
...
764
608
479
556
737
No.
147
96
558
502
163
...
457
907
24
187
166

Price
segment 1 segment 2 segment 3 segment 4 segment 5
-4.38
-4.53
-0.13
-4.74
-4.59
-1.30
-5.04
-5.06
-5.21
-5.89
-3.95
-3.08
-4.08
-4.08
-0.24
-2.51
-2.47
-3.75
-2.66
-0.15
-2.13
-2.02
-2.10
-1.94
-4.71
...
...
...
...
...
-0.52
-0.53
-0.52
-0.53
-0.53
-0.63
-0.64
-0.63
-0.63
-0.63
-0.03
-0.03
-0.03
-0.03
-0.04
-0.42
-0.42
-0.43
-0.42
-0.42
-0.57
-0.56
-0.57
-0.56
-0.57
Featuresegment 1 segment 2 segment 3 segment 4 segment 5
1.06
0.96
1.01
0.91
2.71
2.54
0.90
0.92
0.93
0.96
1.02
1.90
2.51
0.99
1.03
1.03
1.11
1.19
2.53
1.12
0.57
0.50
1.73
0.58
0.56
...
...
...
...
...
0.23
0.22
0.21
0.21
0.23
0.03
0.04
0.05
0.03
0.04
0.06
0.05
0.06
0.05
0.06
0.20
0.20
0.20
0.20
0.21
0.05
0.06
0.06
0.06
0.06

Table 5 Relative preferences of purchased product category for five segments
Segment 1
DessertInstant noodle
Dressing
CookieCoffee
Sauce
YoghurtFish sausage
Dry noodle
Frozen noodle

#
7
4
4
4
3
3
3
3
3
3

Segment 4
Instant noodle
Detergent
Milk
Fresh noodle
Japanese tea
CoffeeCoke
DessertYoghurtDry noodle
Fizzy drink

#
6
4
4
4
4
3
3
3
3
3
3

Segment 2
DessertDry noodle
Chocolate
Snak
Snacks made from rice
Instant noodle
Coke
Detergent
Fish sausage
Milk
Soy sauce
Segment 5
YoghurtInstant noodle
Coke
Fish sausage
Milk
Tea
Snacks made from rice
Japanease Tea
Beans paste
Detergent

29

#
6
4
4
4
4
3
3
3
3
3
3
#
8
6
4
4
4
3
3
3
3
3

Segment 3
Yoghurt
Noodle
Chocolate
Coke
Dressing
Detergent
Softener
Fresh noodle

#
6
4
4
4
3
3
3
3

18
*
*
*
*
*

Display
Product No.
18 110 253 318 742
(a) *
*
*
*
(b)
*
*
*
(c) *
*
*
(d) *
*
*
*
(e) *
*
*

Customer

(a)
(b)
(c)
(d)
(e)

Price
Product No.
110 253 318 742
*
*
*
*
*
*
*
*
*
*
*
*
*
*

Customer

CustomerTable 6 Personalized effective marketing variables for individual consumers and products.
FeatureProduct No.
18 110 253 318 742
*
*
*
*
*
*
*
*
*
*
*

(a)
(b)
(c)
(d)
(e)

Squares marked with * indicate that a consumer (row) is likely to purchase a product(column) based on the marketing variable (price, display, or feature)

Table 7 Estimated parameters by VB and MCMC
No.
1
2
3
4
5
...
95
96
97
99
100

Intercept
VB
MCMC
3.29
1.35
3.95
2.38
-0.12
-0.54
3.61
2.67
0.56
0.60
...
...
-0.13
-0.02
-0.67
-0.88
-0.52
-0.86
1.30
0.34
-0.69
-0.81

VB
-4.91
-5.38
-0.81
-5.33
-2.58
...
-0.44
-2.04
-1.80
-3.08
-0.64

Price
MCMC
-3.42
-4.18
-0.26
-4.26
-2.55
...
-0.50
-1.31
-1.28
-1.97
-0.51

Display
VB
MCMC
0.28
0.29
...
...
0.40
0.14
0.26
0.24
0.17
0.12
0.43
0.44

Real customer database (C = 500, I = 100).

30

Feature
VB
MCMC
0.17
0.44
0.44
0.47
0.71
0.96
0.75
0.74
0.10
0.12
...
...
0.31
0.32
0.93
0.98
1.10
1.21
-0.11
-0.24

Figure 1 Marginal distribution of parameter estimates of individual consumers and products31