Co-occurrence Estimation from Aggregated Data with Auxiliary Information
Tomoharu Iwata, Naoki Marumo
NTT Communication Science Laboratories, Kyoto, Japan{tomoharu.iwata.gy, naoki.marumo.ec}@hco.ntt.co.jp

Abstract
Complete co-occurrence data are unavailable in many applications, including purchase records and medical histories, because of their high cost or privacy protection. Even with suchapplications, aggregated data would be available, such as thenumber of purchasers for each item and the number of patients with each disease. We propose a method for estimating
the co-occurrence of items from aggregated data with auxiliary information. For auxiliary information, we use item features that describe the characteristics of each item. Althoughmany methods have been proposed for estimating the cooccurrence given aggregated data, no existing method canuse auxiliary information. We also use records of a small
number of users. With our proposed method, we introduce
latent co-occurrence variables that represent the amount ofco-occurrence for each pair of items. We model a probabilistic generative process of the latent co-occurrence variables
by a multinomial distribution with Dirichlet priors. The parameters of the Dirichlet priors are parameterized with neuralnetworks that take the auxiliary information as input, where
neural networks are shared across different item pairs. The
shared neural networks enable us to learn unknown relationships between auxiliary information and co-occurrence using
the data of multiple items. The latent co-occurrence variables
and the neural network parameters are estimated by maximizing the sum of the likelihood of the latent co-occurrence variables and the likelihood of the small records. We demonstrate
the effectiveness of our proposed method using user-item rating datasets.

1

Introduction

Co-occurrence is the basic and important statistics for analyzing categorical data. For example, recommender systems in e-commerce services suggest items that are likely
to be purchased by the same user (Sarwar et al. 2001). With
text analysis, words that appear in the same document are
clustered to discover topics (Blei, Ng, and Jordan 2003).
With social network analysis, communities are detected using information about common friends (Girvan and Newman 2002). Medical knowledge about complications is essential for treatment. Complete co-occurrence data are unavailable in many applications, such as purchase records and
Copyright c 2020, Association for the Advancement of Artificial
Intelligence (www.aaai.org). All rights reserved.

medical histories, because of their high cost or privacy protection. However, even within such applications, aggregateddata would be available, such as the number of purchasers
for each item and the number of patients for each disease,since they do not contain any privacy information.
In this paper, we propose a method that estimates cooccurrence from aggregated data with auxiliary information.
For auxiliary information, we use item features. In the caseof purchase records, item features might include genres, release dates and descriptions. Since item features do not contain privacy information, they are often open to the public.
Although many methods have been proposed for estimatingco-occurrence counts given aggregated data (Deming andStephan 1940; Causey 1983; Sheldon and Dietterich 2011),
no existing method can use auxiliary information. We also
use the records of a few users that are sampled from total
users. The sampled records can be obtained by giving incentives to some users if they agree to share their records. We
assume that the number of sampled users with their records
is much smaller than the total number of users.
With the proposed method, we introduce latent cooccurrence variables, which represent the amount of cooccurrence for each pair of items, or the number of users
who purchased both items. The joint probability of the cooccurrence count of two items is assumed to follow a multinomial distribution with Dirichlet priors. The parameters of
the Dirichlet priors are modeled by non-linear functions thattake the auxiliary information as input. For non-linear functions, we use permutation invariant neural networks sinceco-occurrence data are invariant to the permutation of two
items. The neural network parameters are shared across different pairs of items. Shared neural networks enable us tolearn unknown relationships between auxiliary informationand the co-occurrence using the data of multiple items.
The latent co-occurrence variables and the neural network parameters are estimated by maximizing the sum of
the likelihood of the latent co-occurrence variables and thelikelihood of the small sampled records. The parameters ofthe multinomial distributions are analytically marginalized
since we use conjugate Dirichlet priors. By relaxing that thelatent co-occurrence variables are non-negative real values
instead of integers, we efficiently maximize the objective

function with stochastic gradient descent methods.

2

Related work

Previous work proposed learning from aggregated data, suchas imputing individual level records (Park and Ghosh 2012;2014) and regression from aggregated data (Goodman 1953;
Freedman et al. 1991; Bhowmik, Ghosh, and Koyejo 2015;2016). In this paper, we focus on estimating co-occurrence
counts, which can be represented by two-by-two contingency tables. Many methods for estimating contingency tables have been proposed especially in statistics (Irelandand Kullback 1968; Smith 1947; Fienberg 1970; Slavkovic2010), such as least squares (Deming and Stephan 1940;
Stephan 1942), maximum likelihood (Causey 1983) andthe Markov chain Monte Carlo (Dobra, Tebaldi, and West2006), as well as in machine learning (Sheldon and Dietterich 2011; Kumar, Sheldon, and Srivastava 2013; Sheldon et al. 2013; Sun, Sheldon, and Kumar 2015; Nguyen
et al. 2016). Most of these methods use the sampled dataof a small population and/or the aggregated counts for each
item. For other information, samples from different populations (Little and Wu 1991), conditional tables (Slavkovic2010) and short and long form questionnaires (Greco 2016)
have been used. However, these existing methods cannot use
auxiliary item features. Although neural collective graphical models (Iwata and Shimizu 2019) estimate people flowfrom aggregated data using spatio-temporal auxiliary information, they are specialized for people flow, and inapplicable to co-occurrence estimation. The proposed method canbe seen as a multi-task learning approach (Caruana 1997) for
simultaneously multiple contingency tables by transferringknowledge across different two-by-two contingency tables
using shared neural networks.

