Message Boards Message Boards

GROUPS:

Fix issues with Machine Precision in NMaximize?

Posted 1 month ago
733 Views
|
21 Replies
|
16 Total Likes
|

I am having daily problems with the new Machine Precision standard from V11.3 and now in V12.

Please suggest a quick fix that will work for all functions/procedures.

I will not be able to use Mathematica in the future if this problem is not fixed.

21 Replies

We can't help you if you don't tell use what the problem is. Just start telling us what's going on.

There are more examples, but here is my latest example with NMaximize[]:

gbdist[n, q, [Delta]_] := ParameterMixtureDistribution[BinomialDistribution[n, p], p [Distributed] BetaDistribution[ q (-1 + 1/[Delta]), ((-1 + q) (-1 + [Delta]))/[Delta]]];

loglikelihood[data, {q, [Delta]_}] := \!( *UnderoverscriptBox[([Sum]), (i = 1), (Length[data])](Log[ PDF[gbdist[data[([)(i, 1)(])], q, [Delta]], data[([)(i, 2)(])]]]))

antdata={{200, 40}, {200, 78}, {200, 97}, {200, 81}, {200, 96}}

max1 = NMaximize[{loglikelihood[antdata1, {q, [Delta]}], 0 < q < 1, 0 < [Delta] < 1}, {q, [Delta]}]

In Mathematica V 11.2: Out[26]= {-23.49, {q -> 0.370439, [Delta] -> 0.073078}}

In Mathematica V 13: General::munfl: 1.9125910^-158 (-4.167710^-216) is too small to represent as a normalized machine number; precision may be lost.

These will not copy-and-paste. I rewrote so they can be used by others.

gbdist[n, q, del_] := ParameterMixtureDistribution[BinomialDistribution[n, p], Distributed[p, BetaDistribution[q (-1 + 1/del), ((-1 + q) (-1 + del))/del]]]

loglikelihood[data, {q, del_}] := Sum[Log[PDF[gbdist[data[[i, 1]], q, del], data[[i, 2]]]], {i, Length[data]}]

antdata = {{200, 40}, {200, 78}, {200, 97}, {200, 81}, {200, 96}}

max1 = NMaximize[{loglikelihood[antdata, {q, del}], 0 < q < 1, 0 < del < 1}, {q, del}]

As noted, this gives messages. Not noted,but more important, is that it returns unevaluated (which of course is bad). As suggested by @MarvinRayBurns, the following variant gives a viable result.

max1 = NMaximize[{loglikelihood[antdata, {q, del}], 0 < q < 1, 0 < del < 1}, {q, del}, WorkingPrecision->20]

A bug report has been filed about this.

Why limit NMaximize to machine precision? You might want to try the option WorkingPrecision. Like

         NMaximize[Cos[x^2 - 3 y] + Sin[x^2 + y^2], {x, y},WorkingPrecision -> 20]

.

Thank you Marvin and Daniel.

As noted by Daniel it now works.

I hope I can use the same modification in similar problems.

Daniel, I understand from your response that in the future Mathematica will automatically increase working precision instead of not evaluating the expression. Is that correct?

Christian

Dear Christian,

What changed in 11.3 is that Mathematica has aligned the workings of machine precision arithmetic to the rest of the industry, in that machine precision underflows no longer trigger an exception and a fall-back to Mathematica's extended precision arithmetic, but instead machine precision underflows flush to zero.

For example, using SciPy:

In [12]: scipy.stats.distributions.norm.pdf(39.0) # PDF flushes to zero
Out[12]: 0.0

In [13]: scipy.stats.distributions.norm.logpdf(39.0) # solution: use logPDF
Out[13]: -761.4189385332047

With this in mind, your example can be easily made to work should you have replaced Log[PDF[dist, x]] with LogLikelihood[dist, {x}].

