4

Why are the following computed limits different (1 by Sage, 0 by Wolfram), and which (if either) is correct?

EDIT: Revised per @Bill's suggestion to increase the numerical precision in Wolfram. (I don't know how to do the same in Sage.) The Wolfram plot strongly suggests that the limit is indeed $0$, and that the issue is entirely about numerical precision.

Sage: (you can cut/paste/execute this code here) 

#in()=
f(x) = exp(-x^2/2)/sqrt(2*pi)
F(x) = (1 + erf(x/sqrt(2)))/2
num1(a,w) = (a+w)*f(a+w) - a*f(a)
num2(a,w) = f(a+w) - f(a) 
den(a,w) = F(a+w) - F(a)
V(a,w) = 1 - num1(a,w)/den(a,w) - (num2(a,w)/den(a,w))^2
assume(w>0); print(limit(V(a,w), a=oo))
plot(V(a,1),a,0,8)

#out()=
1        #<--------- computed limit = 1

enter image description here

Wolfram: (you can execute this code here)

#in()=
f[x_]:=Exp[-x^2/2]/Sqrt[2*Pi]
F[x_]:=(1 + Erf[x/Sqrt[2]])/2 
num1[a_,w_] := (a+w)*f[a+w] - a*f[a]
num2[a_,w_] := f[a+w] - f[a]  
den[a_,w_] := F[a+w] - F[a]
V[a_,w_] := 1 - num1[a,w]/den[a,w] - (num2[a,w]/den[a,w])^2 
Assuming[w>0, Limit[V[a,w], a -> Infinity]]
Plot[V[a, 10], {a, 0, 100}, WorkingPrecision -> 128] 

#out()=
0        (* <--------- computed limit = 0 *)

enter image description here

(This is supposed to compute the limit, as a -> oo, of the variance of a standard normal distribution when truncated to the interval (a,a+w).)

r.e.s.
  • 227
  • 3
  • 11
  • You're better off asking at ask.sagemath.org, which is where the Sage expertise is. – Iguananaut Aug 29 '18 at 12:33
  • 1
    @Iguananaut - Thanks, I've taken your advice by [cross-posting there](https://ask.sagemath.org/question/43517/conflicting-sage-vs-wolfram-evaluation-of-a-limit/). – r.e.s. Aug 29 '18 at 15:31
  • 1
    It looks like you have some extremely large and extremely small values and are seeing catastrophic cancellation. For the Mathematica plot try `Plot[V[a, 1], {a, 0, 8}, WorkingPrecision -> 128]` which gives you a better behaved plot. – Bill Aug 29 '18 at 17:51
  • @Bill - Thank's very much -- it seems clearly a numerical precision issue, and that 0 really is the limit! Now I don't know whether to delete the question? – r.e.s. Aug 29 '18 at 19:08
  • My guess would be that the incorrect symbolic result given by Sage has nothing to do with the numerical instability. It might be a bug in Sage. I'd report it to them. I would not delete. – Szabolcs Sep 02 '18 at 15:05
  • @Szabolcs that is correct, they are two separate issues. The cross-posted question answers both of them to some extent (symbolic is likely a bug in Maxima). – kcrisman Sep 11 '18 at 03:02

0 Answers0