<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="http://www.r-bloggers.com/numerical-pitfalls-in-computing-variance/"></div></div><div><a data-url="http://www.r-bloggers.com/numerical-pitfalls-in-computing-variance/" data-counturl="http://www.r-bloggers.com/numerical-pitfalls-in-computing-variance/" data-count="horizontal" data-text="Numerical pitfalls in computing variance" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://bayesfactor.blogspot.com/2016/05/numerical-pitfalls-in-computing-variance.html"> BayesFactor: Software for Bayesian inference</a></strong>, and kindly contributed to <a rel="nofollow" href="http://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>One of the most common tasks in statistical computing is computation of sample variance. This would seem to be straightforward; there are a number of algebraically equivalent ways of representing the sum of squares (S), such as [ S = sum_{k=1}^n ( x_k – bar{x})^2 ] or [ S = sum_{k=1}^n x_k^2 + frac{1}{n}bar{x}^2 ] and the sample variance is simply (S/(n-1)).</p>
<p>What is straightforward algebraically, however, is sometimes not so straightforward in the floating-point arithmetic used by computers. Computers cannot represent numbers to infinite precision, and arithmetic operations can affect the precision of floating-point numbers in unexpected ways.</p>
<p><a name="more"></a></p>
<p>Consider the numbers .1 and .3/3. These two numbers are equal. However,
<pre>.1 - .3/3<br></pre>
<pre><code>## [1] 0.00000000000000001388<br></code></pre>
<pre><code><br></code></pre>
<p>is not exactly 0, as one would expect it to be (for more, see “Falling into the Floating<br>Point Trap”, Chapter 4 in the <a target="_blank" href="http://www.burns-stat.com/pages/Tutor/R_inferno.pdf">“R Inferno” by Patrick Burns</a>. Multiple ways of computing the variance that are algebraically equivalent do not necessarily yield equal answers in software such as R, and some ways are better than others.</p>
<p>In a series of posts, <a target="_blank" href="https://twitter.com/JohnDCook">John D. Cook</a> shows that the seemingly reasonable, commonly used second method above, which he calls the “sum of squares” method, <a target="_blank" href="http://www.johndcook.com/blog/2008/09/26/comparing-three-methods-of-computing-standard-deviation/">can be extremely unstable in certain circumstances</a>, even giving impossible negative values. He also discusses <a target="_blank" href="http://www.johndcook.com/blog/standard_deviation/">how to compute the sample variance in a numerically stable way using Welford’s method</a>. Both of these posts are well worth reading.</p>
<p>When I read them, I thought two things. First, I was reminded that I use [used; this was written some time ago] the “sum of squares” method in the <a target="_blank" href="https://cran.r-project.org/web/packages/BayesFactor/index.html"><span>BayesFactor</span> package</a>. Secondly, I thought I would not be affected by the problem, because I represent numbers as logarithms internally for numerical stability and ease of division. Logarithms make many things easier: very large and very small numbers become easier to work with; exponentiation becomes multiplication; and multiplication and division become addition and subtraction. The tricky part of working with logarithms is addition and subtraction. If we have two numbers, (exp(a)) and (exp(b)) represented by their logarithms (a) and (b), and we want to know the logarithm of their sum, we can make use of the identities<br>[<br>begin{eqnarray*}<br>log(exp(a) + exp(b)) = a + log(1 + exp(b – a))\<br>log(exp(a) – exp(b)) = a + log(1 – exp(b – a))<br>end{eqnarray*}<br>]<br>Now arithmetic with (a) and (b) is addition and subtraction, and we can use <a target="_blank" href="http://cran.r-project.org/web/packages/Rmpfr/vignettes/log1mexp-note.pdf">accurate floating point approximations</a> of (log(1+exp(x))) and (log(1-exp(x))).<br>
<h4></h4>
<h3></h3>
<h3></h3>
<h3>When logarithms don’t help</h3>
<p>But I wasn’t really sure whether I would be affected by the instability of the “sum of squares” method, so I decided to check. It turns out, representing numbers logarithmically doesn’t necessarily help. In order to demonstrate this easily, I created an R S4 class that eases arithmetic on logarithmically-represented values. First, we load necessary libraries/files: <pre># Install the BayesFactor and devtools packages, if you don't already have them<br><br># Load my S4 class for representing real numbers with logarithms<br># and performing arithmetic on them<br># See the code at https://gist.github.com/richarddmorey/3c77d0065983e31241bff3807482443e<br>devtools::source_gist('3c77d0065983e31241bff3807482443e')<br><br># set random seed so results are reproducible<br>set.seed(2)<br></pre>
<p><a target="_blank" href="https://gist.github.com/richarddmorey/3c77d0065983e31241bff3807482443e">[Click here to view the R file you’ll be sourcing above]</a></p>
<div><span>To see the S4 class in action, we need to generate some numbers that are logarithmically-represented. The variables </span><span>x</span><span> and </span><span>y</span> below are equal to (exp(1)=2.718) and (exp(2)=7.389)<span>, respectively. The </span><span><span>modulo</span></span><span> argument gives the log-represented magnitude of the number, and the </span><span><span>sign</span></span><span> argument gives the sign (with <span>1L</span> meaning the integer representation of 1):</span></div>
<pre>x = logRepresentedReal(modulo = 1, sign = 1L)<br>y = logRepresentedReal(modulo = 2, sign = 1L)<br></pre>
<p>We can add the two numbers together, for instance: <pre>x + y<br></pre>
<pre><code>##  10.10734<br></code></pre>
<p>Although the result does not look logarithmically-represented, we can verify that it is using the <span>str</span> function: <pre>str( x + y )<br></pre>
<pre><code>## Formal class 'logRepresentedReal' [package &quot;.GlobalEnv&quot;] with 2 slots<br>##   ..@ modulo: num 2.31<br>##   ..@ sign  : int 1<br></code></pre>
<p>The result is of class <span>logRepresentedReal</span>, and the <span>modulo</span> slot tells us (log(x+y)). With the arithmetic on the logarithmically-represented numbers defined using the <span>logRepresentedReal</span> class, we can test whether our logarithms help stabilize the estimate of the variance. Following Cook, we will sample values from a uniform distribution, making use of the fact that if (z) has a uniform distribution, then (−log(z)) has an standard exponential distribution: <pre>runif2 = function(n){<br>  # Sample log values from exponential distribution<br>  x = -rexp(n)<br>  # represent all values logarithmically in a list<br>  lapply(x, logRepresentedReal, sign = 1L)<br>}<br><br>n = 100<br>z = runif2(n)<br></pre>
<p>We sampled (n=100) values from a uniform distribution. We can now compute the variance in several ways. The first way is to use the “sum of squares” method on the exponentiated values: <pre># Sum of squares method<br>var.sumsq.exp = function(z)<br>{<br>  n = length(z)<br>  z = sapply(z, as.numeric)<br>  (sum(z^2) - n*mean(z)^2)/(n-1)<br>}<br><br>var.sumsq.exp(z)<br></pre>
<pre><code>## [1] 0.07419988<br></code></pre>
<p>This presents no problem, since our uniformly-distributed values are rather moderate. We now use Welford’s method on the logarithmically-represented values to compute the variance: <pre>var.welford &lt;- function(z){<br>  n = length(z)<br>  M = list()<br>  S = list()<br>  M[[1]] = z[[1]]<br>  S[[1]] = 0<br><br>  for(k in 2:n){<br>    M[[k]] = M[[k-1]] + ( z[[k]] - M[[k-1]] ) / k<br>    S[[k]] = S[[k-1]] + ( z[[k]] - M[[k-1]] ) * ( z[[k]] - M[[k]] )<br>  }<br>  return(S[[n]] / (n - 1))<br>}<br><br>var.welford(z)<br></pre>
<pre><code>##  0.07419988<br></code></pre>
<p>And finally, we can use the “sum of squares” method on the logarithmically-represented values: <pre>var.sumsq = function(z){<br>  n = length(z)<br>  zsqr = sapply(z, function(x) x^2)<br><br>  sumz = 0<br>  sumz2 = 0<br>  for(k in 1:n){<br>    sumz = sumz + z[[k]]<br>    sumz2 = sumz2 + zsqr[[k]] <br>  }<br>  mnz = sumz/n<br>  ssz = sumz2 - n * mnz^2 <br>  return(ssz / (n-1))  <br>}<br><br>var.sumsq(z)<br></pre>
<pre><code>##  0.07419988<br></code></pre>
<p>Again, this presents no problem, since our uniformly-distributed values are moderate. So far, we see no signs of numerical instability, but none are expected. As Cook did in his example, we can add a very large number — in this case, one billion — to all sampled values. This makes the variance quite small compared to the mean, and would be expected to make the sum of squares estimates unstable. <pre>const = 1e9<br>z = sapply(z, function(x) x + const)<br><br>var.sumsq.exp(z)<br></pre>
<pre><code>## [1] -165.4949<br></code></pre>
<pre>var.welford(z)<br></pre>
<pre><code>##  0.07419973<br></code></pre>
<pre>var.sumsq(z)<br></pre>
<pre><code>##  35886<br></code></pre>
<p>Notice that both the sum of square estimates fail; the logarithmic representation used by <span>var.sumsq</span> does not help the numeric stability. The Welford method, on the other hand, yields an accurate value.  <br>
<h3>Conclusion </h3>
<p>There are many ways of combating numerical instability. Representing numbers as logarithms is one important method. Although representing numbers as logarithms is effective at combating numerical instability from some sources, it does not necessarily help in all cases. The superiority of the Welford method of computing variance, even when numbers are represented logarithmically, shows this clearly. <div></div>
<p><div><div><div><div data-width="55" data-show-faces="false" data-send="true" data-layout="box_count" data-height="62" data-href="http://www.r-bloggers.com/numerical-pitfalls-in-computing-variance/"></div></div><div><a data-url="http://www.r-bloggers.com/numerical-pitfalls-in-computing-variance/" data-counturl="http://www.r-bloggers.com/numerical-pitfalls-in-computing-variance/" data-count="vertical" data-text="Numerical pitfalls in computing variance" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://bayesfactor.blogspot.com/2016/05/numerical-pitfalls-in-computing-variance.html"> BayesFactor: Software for Bayesian inference</a></strong>.</div>
<hr>
<a rel="nofollow" href="http://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="http://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="http://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="http://www.r-bloggers.com/?s=data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="http://www.r-bloggers.com/?s=Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="http://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="http://www.r-bloggers.com/?s=map">maps</a>, <a title="Animation in R" rel="nofollow" href="http://www.r-bloggers.com/?s=animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="http://www.r-bloggers.com/?s=RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="http://www.r-bloggers.com/?s=sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="http://www.r-bloggers.com/?s=LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="http://www.r-bloggers.com/?s=SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="http://www.r-bloggers.com/?s=eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="http://www.r-bloggers.com/?s=git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="http://www.r-bloggers.com/?s=hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="http://www.r-bloggers.com/?s=Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="http://www.r-bloggers.com/?s=PCA">PCA</a>, <a title="Time series" rel="nofollow" href="http://www.r-bloggers.com/?s=time+series">time series</a>, <a title="finance trading" rel="nofollow" href="http://www.r-bloggers.com/?s=trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/h3u9oBjH8yw"></a></div>