I am pretty new to maxima, which I installed a few days ago via brew on MacOS ARM. ergo, I do not know whether I found a bug, or a feature, or I am doing something stupid altogether. The following example shows the unexpected error I get:
$ maxima
Maxima 5.46.0 https://maxima.sourceforge.io
using Lisp SBCL 2.3.3
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) load ("distrib");
(%o1) /opt/homebrew/Cellar/maxima/5.46.0_11/share/maxima/5.46.0/share/distrib/\
distrib.mac
(%i2) cdf_binomial(4,4,p);
(%o2) 1
(%i3) subst( [S=3], cdf_binomial(S,4,p));
3 2
(%o3) (1 - p) (p + p + p + 1)
(%i4) subst( [S=4], cdf_binomial(S,4,p));
beta_incomplete: beta_incomplete(0,5,1-p) is undefined.
-- an error. To debug this try: debugmode(true);
(%i5)
Am I doing something wrong here? (PS: My real use will be integrating over p up to 1, which also triggers this problem.) Are there standard workarounds?
Advice appreciated.
The error you're seeing is actually a subtle consequence of the way unevaluated conditional expressions are handled. (An example of an unevaluated conditional is returned by cdf_binomial(S, 4, p).) A simpler example of the bug (I'll create a bug report) is '(if 0 > 0 then 0/0 else inf) which bumps into a division by 0 error even though 0 > 0 is false.
As a workaround, try ev(myexpr, S = 4) where myexpr is assigned the value returned by cdf_binomial(S, 4, p).
This detour via ev is unnecessary if S has an assigned value -- in that case, cdf_binomial returns one case or another, and doesn't return an unevaluated conditional. The need for special handling comes into play when cdf_binomial was called with S not having an assigned value, so cdf_binomial didn't know which branch to take and therefore returned an unevaluated conditional.
EDIT: I also investigating equipping beta_incomplete_regularized with an additional identity for the first argument being zero. Looks like beta_incomplete_regularized(0, n, z) equals 1, at least when n is an integer. I'll look for a reference for that, but in the meantime one can tell Maxima to recognize that identity via tellsimp.
(%i2) matchdeclare (nn, integerp) $
(%i3) matchdeclare (zz, all) $
(%i4) tellsimp (beta_incomplete_regularized (0, nn, zz), 1) $
(%i5) beta_incomplete_regularized (0, 5, z);
(%o5) 1
(%i6) load (distrib);
(%o6) /usr/local/share/maxima/branch_5_46_base_1803_g61b2abdd9/s\
hare/distrib/distrib.mac
(%i7) cdf_binomial (S, 4, p);
(%o7) if S < 0 then 0 elseif S >= 4 then 1
else beta_incomplete_regularized(4 - floor(S), floor(S) + 1,
1 - p)
Now the substitution into the unevaluated conditional is okay because beta_incomplete_regularized doesn't complain.
(%i8) subst (S = 4, %);
(%o8) if 4 < 0 then 0 elseif 4 >= 4 then 1 else 1
The literal inequalities 4 < 0 and 4 >= 4 require an explicit evaluation (a long story which I'll omit) to derive false and true from them.
(%i9) ''%;
(%o9) 1
Now the call to quad_qags succeeds.
(%i10) quad_qags (subst (S = 4, cdf_binomial (S, 4, p)), p, 0, 1);
(%o10) [1.0, 1.110223024625157e-14, 21, 0]