1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
|
"""
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 :math:`(\vect{x}_1, \dots, \vect{x}_n)` according to a parametric density function :math:`p_{\vect{\theta}}` is:
#
# .. math::
# \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 :math:`\mu=4`, :math:`\sigma=1.5`.
sample = ot.Normal(4.0, 1.5).getSample(200)
# %%
# Create the search interval of (:math:`\mu`, :math:`\sigma`) : the constraint is :math:`\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 :math:`\mu` : the first point,
# - for :math:`\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)
# %%
# Retrieve the estimated parameters.
print(distribution.getParameter())
|
