non-parametric estimation of E(X) for a symmetric density

Question

I am trying to estimate the expectation of a symmetric density using a kernel density estimation. For a known density, the MLE of this quantity is $\sum_{i=1}^n f ( \theta - X_i)$, for $n$ realisations of $X=(X_1 ... X_n)$. I'm replacing the density by a non-parametric estimate (gaussian kernel, cross-validated bandwidth selection).

I have the following function to get a point estimate of the density in R :

density_estimate_gaussian <- function(point,data,bandwidth){

  k_values <- rep(0,length(data))
  for(i in 1:length(data)){
    temp <- (point - data[i])/ bandwidth
    k_values[i] <- dnorm(temp, mean = 0, sd = 1) / bandwidth
  }
  return(sum(k_values) / length(data))
}

Using this function, I am computing a point estimate of the log-likelihood with the following function :

log_likelihood <- function(theta,data,bandwidth){
  lik_values <- rep(0,length(data))
  for(i in 1:length(data)){
    lik_values[i] <- log(density_estimate_gaussian(theta - 
    data[i],data,bandwidth))  
  }

  return(sum(lik_values))
}

Lastly, I am using the optim function in R (Brent method) to get my MLE estimation (bandwidth is from a previous KDE with the density(data) R function) :

MLE <- optim(par = 2.2,log_likelihood,data = data$V1, bandwidth = 1.67,
             method = "Brent",lower = 1, upper = 3)

This is where it gets problematic.

My likelihood is very small in a lot of points, which makes my log-likelihood infinitely small, and I can't get a correct MLE estimate of $E(X)$. An histogram suggests it sits around 2.2.

What am I supposed to do ? Should I remove some extreme data points and optimize my log-likelihood in a smaller region to get a correct estimate ?

My biggest issue would be finding the data points I am supposed to remove. When I compute my log-likelihood in all the grid points used by the density function, the log-likelihood is actually -Inf in most of those. Should I remove the data points closest to the grid point where the log-likelihood is -Inf ?

Are there any obvious errors in my R functions ?

Answer 1

So, your nonparametric model for the density is a mixture of normal distributions:$\DeclareMathOperator{\E}{\mathbb E}$ $$ \hat P = \frac1n \sum_{i=1}^n \mathcal N(X_i, \sigma^2) .$$ You're trying to find the best bandwidth $\sigma$ for this model. This is a good thing to do for density estimation.

But, consider the mean of this approximating distribution: \begin{align} \E_{Y \sim \hat P}[Y] &= \int_{-\infty}^\infty y \, \mathrm{d}\hat P(y) \\&= \int_{-\infty}^\infty y \, \left( \frac1n \sum_{i=1}^n \mathcal N(y; X_i, \sigma^2) \right) \mathrm d y \\&= \frac1n \sum_{i=1}^n \int_{-\infty}^\infty y \, \mathcal N(y; X_i, \sigma^2) \mathrm d y \\&= \frac1n \sum_{i=1}^n X_i ,\end{align} the obvious estimator for the mean in the first place.

So, no matter the bandwidth, you're going to get the same estimator for the mean. If you really only care about the mean, then just don't worry about the bandwidth.