Estimating the Major Cluster by Mean-Shift with Updating Kernel

Abstract

The mean-shift method is a convenient mode-seeking method. Using a principle of the sample mean over an analysis window, or kernel, in a data space where samples are distributed with bias toward the densest direction of sample from the kernel center, the mean-shift method is an attempt to seek the densest point of samples, or the sample mode, iteratively. A smaller kernel leads to convergence to a local mode that appears because of statistical fluctuation. A larger kernel leads to estimation of a biased mode affected by other clusters, abnormal values, or outliers if they exist other than in the major cluster. Therefore, optimal selection of the kernel size, which is designated as the bandwidth in many reports of the literature, represents an important problem. As described herein, assuming that the major cluster follows a Gaussian probability density distribution, and, assuming that the outliers do not affect the sample mode of the major cluster, and, by adopting a Gaussian kernel, we propose a new mean-shift by which both the mean vector and covariance matrix of the major cluster are estimated in each iteration. Subsequently, the kernel size and shape are updated adaptively. Numerical experiments indicate that the mean vector, covariance matrix, and the number of samples of the major cluster can be estimated stably. Because the kernel shape can be adjusted not only to an isotropic shape but also to an anisotropic shape according to the sample distribution, the proposed method has higher estimation precision than the general mean-shift.

1. Introduction

When measuring a certain physical quantity, a few abnormal values, hereinafter designated as outliers, are included among the normal measured values, thereby exacerbating measurement noise. Frequently in science and engineering, some effort is necessary to estimate the true statistical parameters of the physical quantity from these measured values and the included outliers. Because the measurement noise generally follows a Gaussian distribution with mean zero, all samples from the major cluster are Gaussian-distributed around the true value. The problem described above is summarized to estimate the parameters of the major cluster, such as the mean, covariance matrix, and the number of samples included in the major cluster.

Because the mean equals the mode in a Gaussian distribution, if the outliers do not affect the sample mode of the major cluster, then the problem above can be replaced by a mode-seeking problem of the major cluster. Fukunaga and Hostetler [1] first proposed the mean-shift method, which was subsequently generalized by Cheng [2]. It is therefore known as a convenient iterative method for mode-seeking. The mean-shift was shown to be equivalent to the method that seeks a local maximum by the steepest gradient algorithm for the probability density distribution estimated using the kernel method [3,4]. Therefore, the bandwidth, which is the size of the used kernel, deeply affects both the estimation accuracy and precision in the mean-shift as well as in kernel density estimation [5].

Usually in kernel density estimation, the bandwidth is determined such that the difference between the true distribution and the estimated distribution is minimized [6,7,8]. In mean-shift, because the normalized norm affects the convergence speed, a method for determining the bandwidth is proposed for the isotropic kernel [9] and anisotropic kernel [10] such that the norm of the mean-shift vector normalized by the bandwidth is maximized. A method for selecting the most stable bandwidth was also proposed [10,11]. Moreover, mean-shift with bandwidth that varies depending on the coordinate in data space was proposed [9,11]. Nevertheless, these methods entail high calculation costs because they require some provisional estimate of the probability density distribution, which is described as the pilot or initial estimate in some reports of the literature. Other theoretical studies of mean-shift, such as convergence, have been further proven. Li [12] proved its convergence by further imposing some commonly acceptable conditions. Ghassabeh [13] modified the mean-shift to guarantee its convergence. Although the mean-shift has been used widely in many applications [14,15,16], the use of bandwidth for mean-shift has been largely ignored in studies reported in the literature.

As described herein, we propose a new mean-shift method by which adopting the multi-dimensional Gaussian kernel, the kernel bandwidth and shape are updated to fit the major cluster size and shape in each iteration with no provisional estimation. We first derive a calculation equation for calculating the variance (or covariance matrix) of a major cluster from the sample variance in the kernel (or the sample covariance matrix in the multi-dimensional case) around the mode. Then, as the update progresses in the mean-shift method, the variance (or covariance matrix) of a major cluster is estimated using this calculation equation. In addition, the kernel bandwidth and shape are adjusted adaptively based on this estimated value. Therefore, we propose the mean-shift method with such an updating kernel. The proposed mean-shift requires no predetermination of the kernel bandwidth as necessitated by the general mean-shift method.

This paper is organized as follows. A general mean-shift method is introduced in Section 2. In Section 3, we propose the new mean-shift method in a one-dimensional case. Numerical experiments are presented to evaluate the proposed mean-shift compared to the general mean-shift method in Section 4. An explanation of applications and conclusion are presented respectively in Section 5 and Section 6. In the appendices, we describe an extension of the general mean-shift method and the proposed mean-shift to a multi-dimensional case.

2. General Mean-Shift Method

2.1. General Mean-Shift Method

Assuming that the major cluster of $N_N$ points follows a Gaussian distribution with mean ${\mu }_N$ and standard deviation ${\sigma }_N$, we are considering the problem of estimating the mean ${\mu }_N$ of the major cluster when a fewer outliers of $N_O$ points exist in the sample of $N=N_N+N_O$ points. If the mode of the sample is not biased from the mean ${\mu }_N$ under the influence of outliers, then the mean ${\mu }_N$ can be estimated as the mode. The mean-shift is a simple and iterative method to estimate the mode of the major cluster. Letting the sample be $x_n,n=1,\cdots ,N$, then the general mean-shift method is realized using the following iterative process:

• Letting the mean ${\mu }_x$ of sample $x_n,n=1,\cdots ,N$ be the initial value of the mean estimator ${\widehat{\mu }}_N$ of major cluster, then

  $${\widehat{\mu }}_N\leftarrow {\mu }_x.$$

  (1)
• Consider a Gaussian distribution $p(x;{\mu }_W,{\sigma }_W)$ with the mean ${\mu }_W$ and standard deviation ${\sigma }_W$ as the kernel function in the value direction. Here, the mean ${\mu }_W$ of kernel function is found by the mean estimator of major cluster

  $${\mu }_W\leftarrow {\widehat{\mu }}_N.$$

  (2)
  The standard deviation ${\sigma }_W$ is assigned to be an appropriate size as discussed later in Section 2.2.
