##

Masanari Kimura, Howard Bondell
{m.kimura<EMAIL_ADDRESS>
School of Mathematics and Statistics
The University of Melbourne

## Density Ratio Estimation via Sampling along Generalized Geodesics on
Statistical Manifolds

Masanari Kimura, Howard Bondell
{m.kimura<EMAIL_ADDRESS>
School of Mathematics and Statistics
The University of Melbourne

###### Abstract

The density ratio of two probability distributions is one of the fundamental
tools in mathematical and computational statistics and machine learning, and
it has a variety of known applications. Therefore, density ratio estimation
from finite samples is a very important task, but it is known to be unstable
when the distributions are distant from each other. One approach to address
this problem is density ratio estimation using incremental mixtures of the two
distributions. We geometrically reinterpret existing methods for density ratio
estimation based on incremental mixtures. We show that these methods can be
regarded as iterating on the Riemannian manifold along a particular curve
between the two probability distributions. Making use of the geometry of the
manifold, we propose to consider incremental density ratio estimation along
generalized geodesics on this manifold. To achieve such a method requires
Monte Carlo sampling along geodesics via transformations of the two
distributions. We show how to implement an iterative algorithm to sample along
these geodesics and show how changing the distances along the geodesic affect
the variance and accuracy of the estimation of the density ratio. Our
experiments demonstrate that the proposed approach outperforms the existing
approaches using incremental mixtures that do not take the geometry of the
manifold into account.

###  1 Introduction

The density ratio of two probability distributions is one of the fundamental
tools in mathematical and computational statistics. For example, its
applications include estimation under the covariate shift assumption [50, 54],
outlier detection [4, 22], variable selection [5, 21, 45], change point
detection [24, 28, 40] and causal inference [43, 47]. For practical purposes,
since the true density ratio is not accessible, density ratio estimation from
finite size samples is required. Therefore, density ratio estimation (DRE) is
a very important task that has attracted significant interest [26, 57, 63].

A key challenge in DRE is the difficulty of the estimation when the source and
target distributions are far apart. One approach to address this problem is
DRE using incremental mixtures of two distributions. Rhodes et al. [48] have
shown the effectiveness of a method that creates $T>0$ bridging distributions
connecting the source and target distributions and then shrinks the gap
between the two distributions by utilizing a divide-and-conquer framework. The
method is based on the simple fact that a chain of density ratios using
arbitrary bridge distributions restores the original density ratio. They
employ a linear mixture of source and target distributions as the bridge
distributions. In addition, Choi et al. [12] extended this method to use
infinitely continuum bridge distributions. We refer to the estimator that
these methods produce as the Incremental Mixture Density Ratio Estimator
(IMDRE).

We investigate the family of IMDRE within the framework of information
geometry [2, 3], which allows for discussion on manifolds created by a set of
probability distributions. From the viewpoint of information geometry, it is
known that the set of probability distributions constitutes a Riemannian
manifold, which is a non-Euclidean space, and that the parameter space of the
probability distributions plays the role of a coordinate system. With this
geometric tool, we can reinterpret IMDRE with linear mixtures as a sequential
DRE along a specific curve called the $m$-geodesic. A geodesic here is defined
as a straight path on a manifold when a certain coordinate system is used.
This means that different coordinate systems are equipped with different
geodesics. We propose IMDRE with generalized geodesics called
$\alpha$-geodesics. We call this extension as Generalized IMDRE (GIMDRE). To
realize GIMDRE, sampling along the $\alpha$-geodesics connecting the source
and target distributions is required. However, since the classes of the two
probability distributions and their parameters are unknown, direct sampling is
a non-trivial problem. To address this issue, by utilizing the Monte Carlo
method with importance sampling [16, 61], we show that sampling according to
$\alpha$-geodesics can be achieved by giving the density ratio of the source
and target distributions. The next problem here is that the density ratio
estimation by GIMDRE relies on sampling along the $\alpha$-geodesic, while the
sampling depends on the density ratio. To resolve this interdependent
deadlock, we develop an optimization algorithm that alternately estimates the
density ratio and updates the importance weights for sampling. We design
numerical experiments to demonstrate the usefulness of the proposed method and
its behavior. To summarize, our contributions are as follows.

  * •

We geometrically reinterpret IMDRE through the lens of information geometry,
and we show that IMDRE with linear mixtures can be viewed as a sequential DRE
along a special curve called the $m$-geodesic. Furthermore, asymptotic
analysis shows that certain geodesics appear in the evaluation of the DRE
(Section 3, Theorems 1, 2 and Corollary 1).

  * •

We consider extending IMDRE with generalized geodesics called
$\alpha$-geodesics. We refer to this extension as GIMDRE, and we show that the
algorithm to obtain GIMDRE can be achieved by alternating optimization of
density ratio estimation and updating the importance weighting for sampling
based on the Monte Carlo procedure (Section 4 and Figure 1).

  * •

We design numerical experiments, including permutation test, an important task
in statistics, to investigate the behavior of our algorithm and demonstrate
its effectiveness (Section 5, Figures 2, 3, 4, and Tables 1, 2, 3).

Our study provides insights into the connections between statistical
procedures and differential geometry. We show that the use of
$\alpha$-geodesics for the estimation of density ratios improves upon both the
traditional IMDRE and the direct estimation. See the appendix for technical
details such as proofs of theorems and experimental setups, as well as
additional results including discussion on effective sample size, analysis of
variance, and additional numerical experiments.

###  2 Background and Preliminary

In this section, we first provide the background and preliminary knowledge
that is required throughout this study. See Appendix D for related literature.

####  2.1 Problem Setting

Let $\mathcal{X}\subset\mathbb{R}^{d}$ be the $d$-dimensional data domain with
$d\in\mathbb{N}$. Suppose that we are given independent and identically
distributed (i.i.d.) samples $\\{\bm{x}_{i}^{s}\\}^{n_{s}}_{i=1}$ and
$\\{\bm{x}_{i}^{t}\\}^{n_{t}}_{i=1}$ of sizes $n_{s},n_{t}\in\mathbb{N}$, from
$p_{s}(\bm{x})$ and $p_{t}(\bm{x})$, the source and target distributions,
respectively. The goal is to estimate the density ratio

$\displaystyle r(\bm{x})\coloneqq p_{s}(\bm{x})/p_{t}(\bm{x}),$ (1)

where we assume that $p_{t}(\bm{x})$ is strictly positive over the domain
$\mathcal{X}$. Let $\hat{r}(\bm{x})$ be the estimator of $r(\bm{x})$, and what
we want to achieve is to get $\hat{r}(\bm{x})$ and $r(\bm{x})$ as close as
possible.

####  2.2 Incremental Mixture Density Ratio Estimation

The IMDRE family includes Telescoping DRE [48] and $\infty$-DRE [12], and its
formulation is based on a density ratio chain with arbitrary bridge
distributions as follows.

$\displaystyle
r(\bm{x})=\frac{p_{s}(\bm{x})}{p_{t}(\bm{x})}=\frac{p_{s}(\bm{x})}{p_{1}(\bm{x})}\frac{p_{1}(\bm{x})}{p_{2}(\bm{x})}\cdots\frac{p_{m-1}(\bm{x})}{p_{m}(\bm{x})}\frac{p_{m}(\bm{x})}{p_{t}(\bm{x})},$
(2)

where $p_{1}(\bm{x}),\dots,p_{m}(\bm{x})$ are $m$ bridge distributions. They
consider the linear mixtures to gradually transport sample
$\\{\bm{x}_{i}^{s}\\}^{n_{s}}_{i=1}$ from the source distribution
$p_{s}(\bm{x})$ to sample $\\{\bm{x}_{i}^{t}\\}^{n_{t}}_{i=1}$ from the target
distribution $p_{t}(\bm{x})$ as
$\bm{x}_{i}^{k}=(1-\lambda)\bm{x}_{i}^{s}+\lambda\bm{x}_{i}^{t}$ with
$\lambda\in(0,1)$. Note that Rhodes et al. [48] use
$\lambda=\sqrt{1-a^{2}_{k}}$ where $a_{k}$ form an increasing sequence from
$0$ to $1$. For each successive density ratio, once the samples are taken,
standard approaches can be used to estimate the density ratio, such as kernel
logistic regression. Their experimental results report that such a framework
yields good density ratio estimates.

###  3 Geometric Reinterpretation of IMDRE

To begin with, we consider the statistical manifolds constructed by a set of
probability distributions in order to investigate the geometric behavior of
IMDRE. Let $\mathcal{M}=\\{p(\cdot;\bm{\theta})\mid\bm{\theta}\in\Theta\\}$ be
a statistical model parametrized by $\bm{\theta}\in\Theta$, where
$\Theta\subset\mathbb{R}^{b}$ is a parameter space. Let $\bm{G}=(g_{ij})$ be
the Fisher information matrix defined as
$g_{ij}=\mathbb{E}[\partial_{i}\ell(\bm{\theta})\partial_{j}\ell(\bm{\theta})]$,
where $\ell(\bm{\theta})=\ln p(\bm{x};\bm{\theta})$ and
$\partial_{i}=\partial/\partial\theta_{i}$. Here, it is known that
$\mathcal{M}$ is a Riemannian manifold with $\Theta$ as the coordinate system
and $\bm{G}$ as coefficients of the metric. A coordinate system is intuitively
like location information for determining a point on a manifold, which can be
understood from the fact that determining the parameters determines the
corresponding probability distribution. The metric enables us to measure
distances and angles on manifolds and to determine connections, and its
requirements are to be symmetric, positive definite, and non-degenerate. Under
appropriate regularity conditions, the Fisher information matrix satisfies
these requirements. We avoid using geometry notations more than necessary and
require only knowledge of basic covariant tensors. However, additional
knowledge of differential geometry is included in Appendix B.2 to make it
self-contained. For more details, refer to textbooks on differential geometry
[19, 34] and Riemannian geometry [17, 33, 36, 49]. Also see Appendix B.1 for
proofs of the theoretical results in this section.

Consider a smooth curve $\bm{\theta}(\lambda)\colon[0,1]\to\Theta$ in the
parameter space $\Theta$ with $\lambda\in[0,1]$, and a curve
$\gamma(\lambda)\colon[0,1]\to\mathcal{M}$ on the statistical manifold
$\mathcal{M}$ is defined as

$\displaystyle\gamma(\lambda)\coloneqq
p(\bm{x};\bm{\theta}(\lambda)),\quad\forall\lambda\in[0,1].$ (3)

The velocity of the curve $\gamma(\lambda)$ is given by
$\dot{\gamma}(\lambda)=(d/d\lambda)p(\bm{x};\bm{\theta}(\lambda))$. A curve
$\gamma$ is called $\nabla^{(\alpha)}$-autoparallel if $\dot{\gamma}(\lambda)$
is parallel transported along $\gamma(\lambda)$, that is, the acceleration
with respect to the $\nabla^{(\alpha)}$-connection vanishes

$\displaystyle\nabla^{(\alpha)}_{\dot{\gamma}(\lambda)}\dot{\gamma}(\lambda)=0,\quad\forall\lambda\in[0,1].$
(4)

Here, let $\mathcal{X}(\mathcal{M})$ be a set of vector fields on
$\mathcal{M}$, and $\nabla^{(\alpha)}$-connection is defined as

$\displaystyle\nabla^{(\alpha)}_{X}\coloneqq\frac{1+\alpha}{2}\nabla^{(1)}_{X}+\frac{1-\alpha}{2}\nabla^{(-1)}_{X},\quad\forall
X\in\mathcal{X}(\mathcal{M}),$ (5)

where $\nabla^{(1)}$ and $\nabla^{(-1)}$ are operators satisfying

$\displaystyle g(\nabla^{(1)}_{\partial_{i}}\partial_{j},\partial_{k})$
$\displaystyle=\mathbb{E}\left[(\partial_{i}\partial_{j}\ell)(\partial_{k}\ell)\right],$
(6) $\displaystyle g(\nabla^{(-1)}_{\partial_{i}}\partial_{j},\partial_{k})$
$\displaystyle=\mathbb{E}\left[(\partial_{i}\partial_{j}\ell+\partial_{i}\ell\partial_{j}\ell)(\partial_{k}\ell)\right],$
(7)

and $2$-covariant tensor $g(X,Y)$ is a Riemannian metric on $\mathcal{M}$ as

$\displaystyle g(X,Y)=\sum^{b}_{i,j=1}g_{ij}v_{i}w_{j},\quad
X=\sum^{b}_{i=1}v_{i}\partial_{i},\ Y=\sum^{n}_{i=1}w_{i}\partial_{i}.$ (8)

The operation
$\nabla^{(\cdot)}\colon\mathcal{X}(\mathcal{M})\times\mathcal{X}(\mathcal{M})\to\mathcal{X}(\mathcal{M})$
is called a linear connection, and see Appendix B.2 for the details. If the
curve in Eq. (3) is $\nabla^{(\alpha)}$-autoparallel at some
$\alpha\in\mathbb{R}$, it is called $\alpha$-geodesic. The explicit form of
$\alpha$-geodesics connecting two probability distributions is given as
follows.

######  Definition 1 ($\alpha$-geodesics [2]).

Let $p,q\in\mathcal{M}$ be two probability distributions. The
$\alpha$-geodesic $\gamma^{(\alpha)}(\lambda)\colon[0,1]\to\mathcal{M}$
connecting $p(\bm{x})$ and $q(\bm{x})$ is defined as

$\displaystyle\gamma^{(\alpha)}(\lambda)=\begin{cases}\left\\{(1-\lambda)p(\bm{x})^{\frac{1-\alpha}{2}}+\lambda
q(\bm{x})^{\frac{1-\alpha}{2}}\right\\}^{\frac{2}{1-\alpha}},&(\alpha\neq
1),\\\ \exp\left\\{(1-\lambda)\ln p(\bm{x})+\lambda\ln
q(\bm{x})\right\\},&(\alpha=1).\end{cases}$ (9)

ignoring the normalization factor.

Given an $\alpha\in\mathbb{R}$, the $\alpha$-geodesics are known to be the
minimizer of the following $\alpha$-divergence $D_{\alpha}[p\|q]$, that is,
for a fixed $\lambda\in[0,1]$,
$\gamma^{(\alpha)}(\lambda)=\operatorname*{arg\,min}_{\pi\in\mathcal{M}}((1-\lambda)D_{\alpha}[p\|\pi]+\lambda
D_{\alpha}[q\|\pi])$.

######  Definition 2 ($\alpha$-divergence [1]).

Let $\alpha\in\mathbb{R}$. For two probability distributions $p$ and $q$, the
$\alpha$-divergence
$D_{\alpha}\colon\mathcal{M}\times\mathcal{M}\to[0,+\infty)$ is defined as

$\displaystyle D_{\alpha}[p\|q]=\frac{1}{\alpha(\alpha-1)}\left(1-\int
p(\bm{x})^{\alpha}q(\bm{x})^{1-\alpha}d\bm{x}\right).$ (10)

Here, consider the following IMDRE with linear mixtures as bridge
distributions.

$\displaystyle
r(\bm{x})=\frac{p_{s}(\bm{x})}{p_{t}(\bm{x})}=\frac{p_{s}(\bm{x})}{(1-\lambda_{1})p_{s}(\bm{x})+\lambda_{1}p_{t}(\bm{x})}\cdots\frac{(1-\lambda_{m})p_{s}(\bm{x})+\lambda_{m}p_{t}(\bm{x})}{p_{t}(\bm{x})},$
(11)

where $\lambda_{1},\dots,\lambda_{m}\in(0,1)$ is a non-decreasing sequence.
For such an IMDRE, a simple calculation yields the following.

######  Theorem 1.

The IMDRE given by Eq. (11) can be regarded as a DRE along a geodesic with
$\alpha=-1$. Additionally, the curve to which the bridge distributions of this
IMDRE belong is the minimizer of the KL-divergence between $p_{s}$ and
$p_{t}$.

Recall that $\pi^{*}$ is a minimizer of KL-divergence
$D_{\mathrm{KL}}\colon\mathcal{M}\times\mathcal{M}\to[0,+\infty)$ with respect
to $\lambda$, $p_{s}$ and $p_{t}$ if it satisfies
$\pi^{*}=\operatorname*{arg\,min}_{\pi\in\mathcal{M}}(1-\lambda)D_{\mathrm{KL}}[p_{s}\|\pi]+\lambda
D_{\mathrm{KL}}[p_{t}\|\pi]$. Here, the geodesics of $\alpha=-1$ are specially
called $m$-geodesics ($m$ stands for mixture). Similarly, the geodesics of
$\alpha=+1$ are specially called $e$-geodesics ($e$ stands for exponential).

Now consider the evaluation of the density ratios obtained. Given that the Eq.
(1) can be transformed as in $p_{s}(\bm{x})=r(\bm{x})p_{t}(\bm{x})$, a well-
estimated density ratio should minimize

$\displaystyle D_{\mathrm{KL}}[p_{s}\|\hat{r}\cdot
p_{t}]\coloneqq\int_{\mathcal{X}}p_{s}(\bm{x})\log\frac{p_{s}(\bm{x})}{\hat{r}(\bm{x})p_{t}(\bm{x})}dx-1+\int_{\mathcal{X}}p_{t}(\bm{x})\hat{r}(\bm{x})d\bm{x},$
(12)

where $D_{\mathrm{KL}}[p\|q]$ is the unnormalized version of KL-divergence
between $p$ and $q$. We can then show the following theorem.

######  Theorem 2.

The evaluation along the m-geodesics of any density ratio estimator
$\hat{r}(\bm{x})$ is asymptotically evaluated along the $\alpha$-geodesics
with $\alpha=2$.

Furthermore, if the parameter $\alpha$ that determines the geodesic is well
determined, the divergence that appears can be bounded above by the Lipschitz
continuity. In fact, KL-divergence is not Lipschitz continuous since $x\to 0$
and $\ln x\to+\infty$, but the $\alpha$-divergence in the case of $\alpha=0$
satisfies this. Utilizing this, we obtain the following results from Sugiyama
et al. [55].

######  Corollary 1.

Under the appropriate assumptions, the convergence bound of a density ratio
estimator $\hat{r}$ evaluated along $\alpha$-geodesic with $\alpha=0$ is
$O_{p}(n^{-1/(2+\zeta_{n})}+\sqrt{\zeta})$, where $\zeta_{n}$ is a variable
depending on the sample size $n$. and $\zeta$ is some constant.

These results suggest the following:

  * •

Choosing how to construct the bridge distributions in IMDRE is equivalent to
choosing what curve on the statistical manifold to follow for density ratio
estimation.

  * •

The Behavior of the density ratio estimator depends on which curve to
construct the bridge distributions along. In fact, Theorem 2 shows that in
asymptotic situations geodesics with $\alpha=2$ and $\alpha=0$ appear.

These suggestions naturally motivate the development of IMDRE along general
$\alpha$-geodesics.

###  4 Optimization Algorithm for GIMDRE along Generalized Geodesics

Figure 1: Overview of alternating optimization algorithms for GIMDRE along
$\alpha$-geodesics.
$\gamma^{(\alpha)}(\lambda_{1}),\dots,\gamma^{(\alpha)}(\lambda_{m})$ are $m$
bridge distributions on the $\alpha$-geodesic parameterized by
$\lambda\in[0,1]$, respectively.

We now construct IMDRE along general $\alpha$-geodesics, Generalized IMDRE
(GIMDRE). To do so along the geodesics, we propose to utilize the following
importance sampling technique [61]: consider the expectation

$\displaystyle
I\coloneqq\mathbb{E}_{\gamma^{(\alpha)}(\lambda)(\bm{x})}\left[g(\bm{x})\right]=\int_{\mathcal{X}}g(\bm{x})\gamma^{(\alpha)}(\lambda)(\bm{x})d\bm{x},$
(13)

of any function $g(\bm{x})$ of a random variable
$\bm{x}\sim\gamma^{(\alpha)}(\lambda)(\bm{x})$. Let
$\\{\bm{x}_{1},\dots,\bm{x}_{T}\\}$ be the independent instances of sample
size $T$, and the standard Monte Carlo integration is

$\displaystyle\hat{I}\coloneqq\frac{1}{T}\sum^{N}_{i=1}g(\bm{x}_{i}).$ (14)

From the law of large numbers, $\hat{I}$ converges to $I$ as $T\to\infty$.
Now, since it is difficult to sample directly from
$\gamma^{(\alpha)}(\lambda)(\bm{x})$, considering another distribution
$\pi(\bm{x})$, which allows for easy sampling, and we have

$\displaystyle
I=\int_{\mathcal{X}}g(\bm{x})\gamma^{(\alpha)}(\lambda)(\bm{x})d\bm{x}=\int_{\mathcal{X}}\frac{\gamma^{(\alpha)}(\lambda)(\bm{x})}{\pi(\bm{x})}g(\bm{x})\pi(\bm{x})d\bm{x}=\mathbb{E}_{\pi(\bm{x})}\left[\frac{\gamma^{(\alpha)}(\lambda)(\bm{x})}{\pi(\bm{x})}g(\bm{x})\right].$
(15)

Then, let $w(\bm{x})=\gamma^{(\alpha)}(\lambda)(\bm{x})/\pi(\bm{x})$ and we
have

$\displaystyle\hat{I}=\frac{1}{T}\sum^{T}_{i=1}w(\bm{x}^{*}_{i})g(\bm{x}^{*}_{i})=\frac{1}{T}\sum^{T}_{i=1}\frac{\gamma^{(\alpha)}(\lambda)(\bm{x}^{*}_{i})}{\pi(\bm{x}^{*}_{i})}g(\bm{x}^{*}_{i}).$
(16)

where $\bm{x}^{*}_{i}\sim\pi(\bm{x})$. Thus, the original problem falls into
the problem of approximating the importance weighting $w(\bm{x})$. Note that
the normalization constant of the $\alpha$-geodesic is unknown, and we can
consider the self-normalized importance sampling technique. Then, let
$\gamma^{(\alpha)}(\lambda)(\bm{x})=q(\bm{x})/(\int_{\mathcal{X}}q(\bm{x})d\bm{x})\propto
q(\bm{x})$ for $q\in\mathcal{M}$, and rewrite Eq. (15) as
$I=\int_{\mathcal{X}}g(\bm{x})\frac{q(\bm{x})}{\int_{\mathcal{X}}q(x)d\bm{x}}=\frac{\int_{\mathcal{X}}g(\bm{x})q(\bm{x})}{\int_{\mathcal{X}}q(\bm{x})d\bm{x}}$.
Furthermore, if we assume that we are sampling from $p_{s}(\bm{x})$, we have

$\displaystyle
I=\frac{\int_{\mathcal{X}}g(\bm{x})\frac{q(\bm{x})}{p_{s}(\bm{x})}p_{s}(\bm{x})d\bm{x}}{\int_{\mathcal{X}}\frac{q(\bm{x})}{p_{s}(\bm{x})}p_{s}(\bm{x})d\bm{x}}=\frac{\mathbb{E}_{p_{s}(\bm{x})}\left[g(\bm{x})\frac{q(\bm{x})}{p_{s}(\bm{x})}\right]}{\mathbb{E}_{p_{s}(\bm{x})}\left[\frac{q(\bm{x})}{p_{s}(\bm{x})}\right]}.$
(17)

Applying importance sampling to the denominator and numerator of Eq. (17)
yields the following estimator.

$\displaystyle\hat{I}=\frac{\frac{1}{T}\sum^{T}_{i=1}g(\bm{x}^{*}_{i})w(\bm{x}^{*}_{i})}{\frac{1}{T}\sum^{T}_{i=1}w(\bm{x}^{*}_{i})}=\sum^{T}_{i=1}\frac{w(\bm{x}^{*}_{i})}{\sum^{T}_{j=1}w(\bm{x}^{*}_{j})}g(\bm{x}^{*}_{i})=\sum^{T}_{i=1}\frac{\frac{\gamma^{(\alpha)}(\lambda)(\bm{x}^{*}_{i})}{\pi(\bm{x}^{*}_{i})}}{\sum^{T}_{j=1}\frac{\gamma^{(\alpha)}(\lambda)(\bm{x}^{*}_{j})}{\pi(\bm{x}^{*}_{j})}}g(\bm{x}^{*}_{i}).$
(18)

In particular, suppose that the proxy distribution $\pi(\bm{x})=p_{s}(\bm{x})$
for importance sampling, sampling from the $\alpha$-geodesic connecting
$p_{s}(\bm{x})$ and $p_{t}(\bm{x})$ can be accomplished by the following
$w(\bm{x}^{*})$.

$\displaystyle
w_{\lambda}(\bm{x}^{*}_{i})=\frac{\gamma^{(\alpha)}(\lambda)(\bm{x}^{*}_{i})}{p_{s}(\bm{x}^{*}_{i})}=\left\\{1-\lambda+\lambda\left(\frac{p_{t}(\bm{x}^{*}_{i})}{p_{s}(\bm{x}^{*}_{i})}\right)^{\frac{1-\alpha}{2}}\right\\}^{\frac{2}{1-\alpha}}=\left\\{1-\lambda+\lambda\left(\frac{1}{r(\bm{x}^{*}_{i})}\right)^{\frac{1-\alpha}{2}}\right\\}^{\frac{2}{1-\alpha}}$

The key trick to obtain the density ratio
$\gamma^{(\alpha)}(\lambda_{i})(\bm{x})/\gamma^{(\alpha)}(\lambda_{j})(\bm{x})$
is the following cancellation out of the proxy distribution
$\pi(\bm{x})=p_{s}(\bm{x})$ using the ratio of the weighting functions
$w_{\lambda_{i}}(\bm{x})/w_{\lambda_{j}}(\bm{x})$.

$\displaystyle\frac{w_{\lambda_{i}}(\bm{x})}{w_{\lambda_{j}}(\bm{x})}=\frac{\gamma^{(\alpha)}(\lambda_{i})(\bm{x})/p_{s}(\bm{x})}{\gamma^{(\alpha)}(\lambda_{j})(\bm{x})/p_{s}(\bm{x})}=\frac{\gamma^{(\alpha)}(\lambda_{i})(\bm{x})}{\gamma^{(\alpha)}(\lambda_{j})(\bm{x})}.$
(19)

As a result, the ratio of the two distributions on the geodesic can be
rewritten in a form that depends only on the density ratio $r(\bm{x})$. Thus,
if a weighting function is obtained for each $\lambda$, the density ratios
between distributions on the geodesic can be recovered. Note that the case
$\pi(\bm{x})=p_{t}(\bm{x})$ can be computed almost similarly. For example, if
the function $g(\bm{x})$ is a loss function for a probabilistic classifier
$f(\bm{x})$, as in logistic regression, with $y=+1$ for $\bm{x}\sim p_{s}$ and
$y=-1$ for $\bm{x}\sim p_{t}$, the output learned by such a Monte Carlo
integral approximation can be used for density ratio estimation.

However, for the importance sampling, it is necessary to obtain an importance
weighting function $w_{\lambda}(\bm{x})$ that depends on the estimate of the
density ratio $r(\bm{x})$. In other words, the density ratio estimator and the
weighting function for sampling are interdependent and deadlocked.

Table 1: Evaluation results at different step sizes $m$ with $\alpha=3$. The
evaluation metric is the mean absolute error (MAE) between $\hat{r}$ and $r$.
The means and standard deviations of 10 trials are reported.

$m=10$ | $m=20$ | $m=30$ | $m=40$ | $m=50$ | $m=70$ | $m=100$
---|---|---|---|---|---|---
$35.41(\pm 3.55)$ | $34.80(\pm 3.09)$ | $34.61(\pm 2.60)$ | $34.32(\pm 2.54)$ | $34.26(\pm 2.47)$ | $34.18(\pm 2.39)$ | $34.13(\pm 2.34)$

Figure 2: Illustrative examples of GIMDRE convergence. Left panel: an example
of sample transport from $p_{s}$ to $p_{t}$ along $\alpha$-geodesic. Right
panel: examples of $\alpha$-geodesics connecting $p_{s}(\bm{x})$ and
$p_{t}(\bm{x})$ in two dimensional case. Here, we assume that
$p_{s}(\bm{x})=(0.1,0.1)$ and $p_{t}(\bm{x})=(0.9,0.9).$

To address this problem, we develop the following alternating optimization
algorithm.

  * i)

The density ratio estimator $\hat{r}(\bm{x})$ is estimated from the samples
$\\{\bm{x}_{i}\\}^{n_{s}}_{i=1}$ and $\\{\bm{x}_{i}\\}^{n_{t}}_{i=1}$.

  * ii)

