On 01/01/2015 2:43 PM, Mike Miller wrote: > On Thu, 1 Jan 2015, Duncan Murdoch wrote: > >> On 01/01/2015 1:21 PM, Mike Miller wrote: >>> On Thu, 1 Jan 2015, Duncan Murdoch wrote: >>> >>>> On 31/12/2014 8:44 PM, David Winsemius wrote: >>>>> >>>>> On Dec 31, 2014, at 3:24 PM, Mike Miller wrote: >>>>> >>>>>> This is probably a FAQ, and I don't really have a question about it, but >>>>>> I just ran across this in something I was working on: >>>>>> >>>>>>> as.integer(1000*1.003) >>>>>> [1] 1002 >>>>>> >>>>>> I didn't expect it, but maybe I should have. I guess it's about the >>>>>> machine precision added to the fact that as.integer always rounds down: >>>>>> >>>>>> >>>>>>> as.integer(1000*1.003 + 255 * .Machine$double.eps) >>>>>> [1] 1002 >>>>>> >>>>>>> as.integer(1000*1.003 + 256 * .Machine$double.eps) >>>>>> [1] 1003 >>>>>> >>>>>> >>>>>> This does it right... >>>>>> >>>>>>> as.integer( round( 1000*1.003 ) ) >>>>>> [1] 1003 >>>>>> >>>>>> ...but this seems to always give the same answer and it is a little >>>>>> faster in my application: >>>>>> >>>>>>> as.integer( 1000*1.003 + .1 ) >>>>>> [1] 1003 >>>>>> >>>>>> >>>>>> FYI - I'm reading in a long vector of numbers from a text file with no >>>>>> more than three digits to the right of the decimal. I'm converting them >>>>>> to integers and saving them in binary format. >>>>>> >>>>> >>>>> So just add 0.0001 or even .0000001 to all of them and coerce to integer. >>>> >>>> I don't think the original problem was stated clearly, so I'm not sure >>>> whether this is a solution, but it looks wrong to me. If you want to >>>> round to the nearest integer, why not use round() (without the >>>> as.integer afterwards)? Or if you really do want an integer, why add >>>> 0.1 or 0.0001, why not add 0.5 before calling as.integer()? This is the >>>> classical way to implement round(). >>>> >>>> To state the problem clearly, I'd like to know what result is expected >>>> for any real number x. Since R's numeric type only approximates the >>>> real numbers we might not be able to get a perfect match, but at least >>>> we could quantify how close we get. Or is the input really character >>>> data? The original post mentioned reading numbers from a text file. >>> >>> >>> Maybe you'd like to know what I'm really doing. I have 1600 text files >>> each with up to 16,000 lines with 3100 numbers per line, delimited by a >>> single space. The numbers are between 0 and 2, inclusive, and they have >>> up to three digits to the right of the decimal. Every possible value in >>> that range will occur in the data. Some examples numbers: 0 1 2 0.325 >>> 1.12 1.9. I want to multiply by 1000 and store them as 16-bit integers >>> (uint16). >>> >>> I've been reading in the data like so: >>> >>>> data <- scan( file=FILE, what=double(), nmax=3100*16000) >>> >>> At first I tried making the integers like so: >>> >>>> ptm <- proc.time() ; ints <- as.integer( 1000 * data ) ; proc.time()-ptm >>> user system elapsed >>> 0.187 0.387 0.574 >>> >>> I decided I should compare with the result I got using round(): >>> >>>> ptm <- proc.time() ; ints2 <- as.integer( round( 1000 * data ) ) ; >>>> proc.time()-ptm >>> user system elapsed >>> 1.595 0.757 2.352 >>> >>> It is a curious fact that only a few of the values from 0 to 2000 disagree >>> between the two methods: >>> >>>> table( ints2[ ints2 != ints ] ) >>> >>> 1001 1003 1005 1007 1009 1011 1013 1015 1017 1019 1021 1023 >>> 35651 27020 15993 11505 8967 7549 6885 6064 5512 4828 4533 4112 >>> >>> I understand that it's all about the problem of representing digital >>> numbers in binary, but I still find some of the results a little >>> surprising, like that list of numbers from the table() output. For >>> another example: >>> >>>> 1000+3 - 1000*(1+3/1000) >>> [1] 1.136868e-13 >>> >>>> 3 - 1000*(0+3/1000) >>> [1] 0 >>> >>>> 2000+3 - 1000*(2+3/1000) >>> [1] 0 >>> >>> See what I mean? So there is something special about the numbers around >>> 1000. >> >> I think it's really that there is something special about the numbers >> near 1, and you're multiplying that by 1000. >> >> Numbers from 1 to just below 2 are stored as their fractional part, with >> 52 bit precision. Some intermediate calculations will store them with >> 64 bit precision. 52 bits gives about 15 or 16 decimal places. >> >> If your number x is close to 3/1000, it is stored as the fractional part >> of 2^9 * x. This gives it an extra 2 or 3 decimal digits of precision, >> so that's why these values are accurate. >> >> If your number x is close to 2.003, it is stored as the fractional part >> of x/2, i.e. with errors like 1.0015 would have. So I would have >> guessed that 2.006 would have the same problems as 1.003, but I thought >> you didn't see that. So I tried it myself, and I do see that: >> >>> 1000+3 - 1000*(1+3/1000) >> [1] 1.136868e-13 >>> 2000+6 - 1000*(2+6/1000) >> [1] 2.273737e-13 >> >> Reading more closely, I see that you didn't test this particular case, >> so there's no contradiction here. >> >> The one thing I couldn't think of an explanation for is why other >> numbers between 1 and 2 don't have the same sorts of problems. So I >> tried the following: >> >> # Set data to 1.000 thru 1.999 >> data <- 1 + 0:999/1000 >> >> # Find the errors >> errors <- 1000 + 0:999 - 1000*data >> >> # Plot them >> plot(data, errors) >> >> The plot doesn't show a uniform distribution, but much more uniform than >> yours: so I think your data doesn't really cover all possible values >> from 0.000 to 1.999. (I get a similar plot if I look at cases where >> ints != ints2 with my data.) > > No, my data do cover all of the values in the range (I checked that by > listing all of them in a file: > > tr ' ' '\n' < FILE | uniq | sort -n | uniq >| uniq.txt > > They are definitely all there -- all 2001 of them. The thing is, your way > of making the numbers is a little different from reading them from a file. > You construct the number through an arithmetic operation and you get the > error: > >> as.integer( 1000 * (1+118/1000) ) > [1] 1117 > > (Your first error-producing value was 1.118.) But when the number is read > froma file, it starts as character and is translated to numeric. So I > start with the arithmetic, convert to character, back to numeric, and then > I do not see the error: > >> as.integer( 1000 * as.numeric( as.character( 1+118/1000 ) ) ) > [1] 1118 > > That reflects what was happening in my work. To see what I'm seeing, you > have to do this: > > data <- 1 + 0:999/1000 > errors <- 1000 + 0:999 - 1000 * as.numeric( as.character( data ) ) > plot(data, errors)
That's really interesting! I guess the differences basically cancel out. Duncan Murdoch > > Or you could cover he full range from 0 to 2: > > data <- 0:2000/1000 > errors <- 0:2000 - 1000 * as.numeric( as.character( data ) ) > plot(data, errors) > > Mike > ______________________________________________ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
