Two equivalent ways of inputting a function involving gamma function yield different result?

[KvS]

Hi all,

I didn't know how to make the title any clearer, I'm sorry. The
problem is as follows (on Sage v 4.6). Running the following piece of
code (BetaLP1 is some class):

f=BetaLP1.getkthMomentAt(1,1)
print 'f:',f
print f(0)

yields the following error:

f: x |--> e^(9/200*x^2 + 4/3*sqrt(pi)*gamma(-x + 1)/gamma(-x - 1/2) +
4/3*sqrt(pi)*gamma(x + 1)/gamma(x - 1/2) + 4/3)
Traceback (most recent call last): return [(el[1],el[0]) for el
in lst]
File "", line 1, in <module>

File "/tmp/tmpJG9qg9/___code___.py", line 16, in <module>
print f(_sage_const_0 )
File "expression.pyx", line 3473, in
sage.symbolic.expression.Expression.__call__ (sage/symbolic/
expression.cpp:15647)
File "/usr/local/Sage/sage-4.6-linux-32bit-ubuntu_10.04_lts-i686-
Linux/local/lib/python2.6/site-packages/sage/symbolic/callable.py",
line 451, in _call_element_
return SR(_the_element.substitute(**d))
File "expression.pyx", line 3324, in
sage.symbolic.expression.Expression.substitute (sage/symbolic/
expression.cpp:15026)
File "pynac.pyx", line 1047, in
sage.symbolic.pynac.py_doublefactorial (sage/symbolic/pynac.cpp:9456)
ValueError: argument must be >= -1

But if I then copy formula for f from the above output and paste it in
a new cell:

f(x)=e^(9/200*x^2 + 4/3*sqrt(pi)*gamma(-x + 1)/gamma(-x - 1/2) +
4/3*sqrt(pi)*gamma(x + 1)/gamma(x - 1/2) + 4/3)
print f(0)

it yields 1. So in that first cell something happens which I don't
understand. How can both f's evaluate different, while their
respective formulas are the same? Is there some difference between the
two that I can't see from their respective formulas?

Many thanks in advance, Kees

[KvS]

In addition to the previous post, if I change the contents of the
first cell to

f=BetaLP1.getkthMomentAt(1,1)
g=symbolic_expression(str(f(x))).function(x)
print g(0)

it does yield 1 as well. Any hints?

Thanks, Kees

[Jason Grout]

Can you post the code for BetaLP1.getkthMomentAt(1,1) so we can try it?

Thanks,

Jason

[KvS]

Here: http://pastebin.com/1eVSXsXW you can find the relevant part of
the class code, and you can call it e.g. with (there are some
restrictions on the parameters namely):

BetaLP1=BetaLP(mu=0,sigma=3/10,alpha1=1,beta1=1,lambda1=5/2,c1=1,alpha2=1,beta2=1,lambda2=5/2,c2=1)

[KvS]

On Dec 23, 6:42 pm, Jason Grout <jason-s...@creativetrax.com> wrote:

plus: thanks in advance, of course!

[Jason Grout]

Yes, something does seem weird. When I modify your getkthMomentAt
function to:

def getkthMomentAt(self,k,t):
charFunc=self.getCharacteristicExponent()
i=sqrt(-1)
f=(exp(-t*charFunc(-i*x))).function(x)
for _ in range(k):
f=f.derivative()
return f

then I get the following odd sequence of results:

sage: mygammaexpr=f.operands()[0].operands()[0].operands()[1].operands()[0]
sage: mygammaexpr
gamma(-x - 1/2)
sage: mygammaexpr(x=0)

Traceback (most recent call last):

...

ValueError: argument must be >= -1

sage: mygammaexpr
gamma(-x - 1/2)
sage: gamma(-x-1/2) - mygammaexpr
0
sage: mygammaexpr(x=0)
-2*sqrt(pi)

Definitely weird. Apparently subtracting mygammaexpr from gamma(-x-1/2)
makes mygammaexpr suddenly realize that it's okay to evaluate it at x=0.

Jason

[Jason Grout]

On 12/23/10 7:22 AM, Jason Grout wrote:
> Definitely weird. Apparently subtracting mygammaexpr from gamma(-x-1/2)
> makes mygammaexpr suddenly realize that it's okay to evaluate it at x=0.

And lest we think the errors are something from Kees's install, those
computations were done on demo.sagenb.org:
http://demo.sagenb.org/home/pub/66/

Jason

[KvS]

Thanks a lot for your efforts Jason! I can do with the workaround, but
will also keep an eye on this thread to see if it'll be fixed. I
assume if some reporting is required you will do that?

Cheers, Kees