Obtained $\hat{r}(\bm{x})$ is used to update the weighting
$w_{\lambda}(\bm{x})$ for the importance sampling.

  * iii)

Using the updated weighting $w_{\lambda}(\bm{x})$, $\hat{r}(\bm{x})$ is
estimated by GIMDRE along the $\alpha$-geodesic.

  * iv)

Iterate steps ii) and iii) over a predetermined number of times.

Figure 1 shows an overview of the above algorithm. The plug-in estimator that
underlies our algorithm can be expected to have different statistical behavior
depending on two parameters: $\alpha$, which determines the shape of the
curve, and $\lambda$, which determines the position on the curve.

####  4.1 Equally Spaced Transitions on Generalized Geodesic

The right panel of Figure 2 visualizes the trajectory approaching from the
source point to the target point along several $\alpha$-geodesics (see Section
5 for more details). In this figure, we have equally spaced transitions of
$\lambda$ from $0$ to $1$ in increments of $0.1$. As another way of selecting
$\lambda$, we can consider selecting equally spaced transitions on the
geodesic. That is, we can choose $\lambda_{1},\dots,\lambda_{m}$ such that the
arc lengths
$\|\gamma^{(\alpha)}(\lambda_{i})\|,\dots,\|\gamma^{(\alpha)}(\lambda_{m})\|$
of the $\alpha$-geodesic at time $\lambda_{1},\dots,\lambda_{m}$ are all
equal. Let $\gamma^{(\alpha)}(t)$ be the $\alpha$-geodesic. The length of this
curve connecting $\gamma^{(\alpha)}(0)=p_{s}$ and $\gamma^{(\alpha)}(1)=p_{t}$
is defined as

