#3928 instability of numerical integration of the lognormal function

[Pavel Murat]

I'd like to report a numerical instability in the maxima integration of the lognormal function. In the example below, small variations of the parameter sigma change the outcome of the integration from fast convergence to failure due to exhausted memory - the diagnostics and the build information
follow below. The nature of failure suggests that the issue may be more general than just the lognormal function and affect other cases. Just to comment, I'm getting maxima as a part of Fedora distribution, so the version I'm using is not the latest one.

-- thanks for the great program, Pavel

/* ------------------------------------------------------ test.max */
sigma: 0.7147 ;  
/* 0.7  : integration OK, 
  0.714: integration slow, 
  0.7147: integration chokes, running out of memory */
mu   : 0.2 ;

f(x,mu,sigma) := 1/(sqrt(2*%pi)*sigma*x)*exp(-(log(x)-mu)^2/(2*sigma^2)) ;
s    : float(integrate(f(x,mu,sigma),x,0,10)) ;
/* --------------------------------------------------------- end test.max */

1) here is the output of the failed integration:

> cat test.max | maxima 
Maxima 5.43.2 http://maxima.sourceforge.net
using Lisp SBCL 2.0.1-5.fc34
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) 
(%o1)                               0.7147
(%o2)                                 0.2
(%i3) 
                                                                    2
                                     1               - (log(x) - mu)
(%o3)   f(x, mu, sigma) := (-------------------) exp(----------------)
                            sqrt(2 %pi) sigma x                 2
                                                         2 sigma
(%i4) 
rat: replaced -0.9788641882517151 by -50000000/51079609 = -0.9788641882517151

rat: replaced -0.2 by -1/5 = -0.2

rat: replaced -0.9788641882517151 by -50000000/51079609 = -0.9788641882517151

rat: replaced -0.2 by -1/5 = -0.2

rat: replaced -0.9788641882517151 by -50000000/51079609 = -0.9788641882517151

rat: replaced -0.2 by -1/5 = -0.2
Heap exhausted during garbage collection: 0 bytes available, 16 requested.
Gen  Boxed   Code    Raw  LgBox LgCode  LgRaw  Pin       Alloc     Waste        Trig      WP GCs Mem-age
 0    3284      0      0      0      0      0    3   107544592     65520    64542458    3284   1  0.0000
 1    4103      0      0      0      0      0    4   134398832     48272    57743994    4103   1  0.3498
 2    2236      0      0      0      0      0    6    73175408     93840    65586186    2236   1  0.7496
 3    8363      0      0      0      0      0    3   273980928     57856   284718346    8363   1  0.0000
 4   13698      0     32      0      0      0   20   449749632    155008     2000000   13730   0  0.0000
 5       0      0      0      0      0      0    0           0         0     2000000       0   0  0.0000
 6     658      2    295     62      0     35    0    33786480    685456     2000000    1052   0  0.0000
           Total bytes allocated    =    1072635872
           Dynamic-space-size bytes =    1073741824
GC control variables:
   *GC-INHIBIT* = true
   *GC-PENDING* = true
   *STOP-FOR-GC-PENDING* = false
fatal error encountered in SBCL pid 560665(tid 0x7f8855ef9180):
Heap exhausted, game over.

2) and her is the build info:

-------------------------------------------------------------
Maxima version: "5.43.2"
Maxima build date: "2021-07-22 19:02:43"
Host type: "x86_64-redhat-linux-gnu"
Lisp implementation type: "SBCL"
Lisp implementation version: "2.0.1-5.fc34"
User dir: "/home/murat/.maxima"
Temp dir: "/tmp"
Object dir: "/home/murat/.maxima/binary/5_43_2/sbcl/2_0_1_5_fc34"
Frontend: false
-------------------------------------------------------------

[Stavros Macrakis]

integrate performs symbolic, not numeric, integration. Wrapping it in float(...) simply converts the symbolic value to a floating value after first calculating it symbolically. In some cases, floating point exponents can cause problems because they get converted to very messy rationals.

That said, I cannot reproduce your problem in Maxima 5.45.1.

In any case, it is unnecessary to substitute mu and sigma before evaluating the integral. Maxima is perfectly happy to integrate symbolically:

integrate( 1/(sqrt(2*%pi)*sigma*x)*exp(-(log(x)-mu)^2/(2*sigma^2)), x, 0, 10)
    (Is sigma positive or negative? pos)
((sqrt(%pi)*sigma)/sqrt(2)-(sqrt(%pi)*erf((sqrt(2)*mu-sqrt(2)*log(10))/(2*sigma))*sigma)/sqrt(2))/(sqrt(2)*sqrt(%pi)*sigma)

You can then substitute numerical values into that.