ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Fri, 02 Oct 2020 17:42:57 +0200Problem evaluating a very tiny integralhttps://ask.sagemath.org/question/39514/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?Mon, 13 Nov 2017 15:03:58 +0100https://ask.sagemath.org/question/39514/problem-evaluating-a-very-tiny-integral/Comment by mforets for <p>Hi everyone, </p>
<p>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.</p>
<p>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$$</p>
<p>Suppose d = 1, then this is simply the integral </p>
<p>$$ \int_{255.5}^{256.5}\frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}\, dx$$</p>
<p>To evaluate this integral I wrote the following code</p>
<blockquote>
<p>T = RealDistribution('gaussian,1)
print T.cum_distribution_function(256.5)-T.cum_distribution_function(255.5)</p>
</blockquote>
<p>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.</p>
<p>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.</p>
<p>I tried to use the function</p>
<blockquote>
<p>integrate</p>
</blockquote>
<p>to deal with this problem, but something funny and apparently inexplicable happened when I was trying to use it.</p>
<p>To be precise I defined this code:</p>
<pre><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
</code></pre>
<p>And if I let it run, the result I get is</p>
<pre><code>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)
</code></pre>
<p>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</p>
<pre><code>sigma = 2
</code></pre>
<p>with </p>
<pre><code>sigma = VarianceOfvy(2)
</code></pre>
<p>and try to let the whole program run again, Sage is not able anymore to output a result.</p>
<p>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?</p>
https://ask.sagemath.org/question/39514/problem-evaluating-a-very-tiny-integral/?comment=39529#post-id-39529The [error function](http://doc.sagemath.org/html/en/reference/functions/sage/functions/error.html#sage.functions.error.Function_erf) is in Sage and can be useful.Tue, 14 Nov 2017 09:27:51 +0100https://ask.sagemath.org/question/39514/problem-evaluating-a-very-tiny-integral/?comment=39529#post-id-39529Comment by B r u n o for <p>Hi everyone, </p>
<p>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.</p>
<p>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$$</p>
<p>Suppose d = 1, then this is simply the integral </p>
<p>$$ \int_{255.5}^{256.5}\frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}\, dx$$</p>
<p>To evaluate this integral I wrote the following code</p>
<blockquote>
<p>T = RealDistribution('gaussian,1)
print T.cum_distribution_function(256.5)-T.cum_distribution_function(255.5)</p>
</blockquote>
<p>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.</p>
<p>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.</p>
<p>I tried to use the function</p>
<blockquote>
<p>integrate</p>
</blockquote>
<p>to deal with this problem, but something funny and apparently inexplicable happened when I was trying to use it.</p>
<p>To be precise I defined this code:</p>
<pre><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
</code></pre>
<p>And if I let it run, the result I get is</p>
<pre><code>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)
</code></pre>
<p>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</p>
<pre><code>sigma = 2
</code></pre>
<p>with </p>
<pre><code>sigma = VarianceOfvy(2)
</code></pre>
<p>and try to let the whole program run again, Sage is not able anymore to output a result.</p>
<p>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?</p>
https://ask.sagemath.org/question/39514/problem-evaluating-a-very-tiny-integral/?comment=39527#post-id-39527You can enclose the LaTeX code in `$` (or `$$`) to make it readable.Tue, 14 Nov 2017 08:09:53 +0100https://ask.sagemath.org/question/39514/problem-evaluating-a-very-tiny-integral/?comment=39527#post-id-39527Comment by Alessandro Verzasconi for <p>Hi everyone, </p>
<p>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.</p>
<p>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$$</p>
<p>Suppose d = 1, then this is simply the integral </p>
<p>$$ \int_{255.5}^{256.5}\frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}\, dx$$</p>
<p>To evaluate this integral I wrote the following code</p>
<blockquote>
<p>T = RealDistribution('gaussian,1)
print T.cum_distribution_function(256.5)-T.cum_distribution_function(255.5)</p>
</blockquote>
<p>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.</p>
<p>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.</p>
<p>I tried to use the function</p>
<blockquote>
<p>integrate</p>
</blockquote>
<p>to deal with this problem, but something funny and apparently inexplicable happened when I was trying to use it.</p>
<p>To be precise I defined this code:</p>
<pre><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
</code></pre>
<p>And if I let it run, the result I get is</p>
<pre><code>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)
</code></pre>
<p>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</p>
<pre><code>sigma = 2
</code></pre>
<p>with </p>
<pre><code>sigma = VarianceOfvy(2)
</code></pre>
<p>and try to let the whole program run again, Sage is not able anymore to output a result.</p>
<p>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?</p>
https://ask.sagemath.org/question/39514/problem-evaluating-a-very-tiny-integral/?comment=39541#post-id-39541@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?Tue, 14 Nov 2017 16:31:55 +0100https://ask.sagemath.org/question/39514/problem-evaluating-a-very-tiny-integral/?comment=39541#post-id-39541Comment by mforets for <p>Hi everyone, </p>
<p>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.</p>
<p>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$$</p>
<p>Suppose d = 1, then this is simply the integral </p>
<p>$$ \int_{255.5}^{256.5}\frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}\, dx$$</p>
<p>To evaluate this integral I wrote the following code</p>
<blockquote>
<p>T = RealDistribution('gaussian,1)
print T.cum_distribution_function(256.5)-T.cum_distribution_function(255.5)</p>
</blockquote>
<p>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.</p>
<p>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.</p>
<p>I tried to use the function</p>
<blockquote>
<p>integrate</p>
</blockquote>
<p>to deal with this problem, but something funny and apparently inexplicable happened when I was trying to use it.</p>
<p>To be precise I defined this code:</p>
<pre><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
</code></pre>
<p>And if I let it run, the result I get is</p>
<pre><code>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)
</code></pre>
<p>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</p>
<pre><code>sigma = 2
</code></pre>
<p>with </p>
<pre><code>sigma = VarianceOfvy(2)
</code></pre>
<p>and try to let the whole program run again, Sage is not able anymore to output a result.</p>
<p>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?</p>
https://ask.sagemath.org/question/39514/problem-evaluating-a-very-tiny-integral/?comment=39545#post-id-39545Ok, 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
TrueTue, 14 Nov 2017 18:15:21 +0100https://ask.sagemath.org/question/39514/problem-evaluating-a-very-tiny-integral/?comment=39545#post-id-39545Answer by B r u n o for <p>Hi everyone, </p>
<p>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.</p>
<p>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$$</p>
<p>Suppose d = 1, then this is simply the integral </p>
<p>$$ \int_{255.5}^{256.5}\frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}\, dx$$</p>
<p>To evaluate this integral I wrote the following code</p>
<blockquote>
<p>T = RealDistribution('gaussian,1)
print T.cum_distribution_function(256.5)-T.cum_distribution_function(255.5)</p>
</blockquote>
<p>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.</p>
<p>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.</p>
<p>I tried to use the function</p>
<blockquote>
<p>integrate</p>
</blockquote>
<p>to deal with this problem, but something funny and apparently inexplicable happened when I was trying to use it.</p>
<p>To be precise I defined this code:</p>
<pre><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
</code></pre>
<p>And if I let it run, the result I get is</p>
<pre><code>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)
</code></pre>
<p>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</p>
<pre><code>sigma = 2
</code></pre>
<p>with </p>
<pre><code>sigma = VarianceOfvy(2)
</code></pre>
<p>and try to let the whole program run again, Sage is not able anymore to output a result.</p>
<p>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?</p>
https://ask.sagemath.org/question/39514/problem-evaluating-a-very-tiny-integral/?answer=39547#post-id-39547The 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](http://arblib.org/) 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](http://nemocas.org/), written for the language [Julia](https://julialang.org/). 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
Tue, 14 Nov 2017 18:42:29 +0100https://ask.sagemath.org/question/39514/problem-evaluating-a-very-tiny-integral/?answer=39547#post-id-39547Comment by mforets for <p>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 <code>erf</code> is the <em>best</em> (in some sense) exact answer one can give, since the error function is not elementary. </p>
<p>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.</p>
<pre><code>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
</code></pre>
<hr>
<p><strong>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.</strong></p>
<p>Actually, the library <a href="http://arblib.org/">Arb</a> 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 <a href="http://nemocas.org/">Nemo</a>, written for the language <a href="https://julialang.org/">Julia</a>. 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). </p>
<p>The computation I perform is a simplified version of your need:</p>
<p>$$\sum_{c=1}^{d}\int_{256c-0.5}^{256c+0.5}\frac{2}{\sqrt{\pi}}e^{-x^2}dx$$</p>
<p>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$!</p>
<pre><code>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
</code></pre>
https://ask.sagemath.org/question/39514/problem-evaluating-a-very-tiny-integral/?comment=39562#post-id-39562ok, i see. since Arb is already in C there could be some Cython wrapper for this library that can be used from Sage.Wed, 15 Nov 2017 11:26:00 +0100https://ask.sagemath.org/question/39514/problem-evaluating-a-very-tiny-integral/?comment=39562#post-id-39562Comment by mforets for <p>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 <code>erf</code> is the <em>best</em> (in some sense) exact answer one can give, since the error function is not elementary. </p>
<p>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.</p>
<pre><code>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
</code></pre>
<hr>
<p><strong>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.</strong></p>
<p>Actually, the library <a href="http://arblib.org/">Arb</a> 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 <a href="http://nemocas.org/">Nemo</a>, written for the language <a href="https://julialang.org/">Julia</a>. 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). </p>
<p>The computation I perform is a simplified version of your need:</p>
<p>$$\sum_{c=1}^{d}\int_{256c-0.5}^{256c+0.5}\frac{2}{\sqrt{\pi}}e^{-x^2}dx$$</p>
<p>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$!</p>
<pre><code>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
</code></pre>
https://ask.sagemath.org/question/39514/problem-evaluating-a-very-tiny-integral/?comment=39552#post-id-39552thanks for this interesting answer :) have you tried any of julia <--> sage interaction? i attempted [PyCall](https://github.com/JuliaPy/PyCall.jl) but i was getting lots of errors -- that was some time ago though.. for the other way round there exists [pyjulia](https://github.com/JuliaPy/pyjulia) which i haven't tried.Wed, 15 Nov 2017 00:01:06 +0100https://ask.sagemath.org/question/39514/problem-evaluating-a-very-tiny-integral/?comment=39552#post-id-39552Comment by B r u n o for <p>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 <code>erf</code> is the <em>best</em> (in some sense) exact answer one can give, since the error function is not elementary. </p>
<p>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.</p>
<pre><code>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
</code></pre>
<hr>
<p><strong>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.</strong></p>
<p>Actually, the library <a href="http://arblib.org/">Arb</a> 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 <a href="http://nemocas.org/">Nemo</a>, written for the language <a href="https://julialang.org/">Julia</a>. 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). </p>
<p>The computation I perform is a simplified version of your need:</p>
<p>$$\sum_{c=1}^{d}\int_{256c-0.5}^{256c+0.5}\frac{2}{\sqrt{\pi}}e^{-x^2}dx$$</p>
<p>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$!</p>
<pre><code>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
</code></pre>
https://ask.sagemath.org/question/39514/problem-evaluating-a-very-tiny-integral/?comment=39559#post-id-39559I 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.Wed, 15 Nov 2017 09:44:46 +0100https://ask.sagemath.org/question/39514/problem-evaluating-a-very-tiny-integral/?comment=39559#post-id-39559Comment by Alessandro Verzasconi for <p>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 <code>erf</code> is the <em>best</em> (in some sense) exact answer one can give, since the error function is not elementary. </p>
<p>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.</p>
<pre><code>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
</code></pre>
<hr>
<p><strong>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.</strong></p>
<p>Actually, the library <a href="http://arblib.org/">Arb</a> 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 <a href="http://nemocas.org/">Nemo</a>, written for the language <a href="https://julialang.org/">Julia</a>. 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). </p>
<p>The computation I perform is a simplified version of your need:</p>
<p>$$\sum_{c=1}^{d}\int_{256c-0.5}^{256c+0.5}\frac{2}{\sqrt{\pi}}e^{-x^2}dx$$</p>
<p>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$!</p>
<pre><code>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
</code></pre>
https://ask.sagemath.org/question/39514/problem-evaluating-a-very-tiny-integral/?comment=39567#post-id-39567Thank 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.Wed, 15 Nov 2017 14:10:05 +0100https://ask.sagemath.org/question/39514/problem-evaluating-a-very-tiny-integral/?comment=39567#post-id-39567Answer by dan_fulea for <p>Hi everyone, </p>
<p>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.</p>
<p>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$$</p>
<p>Suppose d = 1, then this is simply the integral </p>
<p>$$ \int_{255.5}^{256.5}\frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}\, dx$$</p>
<p>To evaluate this integral I wrote the following code</p>
<blockquote>
<p>T = RealDistribution('gaussian,1)
print T.cum_distribution_function(256.5)-T.cum_distribution_function(255.5)</p>
</blockquote>
<p>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.</p>
<p>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.</p>
<p>I tried to use the function</p>
<blockquote>
<p>integrate</p>
</blockquote>
<p>to deal with this problem, but something funny and apparently inexplicable happened when I was trying to use it.</p>
<p>To be precise I defined this code:</p>
<pre><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
</code></pre>
<p>And if I let it run, the result I get is</p>
<pre><code>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)
</code></pre>
<p>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</p>
<pre><code>sigma = 2
</code></pre>
<p>with </p>
<pre><code>sigma = VarianceOfvy(2)
</code></pre>
<p>and try to let the whole program run again, Sage is not able anymore to output a result.</p>
<p>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?</p>
https://ask.sagemath.org/question/39514/problem-evaluating-a-very-tiny-integral/?answer=39563#post-id-39563
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
Wed, 15 Nov 2017 13:37:08 +0100https://ask.sagemath.org/question/39514/problem-evaluating-a-very-tiny-integral/?answer=39563#post-id-39563Comment by Alessandro Verzasconi for <p>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:</p>
<pre><code>? intnum( x= 255.5, 256.5, exp( -x^2/2 ) / sqrt(2*Pi) )
%18 = 5.8524320023056441353545578190269084186749602193913 E-14179
</code></pre>
<p>Same can be used in sage (as a poor man solution), e.g. as</p>
<pre><code>sage: pari( "intnum( x=255.5, 256.5, exp( -x^2/2 ) / sqrt(2*Pi) )" )
5.85243200230564 E-14179
</code></pre>
<p>or</p>
<pre><code>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
</code></pre>
<p>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\ .$$</p>
<p>The better solution is then:</p>
<pre><code>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()
</code></pre>
<p>This gives:</p>
<pre><code>5.85243200228551e-14179
</code></pre>
<p>For more...</p>
<pre><code>sage: print MyPhi( a,b ).n( digits=100 )
5.852432002305644135354557819026908418674960218659348608849721844978854485293380389274750960850111602e-14179
</code></pre>
https://ask.sagemath.org/question/39514/problem-evaluating-a-very-tiny-integral/?comment=39568#post-id-39568Thank you for your very well written answer!Wed, 15 Nov 2017 14:10:21 +0100https://ask.sagemath.org/question/39514/problem-evaluating-a-very-tiny-integral/?comment=39568#post-id-39568Answer by slelievre for <p>Hi everyone, </p>
<p>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.</p>
<p>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$$</p>
<p>Suppose d = 1, then this is simply the integral </p>
<p>$$ \int_{255.5}^{256.5}\frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}\, dx$$</p>
<p>To evaluate this integral I wrote the following code</p>
<blockquote>
<p>T = RealDistribution('gaussian,1)
print T.cum_distribution_function(256.5)-T.cum_distribution_function(255.5)</p>
</blockquote>
<p>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.</p>
<p>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.</p>
<p>I tried to use the function</p>
<blockquote>
<p>integrate</p>
</blockquote>
<p>to deal with this problem, but something funny and apparently inexplicable happened when I was trying to use it.</p>
<p>To be precise I defined this code:</p>
<pre><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
</code></pre>
<p>And if I let it run, the result I get is</p>
<pre><code>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)
</code></pre>
<p>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</p>
<pre><code>sigma = 2
</code></pre>
<p>with </p>
<pre><code>sigma = VarianceOfvy(2)
</code></pre>
<p>and try to let the whole program run again, Sage is not able anymore to output a result.</p>
<p>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?</p>
https://ask.sagemath.org/question/39514/problem-evaluating-a-very-tiny-integral/?answer=53706#post-id-53706This is to complement the excellent answers
by @dan_fulea and @B_r_u_n_o.
Sage can call Arb via `ComplexBallField` and `RealBallField`.
This allows to get a certified ball around the result.
Taking inspiration from @dan_fulea's answer, we use `erfc`.
Function taking two endpoints and optionally the
precision to use, and returning the desired integral:
def my_phi(a, b, nbits=64):
B = ComplexBallField(nbits)
r = ~AA(2).sqrt()
return ((B(a)*r).erfc() - (B(b)*r).erfc())/2
Example with the default precision:
sage: my_phi(255.5, 266.5)
[5.852432002306e-14179 +/- 3.98e-14192]
Example with increased precision:
sage: my_phi(255.5, 266.5, nbits=128)
[5.85243200230564413535455781902691e-14179 +/- 4.65e-14212]
-----
Note about `erfc` and `RealBallField`.
So far (currently Sage 9.2.beta14), `RealBallField` elements
have no `erfc` method, so we can't use that.
def my_phi(a, b, nbits=64):
B = RealBallField(nbits)
r = ~AA(2).sqrt()
return ((B(a)*r).erfc() - (B(b)*r).erfc())/2
Example:
sage: my_phi(255.5, 266.5)
Traceback (most recent call last)
...
AttributeError: 'sage.rings.real_arb.RealBall' object has no attribute 'erfc'
Fri, 02 Oct 2020 17:42:57 +0200https://ask.sagemath.org/question/39514/problem-evaluating-a-very-tiny-integral/?answer=53706#post-id-53706