Ask Your Question
1

Problem evaluating a very tiny integral

asked 2017-11-13 08:03:58 -0600

updated 2017-11-14 03:53:34 -0600

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 flag offensive close merge delete

Comments

1

You can enclose the LaTeX code in $ (or $$) to make it readable.

B r u n o gravatar imageB r u n o ( 2017-11-14 01:09:53 -0600 )edit

The error function is in Sage and can be useful.

mforets gravatar imagemforets ( 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?

Alessandro Verzasconi gravatar imageAlessandro Verzasconi ( 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
mforets gravatar imagemforets ( 2017-11-14 11:15:21 -0600 )edit

2 answers

Sort by ยป oldest newest most voted
2

answered 2017-11-15 06:37:08 -0600

dan_fulea gravatar image

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
edit flag offensive delete link more

Comments

Thank you for your very well written answer!

Alessandro Verzasconi gravatar imageAlessandro Verzasconi ( 2017-11-15 07:10:21 -0600 )edit
2

answered 2017-11-14 11:42:29 -0600

B r u n o gravatar image

updated 2017-11-14 11:59:22 -0600

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
edit flag offensive delete link more

Comments

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.

mforets gravatar imagemforets ( 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.

B r u n o gravatar imageB r u n o ( 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.

mforets gravatar imagemforets ( 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.

Alessandro Verzasconi gravatar imageAlessandro Verzasconi ( 2017-11-15 07:10:05 -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

1 follower

Stats

Asked: 2017-11-13 08:03:58 -0600

Seen: 68 times

Last updated: Nov 15