A new method of joint nonparametric estimation of probability density and its support

TL;DR

This paper introduces a novel joint nonparametric method for estimating probability density and its support, effectively addressing boundary bias issues in kernel density estimation, including multivariate cases.

Contribution

It proposes a boundary detection technique that eliminates boundary bias in kernel density estimation, extending to multivariate scenarios with an improved estimator.

Findings

Successfully detects boundaries in density estimation

Eliminates boundary bias in univariate and multivariate cases

Provides a more accurate density and support estimation

Abstract

In this paper we propose a new method of joint nonparametric estimation of probability density and its support. As is well known, nonparametric kernel density estimator has "boundary bias problem" when the support of the population density is not the whole real line. To avoid the unknown boundary effects, our estimator detects the boundary, and eliminates the boundary-bias of the estimator simultaneously. Moreover, we refer an extension to a simple multivariate case, and propose an improved estimator free from the unknown boundary bias.

Tables

Table 1. Table 1: Flow of the proposed method (the left boundary l 0 subscript 𝑙 0 l_{0} is known)

0. Assume that the support of the density is from $l_0$ to $u_0$,

where $l_0$ is known, and $u_0$ is unknown but bounded.

Check the assumptions of Theorem 1.

1. Select a boundary bias reduction method,

and choose the kernel and bandwidth, following the method.

${\hat{f}}_u^{†}$ denotes the boundary bias free estimator

of the density function whose support is $l_0$ to $u$

and ${\hat{F}}_u^{†}$ denotes the cumulative distribution estimator.

2. Solve the equation for $u$ :

${\hat{F}}_u^{†}\left(X_{\left(n\right)}\right)-c_n=0$.

3. Set the solution as $\hat{u}$,

and output the boundary-adjusted estimators ${\hat{F}}_{\hat{u}}^{†}$ and ${\hat{f}}_{\hat{u}}^{†}$.

Table 2. Table 2: ISE of the kernel estimators

Beta(1,1)   Naive ^u^[BK] X_(n)^[BK] ^u^[R] X_(n)^[R] n=30 .04414 .06093 .11531 .04626 .04646 n=50 .03823 .03736 .08202 .03444 .03452 n=100 .03161 .02057 .04270 .02369 .02374 n=300 .02340 .00770 .01732 .01454 .01457 Beta(3,1)   Naive ^u^[BK] X_(n)^[BK] ^u^[R] X_(n)^[R] n=30 .23968 .18527 .34458 .21659 .21811 n=50 .21025 .11735 .24758 .17229 .17299 n=100 .17250 .05956 .13326 .12710 .12752 n=300 .11425 .02162 .06087 .07695 .07726 Beta(3,3)   Naive ^u^[BK] X_(n)^[BK] ^u^[R] X_(n)^[R] n=30 .00596 .02917 .00834 .00256 .00256 n=50 .00359 .01594 .00308 .00156 .00156 n=100 .00173 .00520 .00152 .00079 .00079 n=300 .00058 .00185 .00041 .00027 .00027

Table 3. Table 3: Boundary bias of the kernel estimators

Bias    ^f(1)  ^f_^u^[BK](1)  ^f_X_(n)^[BK](1)  ^f_^u^[R](1)  ^f_X_(n)^[R](1) f(1) ¿0    O(1) O(h)  f(1) =f’(1) =0 O(h^2) O(h^2)

Table 4. Table 4: Flow of the proposed method in the multivariate case

0. Assume that the support of the density is of the form,

$supp\left(g\right)=\Xi ={\xi }_1\times \cdots \times {\xi }_d$,

where ${\xi }_j$ is the interval from $l_{j,0}$ to $u_{j,0}$.

Check the assumptions of Theorem 3.

1. Select a boundary bias reduction method,

and choose the kernel and bandwidth as follows.

${\hat{G}}_{{𝒖}_j,j}^{†}$ and ${\hat{g}}_{{𝒖}_j,j}^{†}$ denote the marginal distribution and density estimator,

free from the boundary bias caused by the bounded support ${𝒖}_j=[l_j,u_j]$.

2. Solve the equation for ${𝒖}_j$ for every $j$ :

${\widehat{𝑮}}_{{𝒖}_j,j}^{†}\left({𝑿}_{j,(1,n)}\right)-{𝒄}_n=\mathrm{𝟎}$

(see Section 3).

3. Set the solution as ${\widehat{𝒖}}_j$,

and output the boundary-adjusted estimator ${\hat{G}}_{\hat{\Xi }}^{†}$ and ${\hat{g}}_{\hat{\Xi }}^{†}$,

using a combining nonparametric copula estimator ${\hat{C}}_{n,h}$ (see (14)).

Equations

$$\hat{f}(x):=\frac{1}{nh}\sum _{i=1}^{n}K\left(\frac{x-X_i}{h}\right)$$

$$\hat{F}(x):=\frac{1}{n}\sum _{i=1}^{n}W\left(\frac{x-X_i}{h}\right),$$

$$E[\hat{f}(x)]=\frac{1}{h}{\int }_{-\infty }^{\infty }K\left(\frac{x-y}{h}\right)f(y)={\int }_{-\infty }^{\infty }K\left(z\right)f(x-hz)=f(x)+O(h^2)$$

