[R] Permutations

[Ted Harding]

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.harding at wlandres.net>
Fax-to-email: +44 (0)870 094 0861
Date: 18-Nov-11                                       Time: 21:31:04
------------------------------ XFMail ------------------------------