Access to beta function

[kcrisman]

Although Maxima has the beta function, Sage doesn't:

sage: a, b, c = var('a b c') 
sage: assume(a > 0) 
sage: assume(b > 0) 
sage: x = var('x') 
sage: beta_dist = x**(a-1) * (1 - x)**(b-1) 
sage: c = integral(beta_dist, x, 0, 1) 
sage: c
beta(a, b)
sage: c(a=.5,b=.5)
beta(0.500000000000000, 0.500000000000000)
sage: c(a=.5,b=.5).n()
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)

/Users/karl-dietercrisman/<ipython console> in <module>()

/Users/karl-dietercrisman/Desktop/sage-4.4.2/local/lib/python2.6/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression.n (sage/symbolic/expression.cpp:17042)()

TypeError: cannot evaluate symbolic expression numerically

This is is Ginac, though, and there is even room for defining it in symbolic/expression.pyx . It probably is also included in some of our other libraries, as a standard special function.

---

Apply

• Patches at Implement a symbolic version of the arg function #4498 and Mark random symbolic expression doctests with #random #12507
• attachment: trac_9130-beta_function.3.patch
• attachment: trac_9130-py_float_segfault.take2.patch
• attachment: trac_9130-reviewer.patch
• attachment: trac_9130-random-tests.2.patch

Depends on #4498
Depends on #12507

CC: @benjaminfjones

Component: symbolics

Keywords: special function, pynac, sd35.5 Cernay2012

Author: Karen Kohl, Burcin Erocal, Karl-Dieter Crisman

Reviewer: Benjamin Jones, Burcin Erocal, Karl-Dieter Crisman

Merged: sage-5.0.beta6

Issue created by migration from https://trac.sagemath.org/ticket/9130

[fredrik-johansson]

comment:1

For numerical evaluation, mpmath has beta and also the generalized incomplete beta function for complex arguments. But it's probably easy to do the complete beta function directly with Ginac.

Simplification for rational arguments (beta(0.5,0.5) = pi) would be nice.

Unless someone else wants to work on this, I might have a stab at it within a couple of days.

[kcrisman]

comment:2

I wasn't suggesting which of the several packages in Sage should be used for numerical evaluation, though mpmath did come to mind :-)

I don't think that beta(0.5,0.5) would work, given that

sage: gamma(0.5)
1.77245385090552
sage: gamma(1/2)
sqrt(pi)

but beta(1/2,1/2) becoming pi should work fine once we have a symbolic wrapper (with or without Ginac):

sage: maxima_console()
Maxima 5.21.1 http://maxima.sourceforge.net
using Lisp ECL 10.4.1
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) beta(1/2,1/2);
(%o1)                                 %pi

Please do try to add this! We definitely often get email asking for various special functions both symbolically and numerically. Also, the more examples we have, the easier it is to finish off the rest of them by cut and paste.

[burcin]

comment:3

GiNaC has a beta function, so this can probably be solved simply by wrapping that. See #8864 for an example.

Though I don't know why the beta() function in sage/symbolic/expression.pyx is commented. Maybe there is something I'm missing.

Fredrik, it would be great if you can do this. I'd be happy to answer questions if anything goes wrong.

[kcrisman]

comment:4

Replying to @burcin:

GiNaC has a beta function, so this can probably be solved simply by wrapping that. See #8864 for an example.

Though I don't know why the beta() function in sage/symbolic/expression.pyx is commented. Maybe there is something I'm missing.

I think the same reason the psi and psi2 ones are commented - when those were implemented, they didn't notice that they had been commented earlier. This was probably pretty early in the conversion, maybe when William was dealing with CLN (whatever that is)?

[burcin]

comment:5

Replying to @kcrisman:

Replying to @burcin:

GiNaC has a beta function, so this can probably be solved simply by wrapping that. See #8864 for an example.

Though I don't know why the beta() function in sage/symbolic/expression.pyx is commented. Maybe there is something I'm missing.

I think the same reason the psi and psi2 ones are commented - when those were implemented, they didn't notice that they had been commented earlier. This was probably pretty early in the conversion, maybe when William was dealing with CLN (whatever that is)?

psi() and psi2() were commented because at the time there was no method defined to numerically evaluate those. This is not the case for beta() however. Here is the evalf method (from line 227 of ginac/inifcns_gamma.cpp):

	if (is_exactly_a<numeric>(x) && is_exactly_a<numeric>(y)) {
		try {
			return exp(lgamma(ex_to<numeric>(x))+lgamma(ex_to<numeric>(y))-lgamma(ex_to<numeric>(x+y)));
		} catch (const dunno &e) { }
	}

We'll find out when someone tries this out I suppose.

[sagetrac-ktkohl]

Attachment: trac_9130_beta_function.patch.gz

[sagetrac-ktkohl]

comment:7

Added ginac wrapper for beta function.
There is a problem when one argument is a float.

--Karen

sage: beta(3,2.0)


------------------------------------------------------------
Unhandled SIGSEGV: A segmentation fault occurred in Sage.
This probably occurred because a *compiled* component
of Sage has a bug in it (typically accessing invalid memory)
or is not properly wrapped with _sig_on, _sig_off.
You might want to run Sage under gdb with 'sage -gdb' to debug this.
Sage will now terminate (sorry).
------------------------------------------------------------

[burcin]

Attachment: trac_9130-py_float_segfault.patch.gz

fix segfault in py_float

[burcin]

Author: Karen T. Kohl, Burcin Erocal

[burcin]

comment:8

Hi Karen,

