var() is not equal to sd()^2 in R

Mathematically, it is evident that var = sd^2 (or var^0.5 = sd), but R seems to answer TRUE for the latter and FALSE for the former. To my knowledge, both var() and sd() function use n-1 for degree of freedom. What would be the problem?

data <- c(20, 20, 30, 20, 16, 13, 20)

var(data)^0.5 == sd(data)
var(data) == sd(data)^2

• There is no problem. This is about numerical precision. Use the all.equal function instead to acknowledge the fact that your machine represents real numbers as floating-point numbers with a specific numerical precision. The documentation of all.equal also explains .Machine\$double.eps. Sep 26 at 18:10
• Thanks! Will find the document you mentioned. Sep 26 at 18:31
• Try this: print(x <- 1 + 1e-15, 17); sqrt(x)^2 == x
– whuber
Sep 26 at 19:39
• This is R FAQ 7.31. Sep 27 at 5:43
• A good way to visualize the problem: calculate sqrt(13), write the result on a piece of paper, now square your intermediary result and see if you get 13 again. This is what R (has to) do - it's a calculator, not a symbolic equation solver, it doesn't retain knowledge that a specific intermediary result is supposed to be sqrt(13). Sep 27 at 13:07

0.4^2 == 0.16
0.4 == 0.16^0.5


It's not a statistical issue. It's a rounding issue in floating point arithmetics in the language. The same applies to Python

0.4**2 == 0.16
0.4 == 0.16**0.5


None of these floating point numbers is stored in a computer's memory as two- or three-decimal digits, but rather in a binary format. For instance, 0.16 (8-byte) is stored as 11111111000100011110101110000101000111101011100001010001111011, while 0.4**2 gives us 11111111000100011110101110000101000111101011100001010001111100. You can play with this decimal to binary conversion using the following Python code:

import struct
bin(struct.unpack('!Q',struct.pack('!d',0.16))[0])
bin(struct.unpack('!Q',struct.pack('!d',0.4**2))[0])

• Hey Alex - I understand that the above example from mkt is a issue of floating point, but still cannot quite grasp why the two-digit decimal point example of yours(0.4**2) is not identical with 0.16. Can you please further elaborate? Sep 26 at 20:06
– Alex
Sep 26 at 20:34
• Thanks so much! Sep 26 at 20:38
• @Hwi-youngLee If this answers your question, consider accepting it (or whichever other answer you see as best) by clicking on the green tick mark to the left of the answer.
– mkt
Sep 29 at 7:58

Goldberg, D. (1991). What every computer scientist should know about floating-point arithmetic. ACM computing surveys (CSUR), 23(1), 5-48.

TL;DR you should never use == (or its equivalent in any other programming language) to compare floating point numbers because they represent numbers only approximately. As in your example, two different computations can lead to floating point representations that are different though close enough given the guarantees given by such representation. What you should do instead is to compare them with some precision, e.g. $$|x - y| < \epsilon$$, but this is also not a perfect solution and there are better ones. Happily, in most cases, they would be already implemented for you in the programming language you are using (like the mentioned all.equal in R).

• Dunno why this was downvoted. Just printing sd, one can see that it uses sqrt(var(x)) to perform the calculation. The devil is thus hidden in floating point arithmetic. And one needs to be aware of the problems hidden there. Sep 27 at 2:25
• +1 That paper is a classic. Sep 27 at 4:41

Correct to 14 decimal points, it's a floating point number issue:

options(digits = 20)
> var(data)
[1] 27.476190476190474499
sd(data)^2
[1] 27.476190476190470946


Check the differences between those quantities. They are likely to be small.

Computers are funky about how they store long decimals.

• Thanks, it seems that the difference is indeed very small, but still feeling kind of awkward that it produced FALSE... Sep 26 at 18:11