Conditional Density Estimation
via Least-Squares Density Ratio Estimation

Masashi Sugiyama
Ichiro Takeuchi
Tokyo Institute of Technology & JST Nagoya Institute of Technology
Takafumi Kanamori
Nagoya University

Hirotaka Hachiya
Tokyo Institute of Technology

Abstract

Daisuke Okanohara
The University of Tokyo

the conditional distribution itself. In this paper, weaddress the problem of estimating conditional densities
when x and y are continuous and multi-dimensional.

Estimating the conditional mean of an inputoutput relation is the goal of regression.
However, regression analysis is not suﬃciently informative if the conditional distribution has multi-modality, is highly asymmetric, or contains heteroscedastic noise. In
such scenarios, estimating the conditional
distribution itself would be more useful. In
this paper, we propose a novel method of conditional density estimation. Our basic idea is
to express the conditional density in termsof the ratio of unconditional densities, andthe ratio is directly estimated without going
through density estimation. Experiments using benchmark and robot transition datasets
illustrate the usefulness of the proposed approach.

1

Taiji Suzuki
The University of Tokyo

The mixture density network (MDN) (Bishop, 2006)
models the conditional density by a mixture of parametric densities, where the parameters are estimatedby a neural network. MDN was shown to work well,although its training is time-consuming and only a local optimal solution may be obtained due to the nonconvexity of neural network learning. Similarly, a mixture of Gaussian processes was explored for estimating
the conditional density (Tresp, 2001). The mixture
model is trained in a computationally eﬃcient mannerby an expectation-maximization algorithm (Dempsteret al., 1977). However, since the optimization problem
is non-convex, one may only access to a local optimal
solution in practice.

Introduction

Regression is aimed at estimating the conditional meanof output y given input x. When the conditional
density p(y|x) is unimodal and symmetric, regression
would be suﬃcient for analyzing the input-output dependency. However, estimating the conditional mean
may not be suﬃciently informative, when the conditional distribution possesses multi-modality (e.g., inverse kinematics learning of a robot, see Bishop, 2006)
or a highly skewed proﬁle with heteroscedastic noise(e.g., biomedical data analysis, see Hastie et al., 2001).
In such cases, it would be more informative to estimate
Appearing in Proceedings of the 13th International Conference on Artiﬁcial Intelligence and Statistics (AISTATS)
2010, Chia Laguna Resort, Sardinia, Italy. Volume 9 of
JMLR: W&CP 9. Copyright 2010 by the authors.

781

The
kernelquantile
regression(KQR)
method (Takeuchi et al., 2006; Li et al., 2007)
allows one to predict percentiles of the conditional
distribution. This implies that solving KQR for allpercentiles gives an estimate of the entire conditional
cumulative distribution. KQR is formulated as aconvex optimization problem, and therefore a uniqueglobal solution can be obtained. Furthermore, theentire solution path with respect to the percentile
parameter, which was shown to be piece-wise linear,can be computed eﬃciently (Takeuchi et al., 2009).
However, the range of applications of KQR is limited
to one-dimensional output and solution path tracking
tends to be numerically rather unstable in practice.
In this paper, we propose a new method of conditional density estimation named least-squares conditional density estimation (LS-CDE), which can be applied to multi-dimensional inputs and outputs. The
proposed method is based on the fact that the conditional density can be expressed in terms of unconditional densities as p(y|x) = p(x, y)/p(x). Our key

Conditional Density Estimation via Least-Squares Density Ratio Estimation

idea is that we do not estimate the two densities
p(x, y) and p(x) separately, but we directly estimate
the density ratio p(x, y)/p(x) without going throughdensity estimation. Experiments using benchmark androbot transition datasets show that our method compares favorably with existing methods in terms of theaccuracy and computational eﬃciency.

2

A New Method of Conditional
Density Estimation

In this section, we formulate the problem of conditional density estimation and give a new method.

are basis functions such that ϕ(x, y) ≥ 0b for all(x, y) ∈ DX ×DY . 0b denotes the b-dimensional vector
with all zeros. The inequality for vectors is applied inan element-wise manner.
Note that the number b of basis functions is not necessarily a constant; it can depend on the number n ofsamples. Similarly, the basis functions ϕ(x, y) could
be dependent on the samples {xi , yi }ni=1 . This meansthat kernel models (i.e., b = n and ϕi (x, y) is a kernel
function ‘centered’ at (xi , yi )) are also included in theabove formulation. We explain how the basis functions
ϕ(x, y) are practically chosen in Section 2.5.
2.3