3
3.1

I
U
yi
si
R
πij

Proposed method

Problem formulation

For simplicity, we explain the proposed method with an example of item co-occurrence in user purchase records. However, our proposed method is applicable to any kinds of cooccurrences, such as word co-occurrence in documents anddisease complications in patients.
Assume I items and U users, where each user has purchased some items but we cannot observe their complete
purchase records. We are given marginal counts for each
item, y = {yi }Ii=1 , where yi is the number of users who
have purchased item i. For auxiliary information, we are
given item features S = {si }Ii=1 , and item purchase records∗∗
R = {ru }U U.
u=1 for a small number of users U
D
Here, si ∈ R is a D-dimensional item feature vector
of item i, such as genres, release dates and descriptions,and ru ∈ {0, 1}I is an I-dimensional binary vector, whererui = 1 if user u has purchased item i and rui = 0 otherwise.
Our task is to estimate the joint probability of the occurrence of two items, πij = (πīj̄ , πīj , πij̄ , πij ), for all pairs of
items, i, j ∈ {1, · · · , I}, where πīj̄ = P (ri = 0, rj = 0)
is the probability that a user purchases neither items i nor j,
πīj = P (ri = 0, rj = 1) is the probability that a user does

θ̂i
xij
Ψ

Table 1: Our notation.
number of items
number of users
number of users who have purchased item i
auxiliary feature vector of item ipurchase records of a few sampled users
probability that a user purchases
both items i and jempirical probability that item i is purchased
latent number of users who have purchased
both items i and jparameters of neural networks

Table 2: Co-occurrence count matrix xij between items iand j, number of purchasers yi and yj , and total number ofusers.
Item i \ Item j Not purchased Purchased
Total
Not purchased
xīj̄
xīj
U − yi
Purchased
xij̄
xij
yi
Total
U − yjyj
U
not purchase item i but buys item j, πij̄ = P (ri = 1, rj =0) is the probability that a user purchases item i but not itemj, πij = P (ri = 1, rj = 1) is
Pthe probability that a userpurchases both items i and j, i0 ∈{i,ī},j 0 ∈{j,j̄} πi0 j 0 = 1,
and πi0 j 0 ≥ 0. Table 1 summarizes our notation.

3.2

Model

We introduce latent co-occurrence variables xij =(xīj̄ , xīj , xij̄ , xij ), where xīj̄ is the number of users who
have purchased neither items i nor j, xīj is the number ofusers who have purchased item j but not i, xij̄ is the number of users who have purchased item i but not j, and xij isthe number of users who have purchased both items i and j.
Table 2 shows a co-occurrence count matrix that represents
the relationship among latent co-occurrence variables xij ,item marginal counts yi , yj and the total number of users U .
Given xij , other latent co-occurrence variables xīj̄ , xīj andxij̄ are determined using yi , yj and U as follows,xīj = yj − xij ,xij̄ = yi − xij ,xīj̄ = U − yi − yj + xij .

(1)

Therefore, we do not need to estimate xīj̄ , xīj , xij̄ if xij isestimated. Using Eq.(1) and the non-negativity of all the latent co-occurrence variables, xīj̄ , xīj , xij̄ , xij , the following
constraint on xij is derived,max(0, yi + yj − U ) ≤ xij ≤ min(yi , yj ),(2)which is known as Fréchet inequalities. Range,min(yi , yj ) − max(0, yi + yj − U ), indicates information about co-occurrence xij of aggregated data yi andyj . When the range is narrow, the aggregated data are
informative for estimating co-occurrence. However, since
the range is not generally narrow, the co-occurrence cannot
be determined only from the aggregated data.

We assume that the probability of latent co-occurrence
variables xij is given by the following multinomial distribution with parameters πij ,
Y
U!
x0 0
p(xij |πij ) = P
πi0ij 0j .
0j0 !
x
0
0
i
i ∈{i,ī},j ∈{j,j̄}
0
0

α
si
sj

neural
net

θ̂i
βij

θ̂j

U
πij

xij

yj

i ∈{i,ī},j ∈{j,j̄}(3)Note that a multinomial distribution is derived from the Poisson distributions for each cell when total count U is fixed.
We assume the following Dirichlet distribution for the priorof co-occurrence probabilities πij ,P
Γ( i0 ∈{i,ī},j 0 ∈{j,j̄} βi0 j 0 )p(πij |βij ) = Q
0 0
i0 ∈{i,ī},j 0 ∈{j,j̄} Γ(βi j )Y
βi0 j 0 −1×
π i0 j 0,(4)
i0 ∈{i,ī},j 0 ∈{j,j̄}

