Note

Go to the end to download the full example code.

Fit a distribution by maximum likelihoodΒΆ

In this example we are going to estimate the parameters of a parametric by generic numerical optimization of the likelihood.

The likelihood of a sample (\vect{x}_1, \dots, \vect{x}_n) according to a parametric density function p_{\vect{\theta}} is:

\ell(\vect{x}_1, \dots, \vect{x}_n,\vect{\theta}) = \prod_{i=1}^n p_{\vect{\theta}}(\vect{x}_i)

import openturns as ot
import math as m

ot.Log.Show(ot.Log.NONE)

Create data from the PDF of the Normal distribution with \mu=4, \sigma=1.5.

sample = ot.Normal(4.0, 1.5).getSample(200)

Create the search interval of (\mu, \sigma) : the constraint is \sigma>0.

lowerBound = [-1.0, 1.0e-4]
upperBound = [-1.0, -1.0]
finiteLowerBound = [False, True]
finiteUpperBound = [False, False]
bounds = ot.Interval(lowerBound, upperBound, finiteLowerBound, finiteUpperBound)

Create the starting point of the research:

  • for \mu : the first point,

  • for \sigma : a value evaluated from the two first points.

mu0 = sample[0][0]
sigma0 = m.sqrt((sample[1][0] - sample[0][0]) * (sample[1][0] - sample[0][0]))
startingPoint = [mu0, sigma0]

Create the estimator from a parametric PDF.

pdf = ot.Normal()
factory = ot.MaximumLikelihoodFactory(pdf)
factory.setOptimizationBounds(bounds)

Set the starting point via the solver.

solver = factory.getOptimizationAlgorithm()
solver.setStartingPoint(startingPoint)
factory.setOptimizationAlgorithm(solver)

Estimate the parametric model.

distribution = factory.build(sample)
str(distribution)

'Normal(mu = 3.93644, sigma = 1.43033)'

Retrieve the estimated parameters.

print(distribution.getParameter())

[3.93644,1.43033]