# Estimating parameters of a bounded distribution with MLE? **URL:** \ **Category:** Methodology\ **Created:** [January 12, 2021, 11:42am UTC](https://openturns.discourse.group/t/estimating-parameters-of-a-bounded-distribution-with-mle/83 "2021-01-12T11:42:24Z")\ **Posts on this page:** 2\ **Page:** 1
**Author:** ![MichaelBaudin](https://yyz2.discourse-cdn.com/free1/user_avatar/openturns.discourse.group/michaelbaudin/32/9_2.png) [@MichaelBaudin](https://openturns.discourse.group/u/MichaelBaudin)\ **Post date:** [January 12, 2021, 11:42am UTC](https://openturns.discourse.group/t/estimating-parameters-of-a-bounded-distribution-with-mle/83/1 "2021-01-12T11:42:24Z")
Hi ! I want to estimate the parameters of a `Beta` distribution from a sample based on maximum likelihood estimation (MLE). This distribution has four parameters : \alpha, \beta, a and b. The parameters a and b are bound parameters. In my case, these parameters are known a=0 and b=1, because my variable is a physical quantity between 0 and 1 (it is defined as a ratio where the denominator is greater than the numerator). The thing is that the sample contains an extreme realization “1”. This can be computed this way: ```auto import openturns as ot ot.RandomGenerator.SetSeed(1976) # Why not? # Get sample data = [0.2, 0.3, 0.5, 1.0] # Notice the "1.0" sample = ot.Sample([[x] for x in data]) # Fit : estimate all parameters (the wrong way) factory = ot.MaximumLikelihoodFactory(ot.Beta()) factory.setKnownParameter([0.0, 1.0], [2, 3]) beta = factory.build(sample) print(beta) ``` It produces: ```auto Beta(alpha = 2, beta = 2, a = 0, b = 1) ``` In other words, the MLE failed: the optimization algorithm failed to maximize the likelihood. Notice that **no exception was raised** : it seems to me that it should be so. The reason for the issue is that the sample contains the value “1.0” (the fourth value in the sample). At this point, the log-PDF is -infinity, because the PDF is zero: ```auto logpdf = ot.Beta(2.0, 2.0, 0.0, 1.0).computeLogPDF(sample) print("logpdf=", logpdf) likelihood = logpdf.computeMean() print("likelihood=", likelihood) ``` The previous script prints: ```auto logpdf= 0 : [-0.040822] 1 : [0.231112] 2 : [0.405465] 3 : [-1.79769e+308] likelihood= [-4.49423e+307] ``` Indeed, since OT1.16, the log-PDF can go down to -1.e307, but not -INF. This is a good consequence of the highly discussed PR [https://github.com/openturns/openturns/pull/1566](https://github.com/openturns/openturns/pull/1566) In OT 1.15, the likelihood is approximately equal to -708, which was a numerically wrong answer. There are several ways to fix the issue. One solution would be to drop the “1.0” value in the sample. It makes sense, since this realization does not match the distribution. Another solution is to slightly increase the bounds of the `Beta`: ```auto delta = 0.01 logpdf = ot.Beta(2.0, 2.0, 0.0 - delta, 1.0 + delta).computeLogPDF(sample) likelihood = logpdf.computeMean() ``` The previous script prints: ```auto logpdf= 0 : [-0.0390172] 1 : [0.218678] 2 : [0.385662] 3 : [-2.86287] likelihood= [-0.574386] ``` Now the likelihood is finite. Therefore, one way to estimate the parameters of the `Beta` is: ```auto distribution = ot.Beta(1.0, 1.0, 0.0 - delta, 1.0 + delta) factory = ot.MaximumLikelihoodFactory(distribution) factory.setKnownParameter([0.0 - delta, 1.0 + delta], [2, 3]) beta = factory.build(sample) ``` However, this way is not completely satisfactory, because it depends on the \delta=0.01 parameter. What if \delta \rightarrow 0? Is there a satisfactory way to perform this? It sounds to me that my problem is wrongly established, in the sense that the sample does not match the distribution. As a side comment, if the likelihood was “-INF” in place of “-1.e308”, it would seem to me that the code would be more compliant to the IEEE754 standard (but this does not change the mathematical issue). Best regards, Michaël PS-1 This situation is very similar to estimating the parameters of a `Triangular` or `Trapezoidal` distributions from a sample using the MLE. PS-2 A similar situation would occur with MCMC simulation, as pointed by @MerlinKeller and @josephmure. ---
**Author:** ![regislebrun](https://yyz2.discourse-cdn.com/free1/user_avatar/openturns.discourse.group/regislebrun/32/154_2.png) [@regislebrun](https://openturns.discourse.group/u/regislebrun)\ **Post date:** [January 12, 2021, 10:48pm UTC](https://openturns.discourse.group/t/estimating-parameters-of-a-bounded-distribution-with-mle/83/2 "2021-01-12T22:48:45Z")
Once again, IEEE754 is not the Bible and to be compliant with IEEE754 is not the Holy Graal! By the way, we _are_ compliant with IEEE754 as soon as we use a compliant processing unit and a compliant compiler. And switching from max\_scalar to inf would not have solved your actual issue, which is not to print -inf instead of -1.e308, but to get meaningful parameters estimates. Here is my analysis of your problem: - You didn’t set bounds on the unknown parameters. If you activate the logs (add `ot.Log.Show(ot.Log.INFO+ot.Log.WARN)` in your script) you will see that the optimization algorithm tries negative values for the two shape parameters (e.g. [-16.5669,-21.5641]) You can fix it by adding `factory.setOptimizationBounds(ot.Interval([ot.SpecFunc.ScalarEpsilon]*2, [1]*2, [True]*2, [False]*2)) ` before the call to the `build()` method. Note that the actual value of the upper bound is arbitrary and is not taken into account by the algorithm. - You have a point equal to the upper bound of the range. We decided in OT to keep the lower bound of the range of a distribution in the range, but not its upper bound. It was motivated by two facts: - One can exclude any negligible set from the range without modifying a continuous distribution - It solved an issue when one evaluated the PDF of a mixture of two bounded distributions such that the lower bound of one of them is equal to the upper bound of the other. For example, the mixture `ot.Mixture([ot.Uniform(-1,0), ot.Uniform(0,1)])` is equal to `ot.Uniform(-1,1)` with this choice, otherwise the PDF at 0 would evaluate to 1 instead of 0.5 As a result, one of your point is outside of the distribution range, so you cannot expect to get anything meaningful. You can fix it by adding the slightest margin to the upper bound: `factory.setKnownParameter([0.0, 1.0 + ot.SpecFunc.ScalarEpsilon], [2, 3])` With these two modifications you get `Beta(alpha = 0.291413, beta = 0.0860649, a = 0, b = 1)`. Cheers Régis