I am trying to use maximum likelihood for estimating a parameter for a model based on some experimental data. The likelihood function I am trying to maximize calls a number of custom functions I have written. In short, the likelihood function calls a number of custom functions that modify a list of genotypic counts in a generation to predict the counts in the next generation. The likelihood function compares the predicted and actual genotype numbers to calculate a likelihood. If I provide a value for the parameter (x), the likelihood function can be evaluated almost instantly. The likelihood function can also be plotted pretty quickly for a range of values for the parameter. But when I try to run NMaximize on the likelihood function with that parameter, Mathematica seems to take forever without producing anything. The memory usage seems to be absurd (when I ran it on a High Performance Computing cluster node, max memory usage was >90GB).

The structure of the code is-

Function1 := Function[{x,y}, (* a number of operations done on a list of genotypes *)];

Function2 := Function[{x,y}, (* more operations done on the list produced by Function1 *)];




LogL := Function[{x,y}, (* call Function1[x,y], Function2[x,y], ... produce LLvalue; *)];

LoglikelihoodFunc := (LogL[x,y];LLvalue)

I'm trying to maximize with pretty low precision & accuracy

NMaximize[{LoglikelihoodFunc, 1 >= x && x > 0}, {x}, AccuracyGoal -> 3, PrecisionGoal -> 3]

Giving initial range for x does not help.

But if I run

Plot[LoglikelihoodFunc, {x,0,1}]
(* or *)
Plot3D[LoglikelihoodFunc, {x,0,1},  {y,0,1}]

it runs just fine and produces a plot within 10 seconds.

Grateful for any help/suggestions.

  • 1
    $\begingroup$ not related to your problem per se, but just so you know, you can use simply = instead of := when defining f = Function[{x,y}, ...], since Function does not evaluate its body until it's called with arguments (try f = Function[{x}, 1+1]). You can also use the syntax f[x_, y_] := <function body> instead of f = Function[{x,y}, <function body>] for the same effect; not sure which, if either, is faster. $\endgroup$ – thorimur Mar 11 at 22:39
  • $\begingroup$ Have you tried lowering WorkingPrecision in NMaximize via WorkingPrecision -> 3? $\endgroup$ – thorimur Mar 11 at 22:46
  • $\begingroup$ Even if it's long, I think it might help to include the actual code for these functions, otherwise it could be very difficult to guess what's taking the time and memory. $\endgroup$ – thorimur Mar 11 at 22:48
  • $\begingroup$ @thorimur I tried lowering working precision, but it didn't help. I did figure out a solution to the problem though. And it has to do with how the function is defined. If I had seen your first comment earlier, I would have probably figured it out sooner! $\endgroup$ – Ssdmitten Mar 17 at 16:48

I figured out what was going wrong. I think it has to do with how NMaximize (or NMinimize) handles functions that are defined in different ways.

Defining the function to be maximized as

LoglikelihoodFunc := (LogL[x,y];LLvalue)

works fine for evaluating the function, e.g. in Plot. But, when it is called in NMaximize as

NMaximize[{LoglikelihoodFunc, 1 >= x && x > 0}, {x}]

The value of 'x' is not passed on to the other functions inside the LoglikelihoodFunc function. So in the example I have in my original question, LogL was being evaluated without a value for x (I added Print[x]; inside LogL to check this). And I think this was causing NMaximize to attempt optimization on a VERY large non-numeric function. (Not clear to me why that was ramping up the memory usage)

But, defining the function to be maximized in a different way worked.

LoglikelihoodFunc[x_?NumericQ, y_?NumericQ] := (LogL[ x, y]; LogLValue);

NMaximize[{LoglikelihoodFunc[x, y], 0 <= x <= 1, 0 <= y <= 1}, {x, y}]

When defined this way, the values of the parameters were passed on to nested functions.

Surprisingly, adding the ?NumericQ was necessary. Not sure why, but NMaximize fails to pass on the values of the parameters x and y to nested functions if I removed the ?NumericQ condition, even with the 0<=x<=1,0<=y<=1 constraints.

I feel like this restriction about function definition should be part of the documentation of NMaximize (or shouldn't be a restriction).


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.