# Add Compile function in a code to speed up the execution

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.

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];

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

• I don't understand your last comment. Even if we compare only the execution time of your NMaximize to that of FindDistributionParameters, the built-in is twice as fast as your solution, isn't it? (+1) May 20, 2020 at 0:36
• @MarcoB, FindDistributionParameters uses unconstrained FindMaximum by default, so it's fair to compare it to FindMaximum with objective function constructed by hand as above. The magic of PowerExpand made a 3-fold difference! To be absolutely fair, I should have compared reparametrized FindDistributionParameters to reparametrized FindMaximum, but that would only increase the performance gap. May 20, 2020 at 1:56
• @SAAN, oh. I updated the answer. May 20, 2020 at 11:02
• I think you want CompilationOptions -> {"InlineExternalDefinitions" -> True} to avoid call-backs to the kernel to evaluate logPDFfunction. May 20, 2020 at 13:46
• FWIW, I'm using 12.1 on a Mac, 2.8GHz 4-core i7. I've been assuming you know about CompiledFunctionToolsCompilePrint@logLikelihoodC1 and MainEvaluate (but I mention it just in case). Consider omitting CompilationTarget -> "C". May 20, 2020 at 14:45