How to configure the optimization algorithm during maximum likelihood estimation in OpenTURNS?

152 Views Asked by At

I have a Sample and i want to fit the parameters of a Beta distribution with maximum likelihood estimation. Moreover, I want to truncate its parameters into the [0,100] interval. This should be easy with MaximumLikelihoodFactory, but the problem is that the optimization algorithm fails. How may I change the algorithm so that it can succeed?

Here is a simple example, where I generate a sample with size 100 and configure the parameters a and b with the setKnownParameter.

import openturns as ot

# Get sample
beta_true = ot.Beta(3.0, 1.0, 0.0, 100.0)
sample = beta_true.getSample(100)

# Fit
factory = ot.MaximumLikelihoodFactory(ot.Beta())
factory.setKnownParameter([0.0, 100.0], [2, 3])
beta = factory.build(sample)
print(beta)

The previous script produces:

Beta(alpha = 2, beta = 2, a = 0, b = 100)
WRN - Switch to finite difference to compute the gradient at point=[0.130921,-2.18413]
WRN - TNC went to an abnormal point=[nan,nan]

The algorithm surely fails, since the values of alpha and beta are unchanged with respect to their default values.

I do not know why this fails, perhaps because it uses finite difference derivatives. Anyway, I would like to customize the optimization algorithm and see if it can change anything to the result, but I do not know how to do this.

1

There are 1 best solutions below

0
On

The ResourceMap has a key which allows to configure the optimization algorithm. The value is a string which is the name of the default algorithm:

"MaximumLikelihoodFactory-DefaultOptimizationAlgorithm": "TNC"

The code uses this algorithm to perform the maximum likelihood estimation (MLE). But it does not say what values can be set. Actually, the code of the MLE uses the OptimizationAlgorithm.Build static method to create the optimization algorithm. According to the doc, this is the "Name of the algorithm or problem to solve. For example TNC, Cobyla or one of the NLopt solver names.". So I can configure, say M. J. D. Powell's "Cobyla" algorithm:

ot.ResourceMap.SetAsString("MaximumLikelihoodFactory-DefaultOptimizationAlgorithm", "Cobyla")
factory = ot.MaximumLikelihoodFactory(ot.Beta())
factory.setKnownParameter([0.0, 100.0], [2, 3])
beta = factory.build(sample)
print(beta)

The previous script produces:

Beta(alpha = 2.495, beta = 0.842196, a = 0, b = 100)

This shows that the algorithm now performs correctly. I can also use one of NLOPT's algorithms, e.g. "LN_AUGLAG" (this is a "L"ocal algorithm, with "N"o derivatives).