Geometric Estimation of Multivariate Dependency

TL;DR

This paper introduces a scalable geometric estimator for dependency between multivariate samples based on minimal spanning trees, which converges to a new measure called geometric mutual information, avoiding density estimation.

Contribution

It proposes a novel geometric dependency estimator using MSTs that converges to GMI, a divergence measure equivalent to Henze-Penrose divergence, with proven asymptotic properties.

Findings

Estimator converges to geometric mutual information (GMI).

Method is scalable and does not require density estimation.

Experimental results demonstrate advantages over existing methods.

Abstract

This paper proposes a geometric estimator of dependency between a pair of multivariate samples. The proposed estimator of dependency is based on a randomly permuted geometric graph (the minimal spanning tree) over the two multivariate samples. This estimator converges to a quantity that we call the geometric mutual information (GMI), which is equivalent to the Henze-Penrose divergence [1] between the joint distribution of the multivariate samples and the product of the marginals. The GMI has many of the same properties as standard MI but can be estimated from empirical data without density estimation; making it scalable to large datasets. The proposed empirical estimator of GMI is simple to implement, involving the construction of an MST spanning over both the original data and a randomly permuted version of this data. We establish asymptotic convergence of the estimator and convergence rates of the bias and variance for smooth multivariate density functions belonging to a H\"{o}lder class. We demonstrate the advantages of our proposed geometric dependency estimator in a series of experiments.

Figures

Tables

Table 1. Table 1: Comparison between different scenarios of various dimensions and sample sizes in terms of parameter α 𝛼 \alpha . We applied the FR approach to estimate the GMI ( I α subscript 𝐼 𝛼 I_{\alpha} ) with α = 0.2 , 0.5 , 0.8 𝛼 0.2 0.5 0.8 \alpha=0.2,0.5,0.8 . We varied dimension d = 6 , 8 , 10 𝑑 6 8 10 d=6,8,10 and sample size n = 1000 , 1500 , 2000 𝑛 1000 1500 2000 n=1000,1500,2000 in each scenario. We observe that for α = { 0.2 , 0.5 , 0.8 } 𝛼 0.2 0.5 0.8 \alpha=\{0.2,0.5,0.8\} , the GMI estimator provides less MSE when α = 0.5 𝛼 0.5 \alpha=0.5 .

Overview table for different $d$, $n$, and $\alpha$
Experiments	Dimension ($d$)	Sample size ($n$)	GMI ($I_{\alpha }$)	MSE ($\times {10}^{-4}$)	Parameter ($\alpha )$
Scenario 1-1	6	1000	0.0229	12	0.2
Scenario 1-2	6	1000	0.0143	4.7944	0.5*
Scenario 1-3	6	1000	0.0176	6.3867	0.8
Scenario 2-1	8	1500	0.0246	11	0.2
Scenario 2-2	8	1500	0.0074	1.6053	0.5*
Scenario 2-3	8	1500	0.0137	5.3863	0.8
Scenario 3-1	10	2000	0.0074	2.3604	0.2
Scenario 3-2	10	2000	0.0029	0.54180	0.5*
Scenario 3-3	10	2000	0.0262	11	0.8

Equations

$$H_0:f_{XY}=f_X{f}_Y,versus{H}_1:f_{XY}\ne f_X{f}_Y.$$

$$D_{p}(f,g)=\displaystyle\frac{1}{4pq}\left[\int\frac{\big{(}pf(\mathbf{x})-qg(\mathbf{x})\big{)}^{2}}{pf(\mathbf{x})+qg(\mathbf{x})}\;{\rm d}\mathbf{x}-(p-q)^{2}\right],$$

$$\begin{array}[]{l}I_{p}(\mathbf{X};\mathbf{Y})=D_{p}(f_{XY},f_{X}f_{Y})\\ \\ \qquad\qquad=\displaystyle\frac{1}{4pq}\left[\iint\frac{\big{(}pf_{XY}(\mathbf{x},\mathbf{y})-qf_{X}(\mathbf{y})f_{Y}(\mathbf{y})\big{)}^{2}}{pf_{XY}(\mathbf{x},\mathbf{y})+qf_{X}(\mathbf{x})f_{Y}(\mathbf{y})}\;{\rm d}\mathbf{x}\;{\rm d}\mathbf{y}-(p-q)^{2}\right].\end{array}$$

$$I_p(\mathbf{X};\mathbf{Y})=1-A_p(\mathbf{X};\mathbf{Y})=\frac{u_p(\mathbf{X};\mathbf{Y})}{4pq}-\frac{(p-q{)}^2}{4pq},$$

$$\begin{array}[]{cl}A_{p}(\mathbf{X};\mathbf{Y})=\displaystyle\iint\frac{f_{XY}(\mathbf{x},\mathbf{y})f_{X}(\mathbf{x})f_{Y}(\mathbf{y})}{pf_{XY}(\mathbf{x},\mathbf{y})+qf_{X}(\mathbf{x})f_{Y}(\mathbf{y})}\;{\rm d}\mathbf{x}\;{\rm d}\mathbf{y}=\displaystyle{\mathbb{E}}_{XY}\bigg{[}\Big{(}p\;\frac{f_{XY}(\mathbf{X},\mathbf{Y})}{f_{X}(\mathbf{X})f_{Y}(\mathbf{Y})}+q\Big{)}^{-1}\bigg{]},\;\;\hbox{and}\\ \\ \displaystyle u_{p}(\mathbf{X};\mathbf{Y})=1-4pq\;A_{p}(\mathbf{X};\mathbf{Y}).\end{array}$$

$$\tilde{I}_{p}\big{(}\lambda_{1}f_{Y|X}g_{X}+\lambda_{2}f_{Y|X}h_{X}\big{)}\geq\lambda_{1}\tilde{I}_{p}(f_{Y|X}g_{X})+\lambda_{2}\tilde{I}_{p}(f_{Y|X}h_{X}).$$

$$\tilde{I}_{p}\big{(}\lambda_{1}g_{Y|X}f_{X}+\lambda_{2}h_{Y|X}f_{X}\big{)}\leq\lambda_{1}\tilde{I}_{p}(g_{Y|X}f_{X})+\lambda_{2}\tilde{I}_{p}(h_{Y|X}f_{X}).$$

$$\displaystyle\begin{array}[]{l}I_{p}(\mathbf{X};\mathbf{Y}|\mathbf{Z})=\displaystyle{\mathbb{E}}_{\mathbf{Z}}\Big{[}I_{p}(\mathbf{X};\mathbf{Y}|\mathbf{Z}=\mathbf{z})\Big{]},\;\;\;\hbox{where}\\ \\ I_{p}(\mathbf{X};\mathbf{Y}|\mathbf{Z}=\mathbf{z})=\quad 1-\displaystyle\iint\frac{f_{XY|Z}(\mathbf{x},\mathbf{y}|\mathbf{z})f_{X|Z}(\mathbf{x}|\mathbf{z})f_{Y|Z}(\mathbf{y}|\mathbf{z})}{p\;f_{XY|Z}(\mathbf{x},\mathbf{y}|\mathbf{z})+q\;f_{X|Z}(\mathbf{x}|\mathbf{z})f_{Y|Z}(\mathbf{y}|\mathbf{z})}\;{\rm d}\mathbf{x}\;{\rm d}\mathbf{y}.\end{array}$$

$$\delta_{X_{2},Y|X_{1}}=\displaystyle{\mathbb{E}}_{f}\left[\bigg{(}p\;\frac{f_{X_{2}Y|X_{1}}(x_{2},y|x_{1})}{f_{X_{2}|X_{1}}(x_{2}|x_{1})\;f_{Y|X_{1}}(y|x_{1})}+q\bigg{)}^{-1}\right].$$

$$I_p(X_1,X_2;Y)\ge I_p(X_1;Y)-{\delta }_{X_2,Y|X_1}.$$

$$\delta_{X_{i},Y|\mathbf{X}^{i-1}}=\displaystyle{\mathbb{E}}_{f}\bigg{(}p\;\frac{f_{X_{i}Y|X^{i-1}}(x_{i},y|x^{i-1})}{f_{X_{i}|X^{i-1}}(x_{2}|x^{i-1})\;f_{Y|X^{i-1}}(y|x^{i-1})}+q\bigg{)}^{-1},$$

$$I_p(\mathbf{X};Y)\ge I_p(X_1;Y)-\sum _{i=1}^{d-1}{\delta }_{X_i,Y|{\mathbf{X}}^{i-1}}.$$

$$\displaystyle I_{p}(\mathbf{Y};\mathbf{X})\geq I_{p}(\mathbf{Z};\mathbf{X})-\Big{(}p\;{\mathbb{E}}_{XY}\big{[}\delta_{X,Y}\big{]}+q\Big{)}^{-1},$$

$${\delta }_{X,Y}=\int \frac{f_{X|Y}(\mathbf{x}|\mathbf{y})f_{Z|Y}(\mathbf{z}|\mathbf{y})}{f_{X|Z}(\mathbf{x}|\mathbf{z})}\mathrm{d}\mathbf{z}.$$

$${\mathbb{E}}_Z\left[\frac{f(\mathbf{Z}|\mathbf{Y})}{f(\mathbf{Z}|\mathbf{X})}\right]=\infty ,$$

$${\mathbb{E}}_{XY}\big{[}\delta_{X,Y}\big{]}={\mathbb{E}}_{XY}\left(\displaystyle\frac{f(\mathbf{X}|\mathbf{Y})}{f(\mathbf{X})}{\mathbb{E}}_{Z}\left[\frac{f(\mathbf{Z}|\mathbf{Y})}{f(\mathbf{Z}|\mathbf{X})}\right]\right).$$

$$\big{\{}(\mathbf{x}_{i_{k}},\mathbf{y}_{j_{k}}),\;k=1,\dots,n^{\prime\prime},\;\hbox{selected at random from}\;{\mathcal{Z}}^{\prime\prime}_{n^{\prime\prime}}\big{\}}:$$