Let me also reparameterize your distribution somewhat, and work in terms of gbdist[n, hh/(qq + hh), 1/(hh + 1 + qq] rather than gddist[n, q, del], as inequalities on q and del would translate into positivity of hh and qq. Once parameters are fitted, we can easily recover q and del.

gbdist[n_, q_, del_] := 
 ParameterMixtureDistribution[BinomialDistribution[n, p], 
  Distributed[p, 
   BetaDistribution[q (-1 + 1/del), ((-1 + q) (-1 + del))/del]]]

loglikelihood[data_, {qq_, hh_}] := 
 Sum[LogLikelihood[
   gbdist[Part[data, i, 1], hh/(qq + hh), 
    1/(hh + 1 + qq)], {Part[data, i, 2]}], {i, Length[data]}]

 antdata = {{200, 40}, {200, 78}, {200, 97}, {200, 81}, {200, 96}};

One last thing we've got to do is to replace Log[Beta[a, b]] with combination of LogGamma for the same reason of extending the range of applicability of machine precision:

ll = loglikelihood[antdata, {qq, hh}] /. 
   HoldPattern[Log[Beta[a_, b_]]] :> 
    With[{aa = Simplify[a], bb = Simplify[b]}, 
     LogGamma[aa] + LogGamma[bb] - LogGamma[Simplify[aa + bb]]];

Now NMaximize has not issue finding the solution:

In[325]:= {maxVal, params} = 
 NMaximize[{ll, hh > 0 && qq > 0}, {qq, hh}]

Out[325]= {-22.3822, {qq -> 13.5332, hh -> 8.66654}}

In[326]:= {hh/(qq + hh), 1/(hh + 1 + qq)} /. params

Out[326]= {0.390389, 0.0431039}

I hope you will give Mathematica a chance, and continue using it.

Sincerely, Oleksandr Pavlyk

Dear Oleksandr,

I am sure your solution is very clever. But my problem is that I am not able to reproduce something like that the next time I have a similar problem.

I am not a mathematician but a biologist who work with applied mathematics. For 20+ years Mathematica has helped me to develop mathematical models to population biological questions and apply them to data.

I am extremely sad that it now will be more difficult for me take advantage of the coming new developments in Mathematica. I will be forced to live in a V11.2 world.

Please rethink this decision and get back to the philosophy of providing a simple to use software for users who are not mathematicians but use mathematics to solve applied problems.

Christian

I think there is no need to stay in 11.2 world, you just need to remember to ask Mathematica to use extended precision for its calculations.

Your data coincidentally were given as exact numbers (integers), so it was sufficient to simply use WorkingPrecision->20. If the data itself contains machine precision floating point numbers, it will necessary to execute epdata = SetPrecision[data, 20] and remember to use epdata in place of data.

The only important takeaway from my solution that I would recommend for you to apply going forward to replace Log[PDF[dist, x]] with LogLikelihood[dist, {x}].

Oleksandr

with other words, we have regressed into at state of int, float and double.

What used to be so great about Mathematica was that we could focus on the questions and the relevant math, and not have to worry about how the computer actually did it.

Please - MAKE MATHEMATICA GREAT AGAIN

or at least provide us with an option button so that some of us can continue to choose simplicity over speed.

Christian

The purpose of machine numbers is to give you access to the underlying machine arithmetic. Most machines are optimized for double arithmetic, so this gives you access to efficiency. Underflowing to zero is part of machine arithmetic, so trapping underflow and substituting a different kind of number was a troublesome violation of the basic design.

Using machine numbers inappropriately has always been a problem in Mathematica, as in other computing contexts. You're on shaky ground when you use machine numbers as input to symbolic machinery like Solve or Integrate, and the usual numeric hazards are also present. Fortunately, these are easily avoided. For symbolic calculation, use exact numbers, like 3/2 in place of 1.5. For safer, more rigorous numerics, use controlled precision like 1.5`10 .

You're in control here. Underflowing to zero improves that control: rather than arbitrarily switching the arithmetic domain, it keeps your calculation in the machine domain when that's the way you've formulated your problem.

I believe I may too have similar kind of problem in Mathematica 11.3. enter image description here

In my students' old laptop the same input gives a different correct (the one to be expected) result enter image description here

Can you please suggest a fix?

Thanks. Pallavi

The same happened with Mathematica 10 version and so I installed Mathematica 11.3 it too had the same problem. But worked fine in an old laptop.

My students old laptop gives the same output (which is more correct) but with some error message now in Mathematica 11.3. Please see the output: enter image description here

It would be great help if you can suggest a solution to it in my new machine with Mathematica 11.3 installed. Thanks, Pallavi

It would be easier for us to help if you posted code we could copy rather than an image.

I suggest using controlled precision: replace each machine number (numbers with a decimal point, like 1.5) with a controlled precision number like 1.5`20.

Dear John, Thanks for the reply. I tried what you suggested but that didn't help when {p=5} in the sum from p=1 to p=5 (it just removed the error messages), bit it just works fine when the upper limit is changed to p=6. Sharing the mathematica notebook. So the problem remains and it may pop up again in a corresponding higher order calculation.

The attached notebook does not make clear what is the problem when using higher precision. It seems to give a fine numeric result. When I do the substitutions using machine numbers I see essentially the same numeric result, now with warnings that there were (apparently benign) underflows to zero.

The problem is the 2.141507701982923310^363 result for termn5* at some parameter values, it should be something like .999993 (shown on old laptop, for the same set of parameters). How can I get this result on evaluating by: termn5/.{set of parameters}

So what should I do? using `20 didn't help as shown in the notebook shared now. Thanks. Pallavi

I, like Daniel, see similar results from controlled and machine precision. Why do you trust the old laptop's results?

p = {Subscript[T, 1] -> 5.0`20, Subscript[n, 1] -> 10.0`20, 
  Subscript[c, 1] -> 1.4`20, Subscript[T, 2] -> 5.0`20, 
  Subscript[n, 2] -> 10.0`20, Subscript[c, 2] -> 1.4`20, 
  Subscript[T, 3] -> 5.0`20, Subscript[n, 3] -> 10.0`20, 
  Subscript[c, 3] -> 1.4`20, Subscript[T, 4] -> 5.0`20, 
  Subscript[n, 4] -> 10.0`20, Subscript[c, 4] -> 1.4`20, 
  Subscript[T, 5] -> 5.0`20, Subscript[n, 5] -> 10.0`20, 
  Subscript[c, 5] -> 1.4`20, k -> 0.01`20}

Looking at it term by term:

lmn5 = List @@ termn5 /. p

I won't show the list with a lot of tiny terms along with some enormous ones, but:

Position[lmn5, Max[lmn5]]
(* {{55}} *)

So, the 55th term dominates.

termn5[[55]] /. p

(* 2.141507701982923310^363 )

Looking at it symbolically:

termn5[[55]]
(* (E^(3 \[Pi]^2 (16/Subscript[c, 4] + 25/Subscript[c, 
    5])) k^4 Subscript[c, 4] Subscript[c, 5])/(3600 \[Pi]^4) *)

Suppress the evaluation of the exponential by turning E into (undefined) e, then substitute parameters:

termn5[[55]] /. E -> e /. p

(* 5.58925700532813804710^-14 e^867.11524380999365008 )

So, you have an enormous exponential term. The exponent arises in what looks like a simple, stable calculation. Your result doesn't appear to reflect a problem in numerical calculation.

Dear John, Thanks for your answer. Now the file that I shared with you yesterday for

N[termn6 /. {Subscript[T, 1] -> 5.020, Subscript[n, 1] -> 1020, Subscript[c, 1] -> 1.420, Subscript[T, 2] -> 5.020, Subscript[n, 2] -> 1020, Subscript[c, 2] -> 1.420, Subscript[T, 3] -> 5.020, Subscript[n, 3] -> 1020, Subscript[c, 3] -> 1.420, Subscript[T, 4] -> 5.020, Subscript[n, 4] -> 1020, Subscript[c, 4] -> 1.420, Subscript[T, 5] -> 5.020, Subscript[n, 5] -> 1020, Subscript[c, 5] -> 1.420, Subscript[T, 6] -> 5.020, Subscript[n, 6] -> 1020, Subscript[c, 6] -> 1.420, k -> 0.01`20}]

its no more giving 0.999993 but 4.76145665529708310^546 today*

What might have changed in a day? All the results are from mathematica 11.3. I trust 0.999993 because its more of an answer that one would expect physically as it enters into a probability distribution function. Sharing today's file again. Daniel please see to it too termn6 is giving the same unphysical answer today as termn5 (today & yesterday).

Daniel & John, I would be more than grateful if I can retreive my old results (at least for termn6 of yesterday). Thanks, Pallavi

First of all thank you all for your replies. Specifically the term my term analysis. A term as

(E^(3 [Pi]^2 (16/Subscript[c, 4] + 25/Subscript[c, 5])) k^4 Subscript[c, 4] Subscript[c, 5])/(3600 [Pi]^4)

should not occur at first place it should have been mapped as per the rule: E^(X/Subscript[c, a] + Y/Subscript[c, b]) -> 2^(Subscript[n, a]/2) (( E^Subscript[T, a] Sqrt[1 + (4 X )/Subscript[c, a]] Csch[Subscript[T, a] Sqrt[1 + (4 X )/Subscript[c, a]]])/( 2 + (4 X )/Subscript[c, a] + 2 Sqrt[1 + (4 X )/Subscript[c, a]] Coth[Subscript[T, a] Sqrt[1 + (4 X )/Subscript[c, a]]]))^( Subscript[n, a]/2)* 2^( Subscript[n, b]/2) (( E^Subscript[T, b] Sqrt[1 + (4 Y )/Subscript[c, b]] Csch[Subscript[T, b] Sqrt[1 + (4 Y)/Subscript[c, b]]])/( 2 + (4 Y )/Subscript[c, b] + 2 Sqrt[1 + (4 Y )/Subscript[c, b]] Coth[Subscript[T, b] Sqrt[1 + (4 Y )/Subscript[c, b]]]))^( Subscript[n, b]/2)

which didn't happen. So I changed one of the integers like 6 to 6.0 in Subscript[G, p] -> Subscript[c, p]/(6 [Pi]^2), and the mapping went through and I can get the desired result. Thanks again.

What is surprising is one day the mapping works and other day it doesn't. Nevertheless, thanks.

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract