# Problem evaluating a very tiny integral

Hi everyone,

I'm pretty new with sage. I was forced to change from MATLAB to Sage, because I was told Sage does approximate very tiny numbers better as it can work with sqrt(2) as sqrt(2) and not as the rational number approximating it.

Approximations are very important for my problem. I need to evaluate this integral $$\sum_{c=1}^{d}\int_{\min(256c-0.5,\frac{y(y+1)}{2})}^{\min(256c+0.5,\frac{y(y+1)}{2})}\frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}\, dx$$

Suppose d = 1, then this is simply the integral

$$\int_{255.5}^{256.5}\frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}\, dx$$

To evaluate this integral I wrote the following code

T = RealDistribution('gaussian,1) print T.cum_distribution_function(256.5)-T.cum_distribution_function(255.5)

because the integral above is the same as the difference of the distribution function of a standard distributed gaussian random variable between the boundaries of the integral. However, and you can check yourselves if you don't believe it, the result I get with sage is 0.

I guess that this is due to some approximation, which sage does. Indeed the exact value of the integral is pretty close (and with pretty I mean a lot) to 0. My problem is that I need to be able to have the exact value, because the sum of integrals I'm working with and my whole work behind this integral requires me to be very careful with those little tiny numbers.

I tried to use the function

integrate

to deal with this problem, but something funny and apparently inexplicable happened when I was trying to use it.

To be precise I defined this code:

def VarianceOfvy(y):
temp1 = 0
temp2 = 0
for r in range(0,y+1):
for x in range(0,r+1):
temp1 = temp1 + (255/256)^x * 1/256 * (r-x)^2
for r in range(0,y+1):
for x in range(0,r+1):
temp2 = temp2 + ((255/256)^x * 1/256 * (r-x))^2
sigma = temp1 - temp2
return sqrt(sigma)

def Integerd(y):
b = y*(y+1)/2
d = 1
c = 0
while min((c+1)*256-0.5,b) == (c+1)*256-0.5:
d = d+1
c = c+1
return d

def Probabilityvynequiv(y):
var('c')
b = (y*(y+1))/2
sigma = 2
mu = 1
d = Integerd(y)
factor = 1/(integrate(1/sqrt(2*pi)*e^(-(x/sigma)^2/2),x,-oo,(b-mu)/sigma) - integrate(1/sqrt(2*pi)*e^(-(x)^2/2),x,-oo,(-mu)/sigma))
p = sum(factor*1/sigma*integrate(1/sqrt(2*pi)*e^((-x^2)/(2)),x,c*256+0.5,min((c+1)*256-0.5,b)),c,0,d)
return p


And if I let it run, the result I get is

1/2*(erf(255.75*sqrt(2)) - erf(128.25*sqrt(2)) + erf(127.75*sqrt(2)) - erf(0.25*sqrt(2)))/(3*erf(1/4*sqrt(2)) + 1)


which I assume is correct, and at least it tells me that Sage is able to read my code and output a result. If I call the function VarianceOfvy(2), the result I get is 3/65536*sqrt(11133895), which is also correct. Now, if I'm changing the command

sigma = 2


with

sigma = VarianceOfvy(2)


and try to let the whole program run again, Sage is not able anymore to output a result.

I'm really lost and I don't know what to do. Could someone advise me and give me some hints on how to evaluate those tiny integrals, in such a way that I don't loose any approximation?

edit retag close merge delete

1

You can enclose the LaTeX code in $ (or $$) to make it readable. ( 2017-11-14 01:09:53 -0600 )edit The error function is in Sage and can be useful. ( 2017-11-14 02:27:51 -0600 )edit @mforets I don't know what you mean. If you tell Sage to compute the difference between erf(256.5) and erf(255.5), with the code print erf(256.5)-erf(255.5), it tells you that it is 0. Is there a way to make it more precise? ( 2017-11-14 09:31:55 -0600 )edit Ok, i see. Actually with 10000 bits of precision it's still not enough: sage: abs(erf(RealField(10000)(256.5)) - erf(RealField(10000)(255.5))) > 0 False  I don't know; the exponential of -x^2 goes fast to zero, can we estimate how much precision is needed here? With one more order of magnitude of precision, i get someting non-zero (long time!): sage: abs(erf(RealField(100000)(256.5)) - erf(RealField(100000)(255.5))) > 0 True  ( 2017-11-14 11:15:21 -0600 )edit ## 2 answers Sort by » oldest newest most voted In such cases it is better to use the complementary error function. First of all note that pari/gp does a pretty good job to compute the integral$$ J = \int_{255.5}^{256.5}\frac{1}{\sqrt{2\pi}}\;\exp -\frac{x^2}2\, dx $$"as is", without contorsions: ? intnum( x= 255.5, 256.5, exp( -x^2/2 ) / sqrt(2*Pi) ) %18 = 5.8524320023056441353545578190269084186749602193913 E-14179  Same can be used in sage (as a poor man solution), e.g. as sage: pari( "intnum( x=255.5, 256.5, exp( -x^2/2 ) / sqrt(2*Pi) )" ) 5.85243200230564 E-14179  or sage: a, b = 255.5, 256.5 sage: pari( "intnum( x=%s, %s, exp( -x^2/2 ) / sqrt(2*Pi) )" % (a,b) ) 5.85243200230564 E-14179  Recall that the complementary error function is$$ \textrm{erfc}(z) = 1-\textrm{erf}(z) = \frac 2\pi\int_z^\infty\exp(-t^2)\; dt$$(change of variables, t=u/sqrt2...)$$=2\cdot \frac 1{\sqrt{2\pi}}\int_{z/\sqrt 2}^\infty\exp-\frac {u^2}2\; du\ .$$The better solution is then: def MyPhi( a,b ): return ( erfc(a/sqrt(2)) - erfc(b/sqrt(2)) ) / 2 a, b = 255.5, 256.5 print MyPhi( a,b ).n()  This gives: 5.85243200228551e-14179  For more... sage: print MyPhi( a,b ).n( digits=100 ) 5.852432002305644135354557819026908418674960218659348608849721844978854485293380389274750960850111602e-14179  more ## Comments Thank you for your very well written answer! ( 2017-11-15 07:10:21 -0600 )edit The integrals you are trying to compute are expressible in terms of the error function. Note that the expression of your integral in terms of erf is the best (in some sense) exact answer one can give, since the error function is not elementary. SageMath is able to perform the computation if you allow for enough precision (thanks mforest for testing!), but it takes a very long time. An example of use in SageMath explains why. Below, the result is around 3.1\cdot 10^{-28354}, very close to zero, so hard to distinguish from zero. sage: F = RealField(100000) sage: u = F(256.5) sage: v = F(255.5) sage: d = erf(u) - erf(v) sage: d 3.10226267827805307639865203[... removed loooots of digits! ...]95e-28354  My previous answer, as I thought SageMath was not able to perform the computation. Yet, it is still interesting since you can have the same result much faster using directly (or though Nemo) the library Arb. Actually, the library Arb of Fredrik Johansson is able to perform the computation. The simplest way (to my taste, you may disagree) to use this C library is through the CAS Nemo, written for the language Julia. I give you below an example of use in Nemo (but you can of course write an equivalent C code directly using Arb if you prefer). The computation I perform is a simplified version of your need:$$\sum_{c=1}^{d}\int_{256c-0.5}^{256c+0.5}\frac{2}{\sqrt{\pi}}e^{-x^2}dx$$The code below shows that for$d=50$, the value of the sum is something like$3.1\cdot 10^{-28354}$. The closeness to zero explains why one needs so much precision to be able to distinguish from$0\$!

julia> using Nemo

Welcome to Nemo version 0.5.0

Nemo comes with absolutely no warranty whatsoever

julia> function s(d, prec = 1000)       # using default precision 1000 bits
F = ComplexField(prec)        # For some reason, erf does not work over RealField
sum = F(0)
for c = 1:d
u = F(256*c+.5)
v = F(256*c-.5)
sum += erf(u) - erf(v)
end
return sum
end
s (generic function with 2 methods)

julia> s(50) # Not enough precision: get 0 +/- something
[+/- 8.40e-300] + i*0

julia> s(50, 10000) # Still not enough!
[+/- 4.52e-3009] + i*0

julia> s(50, 100000) # Finally!
[3.102262678278053076398652035217490023670965971936356980763793132685150297353942594965955722762805573392153330406480953855731387398023115984622786616977425540305272338247087011441401529284717218316248016664526853584800590168918124704184766482548557105353651652368296375023569553822793725443102971750875855897752226561077609876866162061711562923329377439218548528916354400384468966317197680008340200698258455463251434793207006670012945890581001148258637826042374760859045345143168726986415658022488421604607834339264038625902663546478278814980043653958176479672213608838686103342175667992448852731310595727573672549161581892009023641358983676574601083376637614888507061187244777656151985731856257941976475105038054800870251422661227024669221253196375384625938232368807030347647045465382519932395272845244305245675307763521925808514638809267924250263503220341147139844099847262900841324568603016291331669298541013247881551863448642400104609478241062200239333495268610840375558431182176503376839668955839352645002793527349614716665957443338984442304676833795005909267696213605538509355001746049287714536707616120292584911365481512840133877579846868859068920299888550789717502048145517061355005325552635169752304689718047493686303050935112466379110899771012977020553628991752738579948440664336445330854288254372713613062555003022198940042814571068769124858706877287484803032872382430528486856936500533970245684004943194205532248702178101002288809615625387993450186953306376601302771382662537102238863693773280261194027940865405478664244791943840999854120140399491991893785064480156556621092433860406215704328279792938871986918708868424499764959378290014074093197179867944667683870991035370327938743608633506239623813944138924732150370855135636489875222e-28354 +/- 5.52e-30101] + i*0

more

thanks for this interesting answer :) have you tried any of julia <--> sage interaction? i attempted PyCall but i was getting lots of errors -- that was some time ago though.. for the other way round there exists pyjulia which i haven't tried.

( 2017-11-14 17:01:06 -0600 )edit

I haven't. For the current question, we should be able to expose Arb's functionalities directly in Sage, rather than using mpmath. I guess this would be much faster.

( 2017-11-15 02:44:46 -0600 )edit

ok, i see. since Arb is already in C there could be some Cython wrapper for this library that can be used from Sage.

( 2017-11-15 04:26:00 -0600 )edit
1

Thank you for your answer, I really appreciate it! I knew the integral I posted above, was a very tiny one, it is actually a simpler version of the one I need to evaluate, because I tried to simplify the problem, in order to make it readable. I like both answers, but I accepted the other one, because it allows me to solve my problem using the converse error function, which I found being preciser than the error function itself. I would also like to thank mforets for his/her (?) hints.

( 2017-11-15 07:10:05 -0600 )edit