bug#58929: 29.0.50; Calc: finding roots utpn doesn't work as advertised

[Matt Armstrong <matt@rfc20.org>]

[-- Attachment #1: Type: text/plain, Size: 1130 bytes --]

Mattias Engdegård <mattiase@acm.org> writes:

> Obviously the claim to guarantee a solution for "any initial guess" is
> fanciful.

Yes.  See attached patch for what I wish the docs had said all along
(would have saved me a lot of time).


> Unfortunately straight Newton won't converge for this function (erf,
> essentially) if you stray too far from the origin (ie, the mean) so I
> suggest you always use that as a starting guess no matter what point
> you are trying to find the inverse for.
>
> Another way is to use an interval as starting guess instead of a
> single point; that seems to make Calc switch tactics and persevere.

I wonder if the implementation can be improved.
https://www.gnu.org/software/emacs/manual/html_node/calc/Root-Finding.html
suggests that there are some cases where calc will switch over to the
bisection method.  Perhaps the heuristics can be tweaked?


> One thing that can happen is running out of stack, because many of the
> implementations like math-newton-root are recursive. This should be
> easy to remedy as they all seem to be tail calls.

Good to know.


[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #2: 0002-In-calc.texi-provide-guidance-for-taking-the-inverse.patch --]
[-- Type: text/x-diff, Size: 2134 bytes --]

From ba68ca5671f9772994bc5c0b2b0b9216a750b692 Mon Sep 17 00:00:00 2001
From: Matt Armstrong <matt@rfc20.org>
Date: Tue, 1 Nov 2022 14:51:36 -0700
Subject: [PATCH 2/2] In calc.texi provide guidance for taking the inverse of
 PDFs.

* doc/misc/calc.texi (Probability Distribution Functions): Stop
claiming that the 'a R' function will "always work" for any guess.
Explain the problem with the PDF functions being very flat.  Provide
recommended guess values for each function.  bug#58929
---
 doc/misc/calc.texi | 36 +++++++++++++++++++++++++++++++++---
 1 file changed, 33 insertions(+), 3 deletions(-)

diff --git a/doc/misc/calc.texi b/doc/misc/calc.texi
index 89a340e7343..e6b416a447a 100644
--- a/doc/misc/calc.texi
+++ b/doc/misc/calc.texi
@@ -19456,9 +19456,39 @@ Probability Distribution Functions
 
 While Calc does not provide inverses of the probability distribution
 functions, the @kbd{a R} command can be used to solve for the inverse.
-Since the distribution functions are monotonic, @kbd{a R} is guaranteed
-to be able to find a solution given any initial guess.
-@xref{Numerical Solutions}.
+@xref{Root Finding}.
+
+The distribution functions, while monotonic, are also be very flat,
+which can cause problems for the Newton's root-finding algorithm used
+by @kbd{a R}.  The problem can be avoided by providing a good starting
+guess.  Prefer guesses that occur where the slope of the plotted
+function is not flat.  Alternatively, use an interval range for the
+guess.  The table below provides a recommendation for each of the
+probability distribution functions.
+
+@multitable @columnfractions .25 .25 .25
+@headitem Function
+@tab When solving for
+@tab Recommended initial guess
+@item @samp{utpb(x,n,p)}
+@tab @expr{x}
+@tab @code{0.5}
+@item @samp{utpc(x,c)}
+@tab @expr{x}
+@tab @code{1.0}
+@item @samp{utpf(F,v1,v2)}
+@tab @expr{F}
+@tab @code{1.0}
+@item @samp{utpn(x,m,s)}
+@tab @expr{x}
+@tab @expr{m}
+@item @samp{utpp(n,x)}
+@tab @expr{n}
+@tab @expr{x}
+@item @samp{utpt(t,v)}
+@tab @expr{t}
+@tab @code{[0..1e4]}
+@end multitable
 
 @node Matrix Functions
 @chapter Vector/Matrix Functions
-- 
2.35.1