Parameter estimation without an explicit likelihood function 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).
 A: 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}} *)

A: 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.
