Opened 5 years ago

Closed 5 years ago

#6189 closed defect (fixed)

[with patch, positive review] 'integrate' produces incorrect answer

Reported by: emchristiansen Owned by: burcin
Priority: major Milestone: sage-4.2
Component: calculus Keywords: integrate, integral, incorrect
Cc: Merged in: sage-4.2.alpha0
Authors: Karl-Dieter Crisman Reviewers: Mike Hansen
Report Upstream: Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Description

Affects: x64 Ubuntu 9.04 (Jaunty) Sage 4.0, compiled from source, and updated (sage -upgrade) as of this posting

By simply casting a number using 'n()', I can cause integrate to return a vastly different result. See below:

var = sage.calculus.calculus.var

def NormalPDF(x,mu,sigma):
    return 1/sqrt(2*pi*sigma^2)*exp(-1/(2*sigma**2)*(x-mu)^2)

x = var('x')
mu = var('m')
sigma = var('s')
assume(sigma>0)
child1 = NormalPDF(x,0,1)
child2 = NormalPDF(x,n(0),n(1))
parent = NormalPDF(x,mu,sigma)

# this expansion helps Sage to do the integral
integrand1 = ((parent-child1)^2).expand()
integrand2 = ((parent-child2)^2).expand()

int1 = integrate(integrand1,x,-infinity,infinity)
int2 = integrate(integrand2,x,-infinity,infinity)

print "THIS EXPRESSION:"
print int1
print "\nSHOULD BE VERY SIMILIAR TO THIS EXPRESSION:"
print int2

print "\nMAKING THIS NUMBER:"
print int1.subs({mu:0,sigma:1}).n()
print "\nVERY SIMILAR TO THIS NUMBER:"
print int2.subs({mu:0,sigma:1}).n()

The above produces the output:

THIS EXPRESSION:
1/2*((s + 1)*sqrt(s^2 + 1)*e^(1/2*m^2/(s^2 + 1)) - 2*sqrt(2)*s)*e^(-1/2*m^2/(s^2 + 1))/(sqrt(s^2 + 1)*sqrt(pi)*s)

SHOULD BE VERY SIMILIAR TO THIS EXPRESSION:
1/2*(sqrt(0.5*s^2 + 0.5)*e^(1/2*m^2/(s^2 + 1)) - sqrt(2)*s)*e^(-1/2*m^2/(s^2 + 1))/(sqrt(0.5*s^2 + 0.5)*sqrt(pi)*s)

MAKING THIS NUMBER:
0.000000000000000

VERY SIMILAR TO THIS NUMBER:
-0.116847488627555

Mathematica claims the correct integral is:

\frac{1+\sigma }{2 \sqrt{\pi } \sigma }-\frac{\sqrt{2} e^{-\frac{\mu ^2}{2+2 \sigma ^2}}}{\sqrt{\pi +\pi  \sigma ^2}}

Attachments (1)

trac_6189-num-approx-integral.patch (1.3 KB) - added by kcrisman 5 years ago.
Based on 4.1.2.alpha2

Download all attachments as: .zip

Change History (5)

comment:1 Changed 5 years ago by kcrisman

  • Summary changed from 'integrate' produces incorrect answer to [with patch, needs review] 'integrate' produces incorrect answer

This is now fixed, presumably in the Maxima upgrade. Note that the integral in fact computes without expand(), and in that case there is no 'experimental error'! Attached patch verifies this.

Changed 5 years ago by kcrisman

Based on 4.1.2.alpha2

comment:2 Changed 5 years ago by kcrisman

Despite what it says, actually based on 4.1.2.alpha4.

comment:3 Changed 5 years ago by mhansen

  • Authors set to Karl-Dieter Crisman
  • Merged in set to sage-4.2.alpha0
  • Reviewers set to Mike Hansen
  • Status changed from needs_review to positive_review
  • Summary changed from [with patch, needs review] 'integrate' produces incorrect answer to [with patch, positive review] 'integrate' produces incorrect answer

Looks good to me.

comment:4 Changed 5 years ago by mhansen

  • Milestone set to sage-4.2
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.