Why does scale not make mean of 0?

Question

I am using R and I want to scale some data. The code looks like this:

data <- read.table(file_name, header = TRUE)
rates <- scale(data[8])
mean <- mean(rates)
sd <- sd(rates)

My understanding is that this scale function should scale the data so the mean is 0 and the standard deviation is 1. The standard deviation seems correct but the mean is not 0. What causes this? And what is the solution to making the mean 0? Or am I interpreting something wrong?

Floating point arithmetic isn't exact. https://stats.stackexchange.com/a/525766/22311 The print out says that the mean is less than $2 \times 10^{-16}$ units away from zero. How much closer to zero do you need it to be? — Sycorax, Jan 26 '23 at 20:56
What is the solution for this. It doesnt make sense to me that everywhere it says that scale makes the mean 0 and this function does not do that. — slipperypete, Jan 26 '23 at 21:03
Floating point arithmetic isn't exact & the inexactness of floating point arithmetic can't be fixed. The statement "scale makes the mean zero" means "scale makes the mean a value that is close to 0, relative to the machine precision of floating point representation." You can learn information about floating point arithmetic in its technical standard, IEEE 754. — Sycorax, Jan 26 '23 at 21:05
Have you ever computed $1/3$ on a calculator and then multiplied by $3$ to obtain $0.99999999$? According to the laws of arithmetic, the difference between this answer and the original value of $1$ is zero. The calculator thereby "proves" that $0 - 1 - 0.9999999 = 10^{-8}.$ The computer is just a big calculator and is subject to the same breakdown of mathematical laws. It is important to understand such things so you can use a computer wisely and well. — whuber, Jan 26 '23 at 21:11
You should never be reporting your results out from the R terminal anyway. If you were producing values to include in a paper or table, you'd need to use the ?formatC to specify such things as significant figures, rounding rules, etc. — AdamO, Jan 26 '23 at 21:14

score 4 · Answer 1 · edited Jan 27 '23 at 13:21

As others have noted, this is a programming side-effect more than a statistical irregularity. For my example below, I set the scipen argument in options super high so you can see how many decimals are produced (this argument just shows the scientific notation explicitly so you can see how long the numbers are). An example should make it obvious that even if the mean isn't exactly zero, it is as close as it gets with floating points.

To demonstrate, I have simulated normally distributed data with a mean of 50 and SD of 20 (the values don't really matter but felt like I'd assign something higher than the defaults which won't say much). I also set a random seed to any generic number so it is reproducible.

#### Simulate Data ####
set.seed(123)
options(scipen = 1000000000)
x <- rnorm(n=100,
           mean=50,
           sd=20)

Then I scale the data and check the mean and standard deviation like you did.

#### Scale and Get Mean/SD ####
scale.x <- scale(x)
mean(scale.x)
sd(scale.x)

When running mean(scale.x), you will see the number is ridiculously long:

0.0000000000000001974408

~~You may notice there are exactly 23 numbers to the right of the decimal. This is not by accident. Computers that run on 32-bit computers only store 23 digits in binary. This makes sense.~~ R would have to jam a theoretically infinite number of real numbers into a finite number of bits to get to zero. R instead creates a limited expression that is as close to zero as possible.

Running the function sd(scale.x) gives you exactly what you would expect though:

[1] 1

This is because unlike the mean, you are not rounding off a large list of numbers to get the answer. The standard deviation is just the squared root of the variance of $x$, which is usually going to be a non-zero value that will not have many decimals.

An accessible video to this topic can be found at this link.

Edit

While my overall point about the topic still stands, it appears I misspoke about the number of binary digits represented in this estimation. See the below comment by Whuber for more details.

23 is incorrect. IEEE double precision floating point explicitly stores 51 binary digits and implicitly has a leading 1, being the equivalent of $(51+1)/\log_2(10)\approx 16$ decimal digits in the mantissa. This is multiplied by a power of $2$ to shift the binary point, but you still get no more than 16 significant decimal digits no matter what. — whuber, Jan 27 '23 at 12:32
My mistake. I have edited the answer to include what you said. — Shawn Hemelstrand, Jan 27 '23 at 12:54

Why does scale not make mean of 0?

1 Answers1

Edit