Thank you Dennis for your explanations!
The results you found are the same as mine. with first an infinity result,
followed by NaN. It seems that, when the number becomes too small, R must
round it up to 0. Hence I was wondering if there might a way to increase the
number of decimals for extremal computations. I have tried to use "double"
variables, but this didn't have a bigger success.
There is no doubt about the function and it orignates from a Kendall's plot
in measuring dependencies.
Thank you in advance,
Best Regards,
Nan
On Fri, Aug 20, 2010 at 23:48, Dennis Murphy <djmu...@gmail.com> wrote:
> 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 <nz...@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]]
______________________________________________
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.
