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 -------------------------------------------------------------
integrate
performs symbolic, not numeric, integration. Wrapping it infloat(...)
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:
You can then substitute numerical values into that.
Last edit: Stavros Macrakis 2022-01-27