Estimate 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.

In [1]:
from __future__ import print_function
import openturns as ot
import math as m
In [2]:
# Create data from a gaussian pdf with mu=4, sigma=1.5
sample = ot.Normal(4.0, 1.5).getSample(200)
In [3]:
# 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)
In [4]:
# Create the starting point of the research
# For mu : the first point
# For sigma : a value evaluated from the two first data
mu0 = sample[0][0]
sigma0 = m.sqrt((sample[1][0] - sample[0][0]) * (sample[1][0] - sample[0][0]))
startingPoint = [mu0, sigma0]
print(startingPoint)
[4.912302476828147, 2.811562130153132]
In [5]:
# Create the estimator from a parametric pdf
pdf = ot.Normal()
factory = ot.MaximumLikelihoodFactory(pdf)
factory.setOptimizationBounds(bounds)
In [6]:
# Set the starting point via the solver
solver = factory.getOptimizationAlgorithm()
solver.setStartingPoint(startingPoint)
factory.setOptimizationAlgorithm(solver)
In [7]:
# Estimate the parametric model
distribution = factory.build(sample)
print(distribution)
Normal(mu = 3.94055, sigma = 1.48893)