Ask Your Question
0

stack overflow in computing entropy

asked 2012-09-22 13:42:12 -0600

kayhan gravatar image

Hi,

I am still new to sage, I am trying to compute use sage to do some symbolic computation. Here is my test to compute entropy of the Gaussian distribution:

x, sigma, mu,a = var('x','sigma','mu','a')
q = 1/((2*pi)^(1/2)*sigma)*exp(-(x - mu)^2/(2*sigma^2))
assume(sigma>0)
assume(mu>0)
show(q)
import sympy
integral(q,(x,-oo,oo)) 
integral(-q*log(q),(x,-oo,oo))

First of all, it complains if I don't assume mu is positive. It is right that sigma should be positive but why mu ?! Second, even after assuming mu>0, it produces stack overflow:

WARNING: Output truncated!  
full_output.txt


;;;
;;; Binding stack overflow.
;;; Jumping to the outermost toplevel prompt
;;;


;;;
;;; Binding stack overflow.
;;; Jumping to the outermost toplevel prompt
;;;

Does anyone know why? I am sure I am missing something.

edit retag flag offensive close merge delete

Comments

The semicolons indicate this is a Maxima overflow. (Maxima handles our symbolic integration, and raises the error messages asking for the sign of various parameters.)

kcrisman gravatar imagekcrisman ( 2012-09-22 16:05:08 -0600 )edit

By the way, importing sympy would not do anything here, unless you used sympy.integral or whatever their syntax is. Which could be an interesting option for you (or so I thought until I saw your other question).

kcrisman gravatar imagekcrisman ( 2012-09-22 16:05:32 -0600 )edit

3 answers

Sort by ยป oldest newest most voted
1

answered 2012-09-22 19:32:57 -0600

achrzesz gravatar image

updated 2012-09-22 22:40:49 -0600

Maxima can compute the entropy

sage:maxima_console()

(%i1) display2d:false$                                       

(%i2) q:1/((2*%pi)^(1/2)*sigma)*exp(-(x - mu)^2/(2*sigma^2))$

(%i3) assume(sigma>0)$

(%i4) assume(mu>0)$

(%i5) ii:integrate(-q*log(q),x,minf,inf)$

(%i6) changevar(ii,(x-mu)/(sqrt(2)*sigma)-y,y,x)$

(%i7) ev(%,integrate)$

(%i8) ratsimp(%);

(%o8) (2*log(sigma)+log(%pi)+log(2)+1)/2
(%i9) logcontract(%);

(%o9) (log(2*%pi*sigma^2)+1)/2

and the last expression is equal to:

log(sqrt(2*pi)*sigma)+1/2

edit flag offensive delete link more

Comments

Nice, except that one has to change the variable explicitly. So Maxima can't really do it without extra intervention :(

kcrisman gravatar imagekcrisman ( 2012-09-24 08:01:04 -0600 )edit
0

answered 2012-09-22 16:13:33 -0600

kcrisman gravatar image

updated 2012-09-24 08:00:33 -0600

This isn't an answer, but apparently the problem isn't Maxima per se. (Recall that Maxima is the normal integration engine within Sage.) After the requisite assumptions, I have (in Sage's native Maxima):

sage: maxima_console()
<snip>
(%i1) q:1/((2*%pi)^(1/2)*sigma)*exp(-(x - mu)^2/(2*sigma^2));
                                             2
                                     (x - mu)
                                   - ---------
                                            2
                                     2 sigma
                                 %e
(%o1)                       -----------------------
                            sqrt(2) sqrt(%pi) sigma
(%i2) assume(sigma>0);
(%o2)                             [sigma > 0]
(%i3) assume(mu>0);
(%o3)                              [mu > 0]
(%i4) integrate(q,x,minf,inf);
(%o4)                                  1

(%i5) integrate(-q*log(q),x,minf,inf);
                                                        2
                                 2              (x - mu)
                         (x - mu)             - ---------
                inf    - ---------                     2
               /                2               2 sigma
               [         2 sigma            %e
               I     %e            log(-----------------------) dx
               ]                       sqrt(2) sqrt(%pi) sigma
               /
                minf
(%o5)        - ---------------------------------------------------
                             sqrt(2) sqrt(%pi) sigma
(%i6)

which shouldn't cause an overflow. See the other answer for more details on how to simplify this.

(In the notebook, one can choose "Maxima" from the combo box of systems, or put %maxima at the beginning of a cell, and then just put the input commands above in that cell to evaluate this.)

edit flag offensive delete link more

Comments

where did you define "inf" and "minf" ? did you import a library to get the answer?

kayhan gravatar imagekayhan ( 2012-09-22 16:24:01 -0600 )edit

Sorry, I'll edit to make this clear.

kcrisman gravatar imagekcrisman ( 2012-09-22 16:25:02 -0600 )edit

sorry for my ignorance, but it is still unclear for me. What is maxima_console() ? I am using web-based interface, are you using something else? Do you think it is a bug? You mentioned default mode of integration in sage is maxima, how can I change it?

kayhan gravatar imagekayhan ( 2012-09-22 16:38:30 -0600 )edit

Hmm, that's true, you probably can't do the interactive console in Maxima. However, you could just type in the commands - I'll edit again.

kcrisman gravatar imagekcrisman ( 2012-09-24 07:59:01 -0600 )edit
0

answered 2012-09-23 15:19:05 -0600

kayhan gravatar image

Thank you for your reply.

The answer is correct but just a few points for clarification of my mind:

  • maxima is able to compute the integral correctly but it totally different software and sage can interface with it. However the syntax is totally different.

  • The mystery that why sage cannot solve the integral is not resolved, right?

  • It is not clear why we need to define (mu > 0) even in maxima.

Again thank you for your reply.

edit flag offensive delete link more

Comments

You are correct that the overflow is still not resolved. Sage doesn't just interface with Maxima, though, but USES it to do integrals. (You can also use Sympy and some other options, with the `algorithm` keyword.) Check out http://www.sagemath.org/doc/reference/sage/symbolic/integration/integral.html#sage.symbolic.integration.integral.integral

kcrisman gravatar imagekcrisman ( 2012-09-24 08:02:29 -0600 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

Stats

Asked: 2012-09-22 13:42:12 -0600

Seen: 352 times

Last updated: Sep 24 '12