possible bug in sage 5.11 SR operations...

[Emmanuel Charpentier]

Dear list,

Setup : sage 5.11 on Debian amd64 self-compiled (Debian's compilers and tuned Atlas library) 

trying to put on paper the elementary proof that the convolution of two normals is a normal, I stumbled on what I think to be a disastrous bug. What follows is shown in text mode, but running it in emacs with sage-view or in the notebook with the "typeset" option is much clearer...

var("x,t,mu,sigma")

assume(x,"real",t,"real",mu,"real",sigma,"real",sigma>0)

phi(x,mu,sigma)=e^-(((x-mu)^2)/(2*sigma^2))/sqrt(2*pi*sigma^2)

var("m1,m2,s1,s2")

assume(m1,"real",m2,"real",s1,"real",s2,"real",s1>0,s2>0)

g(x,m1,m2,s1,s2)=integrate(phi(t,m1,s1)*phi(x-t,m2,s2),t,-oo,oo)

prints :

(x, m1, m2, s1, s2) |--> sqrt(1/2)*sqrt(pi)*e^(-1/2*(m1^2 + 2*m1*m2 + m2^2 - 2*(m1 + m2)*x + x^2)/(s1^2 + s2^2))/(sqrt(pi*s1^2)*sqrt(pi*s2^2)*sqrt((s1^2 + s2^2)/(s1^2*s2^2)))

Ugly (in text mode), but exact.

h=g(x,m1,m2,s1,s2) # to get a symbolic expression
print h.denominator()

This prints

sqrt(pi*s1^2)*sqrt(pi*s2^2)

which is simple, but obviously FALSE. And, by the way,

maxima.denom(h).sage()

gives the much more reasonable :

sqrt(2)*sqrt(pi*s1^2)*sqrt(pi*s2^2)*sqrt((s1^2 + s2^2)/(s1^2*s2^2))

which, simplify_radical()ized (legitimate under our assumptions), gives the expected :

sqrt(2)*pi*sqrt(s1^2 + s2^2)

I think that this is a bug (what is used for this kind of expression ? Sympy ?). Has it been already reported ? Should I report it ?

HTH,

                                                                    Emmanuel Charpentier

[Burcin Erocal]

> which is simple, but obviously *FALSE*. And, by the way,

>
> maxima.denom(h).sage()
>
> gives the much more reasonable :
>
> sqrt(2)*sqrt(pi*s1^2)*sqrt(pi*s2^2)*sqrt((s1^2 + s2^2)/(s1^2*s2^2))

You can get a similar response from Sage (or the underlying library
GiNaC/pynac) by setting normalize=False:

sage: var("x,t,mu,sigma")
(x, t, mu, sigma)
sage: assume(x,"real",t,"real",mu,"real",sigma,"real",sigma>0)
sage: var("m1,m2,s1,s2")
(m1, m2, s1, s2)
sage: phi(x,mu,sigma)=e^-(((x-mu)^2)/(2*sigma^2))/sqrt(2*pi*sigma^2)
sage: ex = integrate(phi(t,m1,s1)*phi(x-t,m2,s2),t,-oo,oo)
sage: ex.denominator(normalize=False)
sqrt(pi*s1^2)*sqrt(pi*s2^2)*sqrt((s1^2 + s2^2)/(s1^2*s2^2))

Our .denominator() method returns the denominator of the "normalized"
expression by default:

sage: ex.normalize()
sqrt(1/2)*sqrt(pi)*sqrt(s1^2*s2^2/(s1^2 + s2^2))*e^(-1/2*(m1^2 +
2*m1*m2 + m2^2 - 2*m1*x - 2*m2*x + x^2)/(s1^2 +
s2^2))/(sqrt(pi*s1^2)*sqrt(pi*s2^2))


Cheers,
Burcin

[Emmanuel Charpentier]

Dear Burcin, dear list

Le mardi 27 août 2013 11:16:18 UTC+2, Burcin Erocal a écrit :

On Sat, 24 Aug 2013 14:16:01 -0700 (PDT)
Emmanuel Charpentier <emanuel.c...@gmail.com> wrote:

> Dear list,

 

[ Snip.. ]

 

> h=g(x,m1,m2,s1,s2) # to get a symbolic expression
> print h.denominator()
>
> This prints
>
> sqrt(pi*s1^2)*sqrt(pi*s2^2)
>
> which is simple, but obviously *FALSE*. And, by the way,
>
> maxima.denom(h).sage()
>
> gives the much more reasonable :
>
> sqrt(2)*sqrt(pi*s1^2)*sqrt(pi*s2^2)*sqrt((s1^2 + s2^2)/(s1^2*s2^2))

You can get a similar response from Sage (or the underlying library
GiNaC/pynac) by setting normalize=False:

sage: var("x,t,mu,sigma")                                          
(x, t, mu, sigma)
sage: assume(x,"real",t,"real",mu,"real",sigma,"real",sigma>0)
sage: var("m1,m2,s1,s2")
(m1, m2, s1, s2)
sage: phi(x,mu,sigma)=e^-(((x-mu)^2)/(2*sigma^2))/sqrt(2*pi*sigma^2)
sage: ex = integrate(phi(t,m1,s1)*phi(x-t,m2,s2),t,-oo,oo)
sage: ex.denominator(normalize=False)
sqrt(pi*s1^2)*sqrt(pi*s2^2)*sqrt((s1^2 + s2^2)/(s1^2*s2^2))

Our .denominator() method returns the denominator of the "normalized"
expression by default:

Thank you. I was under the (false) impression that Sage used Maxima's CRE as a canonical representation. Hence two questions :

1) Shouldn't this "normalized form" be used when one asks for h ?

2) Where is this "normalized" form discussed (possibly against CRE) ? I suppose that, since Sage uses Maxima for a lot of the symbolics, changing its internal representation to one different from Maxima's (or, s case might be, deciding to stay with a previous choice)  has been deemed worth of it, and I'd like to understand  the rationale for this decision.

 

Again, thank you very much !

                                                               Emmanuel Charpentier.

[Burcin Erocal]

On Tue, 27 Aug 2013 11:11:35 -0700 (PDT)

Symbolic operations try to be lazy as much as possible. Normalization
is expensive. We don't want to do this until asked.

> 2) Where is this "normalized" form discussed (possibly against CRE) ?

See 5.8.1 here:

http://www.ginac.de/tutorial/Rational-expressions.html

> I suppose that, since Sage uses Maxima for a lot of the symbolics,
> changing its internal representation to one different from Maxima's
> (or, s case might be, deciding to stay with a previous choice) has
> been deemed worth of it, and I'd like to understand the rationale
> for this decision.

The symbolics subsystem in Sage is supposed to allow development of
native functionality in Sage. We are not aiming to maintain a thin
layer around Maxima.

Unfortunately, there are some rough edges where the backward compatible
interface inherited from the maxima based symbolics clashes with the
(IMHO much cleaner) GiNaC interface. I hope this will be resolved in
time.

The interface described in the link above is a nice way to go from
symbolics expressions to algebraic objects, which can then be used in
various algorithms (e.g. Risch, Karr). The denominator() method is
designed to be compatible with this interface.



Cheers,
Burcin