Add Compile function in a code to speed up the execution

Question

I have a code

 ClearAll["Global`*"];
h[a1_, l_] = 
  ProbabilityDistribution[
   a1 l Exp[-l  y0] (1 - Exp[-l  y0])^(a1 - 1), {y0, 0, \[Infinity]}];
ParallelTable[
  data = RandomVariate[h[3, 1.5], 100];
  lnL1[a1_?NumericQ, l_?NumericQ] := 
   Block[{}, LogLikelihood[h[a1, l], data]];
  e2 = NMaximize[{lnL1[a1, l], {a1 > 0}, {l > 0}}, {a1, l}, 
    Method -> "SimulatedAnnealing"], {5}] // AbsoluteTiming

This code works fine, but takes too much time for huge replications. To save time I know Compile function execute the code in C. But I am unable to add `Compile` function in that code.

Answer 1

There is no need for manual compilation. It's possible to write the code in such a way that NMaximize automatically compiles the likelihood function for you. It will also be able to compute the gradient symbolically, which allows for faster convergence and higher precision.

h[a1_, l_] = 
  ProbabilityDistribution[
   a1 l Exp[-l y0] (1 - Exp[-l y0])^(a1 - 1), {y0, 0, \[Infinity]}];
data = RandomVariate[h[3, 1.5], 100];
(*Your code*)

logPDF[a1_, l_, x_] = PowerExpand@Log@Refine[PDF[h[a1, l], x], x > 0]
(*PowerExpand often results in a simpler log-likelihood function*)

logLikelihood = Total@logPDF[a1, l, data]
(*Note how three of its four terms each appear only once*)

NMaximize[{logLikelihood, {a1 > 0, l > 0}}, {a1, l}] // AbsoluteTiming
(*{0.131648,{-94.9113,{a1->3.18393,l->1.59988}}}*)

transformedParameters = logLikelihood /. {a1 -> Exp[loga1], l -> Exp[logl]};
FindMaximum[transformedParameters, {loga1, logl}] // AbsoluteTiming
(*{0.019304,{-94.9113,{loga1->1.15812,logl->0.469928}}}*)
(*Unconstrained optimization is much faster*)

Edit: I forgot Mathematica has a built-in function specifically for this task. It's not as fast, though.

FindDistributionParameters[data, h[a1, l], {{a1, 1}, {l, 1}}] // AbsoluteTiming
(*{0.066339, {a1->3.18393,l->1.59988}}*)

Edit 2: C Compiled Functions

There are basically two approaches: accumulate log-PDFs sample-by-sample, or incrementally compute the likelihood function with vectorized operations on the dataset as a whole. The first approach is cache-friendly, the second utilizes fast SIMD instructions. It's hard to tell which one will work better on the real problem, so I'll implement both here.

logPDFfunction = Function[{a1, l, x}, Evaluate@PowerExpand@Log@Refine[PDF[h[a1, l], x], x > 0]]
(*Workaround for a weird bug preventing unpure functions from being inlined on some systems. Investigation in progress. Thanks to @Michael E2 for pointing this out.*)

logLikelihoodC1 = Compile[{{a1, _Real}, {l, _Real}, {x, _Real, 1}},
   Module[{i, s = 0.},
    Do[
     s += logPDFfunction[a1, l, x[[i]]],
     {i, Length[x]}];
    s],
   CompilationTarget -> "C",
   CompilationOptions -> {"InlineExternalDefinitions" -> True}];

logLikelihoodC2 = Compile[{{a1, _Real}, {l, _Real}, {x, _Real, 1}},
   Total@logPDFfunction[a1, l, x],
   CompilationTarget -> "C",
   CompilationOptions -> {"InlineExternalDefinitions" -> True}];

objective[function_, a1_?NumericQ, l_?NumericQ, data_] := function[a1, l, data];

NMaximize[{objective[logLikelihoodC1, a1, l, data], a1 > 0, l > 0}, {a1, l}, Method -> "SimulatedAnnealing"] // AbsoluteTiming
(*{0.061932, {-73.9552, {a1 -> 3.77946, l -> 1.99648}}}*)

NMaximize[{objective[logLikelihoodC2, a1, l, data], a1 > 0, l > 0}, {a1, l}, Method -> "SimulatedAnnealing"] // AbsoluteTiming
(*{0.055937, {-73.9552, {a1 -> 3.77946, l -> 1.99648}}}*)