where βij = (βīj̄ , βīj , βij̄ , βij ), βi0 j 0 > 0, and Γ(·) is
the Gamma function. We model Dirichlet parameters βij
as a function of auxiliary information si and sj to incorporate them for estimating the co-occurrence probabilities.
By modeling the Dirichlet parameters, we can handle thevariance of the co-occurrence probabilities as well as their
β0 0mean, where the mean is E[πi0 j 0 ] = βiij0j and the varianceP
β 0 0 (βij0 −βi0 j 0 ), and βij0 = i0 j 0 βi0 j 0 . In
is Var[πi0 j 0 ] = iβj2 (βij0 +1)
ij0
particular, we use the following neural network-based models,βīj̄ = α(1 − θ̂i )(1 − θ̂j ) + f0 (si , sj ),βīj = α(1 − θ̂i )θ̂j + f01 (si , sj ),βij̄ = αθ̂i (1 − θ̂j ) + f01 (sj , si ),βij = αθ̂i θ̂j + f1 (si , sj ),(5)where θ̂i = yUi is the empirical marginal occurrence probability of item i, α > 0 is the scalar positive parameter, and
f0 (·), f01 (·), f1 (·) are functions modeled with neural networks. The first term in Eq.(5) corresponds to the probability
when the occurrences of items i and j are independent, andα controls the strength of independence. Since a simple assumption to estimate the joint probabilities is independence
when there is no information about the co-occurrence, we
include this first term. The second term controls the correlation between the occurrence of two items using auxiliary
information. Neural networks enable us to flexibly model the
relationships between the correlation and auxiliary information. Even when items i and j are transposed, the prior probability should be invariant. Therefore, we use f01 (sj , si ) for
βij̄ , where input auxiliary information si and si are transposed from f01 (si , sj ) for βīj . Also, f0 (·) and f1 (·) need to
be invariant to the permutation of inputs si and sj ,f0 (si , sj ) = f0 (sj , si ), f1 (si , sj ) = f1 (sj , si ).
(6)
To satisfy the above invariance, we use the following permutation invariant networks (Zaheer et al. 2017),
f0 (si , sj ) = ρ0 (φ0 (si ) + φ0 (sj )),f1 (si , sj ) = ρ1 (φ1 (si ) + φ1 (sj )),(7)

yiFigure 1: Proposed model: Shaded and unshaded nodes indicate observed and latent variables. Neural networks takeitem feature vectors si and sj as input. Dirichlet parameters
βij are determined by α, empirical marginal count probabilities θ̂i and θ̂j , and neural network output. Co-occurrence
probabilities πij are generated from Dirichlet distribution
with parameters βij , and latent co-occurrence variables xij
are generated from multinomial distribution with parameters
πij . Total count U = xīj̄ + xīj + xij̄ + xij , marginal counts
yi = xij̄ + xij and yj = xīj + xij are determined by latentco-occurrence variables xij = (xīj̄ , xīj , xij̄ , xij ).
where ρ0 (·), φ0 (·), ρ1 (·), φ1 (·) are neural networks. Since
addition is invariant to permutation, neural networks withstructure of Eq.(7) output the same values even when inputs
si and sj are permuted as in Eq.(6). Fig. 1 illustrates the
proposed model.

3.3

Estimation

First, we count the co-occurrences in small purchase records
R, and represent them by X∗ = {x∗ij }Ii,j=1 , where x∗ij =(x∗īj̄ , x∗īj , x∗ij̄ , x∗ij ), and x∗ij is the number of users withrui = 1 and ruj = 1 in R, which are the users who
have purchased both items i and j in the small purchaserecords. Then, we estimate latent co-occurrence variables
X = {xij }Ii,j=1 , parameter α, and parameters Ψ of neural networks f0 (·), f1 (·), f01 (·) by maximizing the weighted
sum of the likelihoods of X and X∗ .
To automatically satisfy the box constraint in Eq.(2), we
parameterize xij by
xij = max(0, yi + yj − U )min(yi , yj ) − max(0, yi + yj − U )+,1 + exp(−x0ij )(8)where −∞ < x0ij < ∞, and optimize x0ij , instead of xij ,without constraints. Also, we parameterized positive parameter α by α = exp(α0 ), where −∞ < α0 < ∞, for unconstrained optimization.
Since the Dirichlet distribution is the conjugate prior of
the multinomial distribution parameters, we can analyticallymarginalize co-occurrence probabilities πij in Eqs.(3) and(4). Then the probability of xij given βij is as follows,P
U !Γ( i0 ,j 0 βi0 j 0 ) Y Γ(xi0 j 0 + βi0 j 0 )P,p(xij |βij ) =Γ(U + i0 ,j 0 βi0 j 0 ) 0 0 xi0 j 0 !Γ(βi0 j 0 )
i ,j

(9)