$$\displaystyle\left(\begin{array}[]{c}\leavevmode\hbox to15.18pt{\vbox to12.78pt{\pgfpicture\makeatletter\hbox{\hskip 7.59064pt\lower-6.388pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-4.25763pt}{-1.25055pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$x_{1}$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\\ \leavevmode\hbox to15.18pt{\vbox to12.78pt{\pgfpicture\makeatletter\hbox{\hskip 7.59064pt\lower-6.388pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-4.25763pt}{-1.25055pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$x_{2}$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\\ \leavevmode\hbox to14.17pt{\vbox to15.67pt{\pgfpicture\makeatletter\hbox{\hskip 7.08301pt\lower-7.83301pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-3.75pt}{-2.5pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$\vdots$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\\ \leavevmode\hbox to17.5pt{\vbox to12.62pt{\pgfpicture\makeatletter\hbox{\hskip 8.7513pt\lower-6.31078pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-5.41829pt}{-1.32777pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$x_{n^{\prime\prime}}$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\end{array}\right)\qquad\left(\begin{array}[]{c}\leavevmode\hbox to14.73pt{\vbox to12.92pt{\pgfpicture\makeatletter\hbox{\hskip 7.3638pt\lower-6.458pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-4.03079pt}{-1.18056pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$y_{1}$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\\ \leavevmode\hbox to14.73pt{\vbox to12.92pt{\pgfpicture\makeatletter\hbox{\hskip 7.3638pt\lower-6.458pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-4.03079pt}{-1.18056pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$y_{2}$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\\ \leavevmode\hbox to14.17pt{\vbox to15.67pt{\pgfpicture\makeatletter\hbox{\hskip 7.08301pt\lower-7.83301pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-3.75pt}{-2.5pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$\vdots$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\\ \leavevmode\hbox to17.05pt{\vbox to12.92pt{\pgfpicture\makeatletter\hbox{\hskip 8.52446pt\lower-6.458pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-5.19145pt}{-1.18056pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$y_{n^{\prime\prime}}$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\end{array}\right)$$

$$1-{\mathfrak{R}}_{n^{\prime },n^{\prime \prime }}({\mathcal{Z}}_{n^{\prime }}^{\prime },{\tilde{\mathcal{Z}}}_{n^{\prime \prime }})\frac{n}{2n^{\prime }{n}^{\prime \prime }}\to I_{\alpha }(\mathbf{X};\mathbf{Y}),a.s.$$

$$\displaystyle\begin{array}[]{l}\Big{\{}g:\big{\|}g(\mathbf{z})-p_{\mathbf{x}}^{\lfloor\eta\rfloor}(\mathbf{z})\big{\|}_{d}\leq K\;g(\mathbf{x})\;\big{\|}\mathbf{x}-\mathbf{z}\big{\|}_{d}^{\eta},\;\;\mathbf{x},\;\mathbf{z}\in\mathcal{X}\Big{\}},\end{array}$$

$$\Big{\{}g:\big{\|}g(\mathbf{z})-p_{\mathbf{x}}^{\lfloor\eta\rfloor}(\mathbf{z})\big{\|}_{d}\leq K^{\prime}\;\big{\|}\mathbf{x}-\mathbf{z}\big{\|}_{d}^{\eta},\;\;\mathbf{x},\;\mathbf{z}\in\mathcal{X}\Big{\}}.$$

$$\begin{array}[]{l}\displaystyle\bigg{|}\frac{{\mathbb{E}}\big{[}\mathfrak{R}_{n^{\prime},n^{\prime\prime}}\big{]}}{n}-2\alpha\beta\iint\frac{f_{\mathbf{X}\mathbf{Y}}(\mathbf{x},\mathbf{y})f_{\mathbf{X}}(\mathbf{x})f_{\mathbf{Y}}(\mathbf{y})}{\alpha f_{\mathbf{X}\mathbf{Y}}(\mathbf{x},\mathbf{y})+\beta f_{\mathbf{X}}(\mathbf{x})f_{\mathbf{Y}}(\mathbf{y})}\;{\rm d}\mathbf{x}{\rm d}\mathbf{y}\bigg{|}\\ \\ \qquad\quad\leq O\bigg{(}\max\Big{\{}n^{-\eta^{2}\big{/}(d(1+\eta))},\;(\beta n)^{-\eta/(1+\eta)},\;c_{d}2^{d}n^{-1}\Big{\}}\bigg{)},\end{array}$$

$$H(\gamma )=-\gamma \log \gamma -(1-\gamma )\log (1-\gamma ).$$

$$Var\left(\frac{{\mathfrak{R}}_{n^{\prime },n^{\prime \prime }}}{n}\right)\le \frac{(1-\alpha )c_d}{n},\alpha =n^{\prime }/n,$$

$$\frac{({ϵ}_{XY}(\mathbf{x},\mathbf{y})+1/(\beta C_n))(1+{ϵ}_{XY}(\mathbf{xy})+1/(\beta C_n))}{(\alpha +\beta {ϵ}_{XY}(\mathbf{x},\mathbf{y}){)}^2},\beta =1-\alpha .$$

$$\underset{\alpha }{\min }\underset{{ϵ}_{XY}}{\max }$$

$$C_{ϵ}^L\le {ϵ}_{XY}\le C_{ϵ}^U,$$

$${\alpha }_0^L\le \alpha \le {\alpha }_0^U,$$

$$\tilde{\Delta }(\alpha ,{ϵ}_{XY}):=D(n,l_n,d,\eta )+\overline{D}(n,l_n,d)C_{XY}^U{\iint }_{{\mathbb{S}}_{XY}}{\tilde{G}}_{{ϵ}_{XY},n}^{\alpha ,\beta }(\mathbf{x},\mathbf{y})\mathrm{d}\mathbf{x}\mathrm{d}\mathbf{y},$$

$$\displaystyle\begin{array}[]{l}D(n,l_{n},d,\eta)=\displaystyle c_{2}l_{n}^{d}n^{-1}+c_{d}2^{d}n^{-1}+c^{\prime}l_{n}^{d}n^{-\eta/d}+\displaystyle cl_{n}^{d}n^{-1/d}+2c_{1}l_{n}^{d-1}n^{1/d-1},\end{array}$$

$$\begin{array}[]{l}\widetilde{D}(n,l_{n},d)=\displaystyle 2+n^{-1}2c^{\prime\prime}\sum\limits_{i=1}^{M}l_{n}\;l_{n}^{d}a_{i}^{-1}+n^{-3/2}2c^{\prime}_{1}\sum\limits_{i=1}^{M}l_{n}\;l_{n}^{d/2}\sqrt{b_{i}}{a_{i}^{2}}\\ \qquad+n^{-1}\displaystyle\sum\limits_{i=1}^{M}2n^{-3/2}l_{n}^{-d/2}\displaystyle\frac{\sqrt{b_{i}}}{a_{i}^{2}}\big{(}na_{i}l_{n}^{d}+n^{2}a_{i}^{2}\big{)}^{1/2}\big{(}nb_{i}l_{n}^{d}+n^{2}b_{i}^{2}\big{)}^{1/2}.\end{array}$$

Full text

\pubvolume

xx \issuenum1 \articlenumber5

\TitleGeometric Estimation of Multivariate Dependency \AuthorSalimeh Yasaei Sekeh and Alfred O. Hero \AuthorNamesSalimeh Yasaei Sekeh and Alfred O. Hero

1 Introduction

Estimation of multivariate dependency has many applications in fields such as information theory, clustering, structure learning, data processing, feature selection, time series prediction and reinforcement learning, see Lewi et al. (2006) and Peng et al. (2002), Moon et al. (2017), Brillinger (2004), Torkkola (2003); Vergara and A. (2014); Peng et al. (2005, 1997), Sorjamaa et al. (2005), and Mohamed and Rezende (2015), respectively. It is difficult to accurately estimate the mutual information in high dimensional settings where the data is multivariate with a Lebesgue continuous density—the setting considered in this paper. An important and ubiquitous measure of dependency is the Shannon Mutual Information (MI), which has seen extensive use across many application domains. However, the estimation of mutual information can often be challenging. In this paper we focus on a measure of MI that we call the Geometric MI (GMI). This MI measure is defined as the asymptotic large sample limit of a randomized minimal spanning tree (MST) statistic spanning the multivariate sample realizations. The GMI is related to a divergence measure called the Henze-Penrose divergence Neemuchwala and Hero (2005), Neemuchwala et al. (2006), and related to the multivariate runs test Henze and Penrose (1999). In Berisha and Hero (2015), Berisha et al. (2016), it was shown that this divergence measure can be used to specify a tighter bound for the Bayes error rate for testing if a random sample comes from one of two distributions as compared to previous divergence-type bounds such as the Bhattacharrya bound. Furthermore, the authors of Yasaei Sekeh et al. (2018) proposed a nonparametric bound on multi-class classification Bayes error rate using a global MST graph.

Let $\mathbf{X}$ and $\mathbf{Y}$ be random variables with unknown joint density $f_{{X}{Y}}$ and marginal densities $f_{{X}}$ and $f_{{Y}}$, respectively, and consider two hypotheses: $H_{0}$, $\mathbf{X}$ and $\mathbf{Y}$ are independent and $H_{1}$, $\mathbf{X}$ and $\mathbf{Y}$ are dependent,

[TABLE]

The GMI is defined as the Henze-Penrose divergence between $f_{XY}$ and $f_{X}f_{Y}$ which can be used as a dependency measure. In this paper we prove that for large sample size the randomized minimal spanning tree (MST) statistic spanning the original multivariate sample realizations and randomly shuffled data set in which for given the first element in each pair the second element is replaced by a randomly selected from all second elements, convergences to the GMI measure. A direct implication of Berisha and Hero (2015), Berisha et al. (2016) is that the GMI provides a tighter bound on the Bayes misclassification rate for the independence test of $H_{0}$ vs $H_{1}$. In this paper we propose an estimator based on a random permutation modification of the Friedman-Rafsky multivariate test statistic and show that under certain conditions the GMI estimator achieves the parametric mean square error (MSE) rate when the joint and marginal densities are bounded and smooth. Importantly our proposed GMI estimator does not require explicit estimation of the joint and marginal densities.

In addition to achieving the optimal theoretical MSE rate, computational complexity is an important challenge addressed by researchers in machine learning and data science. Most plug-in based estimators like kernel density estimator (KDE) or K-nearest neighbour (KNN) estimator require runtime complexity of $O(n^{2})$ which is not suitable for large scale applications. Noshad et al. proposed a graph theoretic direct estimation method based on nearest neighbor ratios (NNR) Noshad et al. (2017). The NNR estimator is based on $k$-NN graph and computationally more tractable than other competing estimators with complexity $O(kn\log n)$. The construction of the minimal spanning tree lies at the heart of the GMI estimator proposed in this paper. Since the GMI estimator is based on the Euclidean MST the dual-tree algorithm by March et al. March et al. (2010) can be applied. This algorithm is based on the construction of Borůvka Borůvka (1926) and implements the Euclidean MST in approximately $O(nlogn)$ time which is faster than the previous estimators. In this paper we experimentally show that for large number of sample sizes the proposed GMI estimator is less computationally complex than the KDE plug-in method.

1.1 Related work

Estimation of mutual information is related to estimation of the information histogram. The most common estimators of MI are based on plug-in density estimation, e.g., using kernel density or kNN density estimators Kraskov et al. (2004), Moon et al. (2017). Motivated by ensemble methods applied to divergence estimation Moon and Hero (2014)-Moon et al. (2016), in Moon et al. (2017) an ensemble method for combining multiple KDE bandwidths was proposed for estimating MI. Under certain smoothness conditions this ensemble MI estimator was shown to achieve parametric convergence rates.

Another class of estimators of multivariate dependency bypasses the difficult density estimation task. This class includes the statistically consistent estimator of Rényi-$\alpha$ and KL mutual information which is based on the asymptotic of the length of the KNN graph, Leonenko et al. (2008) and Gao et al. (2015). The estimator of Pál et al. (2010) builds on the KNN methods for Rényi entropy estimation. As MI increases, the dependencies between random variables increase which results in less smooth densities. In Gao et al. (2015), it has been shown that for large MI the KNN or KDE approaches are ill-suited candidates for estimating MI since the assumption of local uniformity can be violated when they are strong dependencies. To overcome this issue an assumption on the smoothness of the density is required, see Krishnamurthy et al. (2014), Kandasamy et al. (2015), and Moon and Hero (2014)-Moon et al. (2016). For all of these methods, the optimal parametric rate of MSE convergence is achieved when the densities are either $d$, $(d+1)/2$ or $d/2$ times differentiable Singh and Póczos (2016). In this paper, we assume that joint and marginal densities are smooth in the sense that they belong to a strong Hölder continuous classes of densities $\Sigma^{s}_{d}(\eta,K)$, such that the smoothness parameter $\eta$ ranges in $(0,1]$.

A MI measure based on the Pearson divergence was considered in Sugiyama (2012) that is computational efficient and numerically stable. The authors of Costa and Hero (2004) and Pál et al. (2010) used minimal spanning tree generalized nearest-neighbor graph approaches, respectively, to estimate Rényi mutual information. In Moon et al. (2017), a nonparametric mutual information estimator was proposed using a weighted ensemble method with $O(1/n)$ parametric convergence rate. This estimator was based on plug-in density estimation, which is challenging in high dimension.

Our proposed estimator differs from previous methods in the following ways. First, it estimates a different measure of mutual information, the GMI. Second, instead of using the KNN graph the estimator of GMI uses a randomized minimal spanning tree that spans the multivariate realizations. The proposed GMI estimator is motivated by the multivariate runs test of Friedman and Rafsky (FR) Friedman and Rafsky (1979) which is a multivariate generalization of the univariate Smirnov maximum deviation test Smirnov (1939) and the Wald-Wolfowitz Wald and Wolfowitz (1940) runs test in one dimension. We also emphasize that our proposed GMI estimator does not require boundary correction, in contrast to other graph-based estimators, such as, the NNR estimator Noshad et al. (2017), scalable MI estimator Noshad and Hero (2018), or cross match statistic Yasaei Sekeh et al. (2018).

1.2 Contribution

The contribution of this paper has three components

• (1)

We propose a novel non-parametric multivariate dependency measure, referred to as Geometric Mutual Information (GMI), that is based on graph-based divergence estimation. The geometric mutual information is constructed using a minimal spanning tree and is a function of Friedman-Rafsky multivariate test statistic.

• (2)

We establish properties of the proposed dependency measure analogous to those of standard mutual information, such as, convexity, concavity, chain rule, and the data processing inequality.

• (3)

We derive a bound on the MSE rate for the proposed geometric estimator. An advantage of the estimator is that it achieves the optimal MSE rate without explicitly performing boundary correction, which is required for most plug-in estimators.

1.3 Organization

The paper is organized as follows. In Section 2, we define the geometric mutual information and establish some of its mathematical properties. In Subsection 2.2 and 2.3, we introduce a statistically consistent GMI estimator and derive a bound on its mean square error convergence rate. In Section 3 we verify the theory through experiments.

Throughout the paper, we denote statistical expectation by ${\mathbb{E}}$, the variance by abbreviation ${\rm Var}$. Bold face type indicates random vectors. All densities are assumed to be Lebesgue continuous with no atoms.

2 The geometric mutual information (GMI)

In this section, we first review the definition of Henze-Penrose (HP) divergence measure defined by Berisha and Hero in Berisha and Hero (2015) and Henze and Penrose (1999). The Henze-Penrose divergence between densities $f$ and $g$ with domain ${\mathbb{R}}^{d}$ for parameter $p\in(0,1)$ is defined as follows (see Henze and Penrose (1999), Berisha and Hero (2015), and Berisha et al. (2016)):

[TABLE]

where $q=1-p$. This functional is an $f$-divergence Csiszár (1967), also known as an Ali-Silvey distance Ali and Silvey (1996). That is, it satisfies the properties of non-negativity, monotonicity, and joint convexity Berisha et al. (2016). This measure takes values in $[0,1]$ and $D_{p}=0$ if and only if $f=g$ almost surely.

The geometric mutual information measure is defined as follows. Let $f_{X}$, $f_{Y}$, and $f_{XY}$ be the marginal and joint distributions, respectively, of random vectors $\mathbf{X}\in{\mathbb{R}}^{d_{x}}$, $\mathbf{Y}\in{\mathbb{R}}^{d_{y}}$ where $d_{x}$ and $d_{y}$ are positive integers. Then by using (1), a Henze-Penrose generalization of the mutual information between $\mathbf{X}$ and $\mathbf{Y}$, is defined by

[TABLE]

We will show below that $I_{p}(\mathbf{X};\mathbf{Y})$ has a geometric interpretation in terms of the large sample limit of a minimal spanning tree spanning $n$ sample realizations of $\mathbf{X}\cup\mathbf{Y}$. Thus we call $I_{p}(\mathbf{X};\mathbf{Y})$ the geometric mutual information (GMI) between $\mathbf{X}$ and $\mathbf{Y}$. The GMI satisfies similar properties to other definitions of mutual information, such as Shannon and Rényi mutual information. Recalling (3) in Berisha and Hero (2015), an alternative form of $I_{p}$ is given by

[TABLE]

where

[TABLE]

The affinity $A_{p}(\mathbf{X};\mathbf{Y})$ is called the geometric affinity between $\mathbf{X}$ and $\mathbf{Y}$, and it was originally defined in Henze and Penrose (1999). The next subsection of the paper is dedicated to the basic inequalities and properties of the proposed GMI measure (2).

2.1 Properties of the geometric mutual information

In this subsection we establish basic inequalities and properties of the GMI, $I_{p}$, given in (2).

2.1.1 Basic inequalities

The following theorem shows that $I_{p}(\mathbf{X};\mathbf{Y})$ is a concave function in $f_{X}$ and a convex function in $f_{Y|X}$. The proof is given in Appendix A, Subsection 5.1.

{Theorem}

Denote by $\tilde{I}_{p}(f_{XY})$ the GMI $I_{p}(\mathbf{X};\mathbf{Y})$ when $\mathbf{X}\in{\mathbb{R}}^{d_{x}}$ and $\mathbf{Y}\in{\mathbb{R}}^{d_{y}}$ have joint density $f_{XY}$. Then the GMI satisfies

• (i) Concavity in $f_{X}$: Let $f_{Y|X}$ be conditional density of $\mathbf{Y}$ given $\mathbf{X}$ and let $g_{X}$ and $h_{X}$ be densities on ${\mathbb{R}}^{d_{x}}$. For $p\in(0,1)$

[TABLE]

The inequality is strict unless either $\lambda_{1}$ or $\lambda_{2}$ are zero or $h_{X}=g_{X}$.

• (ii) Convexity in $f_{Y|X}$: Let $g_{Y|X}$ and $h_{Y|X}$ be conditional densities of $\mathbf{Y}$ given $\mathbf{X}$ and let $f_{X}$ be marginal density. Then for $p\in(0,1)$

[TABLE]

The equality occurs when either $\lambda_{1}$ or $\lambda_{2}$ are zero or $h_{Y|X}=g_{Y|X}$.

The GMI, $I_{p}(\mathbf{X};\mathbf{Y})$, satisfies a property analogous to the standard chain rule and the data processing inequality Cover and Thomas (1991). For random variables $\mathbf{X}\in{\mathbb{R}}^{d_{x}},\;\mathbf{Y}\in{\mathbb{R}}^{d_{y}}$, and $\mathbf{Z}\in{\mathbb{R}}^{d_{z}}$ we define the conditional GMI measure by

[TABLE]

The next theorem establishes a relation between the joint and conditional geometric mutual information.

{Theorem} For given $d$-dimensional random vector $\mathbf{X}$ with components $X_{1},X_{2},\dots,X_{d}$ and random variable $Y$,

• (i)

For $d=2$, set

[TABLE]

Here $f$ is the joint PDF of random vector $(X_{1},X_{2},Y)$ and note that $0\leq\delta_{X_{2},Y|X_{1}}\leq 1$. Then we have

[TABLE]

Furthermore $\delta_{X_{i},Y|\mathbf{X}^{i-1}}\in[0,1]$.

• (ii)

Let $\mathbf{X}^{i}:=X_{1},X_{2},\dots,X_{i}$. Define a general form of $\delta$ in (i) by

[TABLE]

where $f$ is the join PDF of random vector $(X_{1},\dots,X_{d},Y)$. A general form of (9) is given by

[TABLE]

Note that in the special case $\displaystyle\sum_{i=1}^{d-1}\delta_{X_{i},Y|\mathbf{X}^{i-1}}\geq 1$, the inequality (11) is trivial. The proof of Theorem 2.1.1 is given in Appendix B, subsection 5.2. Next we apply Theorem 2.1.1 and state a result for Markov chain random vectors. Likewise to previous proofs, the proof of Proposition 2.1.1 is provided in Appendices section, Appendix C. {Proposition} Suppose random vectors $\mathbf{X},\mathbf{Y},\mathbf{Z}$ form a Markov chain denoted, $\mathbf{X}\rightarrow\mathbf{Y}\rightarrow\mathbf{Z}$, in the sense that $f_{XYZ}=f_{X|Y}f_{Y|Z}f_{Z}$. Then for $p\in(0,1)$

[TABLE]

where

[TABLE]

Further, if both $\mathbf{X}\rightarrow\mathbf{Y}\rightarrow\mathbf{Z}$ and $\mathbf{X}\rightarrow\mathbf{Z}\rightarrow\mathbf{Y}$ together hold true, we have $I_{p}(\mathbf{Y};\mathbf{X})=I_{p}(\mathbf{Z};\mathbf{X})$.

The inequality in (12) becomes a tighter inequality interpretable as a standard data-processing inequality $I_{p}(\mathbf{Y};\mathbf{X})\geq I_{p}(\mathbf{Z};\mathbf{X})$, when

[TABLE]

since

[TABLE]

2.2 The Friedman-Rafsky Estimator

Let a random sample $\{\mathbf{x}_{i},\mathbf{y}_{i}\}_{i=1}^{n}$ from $f_{XY}(\mathbf{x},\mathbf{y})$ be available. Here we show that the GMI $I_{p}(\mathbf{X};\mathbf{Y})$ can be directly estimated without estimating the densities. The estimator is inspired by the MST construction of Friedman and Rafsky (1979) that provides a consistent estimate of the Henze-Penrose divergence Berisha and Hero (2015), Berisha et al. (2016). We denote by $\mathbf{z}_{i}$ the $i$-th joint sample $\mathbf{x}_{i},\mathbf{y}_{i}$ and by $\mathcal{Z}_{n}$ the sample set $\{\mathbf{z}_{i}\}_{i=1}^{n}$. Divide the sample set $\mathcal{Z}_{n}$ into two subsets $\mathcal{Z}^{\prime}_{n^{\prime}}$ and $\mathcal{Z}^{\prime\prime}_{n^{\prime\prime}}$ with the proportion $\alpha=n^{\prime}/n$ and $\beta=n^{\prime\prime}/n$, where $\alpha+\beta=1$.

Denote by $\widetilde{\mathcal{Z}}_{n^{\prime\prime}}$ the set

[TABLE]

[TABLE]

This means that for each $\mathbf{z}_{ik}=(\mathbf{x}_{ik},\mathbf{y}_{ik})\in\mathcal{Z}^{\prime\prime}_{n^{\prime\prime}}$ given the first element $\mathbf{x}_{ik}$ the second element $\mathbf{y}_{ik}$ is replaced by a randomly selected $y\in\{y_{jk}\}_{j=1}^{n^{\prime\prime}}$. This results in a shuffling of the binary relation relating $y_{ik}$ in $y_{jk}$. The estimator of $I_{p}(\mathbf{X};\mathbf{Y})$ is derived based on the Friedman-Rafsky (FR) multivariate runs test statistic Friedman and Rafsky (1979) on the concatenated data set , $\mathcal{Z}^{\prime}_{n^{\prime}}\cup\widetilde{\mathcal{Z}}_{n^{\prime\prime}}$. The FR test statistic is defined as the number of edges in the MST spanning the concatenated data set that connect a point in $\mathcal{Z}^{\prime}_{n^{\prime}}$ to a point in $\widetilde{\mathcal{Z}}_{n^{\prime\prime}}$. This test statistic is denoted by $\mathfrak{R}_{n^{\prime},n^{\prime\prime}}:=\mathfrak{R}_{n^{\prime},n^{\prime\prime}}(\mathcal{Z}^{\prime}_{n^{\prime}},\widetilde{\mathcal{Z}}_{n^{\prime\prime}})$. Note that since the MST is unique with probability one (under the assumption that all density functions are Lebesgue continuous) then all inter point distances between nodes are distinct. This estimator converges to $I_{p}(\mathbf{X};\mathbf{Y})$ almost surely as $n\rightarrow\infty$. The procedure is summarized in Algorithm 1.

Theorem 2.2 shows that the output in Algorithm 1 estimates the GMI with parameter $p={\alpha}$. The proof is provided in Appendix D, Subsection 5.4. {Theorem} For given proportionality parameter $\alpha\in(0,1)$, choose $n^{\prime}$, $n^{\prime\prime}$ such that $n^{\prime}+n^{\prime\prime}=n$ and as $n\rightarrow\infty$, we have $n^{\prime}/n\rightarrow\alpha$ and $n^{\prime\prime}/n\rightarrow\beta=1-\alpha$. Then

[TABLE]

Note that the asymptotic limit in (14) depends on the proportionality parameter $\alpha$. Later in Subsection 2.4, we discuss on choosing an optimal parameter $\tilde{\alpha}$. In Fig. 1, we show a visualization of the MST constructed over merged independent ($\rho=0$) and highly dependent ($\rho=0.9$) data sets drawn from two dimensional standard Normal distribution with correlation coefficient $\rho$. Notice that the MST connects different colored samples, corresponding to independent and dependent samples, respectively, indicated in green. The total number of green edges is the FR test statistic.

2.3 Convergence Rates

In this subsection we provide the MSE convergence rates in the form of upper bounds on the bias and the variance. This MSE bound is given in terms of the sample size $n$, the dimension $d$, and the proportionality parameter $\alpha$. Deriving convergence rates for mutual information estimators has been of interest in information theory and machine learning Pál et al. (2010), Moon et al. (2017). The rates are typically derived in terms of a smoothness condition on the densities, such as the Hölder condition Wolfgang (1991). Here we assume $f_{X}$, $f_{Y}$ and $f_{XY}$ with support sets $\mathbb{S}_{X}$, $\mathbb{S}_{Y}$, and $\mathbb{S}_{XY}:=\mathbb{S}_{X}\times\mathbb{S}_{Y}$, respectively, are smooth in the sense that they belong to strong Hölder continuous classes of densities $\Sigma^{s}_{d}(\eta,$ K), $0<\eta\leq 1$ Lorentz (1966), Wolfgang (1991): {Definition} (Strong Hölder class): Let $\mathcal{X}\subset{\mathbb{R}}^{d}$ be a compact space. The strong Hölder class $\Sigma_{d}^{s}(\eta,K)$, with $\eta$-Hölder parameter and constant $K$, of functions on $L_{d}$-norm, consists of the functions $g$ that satisfy

[TABLE]

where $p_{\mathbf{x}}^{k}(\mathbf{z})$ is the Taylor polynomial (multinomial) of $g$ of order $k$ expanded about the point $\mathbf{x}$ and $\lfloor\eta\rfloor$ is defined as the greatest integer strictly less than $\eta$. Note that for the standard Hölder class the factor $g(\mathbf{x})$ on the RHS of (16) is omitted.

To find the optimal parameter $\alpha$ we require explicit expressions for the bias and variance. Bounds are provided in Appendix E, Subsection 5.5. To obtain such expressions, we require several assumptions on the absolutely continuous densities $f_{X}$, $f_{Y}$, $f_{XY}$ and support sets $\mathbb{S}_{X}$, $\mathbb{S}_{Y}$, $\mathbb{S}_{XY}$:

• (A.1)

Each of the densities belong to $\displaystyle\Sigma^{s}_{d}(\eta,K)$ class in their support sets with smoothness parameters $\eta$ and smoothness constants $K$.

• (A.2)

The volumes of the support sets are finite, i.e. $0<\mathbb{V}(\mathbb{S}_{X})<\infty,\;\;\;0<\mathbb{V}(\mathbb{S}_{Y})<\infty$.

• (A.3)

All densities are bounded on their support sets, i.e. there exist two sets of constants $C^{L}_{X},C^{L}_{Y},C^{L}_{XY}$ and $C^{U}_{X},C^{U}_{Y},C^{U}_{XY}$ such that $0<C^{L}_{X}\leq f_{X}\leq C^{U}_{X}<\infty$, $0<C^{L}_{Y}\leq f_{Y}\leq C^{U}_{Y}<\infty$ and $0<C^{L}_{XY}\leq f_{XY}\leq C^{U}_{XY}<\infty$.

Remarks:

The assumption (A.3) implies that the strong Hölder class, (A.1), is a subset of standard Hölder class:

[TABLE]

• 2.

If in (A.3) $C^{L}_{X},C^{L}_{Y},C^{L}_{XY}$ and $C^{U}_{X},C^{U}_{Y},C^{U}_{XY}$ are equal i.e. the density functions are uniform, $g(\mathbf{x})=K^{\prime\prime}$, and if $K^{\prime}=K.K^{\prime\prime}$ then two classes (16) and (17) are equal.

• 3.

By using the Strong Hölder class (16) our proposed bias bound in Theorem 2.3 below, becomes tighter. Another bound can be derived for the standard Hölder class (17) but the proposed bias bound becomes weaker, and so in the result for bounding convergence rate.

The following theorems on the bias and variance follow under assumptions (A.1), (A.2), and (A.3): {Theorem} Under two assumptions (A.1) and (A.3), for given $\alpha\in(0,1)$, $\beta=1-\alpha$, the bias of the $\mathfrak{R}_{n^{\prime},n^{\prime\prime}}$ estimator for $d\geq 2$, $0<\eta\leq 1$ satisfies

[TABLE]

where $c_{d}$ is a constant depending only on $d$.

Note that according to Theorem 13 in Robins and Salowe (1994), the constant $c_{d}$ is lower bounded by $\Omega\left(\sqrt{d}2^{n(1-H(\gamma))}\right)$, $\gamma=2^{-d}$ and $H(\gamma)$ is the binary entropy i.e.

[TABLE]

A proof of Theorem 2.3 is given in Appendix E. The next theorem gives an upper bound on the variance of the FR estimator $\mathfrak{R}_{n^{\prime},n^{\prime\prime}}$. The proof of the variance result requires a different approach than the bias bound. The proof uses the Efron-Stein inequality Efron and stein (1981) and is similar to arguments in (Yasaei Sekeh et al., 2018, Appendix C), and is omitted. In Theorem 2.3 we assume that the densities $f_{X}$, $f_{Y}$, and $f_{XY}$ are absolutely continuous and bounded. Note that we can obtain a weaker bound than (18) when we relax (A.1) to the standard Hölder class with a constant coefficient $K^{\prime}$ instead of $Kg(\mathbf{x})$. {Theorem} Given $\alpha\in(0,1)$, the variance of the estimator $\mathfrak{R}_{n^{\prime},n^{\prime\prime}}(\mathcal{Z}^{\prime}_{n^{\prime}},\widetilde{\mathcal{Z}}_{n^{\prime\prime}})$ is bounded by

[TABLE]

where $c_{d}$ is a constant depending only on the dimension $d$.

2.4 Minimax Parameter $\alpha$

A question in Algorithm 1 is the existance of a best proportion $\alpha$ that partitions sample. In this subsection we obtain on an upper bound of the MSE rate and characterize the parameter $\alpha$ that minimizes the upper bound. This optimal $\alpha$ depends on maximize densities, however.

Recall assumptions (A.1), (A.2), and (A.3) in Subsection 2.3. The constant $\alpha$ can be chosen to ensure that the MSE converges rate obtained from the bias and variance rates (18) and (2.3) by selecting $\alpha$ to minimize the maximum MSE, where the maximum is taken over the space of Hölder smooth joint densities $f_{XY}$.

Throughout this subsection we use the following notations:

• •

\epsilon_{XY}:=f_{XY}\big{/}f_{X}f_{Y},

• •

C^{L}_{\epsilon}:=C^{L}_{XY}\big{/}C^{U}_{X}C^{U}_{Y} and C^{U}_{\epsilon}:=C^{U}_{XY}\big{/}C^{L}_{X}C^{L}_{Y},

• •

$C_{n}:=C^{L}_{XY}\;n/2$,

• •

$\alpha^{L}_{0}:=\displaystyle\frac{2}{C_{n}}$ and $\alpha^{U}_{0}:=\min\left\{\displaystyle\frac{1}{4},\displaystyle\frac{1+1/C_{n}}{4+2C^{U}_{\epsilon}},1-n^{\eta/d-1}\right\}$, where $\eta$ is the smoothness parameter,

• •

l_{n}:=\big{\lfloor}n^{\eta/(d^{2}(1+\eta))}\big{\rfloor}.

Now define $\widetilde{{G}}^{\alpha,\beta}_{\epsilon_{XY},n}(\mathbf{x},\mathbf{y})$ by

[TABLE]

Consider the following optimization problem:

[TABLE]

where

[TABLE]

and

[TABLE]

[TABLE]

Note that in (24), $c,c^{\prime},c_{1},c_{2}$ are constants, and $c_{d}$ only depends on the dimension $d$. Also, in (25), $a_{i}$ and $b_{i}$ are constants. Let $\epsilon_{XY}^{*}$ be the optimal $\epsilon_{XY}$ i.e. $\epsilon_{XY}^{*}$ be the solution of the optimization problem (21). Set

[TABLE]

such that $\displaystyle\widetilde{\Delta}(\alpha,\epsilon^{*}_{XY})$ is (22) when $\epsilon_{XY}=\epsilon^{*}_{XY}$. For $\alpha\in[\alpha^{L}_{0},\alpha^{U}_{0}$], the optimal choice of $\epsilon_{XY}$ in terms of maximizing the MSE is $\epsilon_{XY}^{*}=C^{U}_{\epsilon}$ and the saddle point for the parameter $\alpha$, denoted by $\tilde{\alpha}$, is given as follows:

• •

$\tilde{\alpha}=\alpha_{0}^{U}$, if $\Xi(\alpha_{0}^{U})<0$.

• •

$\tilde{\alpha}=\alpha_{0}^{L}$, if $\Xi(\alpha_{0}^{L})>0$.

• •

$\tilde{\alpha}=\Xi^{-1}(0)$, if $\alpha_{0}^{L}\leq\Xi^{-1}(0)\leq\alpha_{0}^{U}$.

Further details are given in Appendix F.

3 Simulation Study

In this section numerical simulations are presented that illustrate the theory in Section 2. We perform multiple experiments to demonstrate the utility of the proposed direct estimator of the HP-divergence in terms of the dimension $d$ and the sample size $n$. Our proposed MST-based estimator of the GMI is compared to other plug-in GMI estimators, in particular the standard KDE estimator of Moon et al. (2017), where the convergence rates of Theorem 2.3 and 2.3 are validated. We use multivariate Normal simulated data in the experiments. In this section, we also discuss the choice of the proportionality parameter $\alpha$ and compare runtime of our proposed FR estimator approach with KDE method.

Here we perform four sets of experiments to illustrate the estimator and the theory derived above. For the first set of experiments the MSE of the FR estimator proposed in Algorithm 1 is shown in Fig. 2-(left). The samples were drawn from $d$-dimensional Normal distribution, with various sample sizes and dimensions $d=6,10,12$. We selected the proportionality parameter $\alpha=0.3$ and computed the MSE in terms of the sample size $n$. We show the log-log plot of MSE when $n$ varies in $[100,1500]$. Note that the empirically optimal proportion $\alpha$ depends on $n$, so to avoid the computational complexity we fixed $\alpha$ for this experiment. The experimental result shown in Fig. 2-(left) validates the theoretical MSE growth rates derived from (18) and (19), i.e., decreasing sublinearly in $n$ and increasing exponentially in $d$.

In Fig. 2-(right), we compare our estimator with the true GMI and the standard KDE Moon et al. (2017). For the KDE approach, we estimated the joint and marginal densities and then plugged them into the proposed expression (3) in Moon et al. (2017). The bandwidth used for the the KDE plug-in estimator was selected by setting $h=n^{-1/(d+1)}$ to minimize the MSE of the plug in estimator. We generated data from the two dimensional Normal distribution with zero mean and covariance matrix

[TABLE]

We varied $\rho$ in range $[0.1,0.9]$ and computed the proposed estimated GMI and the KDE. We used the Monte Carlo method to approximate the integral (2) and compute the GMI measure (labeled True GMI). We see that as $\rho$ increases the estimated GMI using the FR test statistic outperforms the KDE approach (FR curve approaches to True curve and KDE curve diverges from True curve). In this set of experiments $\alpha=0.6$.

Fig. 3 again compares the FR estimator with the standard KDE estimator. In this setting, we draw samples from the multivariate standard Normal distribution with dimensions $d=4$ and $d=12$. In both cases the proportionality parameter $\alpha=0.5$. The left plots in Fig. 3 show the MSE (100 trials) of the GMI estimator implemented with an KDE estimator (with bandwidth as in Fig. 2 i.e. $h=n^{-1/(d+1)}$) for dimensions $d=4,12$ and various sample sizes. For all dimensions and sample sizes the FR estimator also outperforms the plug-in KDE estimator based on the estimated log-log MSE slope given in Fig. 3 (left plots). The right plots in Fig. 3 show the geometric mutual information estimated by KED and FR approaches. The iteration in this experiment is 100 and error bars are standard deviations. We observe that for higher dimension $d=12$ and larger sample size $n$, the KDE estimator approaches to the true GMI slower than the FR estimator. This reflects the power of graph-based (direct) estimators, and in particular our proposed FR estimator.

The comparison between MSEs for various dimension $d$ is shown in Fig. 4 (left). This experiment highlights higher dimension’s impact on our proposed FR estimator for the GMI measure. As expected, for larger sample size $n$, MSE decreases while for higher dimension it increases. In this setting, we have generated samples from standard Normal distribution with size $n\in[10^{2},4\times 10^{3}]$ and $\alpha=0.5$. From Fig. 4 (left) we observe that for larger sample size, MSE curves are ordered based on their corresponding dimensions. Fig. 4 (right) illustrates the MSE vs proportion parameter $\alpha$ when $n=500,10^{4}$ samples are generated from Normal distribution with $\rho=0.7,0.5$. First, following subsection 2.4, we compute the bound $[\alpha_{0}^{L},\alpha_{0}^{U}]$ and then derive optimal $\alpha$ in this range. Therefore, each experiment with different sample size and $\rho$ provides different range $[\alpha_{0}^{L},\alpha_{0}^{U}]$. We observe that the MSE is not necessarily a monotonic function in $\alpha$ and its behavior strongly depends on sample size $n$, $d$, and density functions’ bounds. This is studied extensively in Appendix F. In this set of experiments $\Xi(\alpha_{0}^{L})>0$, therefore following the results in subsection 2.4, we have $\tilde{\alpha}=\alpha_{0}^{L}$. The optimal $\alpha$ is indicated in the Fig. 4 (right).

Next, we consider three scenarios to analyze parameter $\alpha$. In these scenarios the lower bounds $C^{L}_{X},\;C^{L}_{Y},\;C^{L}_{XY}$ and upper bounds $C^{U}_{X},\;C^{U}_{Y},\;C^{U}_{XY}$ are unknown, therefore results in Section 2.4 are not applicable. In these set of experiments we varied $\alpha$ in the range $(0,1)$ to divide our original sample. We generated sample from the multivariate standard Normal distribution in all three scenarios (all features are independent). Therefore the true GMI is zero and in all scenarios the GMI column is compared with zero. In each scenario we fixed dimension $d$ and sample size $n$ and varied $\alpha=0.2,0.5,0.8$. The dimension and sample size in Scenarios 1,2, and 3 are $d=6,8,10$ and $n=1000,1500,2000$, respectively. In Table 1 the last column ($\alpha$) stars the parameter $\alpha\in\{0.2,0.5,0.8\}$ with the minimum MSE and GMI $(I_{\alpha})$ in each scenario. Table 1 shows that in these sets of experiments when $\alpha=0.5$, the GMI estimator has less MSE (i.e. is more accurate) than when $\alpha=0.2$ or $\alpha=0.8$. This experimentally demonstrates that if we split our training data, the proposed Algorithm 1 performs better than $\alpha=0.2\;\hbox{or}\;0.8$ to estimate the GMI measure.

Finally Fig. 5 shows the runtime (KDE - FR methods) as a function of sample size $n$. We vary sample size in the range $[10^{3},10^{4}]$. Observe that for smaller number of samples the KDE method is slightly faster but as $n$ becomes large we see significant relative speedup of the proposed FR method.

4 Conclusion

In this paper we have proposed a new measure of mutual information, called Geometric MI (GMI), that is related to the Henze-Penrose divergence. The GMI can be viewed as dependency measure that is the limit of the Friedman-Rafsky test statistic that depends on the minimal spanning tree over all data points. We established some properties of the GMI in terms of convexity/concavity, chain rule, and analogous to a data processing inequality. A direct estimator of the GMI was introduced that uses random permutations of observed relationships between variables in the multivariate samples. An explicit form for the MSE convergence rate bound was derived that depends on a free parameter called the proportionality parameter. An asymptotically optimal form for this free parameter was given that minimizes the MSE convergence rate. Simulation studies were performed that illustrate and verify the theory.

Acknowledgments

The work presented in this paper was partially supported by ARO grant W911NF-15-1-0479 and DOE grant de-na0002534. The authors would like to thank Brandon Oselio for the useful comments.

5 Appendices

We organize the appendices as the following: Theorem 2.1.1 which establishes convexity/concavity property is proved in Appendix A. Appendix B and C are arranged to establish the inequality (9) and (12) for given $p\in(0,1)$, respectively. In Appendix D, we first prove that the set $\widetilde{\mathcal{Z}}_{n^{\prime\prime}}$ which is randomly generated from original dependent data, contains independent samples asymptotically. Later by using the generated independent sample $\widetilde{\mathcal{Z}}_{n^{\prime\prime}}$ we show that for given $\alpha$ the FR estimator of the GMI given in Algorithm 1 intends to $I_{\alpha}$ approximately. Appendix E is dedicated to the Theorem 2.3. A full discussion on the proportionality parameter ($\alpha$) optimization strategy is provided in Appendix F.

5.1 Appendix A: Theorem 2.1.1

Proof.

The proof is similar to the result for standard (Shannon) mutual information. However, we require the following lemma, proved in analogue manner as the log-sum inequality: {Lemma} For non-negative real numbers $\alpha_{1},\dots,\alpha_{n}$ and $\beta_{1},\dots,\beta_{n}$, given $p\in(0,1)$, $q=1-p$, the following holds

[TABLE]

Notice this follows by using the convex function u(y)=\displaystyle y^{2}\big{/}(p+q\;y) for any $p\in(0,1)$, $q=1-p$, and the Jensen inequality.

Define the shorthand $\int\limits_{\mathbf{x}}$, $\int\limits_{\mathbf{y}}$, and $\int\limits_{\mathbf{x}\mathbf{y}}$ for $\int{\rm d}\mathbf{x}$, $\int{\rm d}\mathbf{y}$ and $\iint\;{\rm d}\mathbf{x}{\rm d}\mathbf{y}$, respectively. To prove part (i) of the Theorem 2.1.1, we represent the LHS of (5) as:

[TABLE]

Furthermore, the RHS of (5) can be rewritten as

[TABLE]

Thus, in order to prove $\hbox{LHS}\geq\hbox{RHS}$, we use the inequality below:

[TABLE]

In Lemma 5.1, let

[TABLE]

and for $i=1,2$,

[TABLE]

Then the claimed assertion (i) is obtained. Part (ii) follows by convexity of $D_{p}$ and the following expression:

[TABLE]

Therefore the claim in (6) is proved. ∎

5.2 Appendix B: Theorem 2.1.1

Proof.

We prove part (i) and the second part (ii) is shown by repetition method. To show (9), we use the inequality below. Given $p\in(0,1)$ and $q=1-p$, we can easily check that for positive $t>q$, $s>q$, such that $t,s\neq 1$:

[TABLE]

This implies

[TABLE]

By substituting

[TABLE]

we get

[TABLE]

Consequently

[TABLE]

Here $f$ is the joint PDF of random vector $(X_{1},X_{2},Y)$. This completes the proof of part (i). ∎

5.3 Appendix C: Proposition 2.1.1

Proof.

Recall the Theorem 2.1.1, part (i). First from $\mathbf{X}\rightarrow\mathbf{Y}\rightarrow\mathbf{Z}$ we have $f_{XYZ}=f_{XY}f_{Z|Y}$ and then by applying the Jensen inequality, we can write

[TABLE]

where

[TABLE]

Now by Markovian property we can immediately simplify the last line of (31) to the RHS in (12).

Furthermore we can easily show that if $\mathbf{X}\rightarrow\mathbf{Z}\rightarrow\mathbf{Y}$, we have $f_{XYZ}=f_{ZX}f_{Y|Z}$ and therefore $I_{p}(\mathbf{Z};\mathbf{X})=I_{p}(\mathbf{X};\mathbf{Y},\mathbf{Z})$. This together with (31) proves that under both conditions $\mathbf{X}\rightarrow\mathbf{Y}\rightarrow\mathbf{Z}$ and $\mathbf{X}\rightarrow\mathbf{Z}\rightarrow\mathbf{Y}$, the equality $I_{p}(\mathbf{X};\mathbf{Y})=I_{p}(\mathbf{Z};\mathbf{X})$ holds true. ∎

5.4 Appendix D: Theorem 2.2

Proof.

In this appendix, we first derive two required Lemmas 5.4 and 5.4 below: {Lemma} Consider random vector $\mathbf{Z}=(\mathbf{X},\mathbf{Y})$ with joint probability density function (pdf) $f_{XY}$. Let $\mathfrak{Z}_{n}=\{\mathbf{z}_{1},\ldots,\mathbf{z}_{n}\}=\{(\mathbf{x}_{i},\mathbf{y}_{i})\}_{i=1}^{n}$ be a set of samples with pdf $f_{XY}$. Let ${\mathfrak{Z}}^{\prime}_{n^{\prime}}$ and ${\mathfrak{Z}}^{\prime\prime}_{n^{\prime\prime}}$ be two distinct subsets of $\mathfrak{Z}_{n}$ such that $n^{\prime}+n^{\prime\prime}=n$ and sample proportion is $\alpha=n^{\prime}/n$ and $\beta=1-\alpha$. Next, let $\widetilde{\mathfrak{Z}}_{n^{\prime\prime}}=\{\widetilde{\mathbf{z}}_{1},\ldots,\widetilde{\mathbf{z}}_{n^{\prime\prime}}\}$ be a set of pairs such that $\widetilde{\mathbf{z}}_{k}=(\mathbf{x}_{i_{k}},\mathbf{y}_{j_{k}})$, $k=1,\dots,n^{\prime\prime}$ are selected at random from ${\mathfrak{Z}}^{\prime}_{n^{\prime\prime}}$. Denote $\widetilde{\mathbf{Z}}=(\widetilde{\mathbf{X}},\widetilde{\mathbf{Y}})$ as the random vector corresponding to samples in $\widetilde{\mathfrak{Z}}_{n^{\prime\prime}}$. Then as $n\rightarrow\infty$ such that $n^{\prime\prime}$ also grows in a linked manner that $\beta\neq 0$ then the distribution of $\widetilde{\mathbf{Z}}$ convergences to $f_{X}\times f_{Y}$ i.e. random vectors $\widetilde{\mathbf{X}}$ and $\widetilde{\mathbf{Y}}$ become independent.

Proof.

Consider two subsets $\mathbb{A},\;\mathbb{B}\subset{\mathbb{R}}^{n}$, then we have

[TABLE]

Here $\mathbf{I}_{\mathbb{A}}$ stands for the indicator function. Note that

[TABLE]

and $\mathbf{X}_{i}$ and $\mathbf{Y}_{j}$, $i\neq j$ are independent, therefore

[TABLE]

this implies that

[TABLE]

On the other hand, we know that $n^{\prime\prime}=\beta\;n$, so we get

[TABLE]

From (33), we observe that when $\beta$ takes larger values the bound becomes tighter. So if $n\rightarrow\infty$ such that $n^{\prime\prime}$ also becomes large enough in a linked manner that $\beta\neq 0$ then the RHS in (33) tends to zero. This implies that $\widetilde{\mathbf{X}}$ and $\widetilde{\mathbf{Y}}$ become independent when $n\rightarrow\infty$. ∎

An immediate result of Lemma 5.4 is the following: {Lemma} For given random vector $\mathbf{Z}_{n}=(\mathbf{X}_{n},\mathbf{Y}_{n})$ from joint density function $f_{XY}$ and with marginal density functions $f_{X}$ and $f_{Y}$ set \widetilde{\mathfrak{Z}}_{n^{\prime\prime}}=\big{\{}\widetilde{\mathbf{z}}_{1},\dots,\widetilde{\mathbf{z}}_{n^{\prime\prime}}\big{\}} be realization of random vector $\widetilde{\mathbf{Z}}$ as in Lemma 5.4 with parameter $\beta=n^{\prime\prime}/n$. Then for given points of $\widetilde{\mathfrak{Z}}_{n^{\prime\prime}}$ at $\widetilde{\mathbf{z}}=(\widetilde{\mathbf{x}},\widetilde{\mathbf{y}})$, we have

[TABLE]

Now, let us get back to our main goal which is the proof of assertion (14). Consider two subsets $\mathbf{Z}^{\prime}_{n^{\prime}}$ and $\widetilde{\mathbf{Z}}_{n^{\prime\prime}}$ as described in Subsection 2.2 . Assume that the components of sample $\widetilde{\mathbf{Z}}_{n^{\prime\prime}}$ follow density function $\widetilde{f}_{\widetilde{X}\widetilde{Y}}$. Therefore by owing to Lemma 5.4 and 5.4, when $n\rightarrow\infty$ then $\widetilde{f}_{\widetilde{X}\widetilde{Y}}\rightarrow f_{X}f_{Y}$. Let $M_{n^{\prime}}$ and $N_{n^{\prime\prime}}$ be Poisson variables with mean $n^{\prime}$ and $n^{\prime\prime}$ independent of one another and $\{{\mathbf{Z}}^{\prime}_{i}\}$ and $\{\widetilde{\mathbf{Z}}_{j}\}$. Assume two Poisson processes \mathfrak{Z}^{\prime}_{n^{\prime}}=\big{\{}{\mathbf{Z}}^{\prime}_{1},\dots,{\mathbf{Z}}^{\prime}_{M_{n^{\prime}}}\big{\}} and \widetilde{\mathfrak{Z}}_{n^{\prime\prime}}=\big{\{}\widetilde{\mathbf{Z}}_{1},\dots,\widetilde{\mathbf{Z}}_{N_{n^{\prime\prime}}}\big{\}}, with the FR statistic $\mathfrak{R}^{\prime}_{n^{\prime},n^{\prime\prime}}$. So by owing to the arguments in Henze and Penrose (1999), Yasaei Sekeh et al. (2018) we shall prove the following:

[TABLE]

This is because of \big{|}\mathfrak{R}^{\prime}_{n^{\prime},n^{\prime\prime}}-\mathfrak{R}_{n^{\prime},n^{\prime\prime}}\big{|}\leq K_{d}\Big{(}|M_{n^{\prime}}-n^{\prime}|+|N_{n^{\prime\prime}}-n^{\prime\prime}|\Big{)}, where $K_{d}$ is a constant defined in Lemma 1, Henze and Penrose (1999) and $n^{\prime}+n^{\prime\prime}=n$, then (n^{\prime}+n^{\prime\prime})^{-1}\;{\mathbb{E}}\big{|}\mathfrak{R}^{\prime}_{n^{\prime},n^{\prime\prime}}-\mathfrak{R}_{n^{\prime},n^{\prime\prime}}\big{|}\rightarrow 0 as $n\rightarrow\infty$. Let $\mathbf{W}_{1}^{n^{\prime},n^{\prime\prime}},\mathbf{W}_{2}^{n^{\prime},n^{\prime\prime}},\dots$ be independent variables with common density

[TABLE]

for $(\mathbf{x},\mathbf{y})\in{\mathbb{R}}^{d}\times{\mathbb{R}}^{d}$. Let $L_{n^{\prime},n^{\prime\prime}}$ be an independent Poisson variable with mean $n^{\prime}+n^{\prime\prime}$. Let $\mathfrak{F}^{\prime}_{n^{\prime},n^{\prime\prime}}=\left\{\mathbf{W}_{1}^{n^{\prime},n^{\prime\prime}},\dots,\mathbf{W}_{L_{n^{\prime},n^{\prime\prime}}}^{n^{\prime},n^{\prime\prime}}\right\}$ a nonhomogeneous Poisson process of rate $n^{\prime}\;f_{XY}+n^{\prime\prime}\;\widetilde{f}_{\widetilde{X}\widetilde{Y}}$. Assign mark $1$ to a point in $\mathfrak{F}^{\prime}_{n^{\prime},n^{\prime\prime}}$ with probability

[TABLE]

and mark 2 otherwise, having the FR test statistic $\widetilde{\mathfrak{R}}_{n^{\prime},n^{\prime\prime}}$ which by the marking theorem Kingman (1993), Henze and Penrose (1999), has a same distribution as $\mathfrak{R}^{\prime}_{n^{\prime},n^{\prime\prime}}$. Given points of $\mathfrak{F}^{\prime}_{n^{\prime},n^{\prime\prime}}$ at $\mathbf{z}^{\prime}=(\mathbf{x}^{\prime},\mathbf{y}^{\prime})$ and $\mathbf{z}^{\prime\prime}=(\mathbf{x}^{\prime\prime},\mathbf{y}^{\prime\prime})$, the probability that they have different marks is given by (35).

[TABLE]

Set

[TABLE]

then

[TABLE]

Now recall (36). We observe that $g_{n^{\prime},n^{\prime\prime}}(\mathbf{z}^{\prime},\mathbf{z}^{\prime\prime})\rightarrow g(\mathbf{z}^{\prime},\mathbf{z}^{\prime\prime})$. Going back to (37), we can write

[TABLE]

Consider the non-Poisson process

[TABLE]

So by the fact that {\mathbb{E}}\Big{[}M_{n^{\prime}}+N_{n^{\prime\prime}}-(n^{\prime}+n^{\prime\prime})\Big{]}=o(n^{\prime}+n^{\prime\prime}), we have

[TABLE]

Introduce

[TABLE]

Then $\phi_{n^{\prime},n^{\prime\prime}}(\mathbf{x},\mathbf{y})\rightarrow\phi(\mathbf{x},\mathbf{y})$ uniformly because of $\displaystyle{n^{\prime}}/{n}\rightarrow\alpha$ and $\displaystyle{n^{\prime\prime}}/{n}\rightarrow\beta$. Thus by using Proposition 1 in Henze and Penrose (1999), we get

[TABLE]

So, we conclude the proof. ∎

5.5 Appendix E: Theorem 2.3

Proof.

We begin by providing a family of bias rate bound for the FR estimator $\mathfrak{R}_{n^{\prime},n^{\prime\prime}}$ in terms of a parameter $l$. Then by plugging the optimal $l$, we prove the bias rate bound given in (18). {Theorem} Let $\mathfrak{R}_{n^{\prime},n^{\prime\prime}}:=\mathfrak{R}(\mathfrak{Z}_{n^{\prime}},\mathfrak{Z}_{n^{\prime\prime}})$ be the FR test statistic. Then a bound on the bias rate of the $\mathfrak{R}_{n^{\prime},n^{\prime\prime}}$ estimator for $0<\eta\leq 1$, $d\geq 2$ is given by

[TABLE]

where $0<\eta\leq 1$ is the Hölder smoothness parameter and $c_{d}$ is a constant depending only on $d$.

Set

[TABLE]

and

[TABLE]

where

[TABLE]

A more explicit form for the bound on the RHS is given below:

[TABLE]

Proof.

Consider two Poisson variables $M_{n^{\prime}}$ and $N_{n^{\prime\prime}}$ with mean $n^{\prime}$ and $n^{\prime\prime}$ respectively and independent of one another and $\{\mathbf{Z}^{\prime}_{i}\}$ and $\{\widetilde{\mathbf{Z}}_{j}\}$. Let $\mathfrak{Z}^{\prime}_{n^{\prime}}$ and $\widetilde{\mathfrak{Z}}_{n^{\prime\prime}}$ be the Poisson processes $\{\mathbf{Z}^{\prime}_{1},\dots,\mathbf{Z}^{\prime}_{M_{n^{\prime}}}\}$ and $\{\widetilde{\mathbf{Z}}_{1},\dots,\widetilde{\mathbf{Z}}_{N_{n^{\prime\prime}}}\}$. Likewise Appendix D, set $\mathfrak{R}^{\prime}_{n^{\prime},n^{\prime\prime}}=\mathfrak{R}(\mathfrak{Z}^{\prime}_{n^{\prime}},\widetilde{\mathfrak{Z}}_{n^{\prime\prime}})$. Applying Lemma 1, and (12) in Henze and Penrose (1999), we can write

[TABLE]

Here $c_{d}$ denotes the largest possible degree of any vertex of the MST in ${\mathbb{R}}^{d}$. Following the arguments in Yasaei Sekeh et al. (2018), we have {\mathbb{E}}\big{[}\big{|}M_{n^{\prime}}-n^{\prime}\big{|}\big{]}=O\big{(}\;{n^{\prime}}^{1/2}\big{)} and {\mathbb{E}}\big{[}\big{|}N_{n^{\prime\prime}}-n^{\prime\prime}\big{|}\big{]}=O\big{(}\;{n^{\prime\prime}}^{1/2}\big{)}. Hence

[TABLE]

Next let $n^{\prime}_{i}$ and $n^{\prime\prime}_{i}$ be independent binomial random variables with marginal densities $B(n^{\prime},a_{i}l^{-d})$ and $B(n^{\prime\prime},b_{i}l^{-d})$ such that $a_{i},b_{i}$ are non-negative constants $a_{i}\leq b_{i}$ and $\displaystyle\sum\limits_{i=1}^{l^{d}}a_{i}l^{-d}=\displaystyle\sum\limits_{i=1}^{l^{d}}b_{i}l^{-d}=1$. Therefore using subadditivity property in Lemma 2.2, Yasaei Sekeh et al. (2018), we can write

[TABLE]

where $M=l^{d}$, and $\eta>0$ stands with Hölder smoothness parameter. Further, for given $n^{\prime}_{i}$, $n^{\prime\prime}_{i}$, let $\mathbf{W}_{1}^{n^{\prime}_{i},n^{\prime\prime}_{i}},\mathbf{W}_{2}^{n^{\prime}_{i},n^{\prime\prime}_{i}},\dots$ be independent variables with common densities for $(\mathbf{x},\mathbf{y})\in{\mathbb{R}}^{d}\times{\mathbb{R}}^{d}$:

[TABLE]

Denote $L_{n^{\prime}_{i},n^{\prime\prime}_{i}}$ be an independent Poisson variable with mean $n^{\prime}_{i}+n^{\prime\prime}_{i}$ and \mathfrak{F}^{\prime}_{n^{\prime}_{i},n^{\prime\prime}_{i}}=\big{\{}\mathbf{W}_{1}^{n^{\prime}_{i},n^{\prime\prime}_{i}},\dots,\mathbf{W}_{L_{n^{\prime}_{i}.n^{\prime\prime}_{i}}}^{n^{\prime}_{i},n^{\prime\prime}_{i}}\big{\}} a non-homogeneous Poisson of rate $n^{\prime}_{i}f_{XY}+n^{\prime\prime}_{i}\widetilde{f}_{\widetilde{X}\widetilde{Y}}$. Let $\mathfrak{F}_{n^{\prime}_{i},n^{\prime\prime}_{i}}$ be the non-Posisson point process \big{\{}\mathbf{W}_{1}^{n^{\prime}_{i},n^{\prime\prime}_{i}},\dots\mathbf{W}_{n^{\prime}_{i}+n^{\prime\prime}_{i}}^{n^{\prime}_{i},n^{\prime\prime}_{i}}\big{\}}. Assign a mark from the set $\{1,2\}$ to each points of $\mathfrak{F}^{\prime}_{n^{\prime}_{i},n^{\prime\prime}_{i}}$. Let $\widetilde{\mathfrak{Z}}^{\prime}_{n^{\prime}_{i}}$ be the sets of points marked 1 with each probability n^{\prime}_{i}f_{XY}(\mathbf{x},\mathbf{y})\big{/}\big{(}n^{\prime}_{i}f_{XY}(\mathbf{x},\mathbf{y})+n^{\prime\prime}_{i}\widetilde{f}_{\widetilde{X}\widetilde{Y}}(\mathbf{x},\mathbf{y})\big{)} and let $\widetilde{\mathfrak{Z}}^{\prime\prime}_{n^{\prime\prime}_{i}}$ be the set points with mark 2. Note that owing to the marking theorem Kingman (1993), $\widetilde{\mathfrak{Z}}^{\prime}_{n^{\prime}_{i}}$ and $\widetilde{\mathfrak{Z}}^{\prime\prime}_{n^{\prime\prime}_{i}}$ are independent Poisson processes with the same distribution as $\mathfrak{Z}^{\prime}_{n^{\prime}_{i}}$ and $\widetilde{\mathfrak{Z}}_{n^{\prime\prime}_{i}}$ , respectively. Considering $\widetilde{\mathfrak{R}}^{\prime}_{n^{\prime}_{i},n^{\prime\prime}_{i}}$ as FR test statistic on nodes in $\widetilde{\mathfrak{Z}}^{\prime}_{n^{\prime}_{i}}\cup\widetilde{\mathfrak{Z}}^{\prime\prime}_{n^{\prime\prime}_{i}}$, we have

[TABLE]

By the fact that {\mathbb{E}}\big{[}|M_{n^{\prime}}+N_{n^{\prime\prime}}-n^{\prime}-n^{\prime\prime}|\big{]}=O((n^{\prime}+n^{\prime\prime})^{1/2}), we have

[TABLE]

Here $\mathbf{z}^{\prime}=(\mathbf{x}^{\prime},\mathbf{y}^{\prime})$, $\mathbf{z}^{\prime\prime}=(\mathbf{x}^{\prime\prime},\mathbf{y}^{\prime\prime})$, and $P_{n^{\prime}_{i},n^{\prime\prime}_{i}}(\mathbf{z}^{\prime},\mathbf{z}^{\prime\prime})$ is given in below:

[TABLE]

Next set

[TABLE]

By owing to the Lemma B.6 in Yasaei Sekeh et al. (2018) and applying the analogous arguments, we can write the expression in (51):

[TABLE]

where

[TABLE]

Going back to Lemma 5.4, we know that

[TABLE]

Therefore the first term in RHS of (51) turns to be less and equal than

[TABLE]

and the second term is less and equal than

[TABLE]

where

[TABLE]

Recall the definition of the dual MST and FR statistic denoted by $\mathfrak{R}^{*}_{n^{\prime},n^{\prime\prime}}$ from Yasaei Sekeh et al. (2018):

Definition: (Dual MST, ${\rm MST}^{*}$ and dual FR statistic $\mathfrak{R}_{m,n}^{*}$) Let $\mathbb{F}_{i}$ be the set of corner points of the subsection $Q_{i}$ for $1\leq i\leq l^{d}$. Then we define ${\rm MST}^{*}(\mathfrak{X}_{m}\cup\mathfrak{Y}_{n}\cap Q_{i})$ as the boundary MST graph of partition $Q_{i}$ Yukich (1998), which contains $\mathfrak{X}_{m}$ and $\mathfrak{Y}_{n}$ points falling inside the section $Q_{i}$ and those corner points in $\mathbb{F}_{i}$ which minimize total MST length. Notice it is allowed to connect the MSTs in $Q_{i}$ and $Q_{j}$ through points strictly contained in $Q_{i}$ and $Q_{j}$ and corner points are taking into account under condition of minimizing total MST length. Another word, the dual MST can connect the points in $Q_{i}\cup Q_{j}$ by direct edges to pair to another point in $Q_{i}\cup Q_{j}$ or the corner the corner points (we assume that all corner points are connected) in order to minimize the total length. To clarify this, assume that there are two points in $Q_{i}\cup Q_{j}$, then the dual MST consists of the two edges connecting these points to the corner if they are closed to a corner point otherwise dual MST consists of an edge connecting one to another. Further, $\mathfrak{R}_{m,n}^{*}(\mathfrak{X}_{m},\mathfrak{Y}_{n}\cap Q_{i})$ is defined as the number of edges in ${\rm MST}^{*}$ graph connecting nodes from different samples and number of edges connecting to the corner points. Note that the edges connected to the corner nodes (regardless of the type of points) are always counted in dual FR test statistic $\mathfrak{R}^{*}_{m,n}$.

See Yasaei Sekeh et al. (2018) for more detail about the dual MST and dual FR test statistic. Similar discussion as above and in Yasaei Sekeh et al. (2018), consider the Poisson processes samples and the FR test statistic under the union of samples, denoted by ${\mathfrak{R}^{\prime}}^{*}_{n^{\prime},n^{\prime\prime}}$, and superadditivity of dual $\mathfrak{R}^{*}_{n^{\prime},n^{\prime\prime}}$, we have

[TABLE]

The first term of RHS in (54) is greater and equal than

[TABLE]

Furthermore,

[TABLE]

where $c_{d}$ is the largest possible degree of any vertex of the MST in ${\mathbb{R}}^{d}$ as before. Consequently, we have

[TABLE]

where $\mathcal{B}$ is defined in (44) and $\mathcal{A}_{f,n}^{\beta,\alpha}(\mathbf{x},\mathbf{y})$ has been introduced in (42). The last line in (57) is implied from the fact that

[TABLE]

Here $\mathcal{A}_{f,n}^{\beta,n^{\prime}/n}(\mathbf{x},\mathbf{y})$ is given as (42) by substituting $n^{\prime}/n$ in $\alpha$ such that $\beta=1-\alpha$. Hence, the proof of Theorem 5.5 is completed. ∎

Going back to the proof of (18), without loss of generality assume that $(n)l^{-d}>1$. In the range $d\geq 2$ and $0<\eta\leq 1$. We select $l$ as a function of $n$ and $\beta$ to be the sequence increasing in $n$ which minimizes the maximum of these rates:

[TABLE]

The solution $l=l(n,\beta)$ occurs when $l^{d}(n)^{-\eta/d}=l^{-\eta d}$, or equivalently $l=\lfloor(n)^{\eta/(d^{2}(\eta+1))}\rfloor$ and also $l^{d}\beta^{-1}n^{-1}=l^{-\eta d}$ which implies $l=\lfloor(\beta\;n)^{1/(d(1+\eta))}\rfloor$. Substitute this in $l$ in the bound (41), the RHS expression in (18) for $d\geq 2$ is derived. ∎

5.6 Appendix F

Our main goal in subsection 2.4 is to find proportion $\alpha$ such that the parametric MSE rate depending on the joint density $f_{XY}$ and marginal densities $f_{X},f_{Y}$ is minimized. Recalling the explicit bias bound in (44), it can be seen that this function is a complicated function of $f_{XY}$, $f_{X}f_{Y}$ and $\alpha$. By rearranging, we first find an upper bound for $\Delta$ in (44), denoted by $\overline{\Delta}$, as follows:

[TABLE]

where l_{n}:=\big{\lfloor}n^{\eta/(d^{2}(1+\eta))}\big{\rfloor}. From Appendix E we know that optimal $l$ is given in (59). One can check that for $\alpha\leq 1-n^{(\eta/d)-1}$, optimal l=\big{\lfloor}n^{\eta/(d^{2}(1+\eta))}\big{\rfloor} provides a tighter bound. In (61), the constants $D$ and $\overline{D}$ are

[TABLE]

[TABLE]

And the function ${G}^{\alpha,\beta}_{f,n}(\mathbf{x},\mathbf{y})$ is given as the following:

[TABLE]

where $\delta_{f}$ is given in (43). After all still the expression (65) is complicated to optimize therefore we use the fact that $0\leq\alpha,\beta\leq 1$ and bound the function $\mathbf{G}^{\alpha,\beta}_{f,n}(\mathbf{x},\mathbf{y})$. Define set $\Gamma$ as

[TABLE]

where

[TABLE]

Here $K$ is the smoothness constant. Notice the set $\Gamma$ is a convex set. We bound $\overline{\Delta}$ by

[TABLE]

Set $C_{n}=C^{L}_{XY}\;n/2$,

[TABLE]

We simplify (68) by

[TABLE]

Under condition

[TABLE]

$\widetilde{{G}}^{\alpha,\beta}_{n}(\epsilon_{XY})$ is an increasing function in $\epsilon$. Furthermore, for $\alpha\leq\displaystyle\frac{1}{4}$ and

[TABLE]

the function $\widetilde{{G}}^{\alpha,\beta}_{n}(\epsilon_{XY})$ is strictly concave. Next, to find optimal $\alpha$ we consider the following optimization problem:

[TABLE]

here \epsilon_{XY}=f_{XY}\big{/}f_{X}f_{Y}, C^{U}_{\epsilon}=C^{U}_{XY}\big{/}C^{L}_{X}C^{L}_{Y} and C^{L}_{\epsilon}=C^{L}_{XY}\big{/}C^{U}_{X}C^{U}_{Y}, such that $C^{L}_{\epsilon}\leq 1$. We know that under conditions (69) and (70), the function $\widetilde{{G}}^{\alpha,\beta}_{n}$ is strictly concave and increasing in $\epsilon_{XY}$.

We first solve the optimization problem:

[TABLE]

where

[TABLE]

The Lagrangian for this problem is

[TABLE]

In this case the optimum $\epsilon^{*}_{XY}$ is bounded, $\theta^{L}_{\epsilon}(\alpha)\leq\epsilon^{*}_{XY}\leq\theta^{U}_{\epsilon}(\alpha),$ and Lagrangian multiplier $\lambda^{*}_{1},\;\lambda^{*}_{2}\geq 0$ such that

[TABLE]

Set $G^{\prime}_{n}(\epsilon_{XY})=\displaystyle\frac{{\rm d}}{{\rm d}\epsilon_{XY}}\widetilde{{G}}^{\alpha,\beta}_{n}(\epsilon_{XY}).$ In view of the concavity of $\widetilde{{G}}^{\alpha,\beta}_{\epsilon_{XY},n}$ and Lemma 1, page 227 in Luenberger (1969), maximizing $L(\epsilon_{XY},\lambda^{*}_{1},\lambda^{*}_{2})$ over $\epsilon_{XY}$ is equivalent to

[TABLE]

for all $\theta_{\epsilon}^{L}(\alpha)\leq\epsilon^{*}_{XY}\leq\theta_{\epsilon}^{U}(\alpha)$, and

[TABLE]

Denote ${G^{\prime}}^{-1}_{n}$ the inverse function of $G^{\prime}_{n}$. Since $G^{\prime}_{n}$ is strictly decreasing in $\epsilon^{*}_{XY}$ (this is because $\tilde{G}^{\alpha,\beta}_{n}(\epsilon_{XY})$ is strictly concave, then ${G^{\prime}}^{-1}_{n}$ is continues and strictly decreasing in $\epsilon^{*}_{XY}$). From (75) and (76), we see immediately that on any interval $\theta^{L}_{\epsilon}(\alpha)\leq\epsilon^{*}_{XY}\leq\theta^{U}_{\epsilon}(\alpha)$, we have $\epsilon^{*}_{XY}={G^{\prime}}^{-1}_{n}(\lambda^{*}_{1}-\lambda^{*}_{2})$. We can write then

[TABLE]

and $\lambda^{*}_{1},\lambda^{*}_{2}\geq 0$. Next, we intend to find the solution of

[TABLE]

[TABLE]

The function $\overline{{G}}^{\alpha,\beta}_{n}(\lambda_{1},\lambda_{2})$ is increasing in $\lambda_{1}$ and $\lambda_{2}$, therefore it takes its minimum at (\lambda^{*}_{1},\lambda^{*}_{2})=(G^{\prime}_{n}\big{(}{\theta^{U}_{\epsilon}}(\alpha)\big{)},0). This implies that $\epsilon^{*}_{XY}=\theta^{U}_{\epsilon}(\alpha)$. We continue by going back to our primary minimization over $\alpha$:

[TABLE]

where $\alpha_{0}^{L}=\displaystyle\frac{2}{C_{n}}$ and \alpha_{0}^{U}=\min\Big{\{}\displaystyle\frac{1}{4},1-n^{\eta/d-1}\Big{\}}. We know that $\displaystyle\frac{1}{4}\leq\frac{1}{3}+\displaystyle\frac{2}{3C_{n}}$ and $\displaystyle\frac{1}{4}\leq\frac{1}{2}+\frac{1}{2C_{n}}$, therefore the condition below

[TABLE]

implies previous conditions on $\alpha$. Since the objective function (78) is a complicated function in $\alpha$, it is not feasible to determine whether it is a convex function in $\alpha$. For this reason let us solve the optimization problem in (78) in a special case when $C^{U}_{\epsilon}\leq\theta^{U}(\alpha)$. This implies $\epsilon^{*}_{XY}=C^{U}_{\epsilon}$. Under assumption $C^{U}_{\epsilon}$ the objective function in (78) is convex in $\alpha$. Also the case $C^{U}_{\epsilon}\leq\theta^{U}(\alpha)$ is equivalent to $\alpha\leq\displaystyle\frac{1+1/C_{n}}{4+2C^{U}_{\epsilon}}$. Therefore in the optimization problem we have constraint

[TABLE]

We know that $\widetilde{\Delta}(\alpha,\epsilon^{*}_{XY})+c_{d}(1-\alpha)n^{-1}$ is convex over $\alpha\in[\alpha^{L}_{0},\alpha^{U}_{0}]$. So, the problem becomes ordinary convex optimization problem. Let $\tilde{\alpha}$, $\tilde{\lambda}_{1}$ and $\tilde{\lambda}_{2}$ be any points that satisfy the KKT conditions for this problem:

[TABLE]

Recall $\Xi(\alpha)$ from (26):

[TABLE]

where $\widetilde{\Delta}$ is given in (66). So the last condition in (84) becomes $\Xi(\tilde{\alpha})=\tilde{\lambda}_{1}-\tilde{\lambda}_{2}$. We then have

[TABLE]

where $\Xi^{-1}$ is inverse function of $\Xi$. Since $\alpha_{0}^{L}\neq\alpha_{0}^{U}$, at least one of $\tilde{\lambda}_{1}$ or $\tilde{\lambda}_{2}$ should be zero:

• •

$\tilde{\lambda}_{1}=0$, $\tilde{\lambda}_{2}\neq 0$. Then $\tilde{\alpha}=\alpha_{0}^{U}$ and implies $\tilde{\lambda}_{2}=-\Xi(\alpha_{0}^{U})$. Since $\tilde{\lambda}_{2}>0$, so this leads to $\Xi(\alpha_{0}^{U})<0$.

• •

$\tilde{\lambda}_{2}=0$, $\tilde{\lambda}_{1}\neq 0$. Then $\tilde{\alpha}=\alpha_{0}^{L}$ and implies $\tilde{\lambda}_{1}=\Xi(\alpha_{0}^{L})$. We know that $\tilde{\lambda}_{1}>0$, hence $\Xi(\alpha_{0}^{L})>0$.

• •

$\tilde{\lambda}_{1}=0$, $\tilde{\lambda}_{2}=0$. Then $\tilde{\alpha}=\Xi^{-1}(0)$ and so $\alpha_{0}^{L}\leq\Xi^{-1}(0)\leq\alpha_{0}^{U}$.

Consequently, by following the behavior of $\Xi(\alpha)$ with respect to $\alpha^{L}_{0}$ and $\alpha^{U}_{0}$, we are able to detect optimal $\tilde{\alpha}$, $\tilde{\lambda}_{1}$ and $\tilde{\lambda}_{2}$. For instance, if $\Xi(\alpha)$ is positive for all $\alpha\in[\alpha^{L}_{0},\alpha^{U}_{0}]$ then we conclude that $\tilde{\alpha}=\alpha^{L}_{0}$.

\reftitle

References