Finding parameters of a normal distribution which maximize the difference between two likelihood functions

Question

Consider this simple example where I have 2d data that looks like this:

I'm trying to find the parameters of a Normal distribution that maximizes the difference between two likelihood functions $$L_1-L_2=\prod^{n_i}_{i=1}f(x_i|\mu,\sigma^2)-\prod^{n_k}_{k=1}f(x_k|\mu,\sigma^2).$$ Here, the $X_i$ may be the black points and the $X_k$ may be the red points. I'm looking to find $\theta=\{\mu, \sigma^2\}$ that makes the joint density of the black points as large as possible while keeping the joint density of the red points as small as possible.

I already tried another approach by instead maximizing $$\frac{L_1}{L_2}=\frac{\prod^{n_i}_{i=1}f(x_i|\mu,\sigma^2)}{\prod^{n_k}_{k=1}f(x_k|\mu,\sigma^2)}$$ which led me to $$\hat\mu=\frac{\sum^{n_i}x_i-\sum^{n_k}x_k}{n_i-n_k}$$ and $$\hat\sigma^2=\frac{\sum^{n_i}(x_i-\mu)^2-\sum^{n_k}(x_k-\mu)^2}{n_i-n_k}.$$ However, $\hat\sigma^2$ can quickly become negative, which confused me. (I worked with the univariate case to not further complicate things).

So now I'm trying to maximize the difference between the two likelihood functions stated in the beginning. My question is: Is there any mathematical "trick" that simplifies this problem? For the derivation of the standard ML estimators taking logs simplifies the problem a lot but in this case, it does not make much sense given that I no longer have a pure product anymore.

Any hint towards related questions/works/papers is also appreciated!

Answer 1

Here is a simplified example (one-dimensional) of how to do this numerically.

# Create fake data
    x1 <- rnorm(10,2,2)
    x2 <- c(0, 1, 1.5)
    
# Returns the difference in the log likelihood functions;
# equivalent to the ratio of the likelihoods for min/max purposes

    foo <- function(theta) {
      mu <- theta[1]
      sigma <- exp(theta[2])
      
      rll <- sum(dnorm(x1, mu, sigma, log=TRUE)) - 
               sum(dnorm(x2, mu, sigma, log=TRUE))
      rll
    }

# Maximize the ratio; the starting values are the mean and log(std. dev.) of x1     

    optim(c(mean(x1), log(sd(x1))), foo, control=list(fnscale=-1))

When we run it, we get:

> optim(c(mean(x1), log(sd(x1))), foo, control=list(fnscale=-1))
$par
[1] 2.1880969 0.6802588

which converts to $\hat{\mu} = 2.188$, $\hat{\sigma} = 1.974$.

Comparing to the mean and std. deviation of x1 shows that the parameter estimates are shifted away from where x2 is located:

> mean(x1)
[1] 1.781837
> sd(x1)
[1] 1.8948

In a way, for this data, this isn't surprising; x2 is shifted somewhat from the location of x1, so maximizing x1 will result in a relatively low likelihood for the x2 data to start with. Since there's a lot more data in x1 than x2, it's more important to keep the likelihood associated with each observation in x1 high than to reduce the likelihood associated with each observation in x2. Combining these two ideas allows us to conclude that, at least for this example, we won't shift too far away from the MLE for x1.