Saddlepoint approximations to P-values for comparison of density estimates

Summary

This article is concerned with computing approximate p-values for the maximum of the absolute difference between kernel density estimates. The approximations are based on treating the process of local extrema of the differences as a nonhomogeneous Poisson Process and estimating the corresponding local intensity function. The process of local extrema is characterized by the intensity function, which determines the rate of local extrema above a given threshold. A key idea of this article is to provide methods for more accurate estimation of the intensity function by using saddlepoint approximations for the joint density of the difference between kernel density estimates and using the first and second derivative of the difference. In this article, saddlepoint approximations are compared to gaussian approximations. Simulation results from saddlepoint approximations show consistently better agreement between empirical p-value and predetermined value with various bandwidths of kernel density estimates.

Appendix

A-1.

$$\begin{array}{*{20}{c}} {E\Delta \left( t \right) = E\left\{ {\frac{1}{n}\sum\limits_{i = 1}^n {w\left( {{x_i} - t} \right) - \frac{1}{m}\sum\limits_{j = 1}^m {w\left( {{y_j} - t} \right)} } } \right\}} \\\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; { = \frac{1}{n} \cdot n \cdot E\left( {w\left( {{x_i} - t} \right)} \right) - \frac{1}{m} \cdot m \cdot E\left( {w\left( {{y_j} - t} \right)} \right)} \\\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; { = \frac{1}{n} \cdot n \cdot E\left( {w\left( {{u_i}} \right)} \right) - \frac{1}{m} \cdot E \cdot \left( {w\left( {{u_i}} \right)} \right) = \mu - \mu = 0} \end{array}$$

Since u₁,u₂,...,u_n+m all have the same distribution under the null hypothesis, with mean

$$\mu = \frac{1}{{n + m}}\sum\limits_{i = 1}^{n + m} {w\left( {{u_i}} \right)} $$

$$\begin{array}{*{20}{c}} {E\Delta '\left( t \right) = E\frac{\partial }{{\partial t}}\left\{ {\frac{1}{n}\sum\limits_{i = 1}^n {w\left( {{x_i} - t} \right) - \frac{1}{m}\sum\limits_{j = 1}^m {w\left( {{y_j} - 1} \right)} } } \right\}} \\\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\ { = \frac{1}{n} \cdot n \cdot E\left( {w'\left( {{x_i} - t} \right)} \right) - \frac{1}{m} \cdot m \cdot E\left( {w'\left( {{y_j} - t} \right)} \right)} \\\:\:\:\:\:\:\:\:\: { = \frac{1}{{n + m}}\sum\limits_{i = 1}^{n + m} {w'\left( {{u_i}} \right) - \frac{1}{{n + m}}\sum\limits_{i = 1}^{n + m} {w'\left( {{u_i}} \right) = 0} } } \end{array}$$

In a similar way, E∆″(t) = 0 can be proved.

A-2.

$$\begin{array}{*{20}{l}} {\operatorname{var} \Delta \left( t \right) = E{\Delta ^2}\left( t \right) = E{{\left\{ {\frac{1}{n}\sum\limits_{i = 1}^n {w\left( {{x_i} - t} \right) - \frac{1}{m}\sum\limits_{j = 1}^m {w\left( {{y_j} - t} \right)} } } \right\}}^2}} \\\;\;\;\;\;\;\;\;\;\;\;\;\;\; { = E\left\{ {\frac{1}{{{n^2}}}{{\left( {\sum\limits_{i = 1}^n {w\left( {{x_i} - t} \right)} } \right)}^2} + \frac{1}{{{m^2}}}{{\left( {\sum\limits_{j = 1}^m {w\left( {{y_j} - t} \right)} } \right)}^2} - \frac{2}{{nm}}\sum\limits_{i = 1}^n {w\left( {{x_i} - t} \right)\sum\limits_{j = 1}^m {w\left( {{y_j} - t} \right)} } } \right\}} \\\;\;\;\;\;\;\;\;\;\;\;\;\;\; { = E\left\{ {\frac{1}{{{n^2}}} + \left[ {\sum\limits_{i = 1}^n {w{{\left( {{x_i} - t} \right)}^2} + \sum\limits_{i = 1}^n {\sum\limits_{j \ne i} {w\left( {{x_i} - t} \right)w\left( {{x_j} - t} \right)} } } } \right] + \frac{1}{{{m^2}}}\left[ {\sum\limits_{j = 1}^m {w{{\left( {{y_j} - t} \right)}^2}} + \sum\limits_{j = 1}^m {\sum\limits_{i \ne j} {w\left( {{y_j} - t} \right)w\left( {{y_i} - t} \right)} } } \right] - \frac{2}{{nm}}\sum\limits_{i = 1}^n {w\left( {{x_i} - t} \right)\sum\limits_{j = 1}^m {w\left( {{y_j} - t} \right)} } } \right\}} \\\;\;\;\;\;\;\;\;\;\;\;\;\;\; { = \frac{1}{{{n^2}}}nEw{{\left( {{u_i}} \right)}^2} + \frac{1}{{{n^2}}}n\left( {n - 1} \right)Ew\left( {{u_i}} \right)Ew\left( {{u_j}} \right) + \frac{1}{{{m^2}}}mEw{{\left( {{u_i}} \right)}^2} + \frac{{m\left( {m - 1} \right)}}{{{m^2}}} \times Ew\left( {{u_i}} \right)Ew\left( {{u_j}} \right) - \frac{2}{{nm}}nmEw\left( {{u_i}} \right)Ew\left( {{u_j}} \right)\quad \quad \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left( {i \ne j} \right)} \end{array}$$

Since all realizations of u are equally likely,

$$E\frac{n}{{i = 1}}w\left( {{u_i}} \right) = nEw\left( {{u_i}} \right)$$

$$ = \left( {\frac{1}{n} + \frac{1}{m}} \right)\frac{1}{{n + m}}\sum\limits_{i = 1}^{n + m} {w{{\left( {{u_i}} \right)}^2}} + \left( {\frac{{n - 1}}{n} + \frac{{m - 1}}{m} - 2} \right)\frac{1}{{n + m}}\frac{1}{{n + m - 1}} \times \sum\limits_{i = 1}^{n + m} {\sum\limits_{j \ne i} {w\left( {{u_i}} \right)w\left( {{u_j}} \right)} } $$

Since

$$Ew{\left( {{u_i}} \right)^2} = \frac{1}{{n + m}}\sum\limits_{i = 1}^{n + m} {w{{\left( {{u_i}} \right)}^2}} $$

, After a little more arithmetic manipulation, using the relation

$$\sum\limits_{i = 1}^n w \left( {{x_i} - t} \right){\sum\limits_{j = 1}^n {w\left( {{x_j} - t} \right) = \sum\limits_{i = 1}^n w \left( {{x_i} - t} \right)} ^2} + \sum\limits_{i = 1}^n {\sum\limits_{j \ne i} {w\left( {{x_i} - t} \right)w\left( {{x_j} - t} \right)} } $$

$\operatorname{var} \Delta \left( t \right)$ is given by

$$\left( {\frac{1}{n} + \frac{1}{m}} \right)\left\{ {\frac{1}{{n + m}}\sum\limits_{i = 1}^{n + m} {w{{\left( {{u_i}} \right)}^2} - {{\left( {\frac{{\sum\limits_{i = 1}^{n + m} {w\left( {{u_i}} \right)} }}{{n + m}}} \right)}^2}} } \right\}\frac{{n + m}}{{n + m - 1}}.$$

The C code including the module for Appendix A1 and A2 is available upon request.