Parameter estimation without an explicit likelihood function

Question

I have a parametric model, some data $y$, and I would like to find a maximum likelihood estimate for the model parameters $\theta$. My usual approach would be to write down the likelihood function $\mathcal{L}(\theta\,|\,y)$ and use numerical optimization to find the optimal parameters $\theta_\mathrm{MLE}$.

However, I am having difficulty in analytically constructing $\mathcal{L}$. The PDF of $Y$ involves the sum of two random variables with different distributions, and I'm unable to evaluate the convolution integral analytically.

On the other hand, given $\theta$ I can very easily generate independent samples of $Y$, and I have a lot of computing power at my disposal. Is there a way to use these two things to obtain $\theta_\mathrm{MLE}$ without explicitly constructing the likelihood function?

For example, I could imagine fiddling with $\theta$ until the distribution of the simulated sample somehow "looks like" the observed data. With an appropriate loss function, I can run stochastic gradient descent to arrive at a $\theta^\star$ that is optimal in some sense. Is there a more rigorous formulation of this approach? Is this actually equivalent to finding the MLE?

---

In case it's relevant, here is the model that I'm working with. I observe a random variable $Y$ where: $$ \begin{align*} Y &= X + N \\ N &\sim \mathrm{Normal}(0,\sigma^2) \end{align*} $$ and $X$ has a probability $p$ of following a power law distribution, and probability $(1-p)$ of $X=0$, yielding a probability density $$ f_X(x) = \frac{p}{2} \frac{\alpha-1}{(|x|+1)^\alpha} + (1-p)\,\delta(x). $$ So the dimensionality is very low (univariate random variable, 3 parameters) and my data is of a reasonable size (a few hundred thousand samples).

Answer 1

When one can't get a nice closed-form log likelihood, then numerically estimating the log likelihood can sometimes work.

First I'd change the structure of the mixture such that $Y=N$ with probability $p$ and $Y=N+X$ with probability $1-p$. That way we can can numerically determine the likelihood for the sum of $N$ and $X$.

Here is some Mathematica code to do so:

(* Define probability distributions *)
normal = NormalDistribution[0, sigma];
power = ProbabilityDistribution[(alpha - 1)/(2 (Abs[x] + 1)^alpha), 
  {x, -Infinity, Infinity}];
sum = TransformedDistribution[x + n, {n \[Distributed] normal, x \[Distributed] power}];
mixture = MixtureDistribution[{p, 1 - p}, {normal, sum}];

(* Function to calculate log of the likelihood *)
logL[y_, sigma_?NumericQ, alpha_?NumericQ, p_?NumericQ] := 
 Sum[Log[p PDF[NormalDistribution[0, sigma], 
   y[[i]]] + (1 - p) NIntegrate[(alpha - 1)/(2 (Abs[x] + 1)^alpha)
   PDF[NormalDistribution[0, sigma], y[[i]] - x], {x, -Infinity, Infinity}]],
   {i, Length[y]}]

(* Generate a random sample *)
parms = {sigma -> 1, alpha -> 3, p -> 0.8};
SeedRandom[12345];
y = RandomVariate[mixture /. parms, 50];

(* Get maximum likelihood estimates *)
mle = FindMaximum[{logL[y, sigma, alpha, p], sigma > 0 && alpha > 1 && 0 <= p <= 1},
  {{sigma, 1}, {alpha, 3}, {p, 0.8}}]
(* {-73.574, {sigma -> 1.05395, alpha -> 10.2403, p -> 0.999983}} *)

Answer 2

You may use numerical integration to efficiently and accurately compute the likelihood function of the parameters, see for example the Gauss-Hermite quadrature: https://en.wikipedia.org/wiki/Gauss%E2%80%93Hermite_quadrature

Your likelihood is given by a marginalization integral \begin{equation} p(Y) = \int P(Y,X) dX = \int P(Y|X) P(X) dX \end{equation} For the assumed model \begin{equation} P(Y|X) = \frac{1}{\sqrt{2\pi \sigma^2} } \exp\big(-\frac{1}{2\sigma^2} (Y - X)^2\big) \end{equation} and you described $P(X)$ in the edit of your question.

The computations should take only fractions of a second.