Algorithm 1: Estimation procedure with the proposed method.
Input: marginal counts y, total count U , auxiliary
information S, small records R,
hyperparameter λ, batchsize T
Output: estimation of co-occurrence counts X,
neural network parameters Ψ, parameter α
1 Initialize parameters X, Ψ, α;
2 Calculate empirical marginal occurrence probabilityθ̂i = yUi for all items using y;∗
3 Calculate co-occurrence counts for small records X
using R;
4 repeat5
Sample a set of T item pairs Q randomly;
6
Calculate the objective function Eq.(10) for the
sampled item pairs and its gradients Updateparameters X, Ψ, α with a gradient-based
optimization method;
7 until end condition is satisfied;

Table 3: Comparing methods. Methods with check-marks
indicate that they use marginal counts y, small records R,
or auxiliary information S.
IND
EB+A CGM+A
Ours√ ML EB CGM√ EB+CGM
√
√
√
y
√ √
√
√
√
R
√
√
√
S

where i0 ∈ {i, ī} and j 0 ∈ {j, j̄}. This is called the Dirichlet compound multinomial distribution or the multivariatePólya distribution.
The objective function to be maximized is the following
weighted sum of the likelihoods of latent co-occurrence variables X and observed small co-occurrence X∗ ,L(X, Ψ, α) = λ

I
I
X
X

log p(xij |βij )i=1 j=i+1+ (1 − λ)IIXXlog p(x∗ij |βij ),(10)

i=1 j=i+1

where Dirichlet parameters βij are calculated using neuralnetworks by Eq.(5), λ is the hyperparameter, and 0 ≤ λ ≤ 1.
We used fixed hyperparameter λ = 0.5 in our experiments.
By relaxing that the latent co-occurrence variables are nonnegative real values instead of integers, we can efficiently
find a local maximum of the objective function with stochastic gradient-based optimization methods. Algorithm 1 showsthe estimation procedure with the proposed method.

4
4.1

Experiments

Dataand generated occurrence data. For the auxiliary item features, we used style (maki or not), major group (seafood ornot), minor group, heaviness/oiliness in taste, frequency of
consumption, price, and frequency of being sold. We transformed the categorical features, i.e., style, major and minorgroups, into one-hot vectors, and normalized the numericalfeatures, i.e., heaviness, frequency, and price, into a range
from zero to one. The dimensionality of the item feature vectors was D = 20.
The Movielens data were generated from 100,000 ratingsof 943 users of 1,682 movies (Konstan et al. 1997)2 . Whenuser u rated movie i, we assumed that she has purchased theitem rui = 1. For the auxiliary item features, we used therelease date and genres, such as action, comedy, and fantasy.
We normalized the release dates into a range from zero toone, and transformed the genres into vectors by setting thecorresponding element to one when the movie was in thegenre, and zero otherwise. The dimensionality of the itemfeature vectors was D = 20.
With both data sets, we generated item count yi for each
item i using all the users. We randomly sampled the purchase records of ten users for the validation data.

4.2

Evaluation measurement

For the evaluation measurement, we used the following absolute error of the co-occurrence probabilities,
E=

I
I
X
X
2
I(I − 1) i=1 j=i+1 0X|π̂i0 j 0 − πiTRUE|,
0j0i ∈{i,ī},j 0 ∈{j,j̄}(11)where πiTRUE
is the true empirical co-occurrence probabil0j0
ity calculated using all the users, which is not used during
training, and π̂i0 j 0 is the estimated co-occurrence probability.

4.3

Comparing methods

We evaluated the proposed method using sushi and Movielens data sets.
The sushi data were generated from the preferences of
5,000 users among the following ten kinds of sushi: shrimp,
sea eel, tuna, squid, sea urchin, salmon roe, egg, fatty tuna,
tuna roll, and cucumber roll (Kamishima 2003)1 . When user
u ranked sushi item i within the top five among the ten items,
we assumed that the user has purchased the item rui = 1,

We compared the proposed method with the following
seven methods: IND, ML, EB, CGM, EB+CGM, EB+A and
CGM+A. We summarized the characteristics of the comparing methods in Table 3.
The IND method estimates the co-occurrence probability by assuming that items occur independently as follows,π̂īj̄ = (1 − θˆi )(1 − θˆj ), π̂īj = (1 − θˆi )θˆj , π̂ij̄ = θˆi (1 − θˆj ),π̂ij = θˆi θˆj .

1
The original sushi data were obtained from http://www.
kamishima.net/sushi/.

2
The original Movielens data were obtained from https://grouplens.org/datasets/movielens/.

The ML method estimates the probability by the maximum likelihood of the small record data as follows, π̂īj̄ =
x∗
x∗
x∗
x∗
īj̄
īj
ij̄
ij
∗
U ∗ , π̂īj = U ∗ , π̂ij̄ = U ∗ , π̂ij = U ∗ , where U is the total

