Re: multiplication problem - what happens at 46341?

[Brian S Cade <brian_cade@usgs.gov>]

Okay, so if this an issue of integer
vs double precision real storage, why is this same calculation not a problem
in R.  What is done differently between S-Plus and R?   Curious
to avoid future problems.

Brian

Brian S. Cade

U. S. Geological Survey
Fort Collins Science Center
2150 Centre Ave., Bldg. C
Fort Collins, CO  80526-8818

email:  brian_cade@usgs.gov
tel:  970 226-9326

Sundar Dorai-Raj <sundar.dorai-raj@pdf.com>
Sent by: s-news-owner@lists.biostat.wustl.edu

03/03/2006 01:03 PM

To	henrik.parn@bio.ntnu.no
cc	s-news@lists.biostat.wustl.edu
Subject	Re: [S] multiplication problem - what happens at 46341?

Henrik Parn wrote:
> Dear all,
> 
> Why can't calculate 46341*46341 using '*'? I must have missed something
> very basic...
> 
>  > 46340.99999*46340.99999
> [1] 2147488280
>  > 46341*46341
> [1] NA
>  > 4.6341E4
> [1] 46341
>  > 4.6341E4*4.6341E4
> [1] 2147488281
>  > 46341^2
> [1] 2147488281
> 
> My function of for unbiased skew (Sokal & Rohlf 3ed p 115)  involves
a 
> multiplication, and when my samples are 'big' the answer is just 'NA'
> due to the problem described above:
> 
> skew.sokal <- function(x){
>    m3 <- length(x) * sum((x-mean(x))^3)/((length(x)-1)*(length(x)-2))
>    s3 <- sqrt(var(x))^3
>    m3/s3}
> 
>  > x <-  rnorm(46342)
>  > skew.sokal(x)
> [1] 0.0101318
>  > x <-  rnorm(46342)
>  > skew.sokal(x)
> [1] -0.00613662
>  > x <-  rnorm(46343)
>  > skew.sokal(x)
> [1] NA
> 
> Thanks in advance for any suggestion! And sorry if I have overseen
> something basic...
> S-plus 6.2 Build 6713
> WinXP
> 

Hi, Henrik,

You've overlooked how big an integer S-PLUS can store as an integer. 
That's why:

46342.*46342

will work, but

46342*46342

won't. This is becuase the first forces the storage to double.

As for your function, length(x) returns a integer. So, try:

skew.sokal <- function(x) {
   n <- as.numeric(length(x))
   m3 <- n * sum((x - mean(x))^3)/((n - 1) * (n - 2))
   s3 <- sqrt(var(x))^3
   m3/s3
}

set.seed(42)
x <- rnorm(46343)
skew.sokal(x)
--------------------------------------------------------------------
This message was distributed by s-news@lists.biostat.wustl.edu.  To
unsubscribe send e-mail to s-news-request@lists.biostat.wustl.edu with
the BODY of the message:  unsubscribe s-news