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):

enter image description here

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?


1 Answer 1


A simple chi-squared statistic is what you use if the amount of data is constant - a good model should return a value of about N - m, where there's N data and m model parameters. If the amount of data is not constant but is always above say 20, then I believe you should divide the chi-squared value by N-m. This reduced chi-squared statistic can be used to find which model parameters best match the observed values. However, without error budgets, it isn't possible to put meaningful constraints on the model parameters - one can only determine the best fit (which I guess is all right occassionally).

  • $\begingroup$ The number of elements in the observed data set is always the same (N). The number of variable parameters in the model is also always the same (p). What changes is the number of elements in each simulated data set. What is m in your answer above? The number of parameters in the model (p)? The number of elements in each simulated data set (Ni)? $\endgroup$
    – Gabriel
    Oct 7, 2014 at 2:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.