On Sat, 24 Aug 2013 14:16:01 -0700 (PDT)
Emmanuel Charpentier <[email protected]> wrote:
> 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))
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
--
You received this message because you are subscribed to the Google Groups
"sage-support" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To post to this group, send email to [email protected].
Visit this group at http://groups.google.com/group/sage-support.
For more options, visit https://groups.google.com/groups/opt_out.