$\displaystyle\|\gamma^{(\alpha)}\|\coloneqq\int^{1}_{0}\sqrt{\sum_{i,j}g_{ij}\frac{d\gamma^{(\alpha)i}(t)}{dt}\frac{d\gamma^{(\alpha)j}(t)}{dt}},$
(20)

where $\gamma^{(\alpha)i}$ is the $i$-th component of $\gamma^{(\alpha)}$.
Then, we can see that equally spaced transitions
$\lambda_{1},\dots,\lambda_{m}$ on the $\alpha$-geodesic satisfy

$\displaystyle\int^{\lambda_{k}}_{\lambda_{k-1}}\sqrt{\sum_{i,j}g_{ij}\frac{d\gamma^{(\alpha)i}(t)}{dt}\frac{d\gamma^{(\alpha)j}(t)}{dt}}=\frac{\|\gamma^{(\alpha)}\|}{m},\quad
1\leq k\leq m.$ (21)

As can be seen from the right panel of Figure 2, there is a large jump in the
approaching distribution for large $\alpha$. Therefore, adjusting the
transition speed of $\lambda$ may be useful in practical applications.

###  5 Numerical Experiments

In this section, we design numerical experiments to investigate the behavior
of GIMDRE and report the results. Our series of numerical experiments consists
of two parts: i) illustrative examples to clarify GIMDRE behavior, and ii)
experiments on the two-sample test, a major application of density ratio in
statistics. See Appendix C for more details on the experiments.

####  5.1 Illustrative Examples

First, we consider illustrative examples to clarify the behavior of GIMDRE.
The purpose of the experiments in this section is to verify the following:
does the alternating optimization algorithm for GIMDRE converge and give good
estimates, and what geodesics are better to estimate along when varying the
sample size, the number of dimensions, and the distance between the two
distributions. We generate samples
$\\{\bm{x}_{i}^{s}\\}^{n_{s}}_{i=1},\\{\bm{x}_{i}^{t}\\}^{n_{t}}_{i=1}$ of
sample sizes $n_{s}$ and $n_{t}$, respectively, from multivariate Gaussian
distributions
$\mathcal{N}(\bm{\mu}_{s},\bm{\Sigma}_{s}),\mathcal{N}(\bm{\mu}_{t},\bm{\Sigma}_{t})$
with $\bm{\mu}_{s},\bm{\mu}_{t}\in\mathbb{R}^{d}$ and
$\bm{\Sigma}_{s},\bm{\Sigma}_{t}\in\mathbb{R}^{d\times d}$. The values of
these parameters are specified specifically for each experiment.

Figure 3: Behavior of GIMDRE according to $\alpha$ when two distributions are
brought closer together. KLD is the analytically calculated KL-divergence
value

First, for illustrative purposes, we consider the case where $\alpha=3$ and
the samples follow univariate Gaussian distributions as
$p_{s}=\mathcal{N}(8,3)$ and $p_{t}=\mathcal{N}(0,2)$. The sample sizes are
$n_{s}=n_{t}=500$. The base method for DRE is based on a probabilistic
classifier using logistic regression. This setup is a typical example where
classical DRE fails because the two distributions are far apart. That is, it
is known that classical DRE is difficult unless the means of the two
distributions are somewhat close and the variance of the distribution on the
denominator side is sufficiently large compared to that on the numerator side.
In fact, the mean absolute error (MAE) obtained by ordinary DRE with the base
method is $264.07(\pm 116)$. Note that hereafter we will report the mean and
standard deviation of the $10$ trials unless otherwise specified. Table 1
shows the results of GIMDRE evaluations when different step sizes are used.
These results show that even with a small number of steps, there is a
significant improvement compared to the base method, and that increasing the
number of steps further improves the mean and standard deviation of the
evaluation results. The left panel of Figure 2 also shows an illustration of
sample transport using the estimated density ratio $\hat{r}$. This figure
shows that $\hat{r}$ is well able to transport from $p_{s}=\mathcal{N}(8,3)$
to $\mathcal{N}(0,2)$. In addition, the right panel of Figure 2 shows examples
of $\alpha$-geodesics in two dimensional case. In this figure,
$\lambda\in\\{0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0\\}$ for all $\alpha$.
It can be seen that for $\alpha=-1$, one approaches the target distribution at
a uniform rate as the $\lambda$ increases.

Next, Figure 3 shows an illustration of the behavior of GIMDRE according to
$\alpha$ when two distributions are brought together. The parameter
$\mu_{s}\in\\{8,6,4,3,2,1\\}$, and other parameters are fixed as $m=100$,
$n=n_{s}=n_{t}=500$, $p_{s}=\mathcal{N}(\mu_{s},3)$ and
$p_{t}=\mathcal{N}(0,2)$. To quantify the discrepancy between the two
distributions, we analytically compute the KL-divergence. This result shows
that the use of a large $\alpha$ consistently suppresses the variance of the
estimator. In terms of bias, a large $\alpha$ is effective when the KL-
divergence between the two distributions is large, while the effect is
reversed as the distributions get closer together. Tables 2 and 3 show the
evaluation results for different sample sizes and sample dimensions. It can be
seen that the degradation of the estimator with decreasing sample size and
increasing dimensionality is reduced by using a larger $\alpha$.

Table 2: Evaluation results with different $\alpha$ and sample sizes
$n=n_{s}=n_{t}$.

| $n=100$ | $n=200$ | $n=300$ | $n=400$ | $n=500$
---|---|---|---|---|---
$\alpha=-1$ | $180.06(\pm 48.14)$ | $126.48(\pm 25.52)$ | $80.67(\pm 16.16)$ | $64.30(\pm 7.21)$ | $59.88(\pm 4.50)$
$\alpha=3$ | $76.39(\pm 20.77)$ | $54.99(\pm 14.23)$ | $43.60(\pm 7.15)$ | $36.36(\pm 3.59)$ | $34.13(\pm 2.34)$
$\alpha=7$ | $69.14(\pm 18.51)$ | $41.26(\pm 12.75)$ | $35.33(\pm 5.62)$ | $26.11(\pm 2.506)$ | $23.27(\pm 1.86)$

From the above evaluation results, we can make the following suggestions: i)
if the two distributions are close enough, a good estimator can be obtained
with a small $\alpha$, and as the two distributions move away from each other,
the usefulness of using a large $\alpha$ increases. In particular, a large
$\alpha$ consistently leads to a variance reduction in the estimator.

Table 3: Evaluation results with different $\alpha$ and sample dimensions $d$.

| $d=2$ | $d=3$ | $d=4$ | $d=5$ |
---|---|---|---|---|---
$\alpha=-1$ | $260.55(\pm 58.23)$ | $589.62(\pm 103.06)$ | $6030.19(\pm 4845.87)$ | $16337.70(\pm 9379.11)$ |
$\alpha=3$ | $123.83(\pm 18.23)$ | $320.65(\pm 28.31)$ | $504.50(\pm 47.42)$ | $3664.92(\pm 287.02)$ |
$\alpha=7$ | $122.74(\pm 16.17)$ | $331.91(\pm 24.26)$ | $692.35(\pm 44.36)$ | $3792.96(\pm 283.45)$ |

####  5.2 Two-Sample Test based on GIMDRE

In statistics, given two sets of samples, testing whether the probability
distributions behind the samples are equivalent is a fundamental task. This
problem is referred to as the two-sample test or homogeneity test [51, 13,
18]. In our experiments, we utilize the two-sample test based on the
permutation test [15, 56]. Let $X^{s}$ and $X^{t}$ be two samples from $p_{s}$
and $p_{t}$. We first estimate the Pearson divergence using $X^{s}$ and
$X^{t}$ as

$\displaystyle\hat{D}_{\mathrm{PE}}[X^{s}\|X^{t}]\coloneqq\frac{1}{2n_{s}}\sum^{n_{s}}_{i=1}\hat{r}(\bm{x}^{s}_{i})-\frac{1}{n_{t}}\sum^{n_{t}}_{i=1}\hat{r}(\bm{x}^{t}_{i})+\frac{1}{2}.$
(22)

Next, we randomly permute the $\left|X^{s}\cup X^{t}\right|$ samples, and
assign the first $n_{s}$ samples to a set $\tilde{X}$ and the remaining
$n_{t}$ samples to another set $\tilde{X}^{\prime}$. Then, we estimate the
Pearson divergence again using the randomly shuffled samples $\tilde{X}$ and
$\tilde{X}^{\prime}$, and we obtain an estimate
$\hat{D}_{\mathrm{PE}}[\tilde{X}\|\tilde{X}^{\prime}]$. This random shuffling
procedure is repeated $K$ times, and the distribution of
$\hat{D}_{\mathrm{PE}}[\tilde{X}\|\tilde{X}^{\prime}]$ under the null
hypothesis (that is, the two distributions are the same) is constructed. In
this experiment, we set $K=100$. Finally, the $p$-value is approximated by
evaluating the relative ranking of $\hat{D}_{\mathrm{PE}}[X^{s}\|X^{t}]$ in
the distribution of $\hat{D}_{\mathrm{PE}}[\tilde{X}\|\tilde{X}^{\prime}]$.
Suppose that $n_{s}=n_{t}$, and let $F$ be the distribution function of
$\hat{D}_{\mathrm{PE}}[\tilde{X}\|\tilde{X}^{\prime}]$. Let
$\beta\coloneqq\sup\\{t\in\mathbb{R}\mid F(t)\leq 1-\alpha\\}$ be the upper
$100\alpha$-percentile point of $F$. If $p_{s}=p_{t}$, it was shown that
$\mathbb{P}(\hat{D}_{\mathrm{PE}}[X^{s}\|X^{t}]>\beta)\leq\alpha$. That is,
for a given significance level $\alpha$, the probability that
$\hat{D}[X^{s}\|X^{t}]$ exceeds $\beta$ is at most $\alpha$ when
$p_{s}=p_{t}$. Thus, when the null hypothesis is correct, it will be properly
accepted with a specified probability. Let $n=n_{s}=n_{t}=500$, and
$p_{s}=\mathcal{N}(0,1)$. We consider the setups as $\mu_{t}\in\\{0,0.5,1\\}$
and $\sigma_{t}\in\\{1,1.5\\}$. The histogram plots in Figure 4 show the
distribution of estimated
$\hat{D}_{\mathrm{PE}}[\tilde{X}\|\tilde{X}^{\prime}]$, with red crosses
indicating $\hat{D}_{PE}[X^{s}\|X^{t}]$. Also, the line plots show the
transition of the $p$-values with increasing sample size in the different
settings. For all the alphas, when $p_{s}\neq p_{t}$, the $p$-value gets
smaller with each increase in sample size, while if the two distributions are
equal, it always produces a large $p$-value. The figure also shows that the
GIMDRE with large $\alpha$ behaves conservatively.

Figure 4: Experimental results of two-sample test. The histogram plots show
the distribution of estimated
$\hat{D}_{\mathrm{PE}}[\tilde{X}\|\tilde{X}^{\prime}]$, with red crosses
indicating $\hat{D}_{PE}[X^{s}\|X^{t}]$. The line plots show the transition of
the $p$-values with increasing sample size in the different settings.

###  6 Conclusion and Discussion

In this paper, we first geometrically reinterpreted the framework of
incremental density ratio estimation using a mixture between two
distributions. The effectiveness of such a framework for density ratio
estimation has been reported in recent years, and it is very useful to study
this behavior in detail. Our analysis demonstrated that considering what kind
of bridge distributions to create is equivalent to choosing what kind of curve
on a statistical manifold. Such a geometrically intuitive interpretation not
only reveals properties of the algorithm, but also induces natural
generalizations. Next, we considered density ratio estimation along arbitrary
generalized geodesics on the manifold. This procedure requires sampling from
curves on the manifold, called $\alpha$-geodesics, which is non-trivial. This
is because the two probability distributions are generally not given
explicitly. Here we have shown that by using the importance sampling
framework, this sampling can be written in a form that depends on the density
ratio that we originally wanted to obtain. This is just the state of
interdependence, and we have demonstrated that this deadlock can be resolved
with a simple alternating optimization procedure. We have shown that this
strategy works through illustrative examples and experiments on hypothesis
testing, an important task in statistics. Limitlations and broader impacts of
this study are discussed in Appendix A. We believe that this study provides
useful insights from a geometric perspective into a fundamental task in
statistics.