sorry for taking so long to look at this. It seems that I forgot to check for a NULL pointer in py_float(). attachment: trac_9130-py_float_segfault.patch fixes the segfault.

sage: from sage.functions.other import beta
sage: beta(3,2.0)
0.0833333333333333

Will you have time to finish the patch?

You need to add an import statement to sage/functions/all.py so beta is available at the command line. The documentation also needs some care. IIRC there should be an empty line after INPUT, OUTPUT and EXAMPLES. The statement

        It is computed by various libraries within Sage, depending on
        the input type.

is too vague. We should either remove it or explain how GiNaC evaluates this (see beta_eval() and beta_evalf() in inifcns_gamma.cpp).

[sagetrac-ktkohl]

Changed keywords from special function, pynac to special function, pynac, sd35.5

[sagetrac-ktkohl]

Attachment: trac_9130_beta_function.2.patch.gz

[sagetrac-ktkohl]

include beta in random_tests

[sagetrac-ktkohl]

comment:12

Attachment: trac_9130_beta_function.3.patch.gz

Apply all patches in the order they were added:

• trac_9130_beta_function.patch

• trac_9130-py_float_segfault.patch

• trac_9130_beta_function.2.patch

• trac_9130_beta_function.3.patch

[benjaminfjones]

comment:14

I think we discovered that the only complex inputs that the code accepts are ones where one of the parameters is equal to 1. In that case beta(1,x) = 1/x is used to compute the result. Looks like GiNaC can't handle complex inputs at all (or perhaps complex numbers aren't being passed to GiNaC in a way it understands).

On the other hand, mpmath does support evaluation at arbitrary precision complex numbers so that could be a useful enhancement that could take place in a new ticket.

I would change the docstrings to clearly indicate that

1. only real inputs are accepted (for now)
2. beta(1,x) = beta(x,1) = 1/x simplification is automatically applied

[benjaminfjones]

Dependencies: 4498

[benjaminfjones]

Changed dependencies from 4498 to #4498

[benjaminfjones]

Attachment: trac_9130-beta_function.2.patch.gz

[benjaminfjones]

Attachment: trac_9130-random-tests.patch.gz

fixes random tests after rebasing against #4498

[benjaminfjones]

comment:37

I rebased the patches against #4498. This meant moving the changes to sage/symbolic/random_tests.py out of attachment: trac_9130_beta_function.patch and into a new patch attachment: trac_9130-random-tests.patch that reflects the differences in the random tests after #4498 was applied and beta was added.

The tests in random_test.py are getting annoying. Now #9130 must depend on #4498 and in turn my #11888 must depend on #9130, etc. To finish #11143, I need to rebase the changes to random_tests.py to all 3 of these other symbolics patches.

Would it make sense to open a new ticket for changes needed to random_tests.py after all the new symbolic functions are added that will go into sage-5.0? This new ticket would depend on all the other tickets adding new functions. This way we wouldn't have to fix the random tests like this in every ticket that adds a new function, just once per release.

[jdemeyer]

comment:38

Replying to @benjaminfjones:

Would it make sense to open a new ticket for changes needed to random_tests.py after all the new symbolic functions are added that will go into sage-5.0? This new ticket would depend on all the other tickets adding new functions. This way we wouldn't have to fix the random tests like this in every ticket that adds a new function, just once per release.

I think that would be making it more complicated for yourself and for reviewers. It also means that, if just one of those tickets has a problem, none of them can be merged. But I'm not against your suggestion.

[burcin]

comment:39

Replying to @benjaminfjones:

I rebased the patches against #4498. This meant moving the changes to sage/symbolic/random_tests.py out of attachment: trac_9130_beta_function.patch and into a new patch attachment: trac_9130-random-tests.patch that reflects the differences in the random tests after #4498 was applied and beta was added.

Thank you for taking care of this.

The tests in random_test.py are getting annoying. Now #9130 must depend on #4498 and in turn my #11888 must depend on #9130, etc. To finish #11143, I need to rebase the changes to random_tests.py to all 3 of these other symbolics patches.

Would it make sense to open a new ticket for changes needed to random_tests.py after all the new symbolic functions are added that will go into sage-5.0? This new ticket would depend on all the other tickets adding new functions. This way we wouldn't have to fix the random tests like this in every ticket that adds a new function, just once per release.

At this stage I'd be happy to mark that test in random_tests.py with a #random tag so the doctesting framework ignores the output. If you open a new ticket for this and provide a patch, I promise a quick positive review. :)

[benjaminfjones]

comment:40

Replying to @burcin:
Thanks, marking the three "offending" doctests with a random tag is #12507.

[kcrisman]

comment:41

Just checked one last time - yes, everything is fine! A very minor quibble is that the (why here?) added prime pi file could use a few double backticks and one or two other things, but that is immaterial in this saga.

[jdemeyer]

Attachment: trac_9130-random-tests.2.patch.gz

[jdemeyer]

comment:42

Attachment: trac_9130-beta_function.3.patch.gz

Rebased again.

[jdemeyer]

Changed dependencies from #4498 to #4498, #12507

[kcrisman]

comment:43

I think if you would have merged this first, then #12507, there wouldn't have needed to be a rebase. Maybe the authors of #12507 should have put this as a dependency? I think we were all lax with it because it # randomized a doctest...

[jdemeyer]

Merged: sage-5.0.beta6

[fchapoton]

Changed author from Karen T. Kohl, Burcin Erocal, Karl-Dieter Crisman to Karen Kohl, Burcin Erocal, Karl-Dieter Crisman