number of users in the small record data.
The EB method is the empirical Bayesian method thatmaximizes the P
marginalized likelihood of the small record
data, L(B) = i,j log p(x∗ij |βij ), where B = {βij }Ii,j=1
P
is the set of Dirichlet parameters, andi,j denotes
PI PI
.
i=1
j=i+1
The CGM method is the collective graphical model
based method (Kumar, Sheldon, and Srivastava 2013) thatmaximizes the marginalized likelihood of the aggregateddata
P with the constraints on marginal counts, L(X, B) =i,j log p(xij |βij ).
The EB+CGM method is the combination of the
EB and CGM methods, where the weighted sum of
the marginalized likelihood ofPthe aggregated and small
record data, L(X, B) = λ i,j log p(xij |βij ) + (1 −PI PIλ) i=1 j=i+1 log p(x∗ij |βij ), is maximized. We used thefixed weight hyperparameter λ = 0.5 with the EB+CGM
and proposed methods. With the EB, CGM and EB+CGM
methods, Dirichlet parameters βij are directly optimized
without neural networks.
The EB+A method is the empirical Bayesian method thatuses the auxiliary data, where neural networks are trained by
maximizing the marginalized
likelihood of the small record
P∗
data, L(Ψ, α) =logp(xij |βij ), where Dirichlet pai,jrameters βij are parameterized by neural networks in Eq.(5).
The CGM+A is the collective graphical model based
method that uses the auxiliary data, where neural networks
are trained by maximizing the marginalized
likelihood of thePsmall record data, L(X, Ψ, α) = i,j log p(xij |βij ). Here,as with the EB+A method, Dirichlet parameters βij are parameterized by neural networks in Eq.(5).
The EB+A, CGM+A and proposed methods use the auxiliary information, and the other methods do not use it. After
estimating the latent co-occurrence variables, we estimatedx̂
the co-occurrence probability by π̂ij = Uij with the CGM,
EB+CGM, CGM+A and proposed methods, where x̂ij is theestimate of latent co-occurrence variable xij . Since the EB
and EB+A methods did not estimate co-occurrence xij , they
estimated the co-occurrence probability by the expectationβ̂with learned prior π̂ij = P ijβ̂ , where β̂ij is the estii0 j 0

i0 j 0mated Dirichlet parameter.
For neural network ρ0 (·), ρ1 (·), φo (·), φ1 (·), f01 (·) in the
EB+A, CGM+A and proposed methods, we used threelayer feed-forward neural networks with 32 hidden units.
Since βi0 j 0 must be non-negative, we use rectified linear unit
ReLU(x) = max(0, x) at the end of neural networks ρ0 (·),
ρ1 (·) and f01 (·). We optimized using ADAM (Kingma andBa 2015) with learning rate 10−2 , weight decay 10−2 anddropout rate 0.1. For each batch, we randomly sampled 512
item pairs in training. The validation data were used for early
stopping, where the maximum number of training epochs
was 1,000. We implemented all the methods based on PyTorch (Paszke et al. 2017).

4.4

Results

Figure 2(a) shows the absolute error of the co-occurrence
probabilities averaged over 30 experiments with different
numbers of sampled users with records U ∗ , where the number of items was I = 10 with the sushi data, and I = 200
with the Movielens data. Here, items were randomly sampled from all items with the Movielens data. The proposed
method achieved the lowest error for all the cases. The error by the proposed method was statistically lower at the
5% level than that by the other methods by a paired t-test
with all the cases in the sushi data and with all the cases
except for U ∗ = 30 in the Movielens data. This result indicates that the proposed method improved the estimation performance using all of the aggregated marginal count data, asmall number of sampled records, and auxiliary item featuredata. The error decreased as the number of observed users
U ∗ increased with the proposed method, since the estimationof the Dirichlet parameters modeled by the neural networks
became robust with more training data. Although the performance of the ML, EB, and EB+A methods also improved as
the number of observed users increased, their performance
was much worse than the proposed method. This result occurred because they did not use aggregated data information.
In contrast, the proposed method effectively used the aggregated data information in the first term of the objective function in Eq.(10). The reduction of the error of the proposed
method from EB+CGM, which corresponds to the proposed
method without auxiliary information, was 34% and 8% on
Sushi data and Movielens data with U ∗ = 20, respectively.
This result shows the importance to use the auxiliary information for improving performance.
Table 4 shows the average absolute errors for each pair of
items on the sushi data with U ∗ = 20 and I = 10. The IND
method achieved the lowest error with some pairs, where
the user preferences were assumed to be independent. The
performance of the CGM, EB+CGM, and CGM+A methods resembled the IND method, although they were slightly
worse. The proposed method achieved the lowest error withmany pairs, which include item pairs whose preferenceswere not independent. Table 5 shows the rate of item pairsthat each method achieved the best on the Movielens datawith U ∗ = 20 and I = 200. The rate of the the proposed
method was highest.
Figure 2(b) shows the absolute error with different numbers of items I, where the number of sampled users was
U ∗ = 20 with both the sushi and Movielens data, and we
randomly sampled I items from all items. With the sushi
data, the performance of the proposed method was not thebest when the number of items was small. However, its performance improved as the number of items increased; it was
statistically lower at the 5% level than that by the othermethods with I ≥ 6. Since the data for training the neural networks increased with more items, the performance
improved. With the Movielens data, the proposed method
achieved statistically less error than the other methods for
all cases.
Figure 2(c) shows the absolute error with different hyperparameter values λ with the proposed method, where thenumber of sampled users was U ∗ = 20 with both kinds ofSushi