$$V[\hat{f}(x)]=\frac{f(x)}{nh}+O\left(\frac{1}{n}\right),$$

$$E[\hat{F}(x)]=F(x)+O(h^2),V[\hat{F}(x)]=\frac{F(x)}{n}+O\left(\frac{h}{n}\right)$$

$$Bias[\hat{f}(u_0)]$$

$$Bias[\hat{F}(x)]=O(h)\text{uniformly}.$$

$${\hat{f}}^{[R]}(x):=\hat{f}(x)+\hat{f}(2u_0-x)+\hat{f}(2l_0-x)$$

$$=$$

$$Bias[{\hat{f}}^{[R]}(x)]=O(h)=:b_f^{[R]}(x)h\text{(uniformly)},$$

$${\hat{F}}^{[R]}(x):={\int }_{l_0}^x{\hat{f}}^{[R]}(y)dy$$

$$=$$

$$\displaystyle~{}~{}~{}~{}~{}~{}-W\left(\frac{X_{i}+l_{0}-2u_{0}}{h}\right)\bigg{\}}-1$$

$$=$$

$$Bias[{\hat{F}}^{[R]}(x)]=O(h^2)=:b_F^{[R]}(x)h^2\text{(uniformly for x∈R)},$$

$${\hat{F}}^{[BK]}(x)=\{\begin{array}{ll}{\hat{F}}^{[BK,l_0]}(x):=\frac{1}{n}\sum _{i=1}^{n}W\left(\frac{x-X_i}{x-l_0}\right)& (l_0\le x<l_0+h)\\ {\hat{F}}^{[BK,I]}(x):=\hat{F}(x)& (l_0+h\le x<u-h)\\ {\hat{F}}^{[BK,u_0]}(x):=1-\frac{1}{n}\sum _{i=1}^{n}W\left(\frac{X_i-x}{u_0-x}\right)& (u_0-h\le x<u_0)\\ 0,1& (x<l_0),(u_0\le x).\end{array}$$

$$Bias[{\hat{F}}^{[BK]}(x)]=O(h^2)=:b_F^{[BK]}(x)h^2\text{(uniformly)}$$

$$Bias[{\hat{f}}^{[BK]}(x)]=O(h)=:b_f^{[BK]}(x)h\text{(uniformly)}.$$

$$\displaystyle\widehat{\mbox{\boldmath$F$}}_{\mbox{\boldmath$u$}}^{\dagger}(\mbox{\boldmath$X$}_{(1,n)})-\mbox{\boldmath$c$}_{n}=\mbox{\boldmath$0$}$$

$$\displaystyle\widehat{\mbox{\boldmath$F$}}_{\mbox{\boldmath$u$}}^{\dagger}((x_{1},x_{2})^{T}):=\left(\widehat{F}_{\mbox{\boldmath$u$}}^{\dagger}(x_{1}),\widehat{F}_{\mbox{\boldmath$u$}}^{\dagger}(x_{2})\right)^{T}=\left(\int_{-\infty}^{x_{1}}\widehat{f}_{\mbox{\boldmath$u$}}^{\dagger}(z)dz~{},~{}\int_{-\infty}^{x_{2}}\widehat{f}_{\mbox{\boldmath$u$}}^{\dagger}(z)dz\right)^{T}$$

$$\displaystyle\mbox{\boldmath$X$}_{(1,n)}=\left(X_{(1)},X_{(n)}\right)^{T},~{}~{}~{}\mbox{\boldmath$c$}_{n}=\left(\frac{1}{n+1}~{},~{}\frac{n}{n+1}\right)^{T},~{}~{}~{}\mbox{\boldmath$0$}=(0,0)^{T}.$$

$$\widehat{\mbox{\boldmath$F$}}_{{\mbox{\boldmath$u$}}}^{\dagger}(\mbox{\boldmath$X$}_{(1,n)})-\mbox{\boldmath$c$}=\mbox{\boldmath$0$},$$

$$\displaystyle E_{X_{n+1}}\left[\Bigm{(}\widehat{F}_{\mbox{\boldmath$u$}}^{\dagger}(X_{n+1})-F(X_{n+1})\Bigm{)}^{2}I(X_{n+1}\in\mbox{\boldmath$q$}_{n})\Bigm{|}X_{1},\cdots,X_{n}\right]$$

$$=$$

$$0$$

$$\approx$$

$$\mbox{\boldmath$\Psi$}_{\mbox{\boldmath$u$},n}(\mbox{\boldmath$x$})=\frac{1}{n}\sum_{i=1}^{n}\mbox{\boldmath$\psi$}_{i,\mbox{\boldmath$u$},n}(\mbox{\boldmath$x$})$$

$$\mbox{\boldmath$\Psi$}_{\mbox{\boldmath$u$},n}(\mbox{\boldmath$x$}):=(\Psi_{\mbox{\boldmath$u$},n}(x_{1}),\Psi_{\mbox{\boldmath$u$},n}(x_{2}))^{T}$$

$$\mbox{\boldmath$\psi$}_{i,\mbox{\boldmath$u$},n}(\mbox{\boldmath$x$}):=(\psi_{i,\mbox{\boldmath$u$},n}(x_{1}),\psi_{i,\mbox{\boldmath$u$},n}(x_{2}))^{T}.$$

$$\displaystyle\widehat{\mbox{\boldmath$F$}}_{\mbox{\boldmath$u$}}^{\dagger}(\mbox{\boldmath$x$})=\mbox{\boldmath$\Psi$}_{\mbox{\boldmath$u$},n}(\mbox{\boldmath$x$})+\mbox{\boldmath$R$}_{\mbox{\boldmath$u$},n}(\mbox{\boldmath$x$}),$$

$$\mbox{\boldmath$R$}_{\mbox{\boldmath$u$},n}(\mbox{\boldmath$x$})=o_{P}(n^{-1/2}\mbox{\boldmath$1$}),~{}~{}~{}\frac{\partial}{\partial\mbox{\boldmath$u$}}\mbox{\boldmath$R$}_{\mbox{\boldmath$u$},n}(\mbox{\boldmath$x$})=o_{P}(n^{-1/2}\mbox{\boldmath$1$})$$

$$\displaystyle\sup_{\mbox{\boldmath$u$}\in\mbox{\boldmath$\Theta$}}\Biggm{\|}\mbox{\boldmath$\Psi$}_{\mbox{\boldmath$u$},n}(\mbox{\boldmath$u$}_{0})-\mbox{\boldmath$F$}_{\mbox{\boldmath$u$}}(\mbox{\boldmath$u$}_{0})\Biggm{\|}\xrightarrow{p}0,$$

$$\frac{\partial}{\partial x}\Psi_{\mbox{\boldmath$u$}_{0},n}(x_{0})=\frac{\partial}{\partial x}\Psi_{\mbox{\boldmath$u$}_{0},n}(x)\biggm{|}_{x=x_{0}}.$$

Full text

The proposed method detects the boundary and gives the modified density estimator simultaneously. Although it is natural to estimate the support by using the sample maximum (and minimum) and modify the density estimator as in Hall & Park (2002), the new method is numerically superior in the sense of an integrated squared error in the boundary region. Moreover, we discuss an extension to a simple multivariate case and propose a new method for estimating the joint probability density. Using the idea of nonparametric copula estimation, this method combines the marginal densities estimated by the proposed single variable method. The obtained joint density estimator is also boundary bias free.

Keywords: Boundary bias; kernel estimator; integrated squared error; support estimator

1 Introduction

Random phenomena are described uniquely by their probability distribution, and estimations of their distributions give us much information. Nonparametric density estimation not only provides a graphical overview of the shape of a distribution, but also enables us to infer a variety of interesting things. These include estimation of functionals of density and statistical testing, for example, goodness of fit with a parametric model, equality of two distributions, symmetry, multi-modality, and so on. Rosenblatt (1956) has proposed a smooth nonparametric estimator $\widehat{f}$ of the density function $f$, and it has been extensively investigated (Tsybakov (2009)). The kernel cumulative distribution estimator is given by integration of $\widehat{f}$; we will begin by introducing that.

Let $X_{1},X_{2},\cdots,X_{n}$ be independently and identically distributed ($i.i.d.$) random variables with a distribution function $F$, and let $f$ be the density function. The kernel density estimator and distribution estimator are given by

[TABLE]

and

[TABLE]

where $K$ is a symmetric kernel function, $W$ is the integral of $K$ and $h$ is a bandwidth which satisfies $h\rightarrow 0$ and $nh\rightarrow\infty~{}(n\rightarrow\infty)$. We call $\widehat{f}$ the ‘naive density estimator’ and $\widehat{F}$ similarly.

Since the naive kernel estimators are sums of $i.i.d.$ random variables, it is easy to obtain their moments. By making a change of variables and performing a Taylor expansion, if the support of $f$ is all of $\mathbb{R}$, we find that

[TABLE]

and similarly

[TABLE]

[TABLE]

under certain regularity conditions. However, when the support is not all of $\mathbb{R}$, the moments changes, and $\widehat{f}$ loses consistency near the boundary of the support. This situation is known as the ‘boundary bias’ problem. When the support is $(-\infty,u_{0}]$ (i.e. $f(u_{0})>0$), the bias of $\widehat{f}$ at the boundary is

[TABLE]

We call the area $(u_{0}-h,u_{0}]$, where $\widehat{f}$ is biased, the boundary region of the density $f$. In addition, the order of the asymptotic bias of $\widehat{F}$ becomes larger: for $x\in\mathbb{R}$,

[TABLE]

Remark 1

In fact, the boundary bias depends on the left limit $f(u_{0}-)$ rather than $f(u_{0})$. Therefore, it is not a problem of whether $f$ is left or right continuous at $u_{0}$. Hereafter, we will assume that $f(u_{0}-)=f(u_{0})$ and that $f$ is right continuous at the left endpoint $l_{0}$ $($i.e. f(l_{0}+)=f(l_{0})$$) if $l_{0}$ is bounded.

If we know the support of $f$ exactly, we can reduce the bias by using a boundary bias reduction method. The basic reduction methods of $\widehat{f}$ include renormalization (Jones (1993)), reflection (Ćwik & Mielniczuk (1993)), asymmetric kernel (Chen (1999), Chen (2000), etc.) and generalized jackknife (Jones (1993), Terrell & Scott(1980)), the ‘direct’ method for reducing to an arbitrary order (Bearse & Rilstone (2009)), and so on. It seems to be easy to obtain boundary bias free estimators of the cumulative distribution $F$ by integration of the modified kernel density estimators. However, the modified density estimators are not usually $i.i.d.$ sums, and it is often difficult to represent the integral explicitly. That is why it is hard obtain the modified cumulative distribution estimators. Recently, some papers were published that focus on distribution function estimation. The boundary kernel method (Tenreiro(2013)) and generalized reflection method (Kolacek & Karunamuni (2011)) can give boundary bias free estimators of $F$. We should note that all of them control the boundary effect which comes from the ‘known’ support.

When the support is unknown, however, we do not take care with the boundary problem well. This may be because there are almost no papers tackling the unknown boundary effect (except Hall & Park(2002)). In most ‘actual’ cases, the realized values are not large, and so the support must be smaller than $\mathbb{R}$. Therefore, we insist on the necessity of both estimating the support and eliminating the unknown boundary bias appropriately.

To overcome the unknown boundary effect, one can estimate the support and then modify the density estimator which regards the estimated support as true. In fact, Hall & Park(2002) proposed to replace the unknown upper bound of the support with the sample maximum $X_{(n)}$. We call the modified density estimator an ‘$X_{(n)}$-based estimator’. However, it seems obvious that the boundary estimator $X_{(n)}$ is not always numerically best for a density estimation that utilizes a very different boundary bias reduction method. Here, we propose a new method for estimating the probability density and its support simultaneously. The boundary estimator depends on the boundary bias reduction method, which is what we apply, and the proposed density estimator is numerically more accurate in the boundary region. In fact, the new method minimizes a loss function asymptotically.

Section 2 describes some of the basic boundary bias reduction methods for when the support is known and simple, such as $[0,\infty)$ and $[0,1]$. In Section 3, we describe the new method for estimating the population density $f$ which is free from the ‘unknown’ boundary bias and investigate its asymptotic properties. In addition, we confirm that some boundary bias reduction methods satisfy a condition which ensures the proposed estimator works well. In Section 4, we compare the proposed estimator with the naive kernel density estimator $\widehat{f}$ and the $X_{(n)}$-based estimator in the sense of the integrated squared error in the boundary region. When $u_{0}$ is bounded and $f(u_{0})>0$ holds, $X_{(n)}$ has $n$ consistency, while the new boundary estimator has $\sqrt{n}$ consistency. However, the convergence rate does not affect the integrated error asymptotically; in fact, we demonstrate that the new density estimator performs better numerically. Moreover, in section 5, we discuss an extension to the simple multivariate case in which the support is given by an hyper-rectangle. Using the idea of nonparametric copula estimation, we propose a new method for estimating joint density estimator. This method combines the marginals estimated by the single variable method (one-dimensional cases), and the obtained joint density estimator is boundary bias free. We study the proposed method by simulating a number of bivariate distribution estimations. The proofs are given in the appendices.

2 Boundary bias reduction methods

Let us assume that the support of the density function $f$ is known and given by $[l_{0},u_{0}]$. Ćwik & Mielniczuk (1993) discussed a ‘reflection’ method that reduces the boundary bias of the naive kernel density estimator. The estimator is given by

[TABLE]

for $x\in[l_{0},u_{0}]$. The asymptotic bias is as follows: for $x\in\mathbb{R}$,

[TABLE]

where $b_{f}^{[R]}(x)=O(1)$, and the asymptotic order of the variance is the same as that of the naive estimator $\widehat{f}$. Thus, the estimator recovers its consistency in the boundary region. In addition, if $f^{\prime}(u_{0}-)=0$, the bias becomes of order $h^{2}$ uniformly, where $f^{\prime}$ is the derivative of $f$. From the integral of $\widehat{f}^{[R]}(x)$, we can derive the distribution estimator as follows:

[TABLE]

From Ćwik & Mielniczuk (1993, we can see

[TABLE]

where $b_{F}^{[R]}(x)=O(1)$, and the order of the variance is the same as that of the naive estimator.

Tenreiro (2013) proposed the ‘boundary kernel’ method for reducing the bias of the kernel cumulative distribution estimator. The method changes the kernel function in the boundary region. It has a simple form, given as follows:

[TABLE]

The author of that paper shows that for $x\in\mathbb{R}$,

[TABLE]

when the support of $K$ is bounded. We can also derive the following density estimator $\widehat{f}^{[BK]}$ from the derivative of $\widehat{F}^{[BK]}$ and calculate the following bias,

[TABLE]

If $f^{\prime}(u_{0}-)$ is zero, the bias becomes of order $h^{2}$.

3 New joint estimator of probability density and its support in one dimension

Let us assume that both boundaries \mbox{\boldmathu}_{0}=(l_{0},u_{0})^{T} are unknown but that both are bounded. We do not assume that the interval is open, half-open, or closed. Now we put \mbox{\boldmathu}=(l,u)^{T}, and \widehat{f}_{\mbox{\boldmathu}}^{\dagger} denotes the kernel type and boundary bias free estimator of the density whose support is $[l,u]$. We propose the following new method for estimating the probability density and its support simultaneously. Let us define the estimator \mbox{\boldmathu}=\widehat{\mbox{\boldmathu}} of \mbox{\boldmathu}_{0} as the solution of the following function,

[TABLE]

and the density estimator as \widehat{f}_{\widehat{\mbox{\boldmathu}}}^{\dagger} (and the distribution estimator as \widehat{F}_{\widehat{\mbox{\boldmathu}}}^{\dagger}), where

[TABLE]

and

[TABLE]

In this section, all functions in bold-type denote two-dimensional vectors of ‘same’ scalar-valued functions for (e.g., \mbox{\boldmathF}(\mbox{\boldmathx})=(F(x_{1}),F(x_{2}))^{T}). Note that \mbox{\boldmathu}_{0} stands for the true value and $u$ is variable. Intuitively, equation $(1)$ comes from the properties of the maximum and minimum estimation of the uniform distribution on the interval $[0,1]$. This is because \mbox{\boldmathF}((X_{1},X_{2})^{T})\stackrel{{\scriptstyle d}}{{=}}\mbox{\boldmathZ} (in distribution) and \mbox{\boldmathF}(\mbox{\boldmathX}_{(1,n)})\approx\mbox{\boldmathc}_{n} hold, where $Z$ is the two-dimensional uniform random variable on $[0,1]\times[0,1]$. Moreover, $(1+(1/n))Z_{(n)}$ (the sample maximum of the uniform distribution) is known as the minimum-variance unbiased estimator of the boundary value ($u_{0}=1$) of the uniform distribution (the sample minimum is similar).

Remark 2

The solution \mbox{\boldmathu}=\mbox{\boldmathX}_{(1,n)} (that is, \mbox{\boldmathX}_{(1,n)}-based estimator) can be viewed as the solution of the following equation,

[TABLE]

where \mbox{\boldmathc}=(0,1)^{T}.

In fact, the new estimator asymptotically coincides with the minimizer of a local expected predict error. Let us define the following error as the loss function

[TABLE]

where \widehat{F}_{\mbox{\boldmathu}}^{\dagger} is based on the sample $(X_{1},\cdots,X_{n})$ and $I$ is the indicator function $I(A)=1$ (if $A$ occurs), $=0$ (if $A$ fails). Let \mbox{\boldmathq}_{n} be some two areas near both boundaries which contain both $X_{(1)}$ and $X_{(n)}$ as the realized value. Then, replacing $f(x)dx$ by the estimator $dF_{n}(x)$ and using the fact that \mbox{\boldmathF}(\mbox{\boldmathX}_{(1,n)})\approx\mbox{\boldmathc}_{n}, we can see that the minimizer of the error asymptotically coincides with \widehat{\mbox{\boldmathu}} as follows:

[TABLE]

The solution \mbox{\boldmathu}=\widehat{\mbox{\boldmathu}} of $(1)$ is not given as an explicit formula, and the properties depend on the bias reduction method. Next, we state the general properties of the proposed estimator and study some applications.

3.1 Asymptotic properties

To construct the asymptotic properties of the new estimators, we utilize the asymptotic theory of $M$ estimation. To view the solution \mbox{\boldmathu}=\widehat{\mbox{\boldmathu}} as an $M$ estimator asymptotically, we make the following assumptions.

Assumption 1

For all large enough integers $n$, there is $\Theta$ whose interior includes \mbox{\boldmathu}_{0}.

Assumption 2

*There is a function \mbox{\boldmath\Psi}_{\mbox{\boldmathu},n} which satisfies the following for all large enough integers $n$:

$(i)$ \mbox{\boldmath\Psi}_{\mbox{\boldmathu},n} is given by following $i.i.d.$ sum*

[TABLE]

where

[TABLE]

and

[TABLE]

$(ii)$* The following stochastic expansion holds:*

[TABLE]

where the residual \mbox{\boldmathR}_{\mbox{\boldmathu},n} satisfies

[TABLE]

*uniformly for \mbox{\boldmathu}\in\mbox{\boldmath\Theta} and \mbox{\boldmathx}\in\mathbb{R}^{2}, and \mbox{\boldmath1}=(1,1)^{T}.

$(iii)$ The next convergence holds*

[TABLE]

*where \mbox{\boldmathF}_{\mbox{\boldmathu}} is a function such that \mbox{\boldmathF}_{\mbox{\boldmathu}_{0}} equals the underlying

  distribution function $F$.*

Assumption 3

*The following holds for all large enough integers $n$:

$(i)$ \centering{\displaystyle\Psi_{\mbox{\boldmathu}_{0},n}(x)=F(x)+b_{F}^{\dagger}(x,n,h)+O_{P}(n^{-1/2})}\@add@centering

  holds uniformly for $x\in\mathbb{R}$, where $b_{F}^{\dagger}$ satisfies b_{F}^{\dagger}(\mbox{\boldmathu}_{0},n,h)=o(n^{-1/2}\mbox{\boldmath1}).

$(ii)$ \displaystyle\frac{\partial}{\partial x}\Psi_{\mbox{\boldmathu}_{0},n}(x)=f(x)+b_{f}^{\dagger}(x,n,h)+O_{P}((nh)^{-1/2})

  holds uniformly for $x\in\mathbb{R}$, where*

[TABLE]

$(iii)$* For every $l~{}(l=1,2)$,*

[TABLE]

holds uniformly for \mbox{\boldmathu}\in\mbox{\boldmath\Theta} and $x\in\mathbb{R}$ (that is, every component converges some uniformly bounded constants), where

[TABLE]

Assumption 4

*For all large enough integers $n$, the following holds:

$(i)$   E[\mbox{\boldmath\psi}_{i,\mbox{\boldmathu},n}(\mbox{\boldmathu}_{0})]-\mbox{\boldmathc}_{n}=\mbox{\boldmath0}  has the unique solution  \mbox{\boldmathu}=\mbox{\boldmathu}_{n}^{*} which satisfies 

    \mbox{\boldmathu}_{n}^{*}\to\mbox{\boldmathu}_{0}.

$(ii)$   for any $\eta>0$, there is $\kappa_{\eta,\delta_{n},n}=O(1)$ which ensures

    {\displaystyle P\left[\inf_{\mbox{\boldmathu}:\|\mbox{\boldmathu}-\mbox{\boldmathu}_{0}\|>\eta}\|\mbox{\boldmath\Psi}_{\mbox{\boldmathu},n}(\mbox{\boldmathu}_{0})-\mbox{\boldmathc}_{n}\|>\kappa_{\eta,\delta_{n},n}\right]>1-\delta_{n}},

    where $\delta_{n}=o(1)>0$.

$(iii)$   {\displaystyle\left(E\left[\frac{\partial}{\partial\mbox{\boldmathu}}\mbox{\boldmath\psi}_{i,\mbox{\boldmathu}_{n}^{*},n}(X_{i},\mbox{\boldmathu}_{0})\right]\right)^{-1}}   exists, i.e., the matrix is nonsingular.

$(iv)$   E[\|\mbox{\boldmath\psi}_{i,\mbox{\boldmathu},n}(\mbox{\boldmathu}_{0})\|^{2}]  is bounded in some neighborhood of  \mbox{\boldmathu}=\mbox{\boldmathu}_{n}^{*}

 and is continuous at  \mbox{\boldmathu}_{n}^{*}.*

Assumptions 2 and 3 admit the following asymptotic expansion of \widehat{\mbox{\boldmathF}}_{\mbox{\boldmathu}}^{\dagger}(\mbox{\boldmathX}_{(1,n)}) at the $i.i.d.$ sum \mbox{\boldmath\Psi}_{\mbox{\boldmathu},n},

[TABLE]

and the boundary bias reduction of the distribution estimators \widehat{F}_{\mbox{\boldmathu}}^{\dagger} and \widehat{f}_{\mbox{\boldmathu}}^{\dagger}. Assumption 4 keeps the $\sqrt{n}$ consistency of \widehat{\mbox{\boldmathu}}, on the basis of the asymptotic theory of $M$ estimation.

Under these assumptions, we can show the consistency of \widehat{\mbox{\boldmathu}} and the bias reduction of \widehat{F}_{\widehat{\mbox{\boldmathu}}}^{\dagger} and \widehat{f}_{\widehat{\mbox{\boldmathu}}}^{\dagger}.

Theorem 1

Given Assumptions 1 - 4 and that

[TABLE]

we have

[TABLE]

Proof. See the Appendices.

Remark 3

When $f(l_{0})>0$ and $f(u_{0})>0$, \mbox{\boldmathX}_{(1,n)} is an $n$ consistent estimator of \mbox{\boldmathu}_{0} and clearly satisfies the assumptions of Theorem 1.

Remark 4

If we know either $l_{0}$ or $u_{0}$, we can see that Theorem 1 still holds by replacing the all vector values with the scalar ones, for example, \mbox{\boldmathu}_{0} with a scalar $u_{0}$ or $l_{0}$, $u$ with a scalar $u$ or $l$, and \mbox{\boldmathX}_{(1,n)} with $X_{(n)}$ or $X_{(1)}$, and so on. Table 1 shows the main procedure of the proposed method if we know the left support $l_{0}$.

Now, we apply the boundary kernel method and reflection method described above. Since P[\|\mbox{\boldmathX}_{(1,n)}-\mbox{\boldmathu}_{0}\|>h]=o((\sqrt{n}h)^{-1}) holds under the assumptions of Theorem 1, we define the boundary estimator \mbox{\boldmathu}=\widehat{\mbox{\boldmathu}}^{[BK]} as the solution of the following equation,

[TABLE]

where

[TABLE]

The vector-valued equation $(11)$ is divided into two scalar-valued equations, such that one includes only $l$, and the other includes only $u$. We can easily see that $\widehat{F}_{l}^{[BK,l]}(X_{(1)})$ is an increasing function for $l\in(-\infty,X_{(1)}]$ and $\widehat{F}_{u}^{[BK,u]}(X_{(n)})$ is a decreasing function for $u\in[X_{(n)},\infty)$. For any fixed $n(\geq 2)$, we have

[TABLE]

and

[TABLE]

Therefore, we can see that the solution $u=\widehat{u}$ is unique and $l=\widehat{l}$ is also unique. Thus, it is easy to confirm that the minimizer \widehat{\mbox{\boldmathu}}^{[BK]} is unique and satisfies $(11)$. Now, we can deduce the following result.

Corollary 1

Let us assume that $supp(K)=[-1,1]$ and $f^{(1)}$ exists and is continuous. When $(7)$ holds and $h=o(n^{-1/4})$, we have

[TABLE]

By applying the reflection method, we define the boundary estimator \mbox{\boldmathu}=\widehat{\mbox{\boldmathu}}^{[R]} and the distribution estimator \widehat{F}_{\widehat{\mbox{\boldmathu}}}^{[R]} in the same way:

[TABLE]

When $n$ is large enough, equation $(12)$ is divided into two mutually independent parts in the sense of $l$ and $u$ as same as the boundary kernel method. It is easy to see that $\widehat{F}_{l}^{[R]}(X_{(1)})$ is an increasing function for $l\in(-\infty,X_{(1)}]$ and that $\widehat{F}_{u}^{[R]}(X_{(n)})$ is a decreasing function for $u\in[X_{(n)},\infty)$. For any fixed $n$, we have

[TABLE]

and

[TABLE]

The equation

[TABLE]

does not always hold, but we can prove the following asymptotic properties under some assumptions. We can get the following result in a similar manner.

Corollary 2

Let us assume that $f^{(1)}$ exists and is continuous and that $(7)$ holds. Then, if Assumption 1 holds and $h=o(n^{-1/4})$, we have

[TABLE]

4 Simulation study in Case 1

We compared the proposed estimator numerically with the naive estimator and $X_{(n)}$-based estimator. Hall & Park(2002) also discusses bias reduction of $X_{(n)}$ in case $f(u_{0})>0$; however, the modification does not affect the first-order asymptotics of the density estimation since $X_{(n)}$ is $n$ consistent in such case. Therefore, we did not take it into account here. Table 2 shows the averaged values of $ISE$, defined as

[TABLE]

of the density estimator $\widetilde{f}$ in the boundary region $\{x>u_{0}-h\}$. In the table, $Naive$ denotes the $ISE$ of the standard kernel density estimator $\widehat{f}$. $\widehat{u}^{[BK]}$ denotes the $ISE$ of $\widehat{f}_{\widehat{u}}^{[BK]}$ and $X_{(n)}^{[BK]}$ denotes $X_{(n)}$-based estimator whose upper bound of support is given by $X_{(n)}$ using the boundary kernel method. $\widehat{u}^{[R]}$ denotes $\widehat{f}_{\widehat{u}}^{[R]}$ and $X_{(n)}^{[R]}$ is the $X_{(n)}$-based estimator using the boundary kernel method. The number of repetitions was $N=10000$ for all cases. The kernel functions were the Epanechnikov, and the bandwidths were the same and chosen by cross-validation as to which made ‘Naive’ asymptotically best. $\widehat{f}_{\widehat{u}}^{[R]}$, $\widehat{f}$, and $\widehat{f}_{X_{(n)}}^{[R]}$ show the $ISE$ in the same way. The kernel function was a Gaussian.

Remark 5

$f_{Beta(1,1)}(1)>0$, $f_{Beta(3,1)}(1)>0$, $f_{Beta(3,3)}(1)=0$ hold, where $f_{Beta(p,q)}$ denotes the density function of the beta distribution with parameters $(p,q)$. Note that if $f(1)=f^{\prime}(1)=0$, $Bias[\widehat{f}(1)]=O(h^{2})$ holds. Table 3 shows the boundary bias of the kernel estimators.

We can see from these numerical studies that when $f(1)>0$, $\widehat{u}$ converges more slowly than $X_{(n)}$, which has $n$ consistency; however, the proposed $\widehat{f}_{\widehat{u}}^{[\cdot]}$ is mostly better than $\widehat{f}_{X_{(n)}}^{[\cdot]}$. Although the boundary effect is not so great when $f(1)=0$, the proposed estimator is comparable with the other estimators. From the above results, we can claim that the proposed method is better in terms of both the theoretical and numerical local loss, at least when the support seems to be compact. Compared with the reflection method, the boundary kernel method seems to be numerically superior, especially when $f(1)>0$.

5 Extension to a simple multivariate case

5.1 Estimating multivariate joint density with unknown but simple form of support

We want to apply the new method to improve kernel joint probability density estimation, kernel regression and so on; however, it is impossible to apply it to multivariate cases directly. This is because the support of a multivariate joint density is given by a much more general formula and the problem completely changes. In addition, most of the boundary bias reduction methods have only been thoroughly investigated in simple cases, such as $[l_{1},\infty)\times\cdots\times[l_{d},\infty)$, $[0,1]^{d}$ and so on.

Let \mbox{\boldmathX}_{1},\mbox{\boldmathX}_{2},\cdots,\mbox{\boldmathX}_{n} be independently and identically distributed ($i.i.d.$) random variables with distribution function $G$ and density function $g$. \mbox{\boldmathX}_{i}=(X_{1,i},\cdots,X_{d,i})^{T} denotes $d$-dimensional value, and hereafter, we will assume that the support of $g$ is given by an unknown hyper-rectangle $\Xi$.

Assumption 5

There is bounded $l_{1,0}<u_{1,0},\cdots,l_{d,0}<u_{d,0}$ which satisfies

[TABLE]

where $\xi_{j}$ is a bounded and open, half-open, or closed interval from $l_{j,0}$ to $u_{j,0}$.

We can treat the marginal distributions in one dimension, even those whose support is unknown, so let us construct the joint probability density estimator, combining the marginals with the copula $C$ defined as follows:

[TABLE]

where \mbox{\boldmathx}=(x_{1},\cdots,x_{d})^{T}, {\accentset{\diamond}{\mbox{\boldmathG}}}(\mbox{\boldmathx})=(G_{1}(x_{1}),\cdots,G_{d}(x_{d}))^{T} and $G_{j}$ is the marginal distribution of the $j$-th component. The following equivalent form of the copula is known:

[TABLE]

where \mbox{\boldmathv}=(v_{1},\cdots,v_{d})^{T} and the existence and uniqueness is known as Sklar’s theorem. The copula function describes the dependence structure between the marginal variables. Scaillet and Fermanian (2002) proposed the following estimator in the bivariate case ($d=2$):

[TABLE]

where $\widehat{G}_{j}$ is the kernel marginal distribution estimator. In addition, Chen & Huang (2007) discusses reduction of the boundary bias of the copula estimator; however, we can not apply it to the case of $d\geq 3$ because the boundary problem gets rather complicated. Instead, using an appropriate estimator $\widehat{C}_{n,h}$ and the idea of copula, we define the following general form of the joint probability density estimator free from the boundary bias

[TABLE]

where

[TABLE]

\widehat{G}_{j,\widehat{\mbox{\boldmathu}}_{j}}^{\dagger} is the proposed estimator of the marginal distribution $G_{j}$ as given in section 3 and $\widehat{C}_{n,h}$ satisfies the following assumption.

Assumption 6

For every $l_{1}=1,2$ and $l_{2}=1,2$, $\widehat{C}_{n,h}^{(l_{1},l_{2})}$ exists and is continuous. In addition, all of the following hold for any \mbox{\boldmathx}\in\mathbb{R}^{d}:

[TABLE]

Remark 6

The following product-type of the boundary kernel and the reflection estimators of the marginals satisfy Assumption 6 and give an appropriate distribution estimation,

[TABLE]

where n^{-1}\sum_{i=1}^{n}W_{\mbox{\boldmathu}_{j,0}}^{[\cdot]}(x_{j},X_{j,i},n,h)=\widehat{G}_{\mbox{\boldmathu}_{j,0}}^{[\cdot]}(x_{j}) is the marginal distribution estimator of $G_{j}$ and ‘$\cdot$’ is ‘$BK$’ or ‘$R$’.

The critical part of the modification $(14)$ is changing the unknown support (hypercube) to $[0,1]^{d}$. Under the assumption of support, we can remove the boundary bias.

Theorem 2

Let us assume that for every marginal distribution $G_{j}$ $(j=1,\cdots,d)$, Assumptions 1 - 4 hold. In addition, we assume that for every $j$,

[TABLE]

holds, where \mbox{\boldmathX}_{j,(1,n)}=(X_{j,(1)},X_{j,(n)})^{T} and \mbox{\boldmathu}_{j,0}=(l_{j,0},u_{j,0})^{T}. Under Assumptions 5 and 6, we have

[TABLE]

where

[TABLE]

b_{G_{j}}^{\dagger}(\mbox{\boldmathx},n,h)* is the bias of the marginal distribution estimation \widehat{G}_{{\mbox{\boldmathu}}_{j},j}^{\dagger}(\mbox{\boldmathx}) and B_{g}^{\dagger}(\mbox{\boldmathx},n,h) is defined similarly.*

Proof. See the appendices.

Table 4 summarizes the multivariate version of the proposed method.

Acknowledgement

The author thanks Prof. Y. Maesono., Faculty of Mathematics, Kyushu University for his valuable comments.

6 Appendices: Some proofs

Proof of Theorem 1Under the assumptions of Theorem 1, we can prove that the following equation holds by using the asymptotic expansions:

[TABLE]

Therefore, \mbox{\boldmathu}=\widehat{\mbox{\boldmathu}} can be viewed as an $M$ estimator, and we can use asymptotic theory. Using assumptions $(i)$ and $(iv)$ of Assumption 4, we find that

[TABLE]

and

[TABLE]

Combining them and using assumption $(ii)$ of Assumption 4 for any $\eta>0$, we have

[TABLE]

Next, let us expand $\widehat{u}$ around $u_{0}$ in the estimating equation,

[TABLE]

where \widetilde{\mbox{\boldmathu}} is a random variable between \widehat{\mbox{\boldmathu}} and \mbox{\boldmathu}_{0}. Accordingly, we find that

[TABLE]

where

[TABLE]

holds. Therefore, we obtain

[TABLE]

and thus

[TABLE]

Next, expanding \widehat{\mbox{\boldmathu}} around \mbox{\boldmathu}_{0}, we have

[TABLE]

The expectation is given by

[TABLE]

To derive the variance term, we perform another asymptotic expansion as follows:

[TABLE]

Then, we have the following four results:

[TABLE]

From the above results, we can see that

[TABLE]

and that

[TABLE]

In the same way, we can prove that

[TABLE]

Proof of Theorem 2 From Theorem 1, the following asymptotic expansion holds under the assumptions:

[TABLE]

Then, it is easy to see that Theorem 2 follows.