4
$\begingroup$

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.

$\endgroup$
6
$\begingroup$

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}}}*)
$\endgroup$
10
  • 1
    $\begingroup$ 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) $\endgroup$
    – MarcoB
    May 20 '20 at 0:36
  • $\begingroup$ @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. $\endgroup$
    – aooiiii
    May 20 '20 at 1:56
  • 1
    $\begingroup$ @SAAN, oh. I updated the answer. $\endgroup$
    – aooiiii
    May 20 '20 at 11:02
  • 1
    $\begingroup$ I think you want CompilationOptions -> {"InlineExternalDefinitions" -> True} to avoid call-backs to the kernel to evaluate logPDFfunction. $\endgroup$
    – Michael E2
    May 20 '20 at 13:46
  • 1
    $\begingroup$ FWIW, I'm using 12.1 on a Mac, 2.8GHz 4-core i7. I've been assuming you know about CompiledFunctionTools`CompilePrint@logLikelihoodC1 and MainEvaluate (but I mention it just in case). Consider omitting CompilationTarget -> "C". $\endgroup$
    – Michael E2
    May 20 '20 at 14:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.