error

0.8
0.6

0.7
0.6

0.6

0.5

0.4

error

IND
ML
EB
CGM
EB+CGM
EB+A
CGM+A
Proposed

error

1.0

0.4

0.4

0.3
0.2

0.2

0.2

0.1
20

40

60

U*

80

100

2

4

6
I

8

10

0.0

0.2

0.4

0.6

0.8

1.0

0.6

0.8

1.0

λ

Movielens

0.20
0.15

0.12

0.200
0.175

error

error

0.25

0.225

error

IND
ML
EB
CGM
EB+CGM
EB+A
CGM+A
Proposed

0.30

0.150
0.125

0.08

0.100

0.10

0.075

0.05

0.10

0.06

0.050

20

40

60

U*

80

(a) #Users with records U ∗

100

100

150

200
I

250

0.0

300

0.2

0.4
λ

(b) #Items I(c) Hyperparameter λFigure 2: Absolute error of co-occurrence probabilities averaged over 30 experiments with (a) different numbers of sampledusers with records U ∗ , (b) different numbers of items I, and (c) different hyperparameters λ with proposed method. Error bars
show their standard errors.
data, and the number of items was I = 10 with the sushi
data, and I = 200 with the Movielens data. The proposed
method corresponds to the EB+A method when λ = 0, andto the CGM+A method when λ = 1. The best hyperparameter value was λ = 0.4 with the sushi data, and λ = 0.1
with the Movielens data. The co-occurrence data are latent
xij with the first term of the objective function in Eq.(10),but the co-occurrence data are observed x∗ij with the second
term. Therefore, the best hyperparameter values put moreweight on the second term than the first term.
The average computational time with the proposed
method for training using the sushi data with I = 10, Movielens with I = 100, and Movielens with I = 200 were0.12, 5.95 and 24.19 minutes, respectively, on computerswith 2.60GHz CPUs, where we fixed U ∗ = 20. The computational time increased quadratically with the number of
items since we need to consider all pairs of items. The computational time did not depend on the number of users, U or
U ∗.
We also evaluated with Movielens-1M data, which contained 6,040 users. We used the genres as auxiliary information, U ∗ = 10 users with records, and I = 100 items
that were randomly selected from movies rated by more than
300 users. Table 6 shows the results. The proposed method
achieved the lowest error with the Movielens-1M data.

5

Conclusion

We proposed a method for estimating co-occurrence withaggregated marginal count data, the records of a small
number of users, and auxiliary item feature information.
Our proposed method learns unknown relationships between auxiliary information and co-occurrence using neural networks for the parameters of the Dirichlet priors of
joint co-occurrence probabilities. Experiments on real-world
datasets confirmed that the proposed method achieved better estimation performance for co-occurrence than existing
methods. Although our results are encouraging, we must extend our approach in several directions. First, we will apply our framework to general multi-way and/or multivariate
contingency tables; we focused on two-by-two contingency
tables in this paper. Second, we will investigate tuning thehyperparameters using validation data or a Bayesian framework.

