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

  • 2
    $\begingroup$ You just described "likelihood free inference". Check out Approximate Bayes Computation (ABC). Also take a look at ELFI "ELFI: Engine for Likelihood Free Inference" arxiv.org/abs/1708.00707 $\endgroup$ Commented Dec 22, 2017 at 0:01
  • $\begingroup$ @VladislavsDovgalecs Wow, that's exactly what I'm looking for. Is there a good review article or a textbook that you could recommend for more background? $\endgroup$
    – KQS
    Commented Dec 22, 2017 at 0:39
  • $\begingroup$ Unfortunately I can't. I remember reading very few resources about likelihood free inference (mainly from "giants" on twitter) but never actually worked with it. I am still playing with basic probabilistic programming models where one can write the likelihood function. I would like to endorse the PyMC3 project :) $\endgroup$ Commented Dec 22, 2017 at 1:11
  • $\begingroup$ Your problem also reminds me Bayesian optimization using Gaussian Processes. One does not know the functional form the process is described with but can relatively easily compute function values given independent variables. The packages such as BayesOpt can help to find the minimum of such function. At each iteration the algorithm suggests to evaluate new points which you faithfully evaluate and feed back to the algorithm. The process repeats until convergence. These models are actually very simple to use. In your case with 3 parameters, I am willing to bet that it will work nicely. $\endgroup$ Commented Dec 22, 2017 at 1:15
  • 2
    $\begingroup$ In the question title you mention no explicit likelihood yet you seem to be assuming one with your assumption that $N$ is a normal rnadom variable and X is a zero-inflated power law distribution. If you make this assumption the method has a likelihood. The likelihood of $Y$ will be a complicated object but this is not a likelihood free method. $\endgroup$ Commented Dec 22, 2017 at 2:35

2 Answers 2


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};
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}} *)

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.


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.