<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/naive-numerical-sums-in-r/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://ntguardian.wordpress.com/2018/09/10/naive-numerical-sums-in-r/"> R – Curtis Miller's Personal Website</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<h2>Introduction</h2>
<p>The <a rel="nofollow" href="https://en.wikipedia.org/wiki/Kolmogorov%E2%80%93Smirnov_test#Kolmogorov_distribution" target="_blank">Kolmogorov distribution</a> (which I call <img alt="F(x)" title="F(x)" src="https://s0.wp.com/latex.php?latex=F%28x%29&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">) is as follows:</p>
<p><img alt="F(x) = frac{sqrt{2 pi}}{x} sum_{k = 1}^{infty} e^{-(2k - 1)^2 pi^2/(8x^2)}" title="F(x) = frac{sqrt{2 pi}}{x} sum_{k = 1}^{infty} e^{-(2k - 1)^2 pi^2/(8x^2)}" src="https://s0.wp.com/latex.php?latex=F%28x%29+%3D+%5Cfrac%7B%5Csqrt%7B2+%5Cpi%7D%7D%7Bx%7D+%5Csum_%7Bk+%3D+1%7D%5E%7B%5Cinfty%7D+e%5E%7B-%282k+-+1%29%5E2+%5Cpi%5E2%2F%288x%5E2%29%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"></p>
<p><span></span></p>
<p>There is no known simpler form and we have to work with this sum as it is. This is an infinite sum. How can we compute the value of this infinite sum numerically?</p>
<p>Naïvely we can do the following:</p>
<pre>
summand &lt;- function(x, k) sqrt(2 * pi)/x * exp(-(2 * k - 1)^2 * pi^2/(8 * x^2))
# Compute F(1)
sum(summand(1, 1:500))
</pre>
<pre>
[1] 0.7300003
</pre>
<p>In other words, sum up many of the terms and you should be close to the actual infinite sum.</p>
<p>This is a crude approach. The answer is not wrong (numerically) but certainly we should understand why adding up that many terms works. Also, we could have added more terms than necessary… or not enough.</p>
<p>So how can we compute this sum that guarantees some level of precision while at the same time not adding any more terms than necessary? Unfortunately I don’t recall how to do this from my numerical methods classes, but I believe I have found an approach that works well enough for my purposes.</p>
<h2>Geometric Sums</h2>
<p>An infinite sum is defined as the limit of a sequence of finite sums. Let <img alt="S_n = sum_{k = 0}^n a_k" title="S_n = sum_{k = 0}^n a_k" src="https://s0.wp.com/latex.php?latex=S_n+%3D+%5Csum_%7Bk+%3D+0%7D%5En+a_k&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">, where <img alt="a_k" title="a_k" src="https://s0.wp.com/latex.php?latex=a_k&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> is some sequence (we can have <img alt="a_k = frac{1}{k}" title="a_k = frac{1}{k}" src="https://s0.wp.com/latex.php?latex=a_k+%3D+%5Cfrac%7B1%7D%7Bk%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> or <img alt="a_k = 2^{-k}" title="a_k = 2^{-k}" src="https://s0.wp.com/latex.php?latex=a_k+%3D+2%5E%7B-k%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">, for example; in the case of the Kolmogorov distribution, we had <img alt="a_k = frac{sqrt{2 pi}}{x} e^{-(2k - 1)^2 pi^2/(8x^2)}" title="a_k = frac{sqrt{2 pi}}{x} e^{-(2k - 1)^2 pi^2/(8x^2)}" src="https://s0.wp.com/latex.php?latex=a_k+%3D+%5Cfrac%7B%5Csqrt%7B2+%5Cpi%7D%7D%7Bx%7D+e%5E%7B-%282k+-+1%29%5E2+%5Cpi%5E2%2F%288x%5E2%29%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">). Then, by definition, <img alt="S_{infty} = lim_{n to infty} S_n" title="S_{infty} = lim_{n to infty} S_n" src="https://s0.wp.com/latex.php?latex=S_%7B%5Cinfty%7D+%3D+%5Clim_%7Bn+%5Cto+%5Cinfty%7D+S_n&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">.</p>
<p>My attempt to compute this sum simply amounts to trying to find an <img alt="n" title="n" src="https://s0.wp.com/latex.php?latex=n&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> such that the difference between <img alt="S_n" title="S_n" src="https://s0.wp.com/latex.php?latex=S_n&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> and <img alt="S_{infty}" title="S_{infty}" src="https://s0.wp.com/latex.php?latex=S_%7B%5Cinfty%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> is no larger than machine precision: that is, I want <img alt="left| S_n - S_infty right| leq epsilon" title="left| S_n - S_infty right| leq epsilon" src="https://s0.wp.com/latex.php?latex=%5Cleft%7C+S_n+-+S_%5Cinfty+%5Cright%7C+%5Cleq+%5Cepsilon&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> where <img alt="epsilon" title="epsilon" src="https://s0.wp.com/latex.php?latex=%5Cepsilon&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> is the machine precision of the computer.</p>
<p>Since we don’t know what <img alt="S_{infty}" title="S_{infty}" src="https://s0.wp.com/latex.php?latex=S_%7B%5Cinfty%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> is, we can instead decide that machine convergence occurs when <img alt="left| S_n - S_{n + 1}right| leq epsilon" title="left| S_n - S_{n + 1}right| leq epsilon" src="https://s0.wp.com/latex.php?latex=%5Cleft%7C+S_n+-+S_%7Bn+%2B+1%7D%5Cright%7C+%5Cleq+%5Cepsilon&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">; that is, when one summand and the next summand are numerically indistinguishable. Since <img alt="S_{n + 1} - S_n = a_{n + 1}" title="S_{n + 1} - S_n = a_{n + 1}" src="https://s0.wp.com/latex.php?latex=S_%7Bn+%2B+1%7D+-+S_n+%3D+a_%7Bn+%2B+1%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">, this criterion is the same as requiring that <img alt="left| a_{n + 1}right| leq epsilon" title="left| a_{n + 1}right| leq epsilon" src="https://s0.wp.com/latex.php?latex=%5Cleft%7C+a_%7Bn+%2B+1%7D%5Cright%7C+%5Cleq+%5Cepsilon&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">.</p>
<p>Every sum that converges requires the condition <img alt="a_k to 0" title="a_k to 0" src="https://s0.wp.com/latex.php?latex=a_k+%5Cto+0&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">, so this criterion always yields an <img alt="n" title="n" src="https://s0.wp.com/latex.php?latex=n&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> that gives “numerical convergence”. Of course, any Calculus II student who was paying attention in their class can tell you that not all infinite sums with summands going to zero converges, with the classic counterexample being the <a rel="nofollow" href="https://en.wikipedia.org/wiki/Harmonic_series_(mathematics)" target="_blank">Harmonic series</a>. So this approach would claim that sums that diverge are numerically convergent, which is bad. We cannot even expect this method to work in cases where the sum does converge, but it does so slowly (see a later example). However, in some cases this approach may be okay.</p>
<p>Take the case of a geometric series:</p>
<p><img alt="sum_{k = 0}^{infty} rho^k" title="sum_{k = 0}^{infty} rho^k" src="https://s0.wp.com/latex.php?latex=%5Csum_%7Bk+%3D+0%7D%5E%7B%5Cinfty%7D+%5Crho%5Ek&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"></p>
<p>with <img alt="left| rho right|" title="left| rho right|" src="https://s0.wp.com/latex.php?latex=%5Cleft%7C+%5Crho+%5Cright%7C&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> less than 1. These sums converge; in fact, mathematicians consider them as converging quickly. We also have a formula for what the sum is:</p>
<p><img alt="S_{infty} = frac{1}{1 - rho}" title="S_{infty} = frac{1}{1 - rho}" src="https://s0.wp.com/latex.php?latex=S_%7B%5Cinfty%7D+%3D+%5Cfrac%7B1%7D%7B1+-+%5Crho%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"></p>
<p>After some algebra, we can quickly find a rule for determining how many summands we need to attain “numerical convergence”:</p>
<p><img alt="n geq frac{log epsilon}{log left| rho right|} + 2" title="n geq frac{log epsilon}{log left| rho right|} + 2" src="https://s0.wp.com/latex.php?latex=n+%5Cgeq+%5Cfrac%7B%5Clog+%5Cepsilon%7D%7B%5Clog+%5Cleft%7C+%5Crho+%5Cright%7C%7D+%2B+2&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"></p>
<p>We can see that in action with some R examples:</p>
<pre>
.Machine$double.eps  # Numerical accuracy of this system
</pre>
<pre>
[1] 2.220446e-16
</pre>
<pre>
log(.Machine$double.eps)/log(0.5)
</pre>
<pre>
[1] 52
</pre>
<pre>
sum(0.5^(0:53))
</pre>
<pre>
[1] 2
</pre>
<pre>
# 2 is the correct answer, but the interpreter rounds its output; is the answer
# actually 2?
sum(0.5^(0:53)) - 2 == 0
</pre>
<pre>
[1] TRUE
</pre>
<pre>
sum(0.5^(0:52)) - 2 == 0
</pre>
<pre>
[1] FALSE
</pre>
<h2>Caveats</h2>
<p>This method, though, should be met with suspicion. For instance, it will not work for a slowly convergent sum. Take for example <img alt="sum_{k = 1}^{infty} frac{1}{k^2} = frac{pi^2}{6}" title="sum_{k = 1}^{infty} frac{1}{k^2} = frac{pi^2}{6}" src="https://s0.wp.com/latex.php?latex=%5Csum_%7Bk+%3D+1%7D%5E%7B%5Cinfty%7D+%5Cfrac%7B1%7D%7Bk%5E2%7D+%3D+%5Cfrac%7B%5Cpi%5E2%7D%7B6%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">. If you apply the above technique, then “numerical convergence” is achieved for <img alt="n &gt; epsilon^{-1/2}" title="n &gt; epsilon^{-1/2}" src="https://s0.wp.com/latex.php?latex=n+%3E+%5Cepsilon%5E%7B-1%2F2%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">. Not only is that a very large number, it won’t achieve our goal of good numerical accuracy.</p>
<pre>
.Machine$double.eps^(-1/2)
</pre>
<pre>
[1] 67108864
</pre>
<pre>
N &lt;- 67108865
sum((1:N)^(-2))  # This may take a while
</pre>
<pre>
[1] 1.644934
</pre>
<pre>
sum((1:N)^(-2)) - pi^2/6
</pre>
<pre>
[1] -1.490116e-08
</pre>
<p>That difference is much larger than numerical accuracy.</p>
<p>In fact, the technique doesn’t always work for geometric sums either, as demonstrated by these examples.<sup><a rel="nofollow" href="https://ntguardian.wordpress.com/2018/09/10/naive-numerical-sums-in-r/#fn-3431-1" target="_blank">1</a></sup></p>
<pre>
sum(.99^(0:4000)) - 100
</pre>
<pre>
[1] -8.526513e-14
</pre>
<pre>
sum(.999^(0:40000)) - 1000
</pre>
<pre>
[1] -9.094947e-13
</pre>
<h2>Conclusion</h2>
<p>However, while this method cannot guarantee quick convergence or even convergence, I think it’s good enough for the sum I want to compute.</p>
<p>First, the sum converges more quickly than a geometric sum, as the summand decrease at a rate of <img alt="O(rho^{n^2})" title="O(rho^{n^2})" src="https://s0.wp.com/latex.php?latex=O%28%5Crho%5E%7Bn%5E2%7D%29&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> rather than <img alt="O(rho^n)" title="O(rho^n)" src="https://s0.wp.com/latex.php?latex=O%28%5Crho%5En%29&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">. Second, a method trying to attain numerical accuracy would need to be programmed, and if its implementation is written in R, that implementation will likely be much slower than simply using <code>sum()</code>, since the latter is implemented using fast C code. Such an implementation would have to be written from scratch in C++ using a tool such as <strong>Rcpp</strong>. One must wonder whether the tiny numerical efficiency and speed one might potentially gain are worth the work; if <code>x</code> is large it may be best to just round off the CDF at 1.</p>
<p>In the end, using the lessons learned above, I implemented the Kolmogorov distribution in the package I’m writing for my current research project with the code below.</p>
<pre>
pkolmogorov &lt;- function(q, summands = ceiling(q * sqrt(72) + 3/2)) {
  sqrt(2 * pi) * sapply(q, function(x) { if (x > 0) {
    sum(exp(-(2 * (1:summands) - 1)^2 * pi^2/(8 * x^2)))/x
  } else {
    0
  }})
}
pkolmogorov &lt;- Vectorize(pkolmogorov, &quot;q&quot;)
</pre>
<p>Numerical summation, as I mentioned above, is something I know little about, so I’d appreciate any readers with thoughts on this topic (and knowledge of how this is done in R) to share in the comments.</p>
<hr>
<p>I have created a video course published by <a rel="nofollow" href="http://packtpub.com/" target="_blank">Packt Publishing</a> entitled <a rel="nofollow" href="https://ntguardian.wordpress.com/video-courses/training-your-systems-with-python-statistical-modeling/" target="_blank"><em>Training Your Systems with Python Statistical Modelling</em></a>, the third volume in a four-volume set of video courses entitled, <em>Taming Data with Python; Excelling as a Data Analyst</em>. This course discusses how to use Python for machine learning. The course covers classical statistical methods, supervised learning including classification and regression, clustering, dimensionality reduction, and more! The course is peppered with examples demonstrating the techniques and software on real-world data and visuals to explain the concepts presented. Viewers get a hands-on experience using Python for machine learning. If you are starting out using Python for data analysis or know someone who is, please consider <a rel="nofollow" href="https://www.packtpub.com/big-data-and-business-intelligence/training-your-systems-python-statistical-modeling-video" target="_blank">buying my course</a> or at least spreading the word about it. You can buy the course directly or purchase a <a rel="nofollow" href="https://www.packtpub.com/mapt/" target="_blank">subscription to Mapt</a> and watch it there.</p>
<p>If you like my blog and would like to support it, spread the word (if not get a copy yourself)! Also, stay tuned for future courses I publish with Packt at the Video Courses section of my site.</p>
<div>
<hr>
<ol>
<li>
If you’re keeping score at home you’ll notice these sums use many more terms than the rule I describe above suggests, suggesting that the problem is not just that the “stopping” rule is wrong but that adding more terms from the sum won’t lead to greater numerical accuracy since the terms being added are essentially zero. Something else will need to be done, going beyond just adding more terms. <a rel="nofollow" href="https://ntguardian.wordpress.com/2018/09/10/naive-numerical-sums-in-r/#fnref-3431-1" target="_blank"><img alt="↩" src="https://s.w.org/images/core/emoji/2.4/72x72/21a9.png" class="wp-smiley"></a>
</li>
</ol>
</div>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://ntguardian.wordpress.com/2018/09/10/naive-numerical-sums-in-r/"> R – Curtis Miller's Personal Website</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/niC-XWo5UJQ" height="1"></a></div>