References
Bhowmik, A.; Ghosh, J.; and Koyejo, O. 2015. Generalized
linear models for aggregated data. In Artificial Intelligenceand Statistics, 93–101.
Bhowmik, A.; Ghosh, J.; and Koyejo, O. 2016. Sparse
parameter recovery from aggregated data. In International
Conference on Machine Learning, 1090–1099.
Blei, D. M.; Ng, A. Y.; and Jordan, M. I. 2003. LatentTable 4: Absolute error of co-occurrence probabilities averaged over 30 experiments for each pair of items with the sushi data,
where the number of sampled users with records was U ∗ = 20 and number of items was I = 10. Values in bold typeface are
not statistically different (at 5% level) from the best performing method in each row according to a paired t-test. Bottom line
shows number of pairs for which the method achieved the best.
Item i
Item j
IND
ML
EB
CGM
EB+CGM
EB+A
CGM+A
Proposed
shrimpsea eel
0.129
0.285
0.233
0.138
0.120
0.673
0.126
0.027
shrimp
tuna
0.136
0.242
0.219
0.098
0.100
0.726
0.133
0.026
shrimp
squid
0.035
0.285
0.237
0.062
0.082
0.634
0.029
0.165
shrimpsea urchin
0.161
0.289
0.244
0.164
0.144
0.683
0.158
0.026
shrimpsalmon roe
0.190
0.290
0.245
0.207
0.166
0.722
0.186
0.044
shrimp
egg
0.039
0.305
0.265
0.079
0.075
0.617
0.043
0.052
shrimpfatty tuna
0.137
0.257
0.235
0.120
0.124
0.769
0.135
0.101
shrimptuna roll
0.200
0.279
0.231
0.182
0.178
0.658
0.206
0.098
shrimpcucumber roll
0.002
0.242
0.240
0.015
0.028
0.524
0.003
0.039
sea eel
tuna
0.221
0.284
0.243
0.172
0.179
0.760
0.218
0.098
sea eel
squid
0.157
0.316
0.267
0.146
0.124
0.660
0.163
0.024
sea eelsea urchin
0.017
0.316
0.259
0.020
0.049
0.663
0.014
0.147
sea eelsalmon roe
0.126
0.279
0.235
0.137
0.147
0.706
0.123
0.036
sea eel
egg
0.011
0.313
0.280
0.045
0.054
0.587
0.015
0.080
sea eelfatty tuna
0.086
0.271
0.252
0.070
0.081
0.772
0.084
0.039
sea eeltuna roll
0.176
0.334
0.287
0.165
0.154
0.653
0.182
0.048
sea eelcucumber roll
0.071
0.257
0.249
0.085
0.099
0.519
0.073
0.084
tuna
squid
0.141
0.290
0.279
0.188
0.168
0.723
0.145
0.041
tuna
sea urchin
0.203
0.310
0.283
0.154
0.157
0.767
0.201
0.070
tuna
salmon roe
0.157
0.295
0.260
0.115
0.109
0.778
0.154
0.030
tuna
egg
0.130
0.280
0.280
0.163
0.164
0.693
0.133
0.069
tuna
fatty tuna
0.104
0.224
0.251
0.120
0.133
0.745
0.106
0.121
tuna
tuna roll
0.102
0.305
0.265
0.069
0.090
0.668
0.097
0.164
tuna
cucumber roll
0.088
0.242
0.273
0.099
0.108
0.632
0.090
0.110
squid
sea urchin
0.181
0.324
0.276
0.177
0.160
0.650
0.187
0.023
squid
salmon roe
0.205
0.309
0.278
0.188
0.149
0.699
0.210
0.052
squid
egg
0.044
0.300
0.266
0.037
0.056
0.544
0.040
0.064
squid
fatty tuna
0.175
0.260
0.256
0.193
0.206
0.786
0.177
0.094
squid
tuna roll
0.119
0.333
0.272
0.135
0.118
0.598
0.114
0.044
squid
cucumber roll
0.005
0.258
0.257
0.019
0.032
0.446
0.007
0.033
sea urchinsalmon roe
0.194
0.323
0.273
0.191
0.175
0.660
0.197
0.354
sea urchin
egg
0.275
0.308
0.255
0.311
0.218
0.654
0.279
0.154
sea urchinfatty tuna
0.057
0.268
0.254
0.078
0.081
0.733
0.059
0.115
sea urchintuna roll
0.273
0.300
0.254
0.269
0.217
0.661
0.278
0.121
sea urchin
cucumber roll
0.140
0.277
0.270
0.154
0.155
0.537
0.142
0.123
salmon roe
egg
0.170
0.311
0.283
0.212
0.153
0.669
0.174
0.063
salmon roe
fatty tuna
0.007
0.275
0.264
0.024
0.042
0.774
0.009
0.054
salmon roe
tuna roll
0.227
0.293
0.256
0.211
0.165
0.677
0.233
0.094
salmon roe
cucumber roll
0.102
0.257
0.265
0.115
0.127
0.583
0.104
0.104
egg
fatty tuna
0.169
0.266
0.283
0.183
0.193
0.773
0.171
0.110
egg
tuna roll
0.067
0.283
0.243
0.046
0.059
0.530
0.063
0.033
egg
cucumber roll
0.084
0.257
0.274
0.094
0.102
0.377
0.086
0.054
fatty tuna
tuna roll
0.027
0.293
0.261
0.045
0.063
0.736
0.029
0.035
fatty tuna
cucumber roll
0.142
0.230
0.296
0.147
0.151
0.736
0.143
0.135
tuna rollcucumber roll
0.015
0.256
0.251
0.028
0.041
0.428
0.017
0.027
Best
15
0
0
4
4
0
4
30

Table 5: Rate of item pairs that each method achieved the best with the Movielens data, where the number of sampled users
with records was U ∗ = 20 and number of items was I = 200.
IND
ML
EB
CGM
EB+CGM
EB+A
CGM+A
Proposed
0.121
0.031
0.014
0.170
0.142
0.068
0.211
0.242

