On 18-Nov-11 16:03:44, Gyanendra Pokharel wrote:
> Hi all,
> why factorial(150) shows the error out of range in 'gammafn'?
> I have to calculate the number of subset formed by 150 samples
> taking 10 at a time. How is this possible?
> best
Because factorial(150) is nearly 10^263, which is far greater
than any number that R is capable of storing.
I would, perhaps, suggest you work with lgamma(), the logarithm
(to base e) of the Gamma function.
Thus factorial(150) = gamma(151), and
lgamma(151)
# [1] 605.0201
so factorial(150) is close to e^605; to get it as a power of 10:
lgamma(151)*log10(exp(1))
# [1] 262.7569
Hence the "nearly 10^263" above.
If your question means "the number of ways of choosing 10 out
of 150, then, using lgamma(), its log to base e is:
lgamma(151) - lgamma(11) - lgamma(141)
# [1] 34.6954
and its log to base 10 is:
(lgamma(151) - lgamma(11) - lgamma(141))*log10(exp(1))
# [1] 15.06802
so the result is about 10^15 which *is* within R's range.
Hence:
10^( (lgamma(151) - lgamma(11) - lgamma(141))*log10(exp(1)) )
# 1.169554e+15
You can see this (approximately) in an all-digits form by
using print():
X <- 10^( (lgamma(151) - lgamma(11) - lgamma(141))*log10(exp(1)) )
print(X,18)
# [1] 1169554298222353
However, in many cases (such as your example) it will be feasible
to exploit the reduction in numbers of factors if you cancel out
the factors in the larger denominator factorial from the numerator
factorial, so that the number you want is
150*149*148*147*146*145*144*143*142*141/(10*9*8*7*6*5*4*3*2*1)
= (150/10)*(149/9)*(148/8)* ... *(142/2)*(141/1)
for which you can define a function nCr(n,r):
nCr <- function(n,r){
num <- (n:(n-r+1))
den <- (r:1)
prod(num/den)
}
and then:
print(nCr(150,10),18)
# [1] 1169554298222310
which is not quite the same as the result 1169554298222353
obtained above using lgamma(). This is because the previous
result is affected by rounding errors occuring in the
logarithms. The result obtained using the function nCr() is
exact. However, it can fail in its turn if you want one of
the larger possible results, such as nCr(1000,500):
nCr(1500,750)
# [1] Inf
since the result is too large for R to store (the largest on
my machine is about 2*10^308, see below)), and the products
in the function definition accumulate until this number is
exceeded.
You can find out the limits of R storage on your machine
by using the command
.Machine
which gives a long list of the characteristics of the machine
you are using. In particular, you will see:
$double.xmax
[1] 1.797693e+308
is the largest number that R will store; and:
$double.eps
[1] 2.220446e-16
which is the smallest positive number; and:
$double.neg.eps
[1] 1.110223e-16
which is the smallest negative number (note that the latter
is half the former).
They can be individually accessed as:
.Machine$double.xmax
# [1] 1.797693e+308
.Machine$double.eps
# [1] 2.220446e-16
.Machine$double.neg.eps
# [1] 1.110223e-16
Hoping this helps!
Ted.
--------------------------------------------------------------------
E-Mail: (Ted Harding) <ted.hard...@wlandres.net>
Fax-to-email: +44 (0)870 094 0861
Date: 18-Nov-11 Time: 21:31:04
------------------------------ XFMail ------------------------------
______________________________________________
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.