### References

  * Amari [2009] S.-I. Amari.  $\alpha$-divergence is unique, belonging to both $f$-divergence and bregman divergence classes.  _IEEE Transactions on Information Theory_ , 55(11):4925–4931, 2009.
  * Amari [2016] S.-i. Amari.  _Information geometry and its applications_ , volume 194.  Springer, 2016.
  * Amari and Nagaoka [2000] S.-i. Amari and H. Nagaoka.  _Methods of information geometry_ , volume 191.  American Mathematical Soc., 2000.
  * Azmandian et al. [2012a] F. Azmandian, J. G. Dy, J. A. Aslam, and D. R. Kaeli.  Local kernel density ratio-based feature selection for outlier detection.  In _Asian Conference on Machine Learning_ , pages 49–64. PMLR, 2012a.
  * Azmandian et al. [2012b] F. Azmandian, A. Yilmazer, J. G. Dy, J. A. Aslam, and D. R. Kaeli.  Gpu-accelerated feature selection for outlier detection using the local kernel density ratio.  In _2012 IEEE 12th International Conference on Data Mining_ , pages 51–60. IEEE, 2012b.
  * Braga [2014] I. Braga.  A constructive density-ratio approach to mutual information estimation: experiments in feature selection.  _Journal of Information and Data Management_ , 5(1):134–134, 2014.
  * Byrd and Lipton [2019] J. Byrd and Z. Lipton.  What is the effect of importance weighting in deep learning?  In _International conference on machine learning_ , pages 872–881. PMLR, 2019.
  * Chen et al. [2016] X. Chen, M. Monfort, A. Liu, and B. D. Ziebart.  Robust covariate shift regression.  In _Artificial Intelligence and Statistics_ , pages 1270–1279. PMLR, 2016.
  * Chen et al. [2015] Y.-T. Chen, W.-H. Fang, C.-Y. Lee, and K.-W. Cheng.  Abnormal detection in crowded scenes via kernel based direct density ratio estimation.  In _2015 IEEE China Summit and International Conference on Signal and Information Processing (ChinaSIP)_ , pages 230–234. IEEE, 2015.
  * Cheng and Chu [2004] K. F. Cheng and C.-K. Chu.  Semiparametric density estimation under a two-sample density ratio model.  _Bernoulli_ , 10(4):583–604, 2004.
  * Choi et al. [2021] K. Choi, M. Liao, and S. Ermon.  Featurized density ratio estimation.  In _Uncertainty in Artificial Intelligence_ , pages 172–182. PMLR, 2021.
  * Choi et al. [2022] K. Choi, C. Meng, Y. Song, and S. Ermon.  Density ratio estimation via infinitesimal classification.  In _International Conference on Artificial Intelligence and Statistics_ , pages 2552–2573. PMLR, 2022.
  * Conover et al. [1981] W. J. Conover, M. E. Johnson, and M. M. Johnson.  A comparative study of tests for homogeneity of variances, with applications to the outer continental shelf bidding data.  _Technometrics_ , 23(4):351–361, 1981.
  * Dudoit et al. [2003] S. Dudoit, J. P. Shaffer, and J. C. Boldrick.  Multiple hypothesis testing in microarray experiments.  _Statistical Science_ , 18(1):71–103, 2003.
  * Efron and Tibshirani [1993] B. Efron and R. J. Tibshirani.  Permutation tests.  _An introduction to the bootstrap_ , pages 202–219, 1993.
  * Elvira and Martino [2014] V. Elvira and L. Martino.  Advances in importance sampling.  _Wiley StatsRef: Statistics Reference Online_ , pages 1–14, 2014.
  * Gallot et al. [1990] S. Gallot, D. Hulin, J. Lafontaine, et al.  _Riemannian geometry_ , volume 2.  Springer, 1990.
  * Gretton et al. [2006] A. Gretton, K. Borgwardt, M. Rasch, B. Schölkopf, and A. Smola.  A kernel method for the two-sample-problem.  _Advances in neural information processing systems_ , 19, 2006.
  * Guggenheimer [2012] H. W. Guggenheimer.  _Differential geometry_.  Courier Corporation, 2012.
  * Harris et al. [2020] C. R. Harris, K. J. Millman, S. J. Van Der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, et al.  Array programming with numpy.  _Nature_ , 585(7825):357–362, 2020.
  * Heck [2019] D. W. Heck.  A caveat on the savage–dickey density ratio: The case of computing bayes factors for regression parameters.  _British Journal of Mathematical and Statistical Psychology_ , 72(2):316–333, 2019.
  * Hido et al. [2011] S. Hido, Y. Tsuboi, H. Kashima, M. Sugiyama, and T. Kanamori.  Statistical outlier detection using direct density ratio estimation.  _Knowledge and information systems_ , 26:309–336, 2011.
  * Huang et al. [2006] J. Huang, A. Gretton, K. Borgwardt, B. Schölkopf, and A. Smola.  Correcting sample selection bias by unlabeled data.  _Advances in neural information processing systems_ , 19, 2006.
  * Hushchyn and Ustyuzhanin [2021] M. Hushchyn and A. Ustyuzhanin.  Generalization of change-point detection in time series data based on direct density ratio estimation.  _Journal of Computational Science_ , 53:101385, 2021.
  * Kanamori et al. [2008] T. Kanamori, S. Hido, and M. Sugiyama.  Efficient direct density ratio estimation for non-stationarity adaptation and outlier detection.  _Advances in neural information processing systems_ , 21, 2008.
  * Kanamori et al. [2010] T. Kanamori, T. Suzuki, and M. Sugiyama.  Theoretical analysis of density ratio estimation.  _IEICE transactions on fundamentals of electronics, communications and computer sciences_ , 93(4):787–798, 2010.
  * Kanamori et al. [2011] T. Kanamori, T. Suzuki, and M. Sugiyama.  $f$-divergence estimation and two-sample homogeneity test under semiparametric density-ratio models.  _IEEE Transactions on Information Theory_ , 58(2):708–720, 2011.
  * Kawahara and Sugiyama [2009] Y. Kawahara and M. Sugiyama.  Change-point detection in time-series data by direct density-ratio estimation.  In _Proceedings of the 2009 SIAM international conference on data mining_ , pages 389–400. SIAM, 2009.
  * Keziou and Leoni-Aubin [2005] A. Keziou and S. Leoni-Aubin.  Test of homogeneity in semiparametric two-sample density ratio models.  _Comptes Rendus Mathématique_ , 340(12):905–910, 2005.
  * Keziou and Leoni-Aubin [2008] A. Keziou and S. Leoni-Aubin.  On empirical likelihood for semiparametric two-sample density ratio models.  _Journal of Statistical Planning and Inference_ , 138(4):915–928, 2008.
  * Kimura and Hino [2022] M. Kimura and H. Hino.  Information geometrically generalized covariate shift adaptation.  _Neural Computation_ , 34(9):1944–1977, 2022.
  * Kimura and Hino [2024] M. Kimura and H. Hino.  A short survey on importance weighting for machine learning.  _Transactions on Machine Learning Research_ , 2024.  ISSN 2835-8856.  URL https://openreview.net/forum?id=IhXM3g2gxg.  Survey Certification.
  * Klingenberg [1995] W. Klingenberg.  _Riemannian geometry_ , volume 1.  Walter de Gruyter, 1995.
  * Kreyszig [2013] E. Kreyszig.  _Differential geometry_.  Courier Corporation, 2013.
  * LaFond and Watts [2008] R. LaFond and R. L. Watts.  The information role of conservatism.  _The accounting review_ , 83(2):447–478, 2008.
  * Lee [2018] J. M. Lee.  _Introduction to Riemannian manifolds_ , volume 2.  Springer, 2018.
  * Lee and Lee [2018] S. Lee and D. K. Lee.  What is the proper way to apply the multiple comparison test?  _Korean journal of anesthesiology_ , 71(5):353, 2018.
  * Li et al. [2022] K. Li, X. Gao, S. Fu, X. Diao, P. Ye, B. Xue, J. Yu, and Z. Huang.  Robust outlier detection based on the changing rate of directed density ratio.  _Expert Systems with Applications_ , 207:117988, 2022.
  * Liu and Ziebart [2014] A. Liu and B. Ziebart.  Robust classification under sample selection bias.  _Advances in neural information processing systems_ , 27, 2014.
  * Liu et al. [2013] S. Liu, M. Yamada, N. Collier, and M. Sugiyama.  Change-point detection in time-series data by relative density-ratio estimation.  _Neural Networks_ , 43:72–83, 2013.
  * Lu et al. [2020] Y.-h. Lu, A.-y. Liu, M.-j. Jiang, and T. Jiang.  A new two-part test based on density ratio model for zero-inflated continuous distributions.  _Applied Mathematics-A Journal of Chinese Universities_ , 35(2):203–219, 2020.
  * Martin et al. [2023] J. I. S. Martin, S. Mazuelas, and A. Liu.  Double-weighting for covariate shift adaptation.  In _International Conference on Machine Learning_ , pages 30439–30457. PMLR, 2023.
  * Matsushita et al. [2023] Y. Matsushita, T. Otsu, and K. Takahata.  Estimating density ratio of marginals to joint: Applications to causal inference.  _Journal of Business & Economic Statistics_, 41(2):467–481, 2023.
  * Nguyen et al. [2007] X. Nguyen, M. J. Wainwright, and M. Jordan.  Estimating divergence functionals and the likelihood ratio by penalized convex risk minimization.  _Advances in neural information processing systems_ , 20, 2007.
  * Oh et al. [2016] M.-S. Oh, J. Choi, and E. S. Park.  Bayesian variable selection in quantile regression using the savage–dickey density ratio.  _Journal of the Korean Statistical Society_ , 45(3):466–476, 2016.
  * Pedregosa et al. [2011] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, et al.  Scikit-learn: Machine learning in python.  _the Journal of machine Learning research_ , 12:2825–2830, 2011.
  * Reichenheim and Coutinho [2010] M. E. Reichenheim and E. S. Coutinho.  Measures and models for causal inference in cross-sectional studies: arguments for the appropriateness of the prevalence odds ratio and related logistic regression.  _BMC medical research methodology_ , 10:1–12, 2010.
  * Rhodes et al. [2020] B. Rhodes, K. Xu, and M. U. Gutmann.  Telescoping density-ratio estimation.  _Advances in neural information processing systems_ , 33:4905–4916, 2020.
  * Sakai [1996] T. Sakai.  _Riemannian geometry_ , volume 149.  American Mathematical Soc., 1996.
  * Shimodaira [2000] H. Shimodaira.  Improving predictive inference under covariate shift by weighting the log-likelihood function.  _Journal of statistical planning and inference_ , 90(2):227–244, 2000.
  * Stein [1945] C. Stein.  A two-sample test for a linear hypothesis whose power is independent of the variance.  _The Annals of Mathematical Statistics_ , 16(3):243–258, 1945.
  * Storey et al. [2004] J. D. Storey, J. E. Taylor, and D. Siegmund.  Strong control, conservative point estimation and simultaneous conservative consistency of false discovery rates: a unified approach.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 66(1):187–205, 2004.
  * Sugiyama et al. [2007a] M. Sugiyama, M. Krauledat, and K.-R. Müller.  Covariate shift adaptation by importance weighted cross validation.  _Journal of Machine Learning Research_ , 8(5), 2007a.
  * Sugiyama et al. [2007b] M. Sugiyama, S. Nakajima, H. Kashima, P. Buenau, and M. Kawanabe.  Direct importance estimation with model selection and its application to covariate shift adaptation.  _Advances in neural information processing systems_ , 20, 2007b.
  * Sugiyama et al. [2008] M. Sugiyama, T. Suzuki, S. Nakajima, H. Kashima, P. Von Bünau, and M. Kawanabe.  Direct importance estimation for covariate shift adaptation.  _Annals of the Institute of Statistical Mathematics_ , 60:699–746, 2008.
  * Sugiyama et al. [2011] M. Sugiyama, T. Suzuki, Y. Itoh, T. Kanamori, and M. Kimura.  Least-squares two-sample test.  _Neural networks_ , 24(7):735–751, 2011.
  * Sugiyama et al. [2012a] M. Sugiyama, T. Suzuki, and T. Kanamori.  _Density ratio estimation in machine learning_.  Cambridge University Press, 2012a.
  * Sugiyama et al. [2012b] M. Sugiyama, T. Suzuki, and T. Kanamori.  Density-ratio matching under the bregman divergence: a unified framework of density-ratio estimation.  _Annals of the Institute of Statistical Mathematics_ , 64:1009–1044, 2012b.
  * Suzuki and Sugiyama [2010] T. Suzuki and M. Sugiyama.  Sufficient dimension reduction via squared-loss mutual information estimation.  In _Proceedings of the thirteenth international conference on artificial intelligence and statistics_ , pages 804–811. JMLR Workshop and Conference Proceedings, 2010.
  * Suzuki et al. [2008] T. Suzuki, M. Sugiyama, J. Sese, and T. Kanamori.  Approximating mutual information by maximum likelihood density ratio estimation.  In _New challenges for feature selection in data mining and knowledge discovery_ , pages 5–20. PMLR, 2008.
  * Tokdar and Kass [2010] S. T. Tokdar and R. E. Kass.  Importance sampling: a review.  _Wiley Interdisciplinary Reviews: Computational Statistics_ , 2(1):54–60, 2010.
  * Virtanen et al. [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, et al.  Scipy 1.0: fundamental algorithms for scientific computing in python.  _Nature methods_ , 17(3):261–272, 2020.
  * Yamada et al. [2013] M. Yamada, T. Suzuki, T. Kanamori, H. Hachiya, and M. Sugiyama.  Relative density-ratio estimation for robust distribution comparison.  _Neural computation_ , 25(5):1324–1370, 2013.

## Appendix

###  Appendix A Limitations and Broader Impact

In this section, we discuss about both the limitations and the broader impacts
of our study.

####  A.1 Limitations

One limitation is that our algorithm does not include a neural network. Of
course, on the other hand, this is a promising extension direction.
Modification of our algorithm to take advantage of the expressive power of
deep learning could lead to further advances in this research. Deep learning
models have demonstrated remarkable capabilities in capturing complex patterns
and relationships in data, offering a prospect for enhancing the performance
and scope of our algorithm. On the other hand, an existing study reported that
the effect of importance-weighted learning by deep learning decays with the
number of steps [7], so it is not clear whether our method based on importance
sampling is straightforwardly extendable or not. For the application of neural
networks to our research, we need to take into account such negative results
reported by existing studies, and it seems that we need to devote a large
space to this.

####  A.2 Broader Impact

As highlighted in the Introduction section, the density ratio stands out as a
fundamental tool in statistics, playing a pivotal role in various fields.
Therefore, while our algorithm itself does not directly produce social
impacts, its downstream applications may, depending on how it is used. One
example is the observation in Section 5, where the choice of parameters leads
to conservative hypothesis tests. How sensitive a hypothesis test should be
depends on whether it is used in a life-threatening setting, such as the
medical field, or a place where challenging decisions are allowed, such as
marketing [14, 35, 37, 52]. Our numerical experiments suggest that the choice
of parameters can control how sensitive a hypothetical test is constructed,
which could have a sufficient impact on society.

###  Appendix B Technical Details

In this section, we provide supplementary technical details, including proofs
of the theorems and additional discussion of the ESS.

####  B.1 Proofs of Theoretical Results

###### Proof of Theorem 1.

Recall that $\alpha$-geodesics between $p$ and $q$ are given as

$\displaystyle\gamma^{(\alpha)}(\lambda)=\left\\{(1-\lambda)p(\bm{x})^{\frac{1-\alpha}{2}}+\lambda
q(\bm{x})^{\frac{1-\alpha}{2}}\right\\}^{\frac{2}{1-\alpha}}.$ (23)

Substituting $\alpha=-1$, we see that this is an ordinary weighted average.
Also, $\alpha$-divergence becomes KL-divergence with $\alpha=-1$, thus
obtaining the proof. ∎

###### Proof of Theorem 2.

Considering the following Taylor expansion.

$\displaystyle\log\frac{p_{s}(x)}{\hat{r}(x)p_{t}(x)}$
$\displaystyle=\log\frac{r(x)}{\hat{r}(x)}$
$\displaystyle=-\log\frac{\hat{r}(x)}{r(x)}$
$\displaystyle=-\left(\frac{\hat{r}(x)}{r(x)}-1\right)+\frac{1}{2}\left(\frac{\hat{r}(x)}{r(x)}-1\right)^{2}+O_{p}\left(\left|\frac{\hat{r}(x)}{r(x)}-1\right|^{3}\right).$
(24)

Let

$\displaystyle J_{1}$
$\displaystyle=\int_{\mathcal{X}}p_{s}(x)\log\frac{p_{s}(x)}{\hat{r}(x)p_{t}(x)}dx-1,$
(25) $\displaystyle J_{2}$
$\displaystyle=\int_{\mathcal{X}}p_{t}(x)\hat{r}(x)dx.$ (26)

Then,

$\displaystyle J_{1}$
$\displaystyle=\int_{\mathcal{X}}p(x)\log\frac{p_{s}(x)}{\hat{r}(x)p_{t}(x)}dx-1$
$\displaystyle=\int_{\mathcal{X}}\left\\{-\left(\frac{\hat{r}(x)}{r(x)}-1\right)+\frac{1}{2}\left(\frac{\hat{r}(x)}{r(x)}-1\right)^{2}\right\\}p_{s}(x)dx-1+O(\|\hat{r}/r-1\|^{3})$
$\displaystyle=\int_{\mathcal{X}}\left\\{-\frac{\hat{r}(x)}{r(x)}+1+\frac{1}{2}\left(\left(\frac{\hat{r}(x)}{r(x)}\right)^{2}-2\frac{\hat{r}(x)}{r(x)}+1\right)\right\\}p_{s}(x)dx-\int_{\mathcal{X}}p_{s}(x)dx+O(\|\hat{r}/r-1\|^{3})$
$\displaystyle=\int_{\mathcal{X}}\left\\{-\frac{\hat{r}(x)}{r(x)}+\frac{1}{2}\left(\left(\frac{\hat{r}(x)}{r(x)}\right)^{2}-2\frac{\hat{r}(x)}{r(x)}+1\right)\right\\}p_{s}(x)dx+O(\|\hat{r}/r-1\|^{3})$
$\displaystyle=\int_{\mathcal{X}}\left\\{-\frac{\hat{r}(x)}{r(x)}+\frac{1}{2}\frac{\hat{r}(x)^{2}}{r(x)^{2}}-\frac{\hat{r}(x)}{r(x)}+\frac{1}{2}\right\\}p_{s}(x)dx+O(\|\hat{r}/r-1\|^{3})$
$\displaystyle=\int_{\mathcal{X}}\left\\{-2\frac{\hat{r}(x)}{r(x)}+\frac{1}{2}\frac{\hat{r}(x)^{2}}{r(x)^{2}}+\frac{1}{2}\right\\}p_{s}(x)dx+O(\|\hat{r}/r-1\|^{3})$
$\displaystyle=\int_{\mathcal{X}}\left\\{-2\hat{r}(x)+\frac{1}{2}\frac{\hat{r}(x)^{2}}{r(x)}+\frac{1}{2}r(x)\right\\}\frac{1}{r(x)}p_{s}(x)dx+O(\|\hat{r}/r-1\|^{3})$
$\displaystyle=\int_{\mathcal{X}}\left\\{-2\hat{r}(x)+\frac{1}{2}\frac{\hat{r}(x)^{2}p_{t}(x)}{p_{s}(x)}+\frac{1}{2}\frac{p_{s}(x)}{p_{t}(x)}\right\\}p_{t}(x)dx+O(\|\hat{r}/r-1\|^{3}).$
(27)

Here,

$\displaystyle\|\hat{r}/r-1\|\coloneqq\left(\int
p_{s}(x)|\hat{r}(x)/r(x)-1|^{2}dx\right)^{1/2}.$ (28)

We have

$\displaystyle D_{\mathrm{KL}}[p_{s}\|\hat{r}\cdot p_{t}]$
$\displaystyle=J_{1}+J_{2}$
$\displaystyle=\int_{\mathcal{X}}\left\\{-2\hat{r}(x)+\frac{1}{2}\frac{\hat{r}(x)^{2}p_{t}(x)}{p_{s}(x)}+\frac{1}{2}\frac{p_{s}(x)}{p_{t}(x)}\right\\}p_{t}(x)dx+\int_{\mathcal{X}}p_{t}(x)\hat{r}(x)dx$
$\displaystyle=\int_{\mathcal{X}}\left\\{-\hat{r}(x)+\frac{1}{2}\frac{\hat{r}(x)^{2}p_{t}(x)}{p_{s}(x)}+\frac{1}{2}\frac{p_{s}(x)}{p_{t}(x)}\right\\}p_{t}(x)dx$
$\displaystyle=\frac{1}{2}\int_{\mathcal{X}}\left\\{-2\hat{r}(x)+\frac{\hat{r}(x)^{2}p_{t}(x)}{p_{s}(x)}+\frac{p_{s}(x)}{p_{t}(x)}\right\\}p_{t}(x)dx$
$\displaystyle=\frac{1}{2}\int_{\mathcal{X}}\frac{1}{p_{s}(x)}\left\\{-2p_{s}(x)\cdot\hat{r}(x)p_{t}(x)+\hat{r}(x)^{2}p_{t}(x)^{2}+p_{s}(x)^{2}\right\\}dx$
$\displaystyle=\frac{1}{2}\int_{\mathcal{X}}\frac{(p_{s}(x)-\hat{r}(x)p_{t}(x))^{2}}{p_{s}(x)}dx+O(\|\hat{r}/r-1\|^{3})$
$\displaystyle=D_{\mathrm{PE}}[p_{s}\|\hat{r}(x)p_{t}(x)]+O(\|\hat{r}/r-1\|^{3}).$
(29)

Here, $D_{\mathrm{PE}}[\cdot\|\cdot]$ is the Pearson divergence. Hence, the
asymptotic expansion of the unnormalized KL-divergence between $p_{s}(x)$ and
$\hat{r}(x)p_{t}(x)$ is given as

$\displaystyle D_{\mathrm{KL}}[p_{s}\|\hat{r}\cdot
p_{t}]=D_{\mathrm{PE}}[p_{s}\|\hat{r}\cdot p_{t}]+O(n^{-3/2}).$ (30)

We consider the following $\alpha$-divergence.

$\displaystyle D_{\alpha}[p_{s}\|p_{t}]=\frac{1}{\alpha(\alpha-1)}\left(1-\int
p_{s}(x)^{\alpha}p_{t}(x)^{1-\alpha}dx\right).$ (31)

From the simple calculation, we can confirm that Pearson divergence is a
special case of $\alpha$-divergene with $\alpha=2$. Since
$D_{\mathrm{KL}}[p_{s}\|p_{t}]$ and $D_{\alpha}[p_{s}\|p_{t}]$ are minimizers
of $\alpha$-geodesics with $\alpha=-1$ and $\alpha=2$ between $p_{s}$ and
$p_{t}$, we have the proof. ∎

###### Proof of Corollary 1.

We assume that two distributions $p_{s}$ and $p_{t}$ are mutually absolutely
continuous, and satisfy $0<c_{0}\leq dp_{s}/dp_{t}\leq c_{1}$ on the support
of $p_{s}$ and $p_{t}$. We also assume that for any function
$f\colon\to\mathbb{R}$, there exist $\epsilon$ and $\delta$ such that
$\mathbb{E}_{p_{t}}[f(x)]\geq\epsilon,\ \|f\|_{\infty}\leq\delta$. Finally, we
assume that for some constants $0<\zeta<2$ and $M$, $\sup_{q\in\mathcal{M}}\ln
N(\epsilon,F^{M},L_{2}(q))$, where $F$ is a set of finite linear combinations
of integrable functions $f$ with positive coefficients bounded above by $M$,
and $N(\epsilon,F^{M},L_{2}(q))$ is the $\epsilon$-covering number with
$L_{2}$-distance [44]. From the definition of $\zeta_{n}$,

$\displaystyle-\frac{1}{n_{s}}\sum^{n_{2}}_{i=1}\ln\hat{r}(\bm{x}^{s}_{i})\leq-\frac{1}{n_{s}}\sum^{n_{s}}_{i=1}\ln\frac{r(\bm{x}^{s}_{i})}{\frac{1}{n^{t}_{i}}\sum^{n_{t}}_{i=1}r(\bm{x}^{t}_{i})}+\zeta_{n}.$
(32)

By convexity of $-\ln x$,

$\displaystyle-\frac{1}{n_{s}}\sum^{n_{s}}_{i=1}\ln\left(\frac{\hat{r}(\bm{x}^{s}_{i})+r(\bm{x}^{s}_{i})\frac{1}{n^{t}_{i}}\sum^{n_{t}}_{i=1}r(\bm{x}^{t}_{i})}{2r(\bm{x}^{s}_{i})\frac{1}{n^{t}_{i}}\sum^{n_{t}}_{i=1}r(\bm{x}^{t}_{i})}\right)$
$\displaystyle\leq\frac{\zeta_{n}}{2}.$ (33)

Let $\xi_{1}=\hat{r}(\bm{x}^{s}_{i})$ and
$\xi_{2}=2r(\bm{x}^{s}_{i})\frac{1}{n^{t}_{i}}\sum^{n_{t}}_{i=1}r(\bm{x}^{t}_{i})$,
and we have

$\displaystyle(\hat{p}_{t}-p_{t})(\xi_{2}-\xi_{1})-(\hat{p}_{s}-p_{s})\ln\frac{\xi_{2}}{\xi_{1}}-\frac{\zeta_{n}}{2}$
$\displaystyle\leq-
p_{t}(\xi_{2}-\xi_{1})-p_{s}\left(\ln\frac{\xi_{2}}{\xi_{1}}\right)$
$\displaystyle\leq
2p_{s}\left(\sqrt{\frac{\xi_{2}}{\xi_{1}}}-1\right)-p_{t}(\xi_{2}-\xi_{1})$
$\displaystyle=p_{t}\left(2\sqrt{\xi_{2}\xi_{1}}-\xi_{2}-\xi_{1}\right).$ (34)

Here, the last line is the generalized version of Hellinger distance, and it
corresponds to $\alpha$-divergence with $\alpha=0$. The rest of the proof
follows immediately from Sugiyama et al. [55]. ∎

Figure 5: Monotonicity of $\alpha$-geodesics. The left panel shows that the
value of the $\alpha$-geodesic decreases monotonically as $\alpha$ is
increased. Here, $p_{s}(\bm{x})=0.1$, $p_{t}(\bm{x})=0.1$, and $\lambda$ is
set to a constant of $0.5$. The right panel shows that the $\alpha$-geodesic
approaches $p_{s}$ at $\lambda\to 0$ and $p_{t}$ at $\lambda\to 1$ for all
$\alpha$. Figure 6: Examples of $\alpha$-geodesics connecting $p_{s}(\bm{x})$
and $p_{t}(\bm{x})$ in two dimensional case. Here, we assume that
$p_{s}(\bm{x})=(0.1,0.1)$ and $p_{t}(\bm{x})=(0.9,0.9).$

####  B.2 Generalized Geodesics and Information Geometry

In our manuscript, we have tried to deliver to the reader the usefulness of
the geometric interpretation of statistical procedures, while eliminating as
much as possible the notations of differential geometry. However, in order to
make our manuscript more self-contained, in this section we introduce a few
selected geometric concepts that may help in understanding our work.

Let $g_{ij}$ be a Riemannian metric, particularly the Fisher information
matrix for the statistical manifold. Then the most simple connection on the
manifold is defined by the following Christoffel symbols of the first kind.

$\displaystyle\Gamma_{ij,k}\coloneqq\frac{1}{2}\left(\partial_{i}g_{jk}+\partial_{j}g_{ik}-\partial_{k}g_{ij}\right).$
(35)

The Levi-Civita connection $\nabla^{(0)}$ is defined as

$\displaystyle
g(\nabla^{(0)}_{\partial_{i}}\partial_{j},\partial_{k})=\Gamma_{ij,k}.$ (36)

Here, the operation
$\nabla\colon\mathcal{X}(\mathcal{M})\times\mathcal{X}(\mathcal{M})\to\mathcal{X}(\mathcal{M})$
on a differentiable manifold $\mathcal{M}$ is called a linear connection if it
satisfies

  * i)