Table 6: Absolute error of co-occurrence probabilities averaged over 30 experiments with the Movielens-1M data.
IND
ML
EB
CGM
EB+CGM
EB+A
CGM+A
Proposed
0.01690.42190.01670.39520.01640.01540.35970.0147
Dirichlet allocation. Journal of Machine Learning Research3(Jan):993–1022.
Caruana, R. 1997. Multitask learning. Machine Learning
28(1):41–75.
Causey, B. D. 1983. Estimation of proportions for
multinomial contingency tables subject to marginal constraints. Communications in Statistics-Theory and Methods
12(22):2581–2587.
Deming, W. E., and Stephan, F. F. 1940. On a least squares
adjustment of a sampled frequency table when the expectedmarginal totals are known. The Annals of Mathematical
Statistics 11(4):427–444.
Dobra, A.; Tebaldi, C.; and West, M. 2006. Dataaugmentation in multi-way contingency tables with fixedmarginal totals. Journal of Statistical Planning and Inference 136(2):355–372.
Fienberg, S. E. 1970. An iterative procedure for estimationin contingency tables. The Annals of Mathematical Statistics
41(3):907–917.
Freedman, D. A.; Klein, S. P.; Sacks, J.; Smyth, C. A.; andEverett, C. G. 1991. Ecological regression and voting rights.
Evaluation Review 15(6):673–711.
Girvan, M., and Newman, M. E. 2002. Community structure in social and biological networks. Proceedings of the
National Academy of Sciences 99(12):7821–7826.
Goodman, L. A. 1953. Ecological regressions and behavior
of individuals. American Sociological Review.
Greco, F. 2016. Estimation of multi-way tables subject tocoherence constraints. Statistica 76(2):115–125.
Ireland, C. T., and Kullback, S. 1968. Contingency tables
with given marginals. Biometrika 55(1):179–188.
Iwata, T., and Shimizu, H. 2019. Neural collective graphicalmodels for estimating spatio-temporal population flow fromaggregated data. In AAAI Conference on Artificial Intelligence.
Kamishima, T. 2003. Nantonac collaborative filtering: recommendation based on order responses. In Proceedingsof the ninth ACM SIGKDD International Conference on
Knowledge Discovery and Data Mining, 583–588. ACM.
Kingma, D. P., and Ba, J. 2015. ADAM: A method forstochastic optimization. In International Conference on
Learning Representations.
Konstan, J. A.; Miller, B. N.; Maltz, D.; Herlocker, J. L.;
Gordon, L. R.; and Riedl, J. 1997. Grouplens: applying
collaborative filtering to usenet news. Communications of
the ACM 40(3):77–87.
Kumar, A.; Sheldon, D.; and Srivastava, B. 2013. Collective
diffusion over networks: Models and inference. In Proceedings of International Conference on Uncertainty in Artificial
Intelligence.

Little, R. J., and Wu, M.-M. 1991. Models for contingency
tables with known margins when target and sampled populations differ. Journal of the American Statistical Association
86(413):87–95.
Nguyen, T.; Kumar, A.; Lau, H. C.; and Sheldon, D. 2016.
Approximate inference using DC programming for collective graphical models. In Artificial Intelligence and Statistics, 685–693.
Park, Y., and Ghosh, J. 2012. A probabilistic imputation
framework for predictive analysis using variably aggregated,
multi-source healthcare data. In Proceedings of the 2nd
ACM SIGHIT International Health Informatics Symposium,
445–454. ACM.
Park, Y., and Ghosh, J. 2014. Ludia: An aggregateconstrained low-rank reconstruction algorithm to leverage
publicly released health data. In Proceedings of the 20th
ACM SIGKDD International Conference on Knowledge
Discovery and Data Mining, 55–64. ACM.
Paszke, A.; Gross, S.; Chintala, S.; Chanan, G.; Yang, E.;
DeVito, Z.; Lin, Z.; Desmaison, A.; Antiga, L.; and Lerer,
A. 2017. Automatic differentiation in PyTorch. In NIPS
Autodiff Workshop.
Sarwar, B. M.; Karypis, G.; Konstan, J. A.; and Riedl, J.
2001. Item-based collaborative filtering recommendation algorithms. WWW 1:285–295.
Sheldon, D. R., and Dietterich, T. G. 2011. Collective graphical models. In Advances in Neural Information Processing
Systems, 1161–1169.
Sheldon, D.; Sun, T.; Kumar, A.; and Dietterich, T. G. 2013.
Approximate inference in collective graphical models. In
Proceedings of the 30th International Conference on Machine Learning.
Slavkovic, A. B. 2010. Partial information releases for confidential contingency table entries: Present and future research
efforts. Journal of Privacy and Confidentiality 1(2).
Smith, J. H. 1947. Estimation of linear functions of cell proportions. The Annals of Mathematical Statistics 18(2):231–
254.
Stephan, F. F. 1942. An iterative method of adjusting sample
frequency tables when expected marginal totals are known.
The Annals of Mathematical Statistics 13(2):166–178.
Sun, T.; Sheldon, D.; and Kumar, A. 2015. Message passing
for collective graphical models. In Proceedings of the 32nd
International Conference on Machine Learning, 853–861.
Zaheer, M.; Kottur, S.; Ravanbakhsh, S.; Poczos, B.;
Salakhutdinov, R. R.; and Smola, A. J. 2017. Deep sets. In
Advances in Neural Information Processing Systems, 3391–
3401.