• Weight $a_n,n=1,\cdots ,N$ for each sample $x_n,n=1,\cdots ,N$ weighted by such a Gaussian kernel is

  $${\displaystyle a_n=\frac{1}{A}p(x_n;{\mu }_W,{\sigma }_W).}$$

  (3)
  However, A in Equation (3) above is a normalization coefficient for which the sum of the weight $a_n$ is equal to 1, as

  $${\displaystyle A=\sum _{k=1}^{N}p(x_k;{\mu }_W,{\sigma }_W).}$$

  (4)
  We use this weight $a_n$ to calculate the sample mean ${\mu }_x$ with $x_n,n=1,\cdots ,N$ as

  $${\displaystyle {\mu }_x=\sum _{n=1}^N{a}_n{x}_n.}$$

  (5)
• The value of mean estimator ${\widehat{\mu }}_N$ of the major cluster is updated by the following equation:

  $${\widehat{\mu }}_N\leftarrow {\mu }_x.$$

  (6)
• If the variation of the value of mean estimator ${\widehat{\mu }}_N$ is equal to or less than the predetermined fixed value, then the update process is terminated. Otherwise, return to 2 and repeat the iteration.

2.2. Shortcomings and Solution of the General Mean-Shift Method

The general mean-shift method estimates the modes of the underlying probability density function. From the definition of a probability density, if the random variable X of N data points $x_i,i=1,2,3,\cdots ,N$ in one-dimensional space R has density f, then

$${\displaystyle f\left(x\right)=\underset{h\to 0}{lim}\frac{1}{2h}P(x-h<X<x+h).}$$

(7)

For any given h (bin bandwidth or kernel bandwidth), we can estimate $P(x-h<X<x+h)$ by the proportion of the sample falling in the interval $(x-h,x+h)$. Thus, a natural estimator $\widehat{f}$ of the density is given by choosing a small h and setting

$${\displaystyle \widehat{f}\left(x\right)=\frac{1}{2h}\frac{N_x}{N}.}$$

(8)

Here, $N_x$ denotes the number of samples falling in the interval $(x-h,x+h)$. To express the estimator more transparently, define the weight function $\omega (x;h)$ by