$\nabla_{X}Y$ is $\mathcal{F}(\mathcal{M})$-linear in $X$,

  * ii)

$\nabla_{X}Y$ is $\mathbb{R}$-linear in $Y$, and

  * iii)

$\nabla$ satisfies the Leibniz rule, that is

$\displaystyle\nabla_{X}(fY)=(Xf)Y+f\nabla_{X}Y,\quad\forall
f\in\mathcal{F}(\mathcal{M}),$ (37)

where $\mathcal{X}(\mathcal{M})$ is the set of vector fields on $\mathcal{M}$
and $\mathcal{F}(\mathcal{M})$ is the set of all differentiable functions on
$\mathcal{M}$. It is worth nothing that the superscript of $\nabla^{(0)}$
denotes a parameter of the connection, and the generalized connection induced
by the $\alpha$-divergence is denoted as $\nabla^{(\alpha)}$. In this case,
$\alpha=0$ corresponds to the Levi-Civita connection. Two other important
special cases of the $\alpha$-connection are the $\nabla^{(1)}$\- and
$\nabla^{(-1)}$-connections, given as follows.

$\displaystyle g(\nabla^{(1)}_{\partial_{i}}\partial_{j},\partial_{k})$
$\displaystyle=\Gamma^{(1)}_{ij,k}\coloneqq\mathbb{E}\left[(\partial_{i}\partial_{j}\ell)(\partial_{k}\ell)\right],$
(38) $\displaystyle g(\nabla^{(-1)}_{\partial_{i}}\partial_{j},\partial_{k})$
$\displaystyle=\Gamma^{(-1)}_{ij,k}\coloneqq\mathbb{E}\left[(\partial_{i}\partial_{j}\ell+\partial_{i}\ell\partial_{j}\ell)(\partial_{k}\ell)\right].$
(39)

