Opened 11 years ago

#9130 closed enhancement

Reported by: Owned by: kcrisman burcin major sage-5.0 symbolics special function, pynac, sd35.5 Cernay2012 benjaminfjones Karen T. Kohl, Burcin Erocal Benjamin Jones N/A

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

comment:1 Changed 11 years ago by fredrik.johansson

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.

comment:2 Changed 11 years ago by kcrisman

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
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.

comment:3 follow-up: ↓ 4 Changed 11 years ago by 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.

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

comment:4 in reply to: ↑ 3 ; follow-up: ↓ 5 Changed 11 years ago by kcrisman

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)?

comment:5 in reply to: ↑ 4 Changed 11 years ago by 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.

comment:7 Changed 11 years ago by ktkohl

• Status changed from new to needs_work

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).
------------------------------------------------------------

```

Changed 10 years ago by burcin

fix segfault in py_float

comment:8 Changed 10 years ago by burcin

• Authors set to Karen T. Kohl, Burcin Erocal
• Keywords beginner removed
• Milestone changed from sage-5.0 to sage-4.7.1

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).

comment:11 Changed 10 years ago by ktkohl

• Status changed from needs_work to needs_review

Changed 10 years ago by ktkohl

include beta in random_tests

comment:12 Changed 10 years ago by ktkohl

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

comment:13 Changed 10 years ago by kcrisman

• Description modified (diff)

comment:14 Changed 10 years ago by benjaminfjones

• Status changed from needs_review to needs_work

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

comment:15 Changed 10 years ago by benjaminfjones

• Reviewers set to Benjamin Jones

comment:16 Changed 10 years ago by ktkohl

• Status changed from needs_work to needs_review

comment:17 Changed 10 years ago by kcrisman

• Description modified (diff)

Changed 10 years ago by benjaminfjones

combined reviewer patch folding up previously uploaded 5 patches

Note: See TracTickets for help on using tickets.