$$\omega (x;h)=\left\{\begin{array}{cc}\frac{1}{2h}\hfill & \left|x\right|<h,\hfill \\ 0\hfill & \mathrm{others}.\hfill \end{array}\right.$$

(9)

The estimator can be expressed as below [17]:

$${\displaystyle \widehat{f}\left(x\right)=\frac{1}{N}\sum _{i=1}^N\omega (x-x_i;h).}$$

(10)

Replace the weight function $\omega$ by a general kernel function $K(x;\sigma )$ with standard deviation $\sigma$, which satisfies the condition

$${\displaystyle {\int }_{-\infty }^{\infty }K\left(x\right)dx=1,}$$

(11)

and the kernel estimator for the probability density function $\widehat{f}\left(x\right)$ at point x can be expressed as

$${\displaystyle \widehat{f}\left(x\right)=\frac{1}{N}\sum _{i=1}^{N}K(x-x_i;\sigma ).}$$

(12)

The general mean-shift is an attempt to ascertain the local modes of density function $\widehat{f}\left(x\right)$, which correspond to the zeros of the gradient ${▽}_x\widehat{f}\left(x\right)=0$. Therefore, the type of kernel function $K(x;\sigma )$ and the kernel bandwidth $\sigma$ both directly affect the performance of general mean-shift method. Fixing the type of kernel function to Gaussian kernel, we specifically examine the influence of the pre-set of the kernel bandwidth in general mean-shift.

To confirm the influence of fixed kernel bandwidth on estimation accuracy in a general mean-shift method, we set various fixed kernel bandwidths in advance. Here, we summarize the numerical and experimentally obtained results for general mean-shift method as discussed in Section 4. Figure 1a presents the bias error between the estimated value in a general mean-shift method and the true value when we select various kernel bandwidths in advance. The horizontal axis shows a selection of different kernel bandwidths. The vertical axes respectively show the bias error between the estimated value for the mean and the true mean value, and the variance of the mean value. While selecting different fixed kernel bandwidths, we estimated the mean of the major cluster, which is distributed as shown in Figure 2 for 1000 trials. Furthermore, we computed the bias errors using the equation described in Section 4.2. Figure 1a shows that, when we enlarge the fixed kernel bandwidth, the mean estimator is more susceptible to outliers. The bias error in general mean-shift method increases. Otherwise, when we decrease the kernel bandwidth, the number of samples involved in the mean estimation decreases. The local mode can easily become the convergence point of the iterative process. In addition, the bias error in general mean-shift method increases. The kernel bandwidth should be set in the range of 0.5–1.5. As shown in Figure 1b, with enlargement of the kernel bandwidth, the estimation variance in general mean-shift method decreases. Therefore, the optimal kernel bandwidth is 1.5. Because the maximum value of these variances is very small and, because it does not exceed 0.06, if we select the kernel bandwidth within this range of 0.5–1.5, we can ensure the unbiasedness and consistency of the mean estimator in general mean-shift method. However, not knowing the true mean of the major cluster beforehand, we cannot calculate the bias error in general mean-shift method. Therefore, we cannot choose the appropriate kernel bandwidth based on the comparison result shown in Figure 1a. Indeed, the proper pre-set of the kernel bandwidth constitutes an important difficulty.

The optimal kernel bandwidth depends on the existence range of outliers, the number of samples belonging to the major cluster and the distribution that the major cluster follows. In the absence of prior knowledge, the kernel bandwidth is often fixed as appropriate to 1/2 the time of the standard deviation of the whole sample when the whole sample contains the major cluster and the outliers in signal processing [18]. For clustering in image processing or other multiple applications, it is still difficult to preset the kernel bandwidth properly in a general mean-shift method. When the kernel bandwidth is inappropriate, the kernel bandwidth becomes a factor that degrades the estimation accuracy.

As follows, based on the general mean-shift method, we propose a method to change the kernel bandwidth adaptively in accordance with simultaneous estimation of the mean (for a multi-dimensional case, the mean vector) and the standard deviation (for a multi-dimensional case, the covariance matrix) of a major cluster at each iteration. We need not set the kernel bandwidth properly in advance.

3. One-Dimensional Mean-Shift with Updating Kernel

3.1. Derivation of Major Cluster Standard Deviation ${\sigma }_N$ from Sample Standard Deviation ${\sigma }_x$

Here, the Gaussian distribution with mean $\mu$ and standard deviation $\sigma$ is represented by $p(x;\mu ,\sigma )$. It is abbreviated as $p(x;\sigma )$ especially for $\mu =0$. We use the two following equations for the two Gaussian distributions:

$${\displaystyle {\int }_{-\infty }^{\infty }p(x;{\sigma }_W)p(x;{\sigma }_N)dx=\frac{1}{\sqrt{2\pi }\sqrt{{\sigma }_W^2+{\sigma }_N^2}},}$$

(13)

$${\displaystyle {\int }_{-\infty }^{\infty }{x}^{2}p(x;{\sigma }_W)p(x;{\sigma }_N)dx=\frac{{\sigma }_W^2{\sigma }_N^2}{\sqrt{2\pi }{({\sigma }_W^2+{\sigma }_N^2)}^{\frac{3}{2}}}.}$$

(14)

We assume that the influence of outliers is small such that the sample mode is not biased from the mean ${\mu }_N$. If the general mean-shift method with the sufficiently small fixed kernel bandwidth decided by the standard deviation of the kernel starts the iteration from an appropriate initial value, then the influence of the outliers on estimation decreases gradually as the estimate converges. Therefore, it is sufficient to consider only the samples from the major cluster $x_n,n=1,\cdots ,N_N$ when the estimate converges to their true value. In addition, the mean ${\mu }_N$ of the major cluster and the mean ${\mu }_W$ of the Gaussian kernel coincide near the convergence point. Even if coordinate transformation is performed so that both are 0, generality is not lost. Therefore, we let ${\mu }_N={\mu }_W=0$ here for analysis. The variance ${\sigma }_x^2$ of the sample $x_n,n=1,\cdots ,N_N$ weighted by $a_n,n=1,\cdots ,N_N$ is

$${\displaystyle {\sigma }_x^2=\sum _{n=1}^{N_N}{a}_n{x}_n^2.}$$

(15)

Weight $a_n$ is a Gaussian kernel given by Equations (3) and (4). In addition, N is replaced by $N_N$.

The expected value of the sample variance ${\sigma }_x^2$ is calculated after substituting Equation (3) into Equation (15) as

$$E\left[{\sigma }_x^2\right]=E\left[{\displaystyle \frac{1}{A}\sum _{n=1}^{N_N}p(x_n;{\sigma }_W)x_n^2}\right].$$

(16)

The variance of ${\displaystyle \frac{1}{A}}$ is sufficiently smaller than the dispersion of other parts. Therefore, it can be approximated to the following equation based on the assumption that the major cluster follows a Gaussian distribution, as

$${\displaystyle E\left[{\sigma }_x^2\right]\simeq \frac{1}{E\left[A\right]}E\left[{\displaystyle \sum _{n=1}^{N_N}p(x_n;{\sigma }_W)x_n^2}\right].}$$

(17)

The approximation is discussed later in Appendix B. Here, we calculate the expected value of A by Equations (4) and (13) as

$$\begin{array}{ccc}\hfill E\left[A\right]& =& E\left[{\displaystyle \sum _{k=1}^{N_N}p(x_k;{\sigma }_W)}\right]\hfill \\ \hfill & =& {\displaystyle \sum _{k=1}^{N_N}E\left[p(x_k;{\sigma }_W)\right]}\hfill \\ \hfill & =& {\displaystyle N_N{\int }_{-\infty }^{\infty }p(x;{\sigma }_W)p(x;{\sigma }_N)dx}\hfill \\ \hfill & =& {\displaystyle \frac{N_N}{\sqrt{2\pi }\sqrt{{\sigma }_W^2+{\sigma }_N^2}}.}\hfill \end{array}$$

(18)

The expected value of other part becomes

$$\begin{array}{l}E\left[{\displaystyle \sum _{n=1}^{N_N}p(x_n;{\sigma }_W)x_n^2}\right]\\ {\displaystyle =\sum _{n=1}^{N_N}E\left[p(x;{\sigma }_W)x^2\right]}\\ {\displaystyle =N_N{\int }_{-\infty }^{\infty }{x}^{2}p(x;{\sigma }_W)p(x;{\sigma }_N)dx}\\ {\displaystyle =\frac{N_N{\sigma }_W^2{\sigma }_N^2}{\sqrt{2\pi }{({\sigma }_W^2+{\sigma }_N^2)}^{3/2}}}\end{array}$$

(19)

according to Equation (14). In other words, after being weighted by a Gaussian kernel with mean 0 and standard deviation ${\sigma }_W$, the expected value of variance ${\sigma }_x^2$ of the sample which follows a Gaussian distribution with mean 0 and standard deviation ${\sigma }_N$ is

$${\displaystyle E\left[{\sigma }_x^2\right]=\frac{{\sigma }_W^2{\sigma }_N^2}{{\sigma }_W^2+{\sigma }_N^2}}$$

(20)

according to Equations (18) and (19). Equation (20) above can be transformed to

$${\displaystyle {\sigma }_N^2=\frac{{\sigma }_W^{2}E\left[{\sigma }_x^2\right]}{{\sigma }_W^2-E\left[{\sigma }_x^2\right]}.}$$

(21)

This expression shows that standard deviation ${\sigma }_N$ can be estimated from the standard deviation ${\sigma }_x$ of the sample, which is weighted using a Gaussian kernel with mean 0 and standard deviation ${\sigma }_W$ as

$${\displaystyle {\widehat{\sigma }}_N=\sqrt{{\displaystyle \frac{{\sigma }_W^2{\sigma }_x^2}{{\sigma }_W^2-{\sigma }_x^2}}}.}$$

(22)

In addition, using Equation (18), the number $N_N$ of samples belonging to the major cluster can be estimated as

$${\widehat{N}}_N=A\sqrt{2\pi }\sqrt{{\sigma }_W^2+{\widehat{\sigma }}_N^2}.$$

(23)

Adaptive change of the standard deviation ${\sigma }_W$ of the kernel related to the estimated value ${\widehat{\sigma }}_N$ of the standard deviation is sufficient for each update because the mean ${\mu }_N$ of the major cluster and the standard deviation ${\sigma }_N$ can also be estimated. Specifically, the standard deviation ${\sigma }_W$ of the kernel is assigned to be r times the estimated value ${\widehat{\sigma }}_N$, although it depends on the existence range of outliers. We designate this r as a scale factor. Regarding appropriate r, we will examine this point in a numerical experiment discussed later.

3.2. Mean-Shift with Updating Kernel

Based on the discussion presented in Section 3.1, at each iteration of the general mean-shift method, the standard deviation ${\sigma }_N$ is estimated simultaneously in addition to the mean value ${\mu }_N$. Therefore, we propose a new mean-shift method that adaptively changes the standard deviation ${\sigma }_W$ of the kernel. The algorithm is summarized as presented below:

• Let the mean ${\mu }_x$ of sample $x_n,n=1,\cdots ,N$ be the initial value of the mean estimator ${\widehat{\mu }}_N$ of the major cluster and let standard deviation ${\sigma }_x$ of this sample be the initial value of the standard deviation estimator ${\sigma }_N$ of the major cluster as

  $$\begin{array}{ccc}\hfill {\widehat{\mu }}_N& \leftarrow & {\mu }_x,\hfill \end{array}$$

  (24)

  $$\begin{array}{ccc}\hfill {\widehat{\sigma }}_N& \leftarrow & {\sigma }_x.\hfill \end{array}$$

  (25)
• Consider a Gaussian distribution $p(x;{\mu }_W,{\sigma }_W)$ with mean ${\mu }_W$ and standard deviation ${\sigma }_W$ as the kernel function in the value direction. Here, the mean ${\mu }_W$ and the standard deviation ${\sigma }_W$ are given respectively by the estimated value ${\widehat{\mu }}_N$ of the mean and the estimated value ${\widehat{\sigma }}_N$ of the standard deviation of the major cluster:

  $$\begin{array}{ccc}\hfill {\mu }_W& \leftarrow & {\widehat{\mu }}_N,\hfill \end{array}$$

  (26)

  $$\begin{array}{ccc}\hfill {\sigma }_W& \leftarrow & r{\widehat{\sigma }}_N.\hfill \end{array}$$

  (27)
  Here, mean ${\mu }_W$ and variance ${\sigma }_W$ of the Gaussian kernel are not estimators, although they change when the kernel updates.
• Weight $a_n,n=1,\cdots ,N$ for each sample $x_n,n=1,\cdots ,N$ weighted by such a Gaussian kernel $p(x;{\mu }_W,{\sigma }_W)$ is calculated using Equations (3) and (4). We use this weight $a_n$ to calculate the sample mean ${\mu }_x$ and standard deviation ${\sigma }_x$ with $x_n,n=1,\cdots ,N$ as shown below:

  $$\begin{array}{ccc}\hfill {\mu }_x& =& {\displaystyle \sum _{n=1}^N{a}_n{x}_n,}\hfill \end{array}$$

  (28)

  $$\begin{array}{ccc}\hfill {\sigma }_x& =& {\displaystyle \sqrt{{\displaystyle \sum _{n=1}^N{a}_n{(x_n-{\mu }_x)}^2}}.}\hfill \end{array}$$

  (29)
• The values of mean estimator ${\widehat{\mu }}_N$, standard deviation estimator ${\widehat{\sigma }}_N$, and number of samples estimator ${\widehat{N}}_N$ of the sample are updated, respectively, by the following equations:

  $$\begin{array}{ccc}\hfill {\widehat{\mu }}_N& \leftarrow & {\mu }_x,\hfill \end{array}$$

  (30)

  $$\begin{array}{ccc}\hfill {\widehat{\sigma }}_N& \leftarrow & {\displaystyle \sqrt{{\displaystyle \frac{{\sigma }_W^2{\sigma }_x^2}{{\sigma }_W^2-{\sigma }_x^2}}},}\hfill \end{array}$$

  (31)

  $$\begin{array}{ccc}\hfill {\widehat{N}}_N& \leftarrow & A\sqrt{2\pi }\sqrt{{\sigma }_W^2+{\widehat{\sigma }}_N^2}.\hfill \end{array}$$

  (32)
• If the variations of the values of these estimators are equal to or less than the predetermined fixed value, then the update process is terminated. Otherwise, return to 2 and repeat the iteration.

4. Numerical Experiment

4.1. Update Process of Mean-Shift with an Updatable Kernel

For the proposed method, we use iteration to confirm the process by which the estimated values of the mean vector, the covariance matrix, and the number of samples converge to true values of the major cluster. Although no restriction is made of the dimension of data to which the proposed method is applicable, to illustrate and explain the distribution of data and update process, two-dimensional data are targeted for analysis. Herein, we obtain the major cluster with $N_N=3000$ points generated in two-dimensional normal distribution with the mean vector ${\mathbf{\mu }}_N={(0,0)}^T$ and variance covariance matrix as

$${\mathit{C}}_N=\left(\begin{array}{cc}3& 2\\ 2& 3\end{array}\right).$$

The outliers with $N_O=200$ points are distributed uniformly within the range of $x_1\in [-2,-1],x_2\in [3,4]$. Figure 2 shows an example of the generated sample in $(x_1,x_2)$ space. Symbol • in the figure represents the coordinates of each point. The points spreading in the central elliptical shape belong to the major cluster. Other points distributed in a square shape on the upper left are outliers. In the figure, the solid ellipse represents a contour line where 99% of the M-dimensional normal distribution defined by the mean vector ${\mathbf{\mu }}_N$ and the covariance matrix ${\mathit{C}}_N$ fall within it. Later, we present the mean vector ${\mathbf{\mu }}_N$ and covariance matrix ${\mathit{C}}_N$, or their estimates.

In general, as discussed in Appendix C.2, the initial estimated value of the mean ${\widehat{\mathbf{\mu }}}_N$ and covariance matrix ${\widehat{\mathit{C}}}_N$ for major cluster can be assigned respectively to the mean and covariance matrix of all samples. However, we set the initial kernel having mean vector ${\widehat{\mathbf{\mu }}}_N={(-2,3)}^T$ and covariance matrix

$${\widehat{\mathit{C}}}_N=\left(\begin{array}{cc}1.25& -0.75\\ -0.75& 1.25\end{array}\right)$$

intentionally to be located and shaped sufficiently apart from the major cluster. To demonstrate how the estimated value converges to the true value with updating, the scale factor is $r=1.0$. The update ends when it satisfies all conditions for which the sum of squares of the change amount ${\widehat{\mathbf{\mu }}}_N$ is 0.01 or fewer, the sum of squares of the change amount of ${\widehat{\mathit{C}}}_N$ is 0.01 or fewer, and the square of the change amount of ${\widehat{N}}_N$ is 30 or less.

As described earlier, the solid ellipse shown in Figure 3 represents the estimated value of mean vector ${\widehat{\mathbf{\mu }}}_N$, covariance matrix ${\widehat{\mathit{C}}}_N$, and number ${\widehat{N}}_N$ of samples for each update in the proposed method. In Figure 3, the estimated values ${\widehat{\mathbf{\mu }}}_N,{\widehat{\mathit{C}}}_N,{\widehat{N}}_N$ are shown to converge to the true values ${\mathbf{\mu }}_N,{\mathit{C}}_N,N_N$ corresponding to Figure 2 as the update progresses, although they start from more or less bad initial values. Here, for the estimated value ${\widehat{N}}_N$, we have accuracy to one decimal place.

4.2. Influence of Kernel Bandwidth on Estimation Accuracy (Unbiasedness)

An exceedingly important property required for estimators is unbiasedness: a property by which the expected value of the estimated value coincides with the true value. If no statistical bias in the estimated value exists, then it represents that the estimation is accurate. Assuming that the parameter is $\theta$, we investigate the unbiasedness of the estimator $\widehat{\theta }$. If parameter $\theta$ is a scalar, then the bias error is the difference $E\left[\widehat{\theta }\right]-\theta$ between the expected value and the true value $\theta$ of the estimator. Otherwise, if parameter $\theta$ is a vector or matrix, then the bias error is the square root $\sqrt{\parallel E\left[\widehat{\mathbf{\theta }}\right]{-\mathbf{\theta }\parallel }^2}$ of the sum of squares over all the elements. It can be evaluated whether the bias error is zero. As explained below, it demonstrates that the initial value of the kernel bandwidth has less influence on the unbiasedness of the estimated value in the proposed method discussed in Appendix C than in the general mean-shift method introduced in Appendix A.

The distributions that major cluster and outliers follow, the numbers of samples $N_N,N_O$, scale factor r, and update ending condition are the same as those described in Section 4.1. The initial estimated value of mean vector ${\widehat{\mathbf{\mu }}}_N$ is the mean vector of all samples. The initial estimated value of covariance matrix is assigned to ${\widehat{\mathit{C}}}_N={\sigma }^2\mathit{I}$. In the general mean-shift method, the covariance matrix of the kernel is ${\mathit{C}}_W={\sigma }^2\mathit{I}$. Under the conditions presented above, the mean vector ${\mathbf{\mu }}_N$, the covariance matrix ${\mathit{C}}_N$, and the number $N_N$ of samples are estimated using the general mean-shift method and the proposed method. In addition, because it is impossible to obtain the expected value in numerical experiments, the expected value is replaced by the average value of the estimated values for 1000 trials that change the random number.

In the proposed method, ${\sigma }^2$ is the initial value of the kernel bandwidth. It corresponds to the pre-set value of the kernel bandwidth in a general mean-shift method. When this ${\sigma }^2$ is changed to various values, the bias errors of the estimated value of the mean vector ${\mathbf{\mu }}_N$, covariance matrix ${\mathit{C}}_N$, and number $N_N$ of samples are calculated. Results are presented respectively in Figure 4a–c. The horizontal axis shows the selection of different kernel bandwidth. The vertical axes respectively show the bias errors for estimators ${\mathbf{\mu }}_N$, ${\mathit{C}}_N$, and $N_N$. In this figure, symbol ◯ corresponds to the proposed method. The symbol △ represents the bias errors in a general mean-shift method. However, because the covariance matrix and number of samples cannot be estimated in a general mean-shift method, only the results obtained using the proposed method are shown in Figure 4b,c. The scale on the vertical axis of the figures is fixed to represent 10% of errors at full scale. In the following figures, the same scale applied to these figures will be used unless specified otherwise.

Figure 4a shows that the bias error increases linearly and that the unbiasedness is lost when the kernel bandwidth ${\sigma }^2$ approximately exceeds the range of 0.5–1.5 represented by symbol ↓ in a general mean-shift method because, as the kernel becomes larger, the outliers fall within the range of the kernel, which greatly affects the mean estimation of the major cluster. For this reason, the proper set of the kernel bandwidth is an important difficulty in a general mean-shift method. However, the kernel bandwidth is adjusted according to the estimated value of covariance matrix of a major cluster at each iteration in the proposed method. Therefore, it is less susceptible to the influence of initial value ${\sigma }^2$. Furthermore, in Figure 4b,c, it is the same situation in the estimations of covariance matrix ${\mathit{C}}_N$ and number $N_N$ of samples.

While maintaining the ratio of the number $N_N$ of samples of major cluster and the number $N_O$ of samples of the outliers to 3000:200 and changing the number $N=N_N+N_O$ of samples from 1000 to 90,000, the variance of each estimate value of the mean vector ${\mathbf{\mu }}_N$, covariance matrix ${\mathit{C}}_N$, and number $N_N$ of samples are obtained using our proposed method, as shown in Figure 5. The horizontal axis shows the selection of different numbers of samples corresponding to the whole samples. The vertical axes respectively represent the bias errors for estimators ${\mathbf{\mu }}_N$, ${\mathit{C}}_N$, and $N_N$. Because the proposed method is independent of the initial value ${\sigma }^2$, the initial value ${\sigma }^2$ is fixed to 1.5. Figure 5 shows that these estimators are unbiased for a finite number of samples.

4.3. Influence of the Scale Factor r Value on Estimation Accuracy

In the proposed method, we need not select the initial value of kernel bandwidth in advance because the kernel bandwidth is changed adaptively. The pre-set of the initial value shows some difficulty in influencing the estimation accuracy. Instead, the problem of optimal setting of the scale factor r occurred. Scale factor r represents the ratio of the kernel bandwidth (standard deviation) to the major cluster width (standard deviation). Therefore, the smaller the scale factor, the smaller the kernel bandwidth (standard deviation) is set with respect to the major cluster width (standard deviation). From the viewpoint of estimation accuracy, the kernel bandwidth (standard deviation) should be sufficiently large but not cover the outliers. In other words, if the outliers exist at the distance from the mode of major clusters more than three times the standard deviation of the major cluster, according to three-sigma rule of thumb, the kernel bandwidth should be the same as the standard deviation of major cluster, which means $r=1$. Otherwise, if there are a certain number of outliers within a standard deviation away from the mode of the major cluster, the kernel bandwidth is expected to be 1/3 of the standard deviation of the major cluster, which means $r=1/3$. If the distribution of the major cluster and the outliers is specified completely, then it is possible to derive the theoretical formula of the optimal scale factor r as a parameter. However, because the purpose is to estimate the distribution of the major cluster and the outliers, then, even if a theoretical formula for scale factor r is derived, it cannot be used for estimation. Derivation of the theoretical formula for scale factor r has no great value. Therefore, as described below, we investigate the influence of the selected value of this scale factor on the estimation accuracy.

The distributions that major cluster and outliers follow, number $N_N,N_O$ of samples, and update ending condition are the same as those in Section 4.1. As shown in Section 4.2, the initial values of the estimated value of mean vector and covariance matrix ${\widehat{\mathbf{\mu }}}_N,{\widehat{\mathit{C}}}_N$ are given, respectively, by the mean vector and covariance matrix of the whole samples. We select scale factor r to be various values and estimate the mean vector ${\mathbf{\mu }}_N$, covariance matrix ${\mathit{C}}_N$, and number $N_N$ of samples using the proposed method. The bias errors of each estimated value is presented in Figure 6a,c. The horizontal axis represents the selection of various scale factors r. The vertical axes respectively represent the bias errors for estimators ${\mathbf{\mu }}_N$, ${\mathit{C}}_N$, and $N_N$.

Figure 6 shows that the bias errors of any estimate increases and that the unbiasedness is lost when scale factor r is selected as a value larger than a certain value because, when the kernel bandwidth increases, it becomes more susceptible to outliers, as with the general mean-shift method shown in Figure 6. However, when scale factor r is selected as a small value, the bias error is increased extremely. The unbiasedness is lost relative to the covariance matrix ${\mathit{C}}_N$ and number $N_N$ of samples, although it is not readily apparent on mean vector ${\mathbf{\mu }}_N$. The reason for this is explainable as presented below.

If we select scale factor r as a value smaller than one, the kernel bandwidth becomes small because of a lack of the practical number of samples that contribute to the estimation. For that reason, the estimation precisions of mean vector ${\mathbf{\mu }}_x$ and standard deviation ${\sigma }_{x,m}$ are deteriorated. The deterioration of this estimation accuracy results from the small number of samples. Consequently, the estimated error has normality, but does not include bias error. As shown in Equation (A32), the estimated value ${\widehat{\mathbf{\mu }}}_N$ of the mean vector is the sample mean vector ${\mathbf{\mu }}_x$. The estimation equation of the standard deviation ${\widehat{\sigma }}_N$ and number ${\widehat{N}}_N$ of samples is a nonlinear function of the sample standard deviation ${\sigma }_{x,m}$, as shown in Equations (A19) and (A20). In general, normality is lost by a nonlinear transformation. Therefore, the estimation errors of both the standard deviation ${\widehat{\sigma }}_N$ and the number ${\widehat{N}}_N$ of samples are converted to the bias errors by the nonlinear transformations, even if the estimation error of the sample standard deviation ${\sigma }_{x,m}$ had normality.

Figure 6 shows that the appropriate value of the scale factor r is in the range of $0.5\le r\le 1.5$, but it depends on the characteristic of the target data. For example, the lower limit increases when the number of samples is small. The upper limit decreases when the outliers approach a major cluster. Comparing the bias error with the general mean-shift indicates that the selection of scale factor r need not be the same as the situation of kernel bandwidth as shown in Figure 4 because the range in which the bias error can be kept low is wide.

4.4. Verification of Consistency

The goodness of the estimator is evaluated by accuracy and precision. Accuracy is evaluated as the bias error, as discussed in Section 4.2, whereas the precision is evaluated by the variance of estimated values. Before comparing the estimation precision of a general mean-shift method with the proposed method, one must confirm the consistency of the estimated values in both methods. Consistency is an important property required for the estimator. It indicates the characteristics by which the variance of the estimated values approaches 0 as the number of samples used for estimation increases.

The distributions that major cluster and outliers follow, in addition to the update ending conditions, are the same as those described in Section 4.1. As shown in Appendix C.2, the initial values ${\widehat{\mathbf{\mu }}}_N,{\widehat{\mathit{C}}}_N$ of the estimate values of the mean vector and covariance matrix are given respectively by mean vector ${\mathbf{\mu }}_x$ and covariance matrix ${\mathit{C}}_x$ of the whole samples. To ensure that the estimator is unbiased, we select the scale factor as $r=1.0$ based on the discussion of the proposed method in Section 4.3, and the kernel as ${\mathit{C}}_W={\sigma }^2\mathit{I},{\sigma }^2=1.5$ based on the discussion for a general mean-shift method in Section 4.2.

While maintaining the ratio of the number $N_N$ of samples of major cluster and the number $N_O$ of samples of the outliers to $3000:200$ and changing the number $N=N_N+N_O$ of samples from 1000 to 90,000, the variance of each estimate value of the mean vector ${\mathbf{\mu }}_N$, covariance matrix ${\mathit{C}}_N$, and number $N_N$ of samples is obtained using both methods. The estimation variance is replaced by the sample variance of each estimate for 1000 trials as the sample number changes. The estimation variances $Var\left[{\widehat{\mathbf{\mu }}}_N\right],Var\left[{\widehat{\mathit{C}}}_N\right],Var\left[{\widehat{N}}_N\right]$ are shown in Figure 7a–c. The horizontal axis shows the logarithm of various numbers of samples ${\widehat{N}}_N$. The vertical axes respectively show logarithms for estimation variances $Var\left[{\widehat{\mathbf{\mu }}}_N\right],Var\left[{\widehat{\mathit{C}}}_N\right],Var\left[{\widehat{N}}_N\right]$. In this figure, symbol ◯ corresponds to the proposed method. Symbol △ represents the general mean-shift method. Because the covariance matrix and number of samples can not be estimated in the general mean-shift method, only the results obtained using the proposed method are presented in Figure 7b,c.

From Figure 7a–c, it is readily apparent that the variance $Var[·]\to 0$ for sample number $N\to \infty$. Therefore, estimators ${\widehat{\mathbf{\mu }}}_N,{\widehat{\mathit{C}}}_N,{\widehat{N}}_N$ have consistency. Figure 7a–c are drawn as a logarithmic graph; the slope should be $-1$ in fact. Therefore, the relation between the sample number of samples and the estimation variance is approximated using a linear polynomial with the slope fixed at $-1$. The straight line represented by the approximate linear polynomial is superimposed by a solid line in these figures. These results demonstrate the validity of the approximation. Here, we simply define the estimation variance as the 0-order coefficient of the approximate linear polynomial or the virtual estimation variance corresponding to sample number $N=1$. Regarding to the estimation variance, we compare the estimation precision of proposed method with the general mean-shift method.

4.5. Estimation Precisions of the Proposed and General Mean-Shift Methods

The distributions that major cluster and outliers follow, the number $N_N,N_O$ of samples, and the update ending condition are the same as those described in Section 4.1. As shown in Appendix C.2, the initial values of the estimated value of mean vector and covariance matrix ${\widehat{\mathbf{\mu }}}_N,{\widehat{\mathit{C}}}_N$ are given, respectively, by the mean vector ${\mathbf{\mu }}_x$ and covariance matrix ${\mathit{C}}_x$ of the whole samples.

We select scale factor r to be various values and use the proposed method to estimate the mean vector ${\mathbf{\mu }}_N$, covariance matrix ${\mathit{C}}_N$, and number $N_N$ of samples. Figure 8a–c respectively present the estimation variances corresponding to the estimated values of the mean vector ${\mathbf{\mu }}_N$, covariance matrix ${\mathit{C}}_N$, and number $N_N$ of samples. The horizontal axis shows the logarithm of various scale factor r. The vertical axes respectively show the estimation variances $Var\left[{\widehat{\mathbf{\mu }}}_N\right],Var\left[{\widehat{\mathit{C}}}_N\right],Var\left[{\widehat{N}}_N\right]$. Similarly, letting the covariance matrix of kernel be ${\mathit{C}}_W={\sigma }^2\mathit{I}$, we estimate the mean vector ${\mathbf{\mu }}_N$ using the general mean-shift method while the kernel bandwidth ${\sigma }^2$ is changed to various values. The estimation variance of estimated value ${\widehat{\mathbf{\mu }}}_N$ is presented in Figure 9.

From Figure 8a–c, the estimation variance of each estimated value of the mean vector ${\mathbf{\mu }}_N$, covariance matrix ${\mathit{C}}_N$, and number $N_N$ of samples decreases with respect to r, monotonically. If r is small, then the kernel bandwidth decreases. The number of substantial points involved in the estimation decreases. Therefore, the estimation precision deteriorates. On one hand, if r is large, then the estimation precision decreases. Because bias error occurs as shown in Figure 6, it is not desirable as an estimator. However, the estimation variance related to general mean-shift method decreases monotonically with respect to kernel bandwidth ${\sigma }^2$, as shown in Figure 9. The reason is exactly the same as in the case of the proposed method.

Finally, the estimation precision of a general mean-shift method and that of the proposed method are compared. Regarding the general mean-shift method, the estimation is unbiased if ${\sigma }^2\le 1.5$, as shown in Figure 4. However, the estimation precision increases as ${\sigma }^2$ becomes larger, as shown in Figure 9. In the general mean-shift method, the optimal selected value of the kernel bandwidth is ${\sigma }^2=1.5$. The estimation variance at kernel bandwidth ${\sigma }^2=1.5$ is read from Figure 9: its value is shown by a horizontal dotted line in Figure 8a. In the proposed method, $0.5\le r\le 1.5$ is the suitable range of the scale factor r. In this range, the estimation variance of the proposed method is half or less than half of that of the general mean-shift method. The proposed method has higher estimation precision than the general mean-shift method that has the optimal kernel bandwidth for the following reason. In the general mean-shift method, the kernel shape is expressed as an isotopic shape because the covariance matrix of the kernel is represented as a diagonal matrix in which all diagonal elements are equal. Otherwise, in the proposed method, the kernel shape can take an arbitrary anisotropic shape because the covariance matrix of the kernel can take an arbitrary matrix that satisfies the condition as a covariance matrix. The practical number of samples that contribute to the estimation can be maximized by adjusting the kernel shape to the distribution of samples.

4.6. Discussion

Numerical experiments in two-dimensions described in Section 4.2, Section 4.3 and Section 4.4 yield results for the major cluster and outlier model shown in Figure 2. The purpose of our numerical experiment is to confirm whether the estimators (mean, covariance, number of sample of the major cluster) in the proposed method are unbiased and consistent without a proper pre-set of kernel bandwidth. If these estimators are consistent unbiased estimators, then the proposed method can achieve accurate estimation of the mean, covariance, and the number of samples of the major cluster. We chose the two-dimensional numerical experiment to observe the dynamic changes of the kernel more intuitively during the iteration. The iteration process is shown in Figure 3. In the numerical experiments described herein, the major cluster follows the Gaussian distribution. If the proposed method performs well on other distributions, then the scope of application of the proposed method can be expected to expand. We discuss the scope of application of the proposed method in two aspects as presented below.

For a one-dimensional signal processing field, the assumption of normality is not regarded as being such a severe strong assumption. Yokota and Ye [19] proposed the radical root, or r-th root, transform of the power spectrum series that follows the chi-square distribution, such that the transformed series follows a quasi-Gaussian distribution. Lotter and Vary [20] proposed a spectral amplitude estimator with a parametric super-Gaussian speech model for approximating the probability density distribution of the real speech spectral amplitudes. In fact, the parametric super-Gaussian distribution can approximate the Rayleigh–Laplace–Gamma distribution or other distributions exactly. Ye and Yokota [21] applied the radical root transformation to the super-Gaussian distributions. Thereby, they confirmed that the super-Gaussian distribution after r-th radical root transformation can be quasi-Gaussian distributed. By radical root transformation [21], the proposed method is applicable for major clusters that follow different distributions other than a Gaussian distribution. However, for clustering in image processing or other multiple dimensional applications, the major cluster following a Gaussian distribution is truly a strong assumption.

In addition to the problem addressed in this paper, many methods exist to solve this problem other than the mean-shift method. They have been discussed as described below. Under the normality assumption, Grubbs’ test [22,23,24] and Thompson Tau test [25] are known as methods for testing whether the sample farthest from the sample mean is an outlier. These tests are applied sequentially from the samples that are outermost from the sample mean, but the number of outliers is only valid at most to several. Moreover, applying the tests to multi-dimensional data are not easy. If the outliers follow a Gaussian distribution and if the number of clusters in which the outliers are distributed is known, then, by applying a Gaussian mixture model [26,27,28], the mean and covariance matrix of major cluster can be estimated easily using the Expectation-maximization(EM) algorithm [29,30]. However, such an assumption cannot be applied generally to the outliers.

In fact, selection of the kernel bandwidth is an important issue that strongly affects the result of the general mean-shift algorithm compared to setting of the kernel type. Therefore, we only used the Gaussian kernel to make a presumption here. However, there are many commonly used kernel functions in addition to the Gaussian kernel, such as the Epanechnikov Kernel, the Uniform Kernel, the Quartic Kernel, and the Triweight Kernel. Application of it to other kernel functions according to the derivation of this article will undoubtedly make this research more comprehensive and general. Such application is expected to be an important part of our future research.

5. Application

Considering a stochastic process $x\left(t\right)$, the short-time Fourier spectrum centering on time t with a suitable window length is denoted as $X(t,f)$. Here, f represents the frequency. Let $X_f\left(t\right)\equiv X(t,f)$ be denoted as the spectrum series if frequency f is fixed. By applying the non-steady-state analysis of the stochastic process, the spectrogram $P(t,f)={\left|X(t,f)\right|}^2$ denotes the power of the short-time Fourier spectrum $X(t,f)$. Because the frequency f is fixed, $P_f\left(t\right)$ will be designated as the power spectrum series.

Yokota and Ye [19] proposed a power spectrum estimation method robust for sudden noise. The method uses the radical root transformation to quasi-Gaussian distribution. The following concludes the process of the noise estimation algorithm proposed by Yokota and Ye [19]:

(1) Obtainpower spectrogram $P(t,f)$ from the noisy signal. We chose a pulse code modulation(PCM) recording of a noisy signal that contains a certain amount of sudden noise for analysis and computes the spectrogram with a Hamming window length of 10 ms achieving a 50% overlap between adjacent frames by short-time Fourier transformation. Figure 10a presents an example of a noisy signal for analysis and the corresponding spectrogram.

(2) Perform the following process for each frequency f:

(2-1) Use the radical root transformation in the power spectrum series $P_f\left(t\right)$ with the transformation parameter $r^{*}=3.314$. Thereby, obtain the new quasi-Gaussian distributed power spectrum series $P_f^{1/r^{*}}\left(t\right)$. Figure 10b portrays a histogram of the power spectrum series at $f=512$ Hz before the transformation.

(2-2) Compute the mode value of transformed power spectrum series $P_f^{1/r^{*}}\left(t\right)$ by kernel density estimation [5]. Then, put the mode value as the corresponding time average value $P_{noise}\left(f\right)$ of the noise power spectrum series. Figure 10c depicts a histogram of the transformed power spectrum series at $f=512$ Hz, the kernel density estimation [5] with proper kernel bandwidth and the major cluster estimation using our proposed method.

(2-3) Compute the time average value $P\left(f\right)$ of the noise power spectrum series from the time average value $P_{noise}\left(f\right)$ as

$$\begin{array}{c}\hfill P\left(f\right)={\left({\displaystyle \frac{P_{noise}\left(f\right)}{\Gamma \left({\displaystyle \frac{r^{*}+1}{r^{*}}}\right)}}\right)}^{r^{*}}.\end{array}$$

(33)

(3) Obtain $P\left(f\right)$ as an estimation of the noise power spectral density.

In the noise estimation algorithm [19], the mode estimation accuracy directly affects the noise estimation result. As Figure 10c shows, kernel density estimation [5] can be replaced by our proposed method for comparison. The proper pre-setting of kernel bandwidth is also important in kernel density estimation [5]. It exhibits a strong influence on the resulting estimate similarly to the general mean-shift method. To illustrate its effects, we obtained the noise power spectrum series from the PCM recording, which is shown in Figure 10a, for analysis. Figure 11 portrays the relation between the kernel bandwidth and kernel density estimation. The histogram shows the true density. The broken curve is under-smoothed because it includes too many spurious data artifacts arising from use of 0.000001 bandwidth, which is too small. The dotted curve is over-smoothed because using 0.0001 bandwidth obscures much of the underlying structure. The solid curve with 0.00003 bandwidth is regarded as optimally smoothed because its density estimate is close to the true density.

To assess the performance of the proposed method, general mean-shift method, and kernel density estimation for a noise estimation algorithm [19], this study uses PCM recordings of air-conditioning noise with some sudden noise, as shown in Figure 10a and without sudden noise, respectively, as test data and the true value. Noisy signal data in PCM recordings are not compressed. They have no power consumption. Figure 12 presents comparison results for noise estimation using the proposed method and kernel estimation. Here, we preset the kernel bandwidth as 0.0001. As Figure 12 shows, in the case in which an inappropriate kernel bandwidth is set in advance, noise estimation using our proposed method closely approximates the true noise, but the estimation accuracy using the kernel estimation is not high.

6. Conclusions

The study described in this paper has addressed the problem of proper pre-setting for the fixed search kernel in a general mean-shift method. To improve the estimation accuracy, a new mean-shift method was proposed in which the mean vector and covariance matrix of the major cluster are estimated at each iteration. Then, the kernel bandwidth and shape are adjusted corresponding to the estimates. In numerical experiments, we compared the estimation accuracy and precision of the proposed method and of the general mean-shift method. The experimentally obtained results demonstrate that the estimation accuracy and precision of the proposed mean-shift are higher than those of a general mean-shift method. Moreover, the proposed mean-shift can estimate the covariance matrix and the number of samples of major clusters effectively and correctly. Neither can be estimated using the general mean-shift method. These results were confirmed through formal experimentation, the results of which indicated the superior performance of our method compared to that of the general mean-shift method.