Note that $\alpha$-divergence with $\alpha=\pm 1$ is the KL-divergence and its
dual. Here, the relationship between $\alpha$-divergence and
$\alpha$-connection is that under $\alpha$-connection, the $\alpha$-geodesic
given by Eq. (10) becomes a straight path, and $\alpha$-divergence becomes its
minimizer. The $\alpha$-connection is given as

$\displaystyle
g(\nabla^{(\alpha)}_{\partial_{i}}\partial_{j},\partial_{k})=\Gamma^{(\alpha)}_{ij,k}\coloneqq\mathbb{E}\left[\left(\partial_{i}\partial_{j}\ell+\frac{1-\alpha}{2}\partial_{i}\ell\partial_{j}\ell\right)\partial_{k}\ell\right],$
(40)

and

$\displaystyle\nabla^{(\alpha)}=\frac{1+\alpha}{2}\nabla^{(1)}+\frac{1-\alpha}{2}\nabla^{(-1)}.$
(41)

For example, let us consider the following Gaussian distribution.

$\displaystyle
p(x;\bm{\theta})=\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\left\\{-\frac{(x-\mu)^{2}}{2\sigma^{2}}\right\\},$
(42)

where $\bm{\theta}=(\theta_{1},\theta_{2})=(\mu,\sigma)$. The Fisher
information matrix is given by

$\displaystyle g_{ij}=\begin{pmatrix}\frac{1}{\sigma^{2}}&0\\\
0&\frac{2}{\sigma^{2}}\end{pmatrix}.$ (43)

Then, from the simple calculation, the Christoffel symbols of first and second
kind are given as follows.

$\displaystyle\Gamma_{11,2}$ $\displaystyle=\frac{1}{\sigma^{3}},$ (44)
$\displaystyle\Gamma_{12,1}$ $\displaystyle=-\frac{1}{\sigma^{3}},$ (45)
$\displaystyle\Gamma_{22,2}$ $\displaystyle=-\frac{2}{\sigma^{3}},$ (46)
$\displaystyle\Gamma^{1}_{ij}$
$\displaystyle=\begin{pmatrix}0&-\frac{1}{\sigma}\\\
-\frac{1}{\sigma}&0\end{pmatrix},$ (47) $\displaystyle\Gamma^{2}_{ij}$
$\displaystyle=\begin{pmatrix}\frac{1}{2\sigma}&0\\\
0&-\frac{1}{\sigma}\end{pmatrix},$ (48)

and we can see that the geodesics are solutions of the following ordinary
differential equations.

$\displaystyle\ddot{\mu}-\frac{2}{\sigma}\dot{\mu}\dot{\sigma}$
$\displaystyle=0,$ (49)
$\displaystyle\ddot{\sigma}+\frac{1}{2\sigma}(\dot{\mu})^{2}-\frac{1}{\sigma}(\dot{\sigma})^{2}$
$\displaystyle=0.$ (50)

We can then have

$\displaystyle\frac{\ddot{\mu}}{\dot{\mu}}=\frac{2\dot{\sigma}}{\sigma}\Leftrightarrow\dot{\mu}=c\sigma^{2},$
(51)

with some constant $c$. In the case of $c\neq 0$ we have

$\displaystyle\sigma(s)^{2}+(\mu(s)-K)^{2}\frac{J}{\sigma^{4}}=c^{2},$ (52)

where $K>0$ and $J>0$ are positive constants, and the geodesics are half-
ellipses, with $\sigma>0$. In the case of $c=0$,

$\displaystyle\sigma(s)=\sigma(0)e^{\sqrt{J}s},$ (53)

and the geodesics are vertical half-lines. For the manifold of Gaussian
distributions, the Christoffel coefficients of first kind are written as

$\displaystyle\Gamma^{(\alpha)}_{11,1}$
$\displaystyle=\Gamma^{(\alpha)}_{21,2}=\Gamma^{(\alpha)}_{12,2}=\Gamma^{(\alpha)}_{22,1}=0,$
(54) $\displaystyle\Gamma^{(\alpha)}_{11,2}$
$\displaystyle=\frac{1-\alpha}{\sigma^{3}},$ (55)
$\displaystyle\Gamma^{(\alpha)}_{12,1}$
$\displaystyle=\Gamma^{(\alpha)}_{21,1}=\frac{1+\alpha}{\sigma^{3}},$ (56)
$\displaystyle\Gamma^{(\alpha)}_{22,2}$
$\displaystyle=\frac{2(1+2\alpha)}{\sigma^{3}},$ (57)

and the Christoffel symbols of second kind are as follows.

$\displaystyle\Gamma^{1(\alpha)}_{ij}$
$\displaystyle=g_{11}\Gamma^{(\alpha)}_{ij,1}+g_{12}\Gamma^{(\alpha)}_{ij,2}=\sigma^{2}\Gamma^{(\alpha)}_{ij,1}$
$\displaystyle=\sigma^{2}\begin{pmatrix}0&-\frac{1+\alpha}{\sigma^{3}}\\\
-\frac{1+\alpha}{\sigma^{3}}&0\end{pmatrix}=\begin{pmatrix}0&-\frac{1+\alpha}{\sigma}\\\
-\frac{1+\alpha}{\sigma}&0\end{pmatrix}.$ (58)

The ordinary differential equations for the $\alpha$-autoparallel curves are
given as

$\displaystyle\ddot{\mu}-\frac{2(1+\alpha)}{\sigma}\dot{\sigma}\dot{\mu}$
$\displaystyle=0,$ (59)
$\displaystyle\ddot{\sigma}+\frac{1-\alpha}{2\sigma}\dot{\mu}^{2}-\frac{1+2\alpha}{\sigma}\dot{\sigma}^{2}$
$\displaystyle=0.$ (60)

The first equation can be transformed as

$\displaystyle\frac{\ddot{\mu}}{\mu}=2(1+\alpha)\frac{\dot{\sigma}}{\sigma}$
$\displaystyle\Leftrightarrow\frac{d}{ds}\ln\dot{\mu}=2(1+\alpha)\frac{d}{ds}\ln\sigma$
$\displaystyle\Leftrightarrow\ln\dot{\mu}=2(1+\alpha)\ln\sigma+c_{0}$
$\displaystyle\Leftrightarrow\dot{\mu}=c\sigma^{2(1+\alpha)},$ (61)

with constant $c$. Then, we have

$\displaystyle\ddot{\sigma}+\frac{1-\alpha}{2\sigma}c^{2}\sigma^{4(1+\alpha)}-\frac{1+2\alpha}{\sigma}\dot{\sigma}^{2}$
$\displaystyle=0$
$\displaystyle\sigma^{k+1}du+\left(\frac{1-\alpha}{2}c^{2}\sigma^{4(\alpha+1)+k}-(1+2\alpha)\sigma^{k}\mu^{2}\right)d\sigma$
$\displaystyle=0,$ (62)

where $u=\dot{\sigma}$. Here, the following must be satisfied.

$\displaystyle\frac{\partial}{\partial\sigma}\sigma^{k+1}u=\frac{\partial}{\partial\mu}\left(\frac{1-\alpha}{2}c^{2}\sigma^{4(\alpha+1)+k}-(1+2\alpha)\sigma^{k}\mu^{2}\right).$
(63)

By considering this condition, we can determine $k+1=-(4\alpha+2)$, and we
have

$\displaystyle
u\sigma^{-(4\alpha+2)}du+\left(\frac{1-\alpha}{2}c^{2}\sigma-(1+2\alpha)u^{2}\sigma^{-(4\alpha+3)}\right)d\sigma=0.$
(64)

Hence, for a constant $E$, we can solve for $\sigma$ as follows.

$\displaystyle\frac{\mu^{2}}{\sigma^{2}(2\alpha+1)}+\frac{1-\alpha}{2}c^{2}\sigma^{2}=E$
$\displaystyle\Leftrightarrow\left(\frac{\dot{\sigma}}{\sigma^{2\alpha+1}}\right)^{2}+\frac{1-\alpha}{2}c^{2}\sigma^{2}=E$
$\displaystyle\Leftrightarrow\left(\frac{\dot{\sigma}}{\sigma^{2\alpha+1}}\right)^{2}=E-\frac{1-\alpha}{2}c^{2}\sigma^{2}$
$\displaystyle\Leftrightarrow\int\frac{d\sigma}{\sigma^{2\alpha+1}\sqrt{E-\frac{1-\alpha}{2}c^{2}\sigma^{2}}}=\pm
s+s_{0}$
$\displaystyle\Leftrightarrow\int\frac{d\sigma}{\sigma^{2\alpha+1}\sqrt{C^{2}-\sigma^{2}}}=(\pm
s+s_{0})\sqrt{\frac{1-\alpha}{2}}c,$ (65)

where

$\displaystyle C=C_{\alpha}=\frac{2E}{c}\frac{1}{1-\alpha}.$ (66)

Let $t=\sigma^{2}$ and $v=\sqrt{C^{2}-t}$, we have

$\displaystyle\int\frac{d\sigma}{\sigma^{2\alpha+1}\sqrt{C^{2}-\sigma^{2}}}$
$\displaystyle=\int\frac{dt}{2\sigma^{2(\alpha+1)}\sqrt{C^{2}-\sigma^{2}}}$
$\displaystyle=\int\frac{dt}{2t^{(\alpha+1)}\sqrt{C^{2}-1}}$
$\displaystyle=\int\frac{-2vdv}{2t^{\alpha+1}v}$
$\displaystyle=-\int\frac{dv}{(C^{2}-v^{2})^{\alpha+1}},$ (67)

and then,

$\displaystyle-\int\frac{dv}{(C^{2}-v^{2})^{\alpha+1}}=(\pm
s+s_{0})\sqrt{\frac{1-\alpha}{2}}c.$ (68)

The $\mu$ is given by

$\displaystyle\mu=c\int\sigma^{2(1+\alpha)}(s)ds.$ (69)

######  B.2.0.1 The Case of $\alpha=-1$

In the case of $\alpha=-1$, we have

$\displaystyle-v-K=(\pm s+s_{0})\sqrt{\frac{1-\alpha}{2}}c,$ (70)

with solution

$\displaystyle\sigma^{2}(s)=C^{2}-\left((\pm
s+s_{0})\sqrt{\frac{1-\alpha}{2}}c+K\right)^{2},$ (71)

for a constant $K$. Then, we have

$\displaystyle\mu(s)=cs+\mu(0).$ (72)

######  B.2.0.2 The case of $\alpha=1/2$

In the case of $\alpha=1/2$, since

$\displaystyle\int\frac{dv}{(C^{2}-v^{2})^{3/2}}=\frac{v}{C^{2}\sqrt{C^{2}-v^{2}}},$
(73)

and we have

$\displaystyle-\frac{v}{C^{2}\sqrt{C^{2}-v^{2}}}=(\pm s+s_{0})\frac{c}{2}+K,$
(74)

with solution

$\displaystyle\sigma(s)=\frac{C}{\sqrt{1+C^{4}\left(\frac{c}{2}(\pm
s+S_{0})+K\right)^{2}}},$ (75)

for a constant $K$. The, we have

$\displaystyle\mu(s)=c\int\sigma^{3}(s)ds.$ (76)

####  B.3 Approximation of Arc Lengths of Generalized Geodesics for Equally
Spaced Transitions

Let $\gamma^{(\alpha)}(t)$ be the $\alpha$-geodesic, and we consider its
Taylor expansion as

$\displaystyle\gamma^{(\alpha)i}(t)=\gamma^{(\alpha)i}(0)+t\frac{d\gamma^{(\alpha)i}}{dt}\Bigr{|}_{t=0}+\frac{t^{2}}{2}\frac{d^{2}\gamma^{(\alpha)i}}{dt^{2}}\Bigr{|}_{t=0}.$
(77)

Here, $\alpha$-geodesics satisfy the following Euler-Lagrange equation.

$\displaystyle\frac{d^{2}\gamma^{(\alpha)i}}{dt^{2}}+\Gamma^{(\alpha)i}_{jk}\frac{d\gamma^{(\alpha)j}}{dt}\frac{d\gamma^{(\alpha)k}}{dt}=0,$
(78)

where

$\displaystyle\Gamma^{(\alpha)}_{ij,k}=\mathbb{E}\left[\left(\partial_{i}\partial_{j}\ell+\frac{1-\alpha}{2}\partial_{i}\ell\partial_{j}\ell\right)\partial_{k}\ell\right],$
(79)

and $\ell=\ln p(x)$. Then, we can rewrite as

$\displaystyle\gamma^{(\alpha)i}(t)=\gamma^{(\alpha)i}(0)+t\frac{d\gamma^{(\alpha)i}(t)}{dt}\Bigr{|}_{t=0}-\frac{t^{2}}{2}\left(\Gamma^{(\alpha)i}_{jk}\frac{d\gamma^{(\alpha)j}}{dt}\frac{d\gamma^{(\alpha)k}}{dt}\right)\Bigr{|}_{t=0}.$
(80)

The length of this curve connecting $\gamma^{(\alpha)}(0)=p_{s}$ and
$\gamma^{(\alpha)}(1)=p_{t}$ is defined as

$\displaystyle\|\gamma^{(\alpha)}\|\coloneqq\int^{1}_{0}\sqrt{\sum_{i,j}g_{ij}\frac{d\gamma^{(\alpha)i}(t)}{dt}\frac{d\gamma^{(\alpha)j}(t)}{dt}}.$
(81)

Finally, we can see that equally spaced transitions
$\lambda_{1},\dots,\lambda_{m}$ on the $\alpha$-geodesic satisfies

$\displaystyle\int^{\lambda_{i}}_{\lambda_{i-1}}\sqrt{\sum_{i,j}g_{ij}\frac{d\gamma^{(\alpha)i}(t)}{dt}\frac{d\gamma^{(\alpha)j}(t)}{dt}}=\int^{\lambda_{j}}_{\lambda_{j-1}}\sqrt{\sum_{i,j}g_{ij}\frac{d\gamma^{(\alpha)i}(t)}{dt}\frac{d\gamma^{(\alpha)j}(t)}{dt}},\quad
1\leq i,j\leq m.$ (82)

Figure 7: Relationship between ESS and the parameters $\alpha$ and $\lambda$.

####  B.4 Additional Discussion on Variance and Effective Sample Size

We provide some insight into the question about the choice of $\alpha$ through
the following analysis.

######  Proposition 1.

For the estimator $\hat{I}$, we can see that the variance
$\mathrm{Var}(\hat{I})$ decreases as $\alpha\to+\infty$.

###### Proof.