2.1

Conditional Density Estimation via
Density Ratio Estimation

Let DX (⊂ RdX ) and DY (⊂ RdY ) be input and output
data domains, where dX and dY are the dimensionality of the data domains, respectively. Let us consider ajoint probability distribution on DX × DY with probability density function p(x, y), and suppose that we are
given n independent and identically distributed (i.i.d.)paired samples of input x and output y:{zi | zi = (xi , yi ) ∈ DX × DY }ni=1 .
The goal is to estimate the conditional density p(y|x)
from the samples {zi }ni=1 . Our primal interest is
the case where both variables x and y are multidimensional and continuous.
A key idea of our proposed approach is to consider theratio of two unconditional densities:p(y|x) =p(x, y):= r(x, y),p(x)where we assume p(x) > 0 for all x ∈ DX . However,naively estimating the two unconditional densities andtaking their ratio can result in large estimation error. In order to avoid this, we propose to estimate thedensity ratio function r(x, y) directly without going
through density estimation of p(x, y) and p(x).
2.2

A Least-squares Approach to Conditional
Density Estimation

We determine the parameter α in the model rbα (x, y)so that the following squared error J0 is minimized:∫∫
12J0 (α) :=(brα (x, y) − r(x, y)) p(x)dxdy.
2This can be expressed as∫∫
1J0 (α) =rbα (x, y)2 p(x)dxdy
2∫∫−rbα (x, y)r(x, y)p(x)dxdy + C
∫∫
( ⊤
)2
1
α ϕ(x, y) p(x)dxdy=
2∫∫−α⊤ ϕ(x, y)p(x, y)dxdy + C, (2)∫∫where C = 21r(x, y)p(x, y)dxdy is a constant andtherefore can be safely ignored. Let us denote the ﬁrst
two terms of Eq.(2) by J:

whereJ(α) := J0 (α) − C =1 ⊤α Hα − h⊤ α,
2

∫

∫∫

Φ(x)p(x)dx, h :=ϕ(x, y)p(x, y)dxdy,∫
Φ(x) := ϕ(x, y)ϕ(x, y)⊤ dy.
(3)
H :=

Linear Density-ratio Model

We model the density ratio function r(x, y) by thefollowing linear model:rbα (x, y) := α⊤ ϕ(x, y),(1)

The matrix H and the vector h included in J(α) contain the expectations over unknown densities p(x) andp(x, y), so we approximate the expectations by sample
averages. Then we have1b ⊤ α,c −hbJ(α):= α⊤ Hα2where ⊤ denotes the transpose of a matrix or a vector,
α = (α1 , α2 , . . . , αb )⊤

where
n
n
∑
∑
b := 1c := 1HΦ(xi ), hϕ(xi , yi ).
n i=1
n i=1

are parameters to be learned from samples, andϕ(x, y) = (ϕ1 (x, y), ϕ2 (x, y), . . . , ϕb (x, y))⊤

782

(4)

Sugiyama, Takeuchi, Suzuki, Kanamori, Hachiya, and Okanohara

Note that the integral over y included in Φ(x) (seeEq.(3)) can be computed in principle since it does
not contain any unknown quantity. As shown in Section 2.5, this integration can be computed analyticallyin our basis function choice.

Let G be a general set of functions on DX × DY . Notethat G corresponds to the span of our model, which
could be non-parametric (i.e., an inﬁnite dimensional
linear space). For a function g (∈ G), let us consider anon-negative function R(g) such that

Now our optimization criterion is summarized as[]b
e := argmin J(α)α+ λα⊤ α ,
(5)

{
[∫
]
}
max supg(x, y)dy , sup [g(x, y)] ≤ R(g).
x

x,y

α∈Rb

