Goodness of fit: observed vs simulated data

I have a set of 2-dimensional "observed" data of sample size N:

$$O = \{(x_1, y_1), (x_2, y_2), ..., (x_N, y_N)\}$$

The hypothesis is that $O$ is a realization of a model $M(v_1, v_2, ..., v_p) = M(\theta)$ with some unknown combination of values for its $p$ parameters.

I can use this model to generate a large number of "simulated" sets, $S_i$:

$$S_i=M(\theta_i) = \{(x'_1, y'_1), (x'_2, y'_2), ..., (x'_{Ni}, y'_{Ni})\}$$

by varying the values of its $p$ parameters.

Each simulated set attempts to replicate the observed set $O$ (Notice that $Ni$, the number of points in each simulated set, is not necessarily equal to $N$, the number of points in the observed set).

I need to find the set of parameter values $\theta^\star$, which generates the simulated set $S_\star=M(\theta^\star)$, such that it produces the best fit with the observed data $O$. This combination of parameters, $\theta^\star$, is thus the one I associate with $O$.

In order to find the best fit simulation-observation I thought of using a maximum likelihood estimator of the form:

$$L = \prod_{i=1}^N \frac{1}{Ni} \sum_{j=1}^{Ni} e^{-\frac{(x_i - x'_j)^2}{2}} e^{-\frac{(y_i - y'_j)^2}{2}}$$

or its logarithmic form (which I minimize):

$$L = -\sum_{i=1}^N \log ( \frac{1}{Ni} \sum_{j=1}^{Ni} e^{-\frac{(x_i - x'_j)^2}{2}} e^{-\frac{(y_i - y'_j)^2}{2}} )$$

(disregarding for the moment Gaussian errors in the observed data set)

I assumed the $1/Ni$ term would take care of those models with a very large number of elements, but this is not happenning. As can be seen in the simple example below, as the simulated sets get larger, the log-likelihood drops:

(Left column: Observed data set. Right column: Log-ikelihood of the observed data set versus various models. In red, top left, is the observed vs observed log-likelihood value, the rest (in blue) are simulated vs observed data set with increasing number of elements in the simulations):

For what I can gather, this is a knonw behaviour of the maximum likelihood estimator.

So my question is: how can I find the best fit simulated set, given my observed data set $O$ and the model $M(\theta)$, and accounting for the fact that the simulations have different number of elements?