The variance of $\hat{I}$ with $p_{\alpha}^{\lambda}$ is calculated as
follows.

$\displaystyle\mathrm{Var}_{\alpha,\lambda}(\hat{I})$
$\displaystyle=\int\left(\frac{\gamma^{(\alpha)}(\lambda)(x)}{p_{s}(x)}g(x)-I\right)^{2}p_{s}(x)dx$
$\displaystyle=\int\left(\frac{\gamma^{(\alpha)}(\lambda)(x)^{2}}{p_{s}(x)^{2}}g(x)^{2}-2\frac{\gamma(\alpha)(\lambda)(x)}{p_{s}(x)}g(x)I+I^{2}\right)p_{s}(x)dx$
$\displaystyle=\int\frac{\gamma^{(\alpha)}(\lambda)(x)^{2}}{p_{s}(x)}g(x)^{2}dx-I^{2}$
$\displaystyle=\int\frac{\gamma^{(\alpha)}(\lambda)(x)}{p_{s}(x)}g(x)^{2}\gamma^{(\alpha)}(\lambda)(x)dx-I^{2}$
$\displaystyle=\mathbb{E}_{\gamma^{(\alpha)}(\lambda)}\left[\frac{\gamma^{(\alpha)}(\lambda)(x)}{p_{s}(x)}g(x)^{2}\right]-I^{2}$

Thus, the variance of the estimator by sampling from the alpha-mixture depends
on the expectation of $g(x)$ scaled by the weight
$w(x)=\gamma^{(\alpha)}(\lambda)(x)/p_{s}(x)$. Here, for all $\alpha$, we have
$\frac{\partial}{\partial\alpha}\frac{\gamma^{(\alpha)}(\lambda)(x)}{p_{s}(x)}\leq
0$, and then we have the proof. ∎

The above proof is based on the monotonicity of $\alpha$-geodesics with
respect to $\alpha$. Therefore, we give a few more observations on this
important key monotonicity. Figure 5 shows the monotonicity of the
$\alpha$-geodesic. We use $p_{s}=0.1$ and $p_{t}=0.9$. The left panel shows
that the value of the $\alpha$-geodesic decreases monotonically as $\alpha$ is
increased. Here, $\lambda$ is set to a constant of $0.5$. Recall that
$\alpha=-1$ is an ordinary weighted average, and we can see that the geodesic
value is $(p_{s}(\bm{x})+p_{t}(\bm{x}))/2$ at $\alpha=-1$ in the figure. The
right panel shows that the $\alpha$-geodesic approaches $p_{s}$ at $\lambda\to
0$ and $p_{t}$ at $\lambda\to 1$ for all $\alpha$. Further, for larger values
of $\alpha$, the speed at which the $\alpha$-geodesic approaches $p_{t}$ is
found to be non-uniform. This behavior, in which the geodesic gradually
approaches the target distribution from the source distribution, can be
expected to affect the stability of the estimator. In addition, Figure 6 shows
the examples of $\alpha$-geodesics connecting $p_{s}(\bm{x})$ and
$p_{t}(\bm{x})$ in two dimensional case. Here, we assume that
$p_{s}(\bm{x})=(0.1,0.1)$ and $p_{t}(\bm{x})=(0.9,0.9).$ In this figure,
$\lambda\in\\{0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0\\}$ for all $\alpha$.
It can be seen that for $\alpha=-1$, one approaches the target distribution at
a uniform rate as the $\lambda$ increases, while for larger values of
$\alpha$, there is an acceleration in approaching the target distribution. The
above results follow from the monotonicity of $\alpha$-geodesics with respect
to $\alpha$. Note that, for arbitrary $\lambda\in[0,1]$,
$\lim_{\alpha\to+\infty}w(\bm{x})=\min\left\\{1,p_{t}(\bm{x})/p_{s}(\bm{x})\right\\}$.
This can be seen as a truncation of the density ratio by $1$ in the case of
$\alpha\to+\infty$.

Finally, the choice of the geodesic parameter $\alpha$ is discussed in more
detail. The appropriate parameter $\alpha$ for sampling from
$\alpha$-geodesics depends on $g(\bm{x})$, the target of the Monte Carlo
integration. Let
$\bar{Z}_{\alpha}\coloneqq\sum^{T}_{i=1}w(\bm{x}^{*})\bm{x}_{i}/\sum^{T}_{j=1}w(\bm{x}_{j})$,
this variance is given by

$\displaystyle\mathrm{Var}(\bar{Z}_{\alpha})=\mathrm{Var}(\bm{x})\frac{\sum^{T}_{i=1}w(\bm{x})^{2}}{\left(\sum^{T}_{i=1}w(\bm{x})\right)^{2}}=\mathrm{Var}(\bm{x})\frac{\sum^{T}_{i=1}\left\\{1-\lambda+\lambda\left(1/\hat{r}(\bm{x}^{s}_{i})\right)^{\frac{1-\alpha}{2}}\right\\}^{\frac{4}{1-\alpha}}}{\left(\sum^{T}_{i=1}\left\\{1-\lambda+\lambda\left(1/\hat{r}(\bm{x}^{s}_{i})\right)^{\frac{1-\alpha}{2}}\right\\}^{\frac{2}{1-\alpha}}\right)^{2}}.$
(84)

On the other hand, the variance of
$\bar{Z}=\frac{1}{T_{e}}\sum^{T_{e}}_{i=1}\bm{x}_{i}$ is
$\mathrm{Var}(\bar{Z})=\mathrm{Var}(\bm{x})/T_{e}$ for some $T_{e}>0$. Then,
let $\mathrm{Var}(\bar{Z}_{\alpha})=\mathrm{Var}(\bar{Z})$ and solve for
$T_{e}$ to obtain the following.

$\displaystyle\hat{T}_{e}=\frac{\left(\sum^{T}_{i=1}w(\bm{x})\right)^{2}}{\sum^{T}_{i=1}w(\bm{x})^{2}}=\frac{\left(\sum^{T}_{i=1}\left\\{1-\lambda+\lambda\left(1/\hat{r}(\bm{x}^{s}_{i})\right)^{\frac{1-\alpha}{2}}\right\\}^{\frac{2}{1-\alpha}}\right)^{2}}{\sum^{T}_{i=1}\left\\{1-\lambda+\lambda\left(1/\hat{r}(\bm{x}^{s}_{i})\right)^{\frac{1-\alpha}{2}}\right\\}^{\frac{4}{1-\alpha}}}.$
(85)

$\hat{T}_{e}$ represents the number of data required to obtain the same
accuracy as estimation by $\bar{Z}_{\alpha}$ when estimating the expected
value of a sample by $\bar{Z}$. It is called the effective sample size (ESS).
ESS holds for $\hat{T}_{e}\leq T$, and the equality holds when w is a
constant. When the value of $\hat{T}_{e}$ is small, estimation by
$\bar{Z}_{\alpha}$ requires more samples than $\bar{Z}$, meaning that
estimation is inefficient. From the monotonicity of the $\alpha$-geodesic with
respect to $\alpha$, we obtain the following.

######  Proposition 2.

The estimation efficiency in terms of ESS of GIMDRE improves with
$\alpha\to+\infty$.

Recall that ESS $\hat{T}_{e}$ is given as

$\displaystyle\hat{T}_{e}=\frac{\left(\sum^{T}_{i=1}w(\bm{x})\right)^{2}}{\sum^{T}_{i=1}w(\bm{x})^{2}}=\frac{\left(\sum^{T}_{i=1}\left\\{1-\lambda+\lambda\left(1/\hat{r}(\bm{x}^{s}_{i})\right)^{\frac{1-\alpha}{2}}\right\\}^{\frac{2}{1-\alpha}}\right)^{2}}{\sum^{T}_{i=1}\left\\{1-\lambda+\lambda\left(1/\hat{r}(\bm{x}^{s}_{i})\right)^{\frac{1-\alpha}{2}}\right\\}^{\frac{4}{1-\alpha}}}.$
(86)

Here, we can rewrite this as follows.

$\displaystyle\hat{T}_{e}$
$\displaystyle=\frac{T}{1+\frac{\frac{1}{T}\sum^{T}_{i=1}\left\\{1-\lambda+\lambda\left(1/\hat{r}(\bm{x}^{s}_{i})\right)^{\frac{1-\alpha}{2}}\right\\}^{\frac{4}{1-\alpha}}-\left(\frac{1}{T}\sum^{T}_{i=1}\left\\{1-\lambda+\lambda\left(1/\hat{r}(\bm{x}^{s}_{i})\right)^{\frac{1-\alpha}{2}}\right\\}^{\frac{2}{1-\alpha}}\right)^{2}}{\left(\frac{1}{T}\sum^{T}_{i=1}\left\\{1-\lambda+\lambda\left(1/\hat{r}(\bm{x}^{s}_{i})\right)^{\frac{1-\alpha}{2}}\right\\}^{\frac{2}{1-\alpha}}\right)^{2}}}$
$\displaystyle=\frac{T}{1+V^{2}},$ (87)

where

$\displaystyle
V\coloneqq\frac{\frac{1}{T}\sum^{T}_{i=1}\left\\{1-\lambda+\lambda\left(1/\hat{r}(\bm{x}^{s}_{i})\right)^{\frac{1-\alpha}{2}}\right\\}^{\frac{4}{1-\alpha}}-\left(\frac{1}{T}\sum^{T}_{i=1}\left\\{1-\lambda+\lambda\left(1/\hat{r}(\bm{x}^{s}_{i})\right)^{\frac{1-\alpha}{2}}\right\\}^{\frac{2}{1-\alpha}}\right)^{2}}{\left(\frac{1}{T}\sum^{T}_{i=1}\left\\{1-\lambda+\lambda\left(1/\hat{r}(\bm{x}^{s}_{i})\right)^{\frac{1-\alpha}{2}}\right\\}^{\frac{2}{1-\alpha}}\right)^{2}}.$
(88)

The value $V$ is called the coefficient of variation, and from the fact that
$V^{2}\geq 0$, we can see that $\hat{T}_{e}\leq T$. Figure 7 shows the
relationship between ESS and the parameters $\alpha$ and $\lambda$. In this
figure, $p_{s}=\mathcal{N}(8,3)$ and $p_{t}=\mathcal{N}(0,2)$. The left panel
of Figure 7 shows the relationship between $\alpha$ and ESS with $\lambda=0.5$
as a constant. We see that the estimation efficiency in terms of the ESS is
improved by using a large $\alpha$. Furthermore, the right panel of Figure 7
shows that the estimation efficiency deteriorates as one moves away from the
source distribution with larger $\lambda$. However, it can be observed that
this degradation can be reduced by using a large $\alpha$. Figures 8 and 9
show the results for the two skewed distributions, Log-normal and Power-law,
and reveal similar observations. Here, the density functions of these two
distributions $p_{\mathrm{LogNormal}}$ and $p_{\mathrm{PowerLaw}}$ are defined
as follows.

$\displaystyle p_{\mathrm{LogNormal}}(x;\mu,\sigma)$
$\displaystyle=\frac{1}{x\sigma\sqrt{2\pi}}\exp\left\\{-\frac{(\ln
x-\mu)^{2}}{2\sigma^{2}}\right\\},\quad x\in(0,+\infty),\
\mu\in(-\infty,+\infty),\ \sigma>0,$ $\displaystyle
p_{\mathrm{PowerLaw}}(x;a)$ $\displaystyle=ax^{a-1},\quad 0\leq 1,\ a>0.$

For Log-normal distributions $p_{s}=\mathrm{LogNormal}(\mu_{s},\sigma_{s})$
and $\mathrm{LogNormal}(p_{t},\sigma_{t})$, we use
$\mu_{s}\in\\{3,4\\},\sigma_{s}=0.5$, $\mu_{t}=0$ and
$\sigma_{t}\in\\{1.5,2,2.5\\}$. Also, for Power-law distributions
$p_{s}=\mathrm{PowerLaw}(a_{s})$ and $p_{t}=\mathrm{PowerLaw}(a_{t})$, we use
$a_{s}=3$ and $a_{t}\in\\{0.05,0.1,0.15,0.2,0.25,0.3\\}$.

Figure 8: Relationship between ESS and the parameters $\alpha$ and $\lambda$
for Log-normal distributions. Figure 9: Relationship between ESS and the
parameters $\alpha$ and $\lambda$ for Power-law distributions.

###  Appendix C Details of Numerical Experiments and Additional Experimental
Results

This appendix provides details of the numerical experiments and additional
experimental results.

####  C.1 Implementation Details

#####  C.1.1 Software

The results of all numerical experiments are obtained by our Python 3.11.0
111https://www.python.org/downloads/release/python-3110/ implementation. We
use numpy [20] for random variable generation and matrix calculations, scipy
[62] for computing the density function of probability distributions,
matplotlib 222https://matplotlib.org/ and seaborn
333https://seaborn.pydata.org/ for visualization of the results. In addition,
scikit-learn [46] is used to implement logistic regression, and default
parameters are used unless otherwise noted.

#####  C.1.2 Computing Resources

All our numerical experiments are performed on a machine with 16 GiB of system
memory and 4 vCPUs

#####  C.1.3 Base Model used for Density Ratio Estimation

In our experiments, we use a kernel logistic regression classifier as the
density ratio estimator. A kernel logistic regression classifier employs a
parametric model of the following form for expressing the class-posterior
probability $p(y\mid\bm{x})$.

$\displaystyle
p(y\mid\bm{x};\bm{\theta})=\frac{1}{1+\exp(-y\psi(\bm{x})^{\top}\bm{\theta})},$
(89)

where $\psi(\bm{x})\colon\mathbb{R}^{d}\to\mathbb{R}^{c}$ is a basis function
vector and $\bm{\theta}\in\mathbb{R}^{c}$ is a parameter vector. In our
experiments, we use the following basis functions corresponding to the
polynomial kernel and B-spline kernel in addition to the linear kernel. The
penalized log-likelihood maximization problem reduces the following
minimization.

$\displaystyle\hat{\bm{\theta}}\coloneqq\operatorname*{arg\,min}_{\bm{\theta}\in\mathbb{R}^{c}}\left\\{\sum^{n}_{i=1}\ln\left(1+\exp\left(-y_{i}\psi(\bm{x}_{i})^{\top}\bm{\theta}\right)\right)+\lambda\bm{\theta}^{\top}\bm{\theta}\right\\},$
(90)

where $\lambda\bm{\theta}^{\top}\bm{\theta}$ is a penalty term included for
regularization process. Thus, the loss function $L(y_{i},\bm{\theta})$ is

$\displaystyle
L(y_{i},\bm{\theta})\coloneqq\ln\left(1+\exp\left(-y_{i}\psi(\bm{x}_{i})^{\top}\bm{\theta}\right)\right).$
(91)

