R help - Aug 2010 - Problem to compute a function with very large numbers

Dear R users,

I have been trying to compute the following function and need it to work
with n=15000, but it would only compute for smaller ns, such as n=1000 and
not above. I was wondering if anyone would have a solution for this problem!
Thank you very much for your kind support!

Sincerely,

Nan

------

Wi <- function(n)    {
    fun <- function(w,i){
        w*(-log(w))*((w-w*log(w))^(i-1))*((1-w+w*log(w))^(n-i));
    }
    out <- as.double(n);

    for(i in 1:n){
        out[i] <- integrate(fun,lower=0,upper=1,i=i)$value;
        out[i] <- out[i]*n*choose(n-1,i-1);
    }
    out
}

	[[alternative HTML version deleted]]

Hi:

I had no problem getting the function to 'run', but it appears the
execution
time starts to rise dramatically after 10000:
> system.time(x <- Wi(100))
   user  system elapsed
   0.03    0.00    0.03
> system.time(x <- Wi(1000))
   user  system elapsed
   0.28    0.00    0.28
> system.time(x <- Wi(10000))
   user  system elapsed
   1.60    0.00    1.61
> system.time(x <- Wi(15000))
   user  system elapsed
   2.56    0.00    2.56
> system.time(x <- Wi(100000))
   user  system elapsed
  57.39    0.72   58.25

I'm using 64-bit R-2.11.1 on a Windows 7 system with 8Gb RAM to play with.
If you're running 32-bit R on a machine with comparatively low memory (say
1Gb or 2Gb), I could see where you might have problems.

There is also a problem with NaN generation (not you!) shortly after n 1000. I
started picking them up at 1100, but didn't test back from
there:
> x1000 <- Wi(1000)
> sum(is.nan(x1000))
[1] 0
> x1100 <- Wi(1100)
> sum(is.nan(x1100))
[1] 232
> x1500 <- Wi(1500)
> sum(is.nan(x1500))
[1] 920
> x2000 <- Wi(2000)
> sum(is.nan(x2000))
[1] 1518
> x5000 <- Wi(5000)
> sum(is.nan(x5000))
[1] 4670
> sum(is.nan(x10000))
[1] 9722
> sum(is.nan(x15000))
[1] 14749

Evidently, numerical precision becomes a problem once n exceeds 1000, which
might explain why you are experiencing problems with getting the function to
run afterward.

I'm kind of wondering about the kernel of the function:  (w * (1 -
log(w)))^(i - 1) * (1 - w * (1 - log(w)))^(n - i). This is a beta kernel, so
it would seem  that you need to have 0 < w * (1 - log(w)) < 1   => 
log(w) <
w < 1 + log(w) for the kernel to be positive-valued. If you look at the
first few values of the output of Wi(1000), you'll see how small they are.
It might be worth looking at the function again just in case.

HTH,
Dennis


On Fri, Aug 20, 2010 at 9:48 AM, Nan Zhao <nzhao@student.ethz.ch> wrote:
> Dear R users,
>
> I have been trying to compute the following function and need it to work
> with n=15000, but it would only compute for smaller ns, such as n=1000 and
> not above. I was wondering if anyone would have a solution for this
> problem!
> Thank you very much for your kind support!
>
> Sincerely,
>
> Nan
>
> ------
>
> Wi <- function(n)    {
>    fun <- function(w,i){
>        w*(-log(w))*((w-w*log(w))^(i-1))*((1-w+w*log(w))^(n-i));
>    }
>    out <- as.double(n);
>
>    for(i in 1:n){
>        out[i] <- integrate(fun,lower=0,upper=1,i=i)$value;
>        out[i] <- out[i]*n*choose(n-1,i-1);
>    }
>    out
> }
>
>        [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help@r-project.org mailing list
> 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.
>
	[[alternative HTML version deleted]]