# Maximum likelihood estimation by simulating the probability density function

Let's say i want to fit a certain exotic distribution using maximum likelihood estimation. However, i have no access to the theoretical probability density function to retrieve the likelihood (it is an exotic distribution). On the other hand, i can use a model to simulate the process that is supposed to give rise to this exotic distribution, and therefore, i can have the simulation return a "simulated" probability density function (pdf) given some parameters.

The question is: can i use this simulated pdf to retrieve the likelihood of my data? I don't see any immediate problem with doing this, if the simulation is based enough samples so that the simulated pdf gets close enough to the theoretical pdf (to which i have no access).

I am mainly wondering whether this approach is correct? Or am i missing something?

The simulated method of moments is a way in econometrics to exploit this ability to produce new samples for a given value of the parameter $$\theta$$.

If $$\mathbf x^\text{obs}$$ is the observed sample and if $$\mathbf x(\theta)$$ is a simulated sample associated with the value $$\theta$$ of the parameter, then a consistent estimate of $$\theta$$ is produced by the minimisation programme $$\arg\min_\theta d(\hat\theta(\mathbf x^\text{obs}),\hat\theta(\mathbf x(\theta)))$$where

• $$d(\cdot,\cdot)$$ is a measure of distance
• $$\hat\theta(\cdot)$$ is a statistic like (i) an estimator of $$\theta$$ (b) a vector of empirical moments, enough to identify $$\theta$$ (iii) an estimator of a parameter of an auxiliary model, complex enough to identify $$\theta$$ [in which case this is called indirect inference] or (iv) a score function.

For instance, if $$\theta$$ is the mean of the elements in the sample, one could estimate $$\theta$$ as $$\arg\min_\theta ||\overline{\mathbf x^\text{obs}}-\overline{\mathbf x(\theta)}||^2$$by making the empirical means as close as possible.

Related methods of inference for intractable likelihoods of generative models are approximate Bayesian computation and synthetic likelihood.

From the outside, it's a little difficult to understand what the difference between the simulated likelihood and the "real" likelihood is in your example.

In order to answer the question, I will assume that you can evaluate the likelihood function but you have some noise on top. In other words, you want to maximize $$f(x)$$, but at best you are able to evaluate $$f(x) + \epsilon$$, where $$\epsilon$$ is a random error term. This also implies you probably can't evaluate something like $$f'(x)$$ either so standard gradient-based methods are out.

If that describes your problem relatively well, you're in luck! Bayesian optimization is a method for addressing exactly that. In a nutshell, we assume $$f$$ is a Gaussian-Process. This allows us to evaluate a few points of $$f(x)$$ or $$f(x) + \epsilon$$ and from this, evaluate the uncertainty of $$f(x)$$ over a grid of unevaluated points. We then pick the points that have promise (i.e., we are not certain they are not a global max), try those points, update uncertainty, repeat.

One thing to note about this method is that the search for "points of promise" over our Gaussian process is particularly slow. This means Bayesian Optimization can be very slow compared to classical methods such as Newton's Method if the evaluation function (and derivatives for Newton's Method) is not too expensive. However, if your target function is expensive to compute anyways (such as cross-validation error) then this cost is irrelevant.