From Bayes’ theorem, the density ratio can be expressed as

$\displaystyle r(\bm{x})$ $\displaystyle=\frac{p_{s}(\bm{x})}{p_{t}(\bm{x})}$
$\displaystyle=\left(\frac{p(y=+1\mid\bm{x})p(\bm{x})}{p(y=+1)}\right)\left(\frac{p(y=-1\mid\bm{x})p(\bm{x})}{p(y=-1)}\right)^{-1}$
$\displaystyle=\frac{p(y=-1)}{p(y=+1)}\frac{p(y=+1\mid\bm{x})}{p(y=-1\mid\bm{x})}.$
(92)

The ratio $p(y=-1)/p(y=+1)$ can be approximated by

$\displaystyle\frac{p(y=-1)}{p(y=+1)}\approx\frac{n_{t}/(n_{t}+n_{s})}{n_{s}/(n_{t}+n_{s})}=\frac{n_{t}}{n_{s}}.$
(93)

A density ratio estimator $\hat{r}$ is then given by

$\displaystyle\hat{r}(\bm{x})$
$\displaystyle=\frac{n_{t}}{n_{s}}\frac{1+\exp\left(\psi(\bm{x})^{\top}\hat{\bm{\theta}}\right)}{1+\exp\left(-\psi(\bm{x})^{\top}\hat{\bm{\theta}}\right)}$
$\displaystyle=\frac{n_{t}}{n_{s}}\frac{\exp\left(\psi(\bm{x})^{\top}\hat{\bm{\theta}}\right)\left\\{\exp\left(-\psi(\bm{x})^{\top}\hat{\bm{\theta}}\right)+1\right\\}}{1+\exp\left(-\psi(\bm{x})^{\top}\hat{\bm{\theta}}\right)}$
$\displaystyle=\frac{n_{t}}{n_{s}}\exp\left(\psi(\bm{x})^{\top}\hat{\bm{\theta}}\right).$
(94)

#####  C.1.4 Estimator of Pearson Divergence

The Pearson divergence $D_{\mathrm{PE}}[p\|q]$ is defined as follows.

$\displaystyle
D_{\mathrm{PE}}[p\|q]\coloneqq\frac{1}{2}\int_{\mathcal{X}}\left(\frac{p(\bm{x})}{q(\bm{x})}-1\right)^{2}q(\bm{x})d\bm{x}.$
(95)

Here, we have

$\displaystyle D_{\mathrm{PE}}[p\|q]$
$\displaystyle=\frac{1}{2}\int_{\mathcal{X}}\left(\frac{p(\bm{x})}{q(\bm{x})}-1\right)^{2}q(\bm{x})d\bm{x}$
$\displaystyle=\frac{1}{2}\int_{\mathcal{X}}\left\\{\left(\frac{p(\bm{x})}{q(\bm{x})}\right)^{2}-2\frac{p(\bm{x})}{q(\bm{x})}+1\right\\}q(\bm{x})d\bm{x}$
$\displaystyle=\frac{1}{2}\int_{\mathcal{X}}\frac{p(\bm{x})^{2}}{q(\bm{x})}-2p(\bm{x})+q(\bm{x})d\bm{x}$
$\displaystyle=\frac{1}{2}\int_{\mathcal{X}}\frac{p(x)^{2}}{q(\bm{x})}d\bm{x}-\int_{\mathcal{X}}p(\bm{x})d\bm{x}+\frac{1}{2}\int_{\mathcal{X}}q(\bm{x})d\bm{x}$
$\displaystyle=\frac{1}{2}\int_{\mathcal{X}}\frac{p(\bm{x})}{q(\bm{x})}p(\bm{x})d\bm{x}-\int_{\mathcal{X}}\frac{p(\bm{x})}{q(\bm{x})}q(\bm{x})d\bm{x}+\frac{1}{2}$
$\displaystyle=\frac{1}{2}\int_{\mathcal{X}}r(\bm{x})p(\bm{x})d\bm{x}-\int_{\mathcal{X}}r(\bm{x})q(\bm{x})d\bm{x}+\frac{1}{2}=\frac{1}{2}\mathbb{E}_{p(\bm{x})}\left[r(\bm{x})\right]-\mathbb{E}_{q(\bm{x})}\left[r(\bm{x})\right]+\frac{1}{2}.$

The estimator of the Pearson divergence $\hat{D}_{\mathrm{PE}}[X^{s}\|X^{t}]$
is

$\displaystyle\hat{D}_{\mathrm{PE}}\left[X^{s}\|X^{t}\right]\coloneqq\frac{1}{2n_{s}}\sum^{n_{s}}_{i=1}\hat{r}(\bm{x}^{s}_{i})-\frac{1}{n_{t}}\sum^{n_{t}}_{i=1}\hat{r}(\bm{x}^{t}_{i})+\frac{1}{2}.$
(96)

Combining with the estimator by Eq. (94),

$\displaystyle\hat{D}_{\mathrm{PE}}\left[X^{s}\|X^{t}\right]$
$\displaystyle=\frac{1}{2n_{s}}\sum^{n_{s}}_{i=1}\frac{n_{t}}{n_{s}}\frac{1+\exp\left(\psi(\bm{x}^{s}_{i})^{\top}\hat{\bm{\theta}}\right)}{1+\exp\left(\\-\psi(\bm{x}^{s}_{i})^{\top}\hat{\bm{\theta}}\right)}-\frac{1}{n_{t}}\sum^{n_{t}}_{i=1}\frac{n_{t}}{n_{s}}\frac{1+\exp\left(\psi(\bm{x}^{t}_{i})^{\top}\hat{\bm{\theta}}\right)}{1+\exp\left(\\-\psi(\bm{x}^{t}_{i})^{\top}\hat{\bm{\theta}}\right)}+\frac{1}{2}$
$\displaystyle=\frac{n_{t}}{2n^{2}_{s}}\sum^{n_{s}}_{i=1}\frac{1+\exp\left(\psi(\bm{x}^{s}_{i})^{\top}\hat{\bm{\theta}}\right)}{1+\exp\left(\\-\psi(\bm{x}^{s}_{i})^{\top}\hat{\bm{\theta}}\right)}-\frac{1}{n_{s}}\sum^{n_{t}}_{i=1}\frac{1+\exp\left(\psi(\bm{x}^{t}_{i})^{\top}\hat{\bm{\theta}}\right)}{1+\exp\left(\\-\psi(\bm{x}^{t}_{i})^{\top}\hat{\bm{\theta}}\right)}+\frac{1}{2}.$
(97)

#####  C.1.5 Analytic Calculation of Divergence

In our numerical experiments, we compute KL-divergence between two Gaussian
distributions $p(x)=\mathcal{N}(\mu_{1},\sigma_{1})$ and
$q(x)=\mathcal{N}(\mu_{2},\sigma_{2})$ analytically. The calculation is
performed as follows.

$\displaystyle D_{\mathrm{KL}}[p\|q]$
$\displaystyle=\int_{\mathcal{X}}\left\\{\ln p(x)-\ln q(x)\right\\}p(x)dx$
$\displaystyle=\int_{\mathcal{X}}\left\\{-\ln\sigma_{1}-\frac{1}{2}\left(\frac{x-\mu_{1}}{\sigma_{1}}\right)^{2}+\ln\sigma_{2}+\frac{1}{2}\left(\frac{x-\mu_{s}}{\sigma^{2}}\right)^{2}\right\\}\frac{1}{\sqrt{2\pi\sigma}}e^{-\frac{1}{2}\left(\frac{x-\mu_{1}}{\sigma_{1}}\right)^{2}}$
$\displaystyle=\int_{\mathcal{X}}\left\\{\ln\frac{\sigma_{2}}{\sigma_{1}}+\frac{1}{2}\left(\left(\frac{x-\mu_{2}}{\sigma_{2}}\right)^{2}-\left(\frac{x-\mu_{1}}{\sigma_{1}}\right)^{2}\right)\right\\}\frac{1}{\sqrt{2\pi\sigma}}e^{-\frac{1}{2}\left(\frac{x-\mu_{1}}{\sigma_{1}}\right)^{2}}dx$
$\displaystyle=\mathbb{E}_{p(x)}\left[\ln\frac{\sigma_{2}}{\sigma_{1}}+\frac{1}{2}\left(\left(\frac{x-\mu_{2}}{\sigma_{2}}\right)^{2}-\left(\frac{x-\mu_{1}}{\sigma_{1}}\right)^{2}\right)\right]$
$\displaystyle=\ln\frac{\sigma_{2}}{\sigma_{1}}+\frac{1}{2\sigma_{2}^{2}}\mathbb{E}_{p(x)}\left[(x-\mu_{2})^{2}\right]-\frac{1}{2\sigma_{1}^{2}}\mathbb{E}_{p(x)}\left[(x-\mu_{1})^{2}\right]$
$\displaystyle=\ln\frac{\sigma_{2}}{\sigma_{1}}+\frac{1}{2\sigma_{2}^{2}}\mathbb{E}_{p(x)}\left[(x-\mu_{2})^{2}\right]-\frac{1}{2}$
$\displaystyle=\ln\frac{\sigma_{2}}{\sigma_{1}}+\frac{1}{2\sigma_{2}^{2}}\left\\{\mathbb{E}_{p(x)}\left[(x-\mu_{1})^{2}\right]+2(\mu_{1}-\mu_{2})\mathbb{E}\left[(x-\mu_{1})\right]+(\mu_{1}-\mu_{2})^{2}\right\\}-\frac{1}{2}$
$\displaystyle=\ln\frac{\sigma_{2}}{\sigma_{1}}+\frac{\sigma_{1}^{2}+(\mu_{1}-\mu_{2})^{2}}{2\sigma_{2}^{2}}-\frac{1}{2}.$
(98)

Table 4: Evaluation results with different $\alpha$ and sample dimensions $d$
with polynomial kernel for Gaussian distributions. DRE refers the density
ratio estimation without incremental mixtures.

| $d=2$ | $d=3$ | $d=4$ | $d=5$ |
---|---|---|---|---|---
$DRE$ | $35.98(\pm 27.33)$ | $1335.85(\pm 1447.51)$ | $5038.17(\pm 2666.30)$ | $48586.57(\pm 19828.46)$ |
$\alpha=-1$ | $30.61(\pm 31.03)$ | $715.52(\pm 143)$ | $723.75(\pm 320.86)$ | $1194.64(\pm 400.47)$ |
$\alpha=3$ | $8.27(\pm 9.55)$ | $63.14(\pm 80.29)$ | $296.16(\pm 339.05)$ | $480.41(\pm 89.07)$ |
$\alpha=7$ | $6.09(\pm 0.02)$ | $19.81(\pm 29.53)$ | $42.52(\pm 126.12)$ | $412.99(\pm 71)$ |

Table 5: Evaluation results with different $\alpha$ and sample dimensions $d$
with spline kernel for Gaussian distributions. DRE refers the density ratio
estimation without incremental mixtures.

| $d=2$ | $d=3$ | $d=4$ | $d=5$ |
---|---|---|---|---|---
$DRE$ | $38.13(\pm 4.04)$ | $485.40(\pm 71.28)$ | $5167.16(\pm 1710.77)$ | $48656.67(\pm 14837.67)$ |
$\alpha=-1$ | $6.09(\pm 0.01)$ | $123.38(\pm 122)$ | $489.74(\pm 409.68)$ | $2767.53(\pm 3596.53)$ |
$\alpha=3$ | $6.10(\pm 0.00)$ | $7.69(\pm 0.01)$ | $22.59(\pm 0.01)$ | $119.83(\pm 0.02)$ |
$\alpha=7$ | $6.10(\pm 0.00)$ | $7.70(\pm 0.01)$ | $22.60(\pm 0.01)$ | $119.85(\pm 0.01)$ |

Figure 10: Evaluation results with different $\alpha$ and sample dimensions
$d$ with polynomial kernel for Gaussian distributions. Figure 11: Evaluation
results with different $\alpha$ and sample dimensions $d$ with spline kernel
for Gaussian distributions.

####  C.2 Additional Experimental Results

In additional experiments, we report results based on density ratio estimation
by nonlinear kernel logistic regression. We use the following two-dimensional
polynomial and cubic spline kernels.

$\displaystyle
K_{\mathrm{polynomial}}(\bm{u},\bm{v})\coloneqq(\bm{u}^{\top}\bm{v}+c)^{2},$
$\displaystyle
K_{\mathrm{spline}}(\bm{u},\bm{v})\coloneqq\|\bm{u}-\bm{v}\|^{3}.$ (99)

Tables 4, 5 and Figures 10, 11 show the results of density ratio estimation
based on nonlinear kernel logistic regression. The settings in these
experiments are the same as in the main manuscript. These results show that
density ratio estimation along geodesics with large $\alpha$ reduces
variances, even when nonlinear kernels are used.

###  Appendix D Related Work

####  D.1 Applications of Density Ratio

The density ratio of the two probability distributions has many known
applications. Shimodaira [50] shows that in supervised learning under the
covariate shift assumption, where the marginal distribution of the input data
differs between training and testing, learning by importance weighting using
density ratios restores consistency. This framework is called Importance
Weighted Empirical Risk Minimization (IWERM), and many variants have been
investigated [8, 31, 32, 39, 42, 53, 63]. Density ratios are also known to be
effective for outlier detection, and many studies have used estimated density
ratios from samples as outlier scores [9, 22, 25, 38]. Another useful
application of the density ratio is the two-sample test, which considers
whether, given two sample pairs, they are generated from the same distribution
[10, 27, 29, 30, 41]. Further applications of density ratios include variable
selection [45], dimensionality reduction [59], causal inference [43], and
estimation of mutual information [6, 60]. Many of them are also discussed
theoretically in the use of density ratios.

####  D.2 Methods of Density Ratio Estimation

Because of this wide range of applications, density ratio estimation has
become one of the most important tasks. The simplest idea is to estimate the
density ratio based on a separate density estimation, but these two-step
approaches are known to be computationally unstable [57]. A more promising
idea is a framework for direct density ratio estimation. One idea is the
moment matching approach which tries to match the moments of $p_{s}(\bm{x})$
and $p_{t}(\bm{x})$ by considering $\hat{r}(\bm{x})$ as a transformation
function [23]. Another well-known approach to density ratio estimation is to
optimize the divergence between one distribution $p_{s}(\bm{x})$ and the
transformed other $\hat{r}(\bm{x})p_{t}(\bm{x})$ [58, 54]. Recently,
algorithms for density ratio estimation based on generative models have also
been proposed [11].