where a regularizer λα⊤ α (λ > 0) is included for stabilization purposes. Taking the derivative of the above
objective function and equating it to zero, we can see
e can be obtained just by solving thethat the solution αbc + λIb )α = h,following system of linear equations. (Hwhere Ib denotes the b-dimensional identity matrix.
e is given analytically as
Thus, the solution αbc + λIb )−1 h.
e = (Hα(6)Since the density ratio function is non-negative by defe as1
inition, we modify the solution αb := max(0b , α),
e
α

(7)

where the ‘max’ operation for vectors is applied inan element-wise manner. Thanks to this rounding-up
b tends to be sparse, which
processing, the solution αcontributes to reducing the computation time in thetest phase.
In order to assure that the obtained density-ratio function is a conditional density, we renormalize the solue, ourtion in the test phase—given a test input point x
ﬁnal solution is given ase) = ∫
pb(y|x = xb ⊤ ϕ(eαx, y).
b ⊤ ϕ(eαx, y)dy(8)

We call the above method Least-Squares Conditional
Density Estimation (LS-CDE). LS-CDE can be regarded as an application of the direct density-ratio
estimation method called the unconstrained LeastSquares Importance Fitting (uLSIF) (Kanamori et al.,2009) to the problem of density ratio estimation.
2.4

Convergence Analysis

Here, we brieﬂy show a non-parametric convergence
rate of the LS-CDE solution.

Then the problem (5) can be generalized as[

1 ∑rb := argmin
2n i=1g∈Gn∫g(xi , y)2 dy

]
n
1∑
2
−
g(xi , yi ) + λn R(g) ,n i=1

where λn is the regularization parameter depending
on n. We assume that the true density ratio functionr(x, y) is contained in G and there exists M (> 0) suchthat R(r) < M . We also assume that there exists γ(0 < γ < 2) such thatγH[] (GM , ϵ, L2 (px × µY )) = O ((M/ϵ) ) ,where GM := {g ∈ G | R(g) ≤ M }. µY is the Lebesgue
measure on DY , px × µY is a product measure of px
and µY , and H[] is the bracketing entropy of GM withrespect to the L2 (px × µY )-norm (van der Vaart &Wellner, 1996).
Under the above assumptions, we have the following
theorem (its proof is omitted since it follows essentiallythe same line as Sugiyama et al., 2008).
Theorem 1 Under the above setting, if λn → 0 and2/(2+γ)λ−1), thenn = o(n∥br − r∥2 = Op (λ1/2n ),where ∥ · ∥2 denotes the L2 (px × µY )-norm and Op
denotes the asymptotic order in probability.
Note that the conditions λn → 0 and λ−1=no(n2/(2+γ) ) intuitively means that λn should converge
to zero as n tends to inﬁnity, but the speed of convergence should not be too fast.

1

A variant of the proposed method would be to include
the positivity constraint α ≥ 0n directly in Eq.(6). Our
preliminary experiments showed that the estimation accuracy of this modiﬁed algorithm turned out to be comparable to Eq.(7), while the constrained version was computationally less eﬃcient than Eq.(7) since we need to use anumerical quadratic program solver for computing the solution. For this reason, we only consider Eq.(7) in the restof this paper.

783

2.5

Basis Function Design

It is straightforward to show that cross-validation isavailable for model selection. A good model may
be chosen by cross-validation, given that a family ofpromising model candidates is prepared. As modelcandidates, we propose to use a Gaussian kernel model:

Conditional Density Estimation via Least-Squares Density Ratio Estimation

for z = (x⊤ , y ⊤ )⊤ ,()∥z − wℓ ∥2
ϕℓ (x, y) = exp −
2σ 2
)
(
)
(
∥x − uℓ ∥2∥y − vℓ ∥2= exp −
exp −
,
2σ 2
2σ 2
(9)
⊤ ⊤ b
where {wℓ | wℓ = (u⊤ℓ , vℓ ) }ℓ=1 are center points⊤ ⊤ nrandomly chosen from {zi | zi = (x⊤i , yi ) }i=1 . We
may use diﬀerent Gaussian widths for x and y. However, for simplicity, we decided to use the common
Gaussian width σ for both x and y under the setting
where the variance of each element of x and y is normalized to one.

An advantage of the above Gaussian kernel model isthat the integrals over y in matrix Φ (see Eq.(3)) andin the normalization factor (see Eq.(8)) can be computed analytically; indeed, a simple calculation yields∫
Φℓ,ℓ′ (x) = ϕℓ (x, y)ϕℓ′ (x, y)dy()√
ξℓ,ℓ′ (x)dY
= ( πσ) exp −
,
4σ 2
(
)
∫
b
∑
√
∥x − uℓ ∥2
⊤
dY
b ϕ(e
α
x, y)dy = ( 2πσ)αbℓ exp −
,
2σ 2
ℓ=1

where ξℓ,ℓ′ (x) := 2∥x−uℓ ∥2 +2∥x−uℓ′ ∥2 +∥vℓ −vℓ′ ∥2 .
In the experiments, we ﬁx the number of basis functions to b = min(100, n), and choose the Gaussian
width σ and the regularization parameter λ by CVfrom σ, λ ∈ {0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10}.

3Discussion

In this section, we discuss the characteristics of existing and proposed methods of conditional density estimation.
3.1

ϵ-neighbor Kernel Density Estimation(ϵ-KDE)For estimating the conditional density p(y|x), ϵneighbor kernel density estimation (ϵ-KDE) employs
the standard kernel density estimator using a subset ofsamples, {yi }i∈Ix,ϵ for some threshold ϵ (≥ 0), where
Ix,ϵ is the set of sample indices such that ∥xi −x∥ ≤ ϵ.
In the case of Gaussian kernels, ϵ-KDE is expressed as∑1
pb(y|x) =N (y; yi , σ 2 IdY ),|Ix,ϵ |i∈Ix,ϵwhere N (y; µ, Σ) denotes the Gaussian density withmean µ and covariance matrix Σ. The threshold ϵ784

and the bandwidth σ may be chosen based on CV.
ϵ-KDE is simple and easy to use, but it may not be
reliable in high-dimensional problems. Slightly moresophisticated variants have been proposed based onweighted kernel density estimation (Fan et al., 1996;
Wolﬀ et al., 1999), but they may still share the same
weakness.
3.2

Mixture Density Network (MDN)

The mixture density network (MDN) models the conditional density by a mixture of parametric densities (Bishop, 2006). In the case of Gaussian densities,
MDN is expressed aspb(y|x) =t∑πℓ (x)N (y; µℓ (x), σℓ2 (x)IdY ),

ℓ=1

where
πℓ (x) denotes the mixing coeﬃcient such that∑tπℓ=1 ℓ (x) = 1 and 0 ≤ πℓ (x) ≤ 1 for all x ∈ DX . All
the parameters {πℓ (x), µℓ (x), σℓ2 (x)}tℓ=1 are learned
as a function of x by a neural network with regularized maximum likelihood estimation. The number t ofGaussian components, the number of hidden units inthe neural network, and the regularization parameter
may be chosen based on CV. MDN has been shown
to work well, although its training is time-consuming
and only a local solution may be obtained due to thenon-convexity of neural network learning.
3.3

Kernel Quantile Regression (KQR)

Kernel quantile regression (KQR) allows one to predictthe 100τ -percentile of conditional distributions for agiven τ (∈ (0, 1)) when y is one-dimensional (Takeuchi
et al., 2006; Li et al., 2007). For the Gaussian kernel
model
n
∑
fbτ (x) =αi,τ ϕi (x) + bτ ,
i=1

(
)
2
i∥
where ϕi (x) = exp − ∥x−x
, the parameters22σ{αi,τ }ni=1 and bτ are learned by


n
n
∑
∑
 ψτ (yi − fbτ (xi ))+λmin
ϕi (xj )αi,τ αj,τ ,n{αi,τ }i=1 ,bτ

i=1

i,j=1

where ψτ (r) denotes the pin-ball loss function deﬁnedby{(1 − τ )|r| (r ≤ 0),ψτ (r) =τ |r|(r > 0).
Thus, solving KQR for all τ ∈ (0, 1) gives an estimateof the entire conditional distribution. The bandwidth
σ and the regularization parameter λ may be chosen
based on CV.

Sugiyama, Takeuchi, Suzuki, Kanamori, Hachiya, and Okanohara

A notable advantage of KQR is that the solution of
KQR is piece-wise linear with respect to τ , so the entire
solution path can be computed eﬃciently (Takeuchi
et al., 2009). This implies that the conditional cumulative distribution can be computed eﬃciently. However, solution path tracking tends to be numericallyrather unstable and the range of applications of KQR
is limited to one-dimensional output y. Furthermore,some heuristic procedure is needed to convert conditional cumulative distributions into conditional densities, which can cause additional estimation errors.

3.4

4.1

Illustrative Examples

Here we illustrate how the proposed LS-CDE methodbehaves using toy datasets.
Let dX = dY = 1. Inputs {xi }ni=1 were independently
drawn from U (−1, 1), where U (a, b) denotes the uniform distribution on (a, b). Outputs {yi }ni=1 were generated by the following heteroscedastic noise model:
yi = sinc(2πxi ) +1
exp(1 − xi ) · εi .
8

We tested the following three diﬀerent distributionsfor {εi }ni=1 :Other Methods of Density Ratio
Estimation

i.i.d.

(a) Gaussian: εi ∼ N (0, 1).
A naive method for estimating the density ratiop(x, y)/p(x) is to ﬁrst approximate the two densities
p(x, y) and p(x) by standard kernel density estimationand then taking the ratio of the estimated densities.
We refer to this method as the ratio of kernel density
estimators (RKDE). As we will show through experiments in the next section, RKDE does not work wellsince taking the ratio of estimated quantities signiﬁcantly magniﬁes the estimation error.
To overcome the above weakness, we decided to directly estimate the density ratio without going throughdensity estimation under the squared-loss. The kernelmean matching method (Huang et al., 2007) and thelogistic regression based method (Qin, 1998; Cheng &Chu, 2004; Bickel et al., 2007) also allow one to directly estimate a density ratio q(x)/q ′ (x). However,the derivation of these methods heavily relies on thefact that the two density functions q(x) and q ′ (x)share the same domain, which is not fulﬁlled in thecurrent setting. For this reason, these methods may
not be employed for conditional density estimation.
Other methods of direct density ratio estimation (Sugiyama et al., 2008; Nguyen et al., 2008) employs the Kullback-Leibler (KL) divergence as the loss
function, instead of the squared-loss. It is possible touse these methods for conditional density estimationin the same way as the proposed method. Indeed,
our preliminary experiments showed that a KL-based
method was comparable to the squared-loss method
in terms of accuracy. However, the KL-based method
was computationally less eﬃcient. For this reason, we
decided to focus on the squared-loss method.

4

i.i.d. 1414
2 N (−1, 9 ) + 2 N (1, 9 ).

(b) Bimodal: εi ∼

i.i.d. 31
3 1
4 N (0, 1) + 4 N ( 2 , 9 ).

(c) Skewed: εi ∼
i.i.d.

‘ ∼ ’ denotes ‘independent and identically distributed’and N (µ, σ 2 ) denotes the Gaussian distribution withmean µ and variance σ 2 . The number of training samples was set to n = 200. The numerical results were depicted in Figure 1, illustrating that LS-CDE well captures heteroscedasticity, bimodality, and asymmetricity.
We have also investigated the experimental performance of LS-CDE using the following real datasets:(d) Bone Mineral Density dataset: Relative
spinal bone mineral density measurements on 485
North American adolescents (Hastie et al., 2001),
having a heteroscedastic asymmetric conditional
distribution.
(e) Old Faithful Geyser dataset: The durationsof 299 eruptions of the Old Faithful Geyser (Weisberg, 1985), having a bimodal conditional distribution.
Figure 2 depicts the experimental results, showing thatheteroscedastic and multi-modal structures were nicely
revealed by LS-CDE.
4.2

Benchmark Datasets

We applied the proposed and existing methods to thebenchmark datasets accompanied with the R package(see Table 1) and evaluate their experimental performance.

Numerical Experiments

In this section, we investigate the experimental performance of the proposed and existing methods.

785

In each dataset, 50% of samples were randomly chosen for conditional density estimation and the rest wasConditional Density Estimation via Least-Squares Density Ratio Estimation

3

0.25
Relative Change in Spinal BMD

Truth
Estiamted

2

Output y

1
0
-1
-2

0.2
0.15
0.1
0.05
0
-0.05

-3

-0.1
-1

-0.5

0
Input x0.518(a) Heteroscedastic Gaussian
3

12

14

16

18 20
Age

22

24

26

28

(a) Bone Mineral Density
6

Truth
Estiamted
Duration time [minutes]

21
Output y

10

0
-1
-2

5
4
3
2
1

-3

0
-1

-0.5

0
Input x

0.5

1

40

(b) Bimodal
3

50

60
70
80
Waiting time [minutes]

90

100

(b) Old Faithful Geyser

Figure 2: Illustrative examples of LS-CDE for real
datasets.

Truth
Estiamted

2

Output y

1
0
-1
-2
-3
-1

-0.5

0
Input x

0.51(c) Skewed

Figure 1: Illustrative examples of LS-CDE for artiﬁcial
datasets.
used for computing the estimation accuracy. The accuracy of a conditional density estimator pb(y|x) was
measured by the negative log-likelihood for test sameples {ezi | zei = (exi , yei )}ni=1:ne

NLL := −1∑
log pb(e
yi |e
xi ).
n
e i=1

(10)

Thus, the smaller the value of NLL is, the betterthe performance of the conditional density estimatorpb(y|x) is.
We compared LS-CDE, ϵ-KDE, MDN, KQR, and

786

RKDE. For model selection, we used CV based on the
log-likelihood. In MDN, CV over three tuning parameters (the number of Gaussian components, the numberof hidden units in the neural network, and the regularization parameter) was unbearably slow, so the number of Gaussian components was ﬁxed to t = 3 andthe other two tuning parameters were chosen by CV.
The experimental results are summarized in Table 1.
ϵ-KDE was computationally very eﬃcient, but it
tended to perform rather poorly. MDN worked well,
but it is computationally highly demanding. KQR
overall performed well and it was computationallyslightly more eﬃcient than LS-CDE. However, its solution path tracking algorithm was numerically rather
unstable and we could not obtain solutions for the ‘engel’ and ‘cpus’ datasets. RKDE did not perform well
for all cases, implying that density ratio estimation via
density estimation is not reliable in practice. Overall,
the proposed LS-CDE was shown to be a promising
method for conditional density estimation in terms of
the accuracy and computational eﬃciency.

Sugiyama, Takeuchi, Suzuki, Kanamori, Hachiya, and Okanohara

Table 1: Experimental results on benchmark datasets (dY = 1). The average and the standard deviation of NLL(see Eq.(10)) over 10 runs are described (smaller is better). The best method in terms of the mean error andcomparable methods according to the two-sided paired t-test at the signiﬁcance level 5% are speciﬁed by boldface. Mean computation time is normalized so that LS-CDE is one.
Datasetcautionftcollinssnow
highwayheights
sniﬀersnowgeeseufcbirthwt
crabs
GAGurine
geysergilgais
topo
BostonHousing
CobarOreengel
mcycleBigMac2003
UN3
cpus
Time

4.3

(n, dX )(50,2)(46,1)(19,11)(687,1)(62,4)(22,2)(117,4)(94,7)(100,6)(157,1)(149,1)(182,8)(26,2)(253,13)(19,2)(117,1)(66,1)(34,9)(62,6)(104,7)

LS-CDE
1.24 ± 0.29
1.48 ± 0.01
1.71 ± 0.41
1.29 ± 0.00
0.69 ± 0.16
0.95 ± 0.10
1.03 ± 0.01
1.43 ± 0.01
-0.07 ± 0.11
0.45 ± 0.04
1.03 ± 0.00
0.73 ± 0.05
0.93 ± 0.02
0.82 ± 0.05
1.58 ± 0.06
0.69 ± 0.04
0.83 ± 0.03
1.32 ± 0.11
1.42 ± 0.12
1.04 ± 0.07
1ϵ-KDE
1.25 ± 0.19
1.53 ± 0.05
2.24 ± 0.64
1.33 ± 0.01
0.96 ± 0.15
1.35 ± 0.17
1.40 ± 0.02
1.48 ± 0.01
0.99 ± 0.09
0.92 ± 0.05
1.11 ± 0.02
1.35 ± 0.03
1.18 ± 0.09
1.03 ± 0.05
1.65 ± 0.09
1.27 ± 0.05
1.25 ± 0.23
1.29 ± 0.14
1.78 ± 0.14
1.01 ± 0.10
0.004

MDN
1.39 ± 0.18
1.48 ± 0.03
7.41 ± 1.22
1.30 ± 0.01
0.72 ± 0.09
2.49 ± 1.02
1.02 ± 0.06
1.46 ± 0.01
-0.70 ± 0.35
0.57 ± 0.15
1.23 ± 0.05
0.10 ± 0.04
2.11 ± 0.46
0.68 ± 0.06
1.63 ± 0.08
0.71 ± 0.16
1.12 ± 0.10
2.64 ± 0.84
1.32 ± 0.08
-2.14 ± 0.13
267

KQR
1.73 ± 0.86
2.11 ± 0.44
5.69 ± 1.69
1.29 ± 0.00
0.68 ± 0.21
2.96 ± 1.13
1.02 ± 0.06
1.58 ± 0.05
-1.03 ± 0.16
0.40 ± 0.08
1.10 ± 0.02
0.45 ± 0.15
2.88 ± 0.85
0.48 ± 0.10
6.33 ± 1.77
N.A.
0.72 ± 0.06
1.35 ± 0.26
1.22 ± 0.13
N.A.
0.755

RKDE
17.11 ± 0.25
46.06 ± 0.78
15.30 ± 0.76
54.79 ± 0.10
26.80 ± 0.58
28.43 ± 1.02
11.10 ± 0.49
15.95 ± 0.53
12.60 ± 0.45
53.43 ± 0.27
53.49 ± 0.38
10.44 ± 0.50
10.80 ± 0.35
17.81 ± 0.25
11.42 ± 0.51
52.83 ± 0.16
48.35 ± 0.79
13.34 ± 0.52
11.43 ± 0.58
15.16 ± 0.72
0.089

Robot Transition Estimation

We further applied the proposed and existing methods
to the problem of robot transition estimation. We used
the pendulum robot and the Khepera robot simulators
illustrated in Figure 3.
The pendulum robot consists of wheels and a pendulum hinged to the body. The state of the pendulum
robot consists of angle θ and angular velocity θ̇ of thependulum. The amount of torque τ applied to thewheels can be controlled, by which the robot can moveleft or right and the state of the pendulum is changed
to θ′ and θ̇′ . The task is to estimate p(θ′ , θ̇′ |θ, θ̇, τ ),the transition probability density from state (θ, θ̇) tostate (θ′ , θ̇′ ) by action τ .
The Khepera robot is equipped with two infra-red sensors and two wheels. The infra-red sensors dL and dR
measure the distance to the left-front and right-front
walls. The speed of left and right wheels vL and vR canbe controlled separately, by which the robot can move
forward/backward and rotate left/right. The task isto estimate p(d′L , d′R |dL , dR , vL , vR ), where d′L and d′R
are the next state.
The state transition of the pendulum robot is highly
stochastic due to slip, friction, or measurement errors with strong heteroscedasticity. Sensory inputs of
the Khepera robot suﬀer from occlusions and containhighly heteroscedastic noise, so the transition probability density may possess multi-modality and het-

787(a) Pendulum robot.

(b) Khepera robot.

Figure 3: Illustration of robots used for experiments.

eroscedasticity. Thus transition estimation of dynamic
robots is a challenging task. Note that transition estimation is highly useful in model-based reinforcement
learning.
For both robots, 100 samples were used for conditional
density estimation and additional 900 samples wereused for computing NLL (see Eq.(10)). The numberof Gaussian components was ﬁxed to t = 3 in MDN,
and all other tuning parameters were chosen by CV
based on the log-likelihood. Experimental results are
summarized in Table 2, showing that LS-CDE is still
useful in this challenging task of robot transition estimation.

Conditional Density Estimation via Least-Squares Density Ratio Estimation

Table 2: Experimental results on robot transition estimation. The average and the standard deviation of NLL(see Eq.(10)) over 10 runs are described (smaller is better). The best method in terms of the mean error andcomparable methods according to the two-sided paired t-test at the signiﬁcance level 5% are speciﬁed by boldface. Mean computation time is normalized so that LS-CDE is one. KQR was not included here since it isapplicable only when dY = 1.
Dataset(n, dX , dY )Pendulum1(100,3,2)
Pendulum2(100,3,2)
Khepera1(100,4,2)
Khepera2(100,4,2)Time

5LS-CDE
1.27 ± 0.05
1.38 ± 0.05
1.69 ± 0.01
1.86 ± 0.01
1ϵ-KDE
2.04 ± 0.10
2.07 ± 0.10
2.07 ± 0.02
2.10 ± 0.01
0.164

Conclusions

We proposed a novel approach to conditional density
estimation called LS-CDE. Our basic idea was to directly estimate the ratio of unconditional density functions without going through density estimation. LSCDE was shown to oﬀer a sparse solution in an analytic
form and therefore is computationally eﬃcient. A nonparametric convergence rate of the LS-CDE algorithm
was also provided. Experiments on benchmark androbot-transition datasets demonstrated the usefulnessof LS-CDE.

MDN
1.44 ± 0.67
1.43 ± 0.58
1.90 ± 0.36
1.92 ± 0.26
1134

RKDE
11.24 ± 0.32
11.24 ± 0.32
11.03 ± 0.03
11.09 ± 0.02
0.431

by unlabeled data. In B. Schölkopf, J. Platt and T. Hoﬀman (Eds.), Advances in neural information processing
systems 19, 601–608. Cambridge, MA: MIT Press.
Kanamori, T., Hido, S., & Sugiyama, M. (2009). A leastsquares approach to direct importance estimation. Journal of Machine Learning Research, 10, 1391–1445.
Li, Y., Liu, Y., & Zhu, J. (2007). Quantile regression in
reproducing kernel Hilbert spaces. Journal of the American Statistical Association, 102, 255–268.

Acknowledgments

Nguyen, X., Wainwright, M., & Jordan, M. (2008). Estimating divergence functionals and the likelihood ratio
by penalized convex risk minimization. In J. C. Platt,
D. Koller, Y. Singer and S. Roweis (Eds.), Advancesin neural information processing systems 20, 1089–1096.
Cambridge, MA: MIT Press.

MS was supported by AOARD, SCAT, and the JST
PRESTO program.

Qin, J. (1998). Inferences for case-control and semiparametric two-sample density ratio models. Biometrika, 85,
619–639.

References
Bickel, S., Brückner, M., & Scheﬀer, T. (2007). Discriminative learning for diﬀering training and test distributions. Proceedings of the 24th International Conference
on Machine Learning (pp. 81–88).
Bishop, C. M. (2006). Pattern recognition and machinelearning. New York, NY, USA: Springer.
Cheng, K. F., & Chu, C. K. (2004). Semiparametric density estimation under a two-sample density ratio model.
Bernoulli, 10, 583–604.
Dempster, A. P., Laird, N. M., & Rubin, D. B. (1977).
Maximum likelihood from incomplete data via the EM
algorithm. Journal of the Royal Statistical Society, series
B, 39, 1–38.
Fan, J., Yao, Q., & Tong, H. (1996). Estimation of conditional densities and sensitivity measures in nonlinear
dynamical systems. Biometrika, 83, 189–206.
Hastie, T., Tibshirani, R., & Friedman, J. (2001). The
elements of statistical learning: Data mining, inference,and prediction. New York: Springer.
Huang, J., Smola, A., Gretton, A., Borgwardt, K. M., &Schölkopf, B. (2007). Correcting sample selection bias

788

Sugiyama, M., Suzuki, T., Nakajima, S., Kashima, H., vonBünau, P., & Kawanabe, M. (2008). Direct importanceestimation for covariate shift adaptation. Annals of theInstitute of Statistical Mathematics, 60, 699–746.
Takeuchi, I., Le, Q. V., Sears, T. D., & Smola, A. J. (2006).
Nonparametric quantile estimation. Journal of Machine
Learning Research, 7, 1231–1264.
Takeuchi, I., Nomura, K., & Kanamori, T. (2009).
Nonparametric conditional density estimation using
piecewise-linear solution path of kernel quantile regression. Neural Computation, 21, 533–559.
Tresp, V. (2001). Mixtures of gaussian processes. Advancesin Neural Information Processing Systems 13 (pp. 654–
660). MIT Press.
van der Vaart, A. W., & Wellner, J. A. (1996). Weakconvergence and empirical processes with applications tostatistics. New York, NY, USA: Springer.
Weisberg, S. (1985). Applied linear regression. New York,NY, USA: John Wiley.
Wolﬀ, R. C. L., Yao, Q., & Hall, P. (1999). Methods for
estimating a conditional distribution function. Journalof the American Statistical Association, 94, 154–163.