Re: [R] Stirling numbers

From: Martin Maechler <maechler_at_stat.math.ethz.ch>
Date: Thu 20 Jul 2006 - 00:19:15 EST

>>>>> "Robin" == Robin Hankin <r.hankin@noc.soton.ac.uk> >>>>> on Wed, 19 Jul 2006 14:25:21 +0100 writes:

    Robin> Hi anyone coded up Stirling numbers in R?

"Sure" ;-)

    Robin> [I need unsigned Stirling numbers of the first kind]

but with my quick search, I can only see those for which I had "2nd kind" :

-o<--o<-------------------------------------------------------------------

##-- Stirling numbers of the 2nd kind
##-- (Abramowitz/Stegun: 24,1,4 (p. 824-5 ; Table 24.4, p.835)

##> S^{(m)}_n = number of ways of partitioning a set of $n$ elements into $m$
##> non-empty subsets

Stirling2 <- function(n,m)
{

    ## Purpose:  Stirling Numbers of the 2-nd kind
    ## 		S^{(m)}_n = number of ways of partitioning a set of
    ##                      $n$ elements into $m$ non-empty subsets
    ## Author: Martin Maechler, Date:  May 28 1992, 23:42
    ## ----------------------------------------------------------------
    ## Abramowitz/Stegun: 24,1,4 (p. 824-5 ; Table 24.4, p.835)
    ## Closed Form : p.824 "C."
    ## ----------------------------------------------------------------

    if (0 > m || m > n) stop("'m' must be in 0..n !")     k <- 0:m
    sig <- rep(c(1,-1)*(-1)^m, length= m+1)# 1 for m=0; -1 1 (m=1)     ## The following gives rounding errors for (25,5) :     ## r <- sum( sig * k^n /(gamma(k+1)*gamma(m+1-k)) )     ga <- gamma(k+1)
    round(sum( sig * k^n /(ga * rev(ga)))) }

options(digits=15)
for (n in 11:31) cat("n =", n," S_n(5) =", Stirling2(n,5), "\n") n <- 25
for(k in c(3,5,10))

    cat(" S_",n,"^(",formatC(k,wid=2),") = ", Stirling2(n,k),"\n",sep = "")

for (n in 1:15)

    cat(formatC(n,w=2),":", sapply(1:n, Stirling2, n = n),"\n")

##  1 : 1
##  2 : 1 1
##  3 : 1 3 1
##  4 : 1 7 6 1
##  5 : 1 15 25 10 1
##  6 : 1 31 90 65 15 1
##  7 : 1 63 301 350 140 21 1
##  8 : 1 127 966 1701 1050 266 28 1
##  9 : 1 255 3025 7770 6951 2646 462 36 1
## 10 : 1 511 9330 34105 42525 22827 5880 750 45 1
## 11 : 1 1023 28501 145750 246730 179487 63987 11880 1155 55 1
## 12 : 1 2047 86526 611501 1379400 1323652 627396 159027 22275 1705 66 1
## 13 : 1 4095 261625 2532530 7508501 9321312 5715424 1899612 359502 39325 2431 78 1
## 14 : 1 8191 788970 10391745 40075035 63436373 49329280 20912320 5135130 752752 66066 3367 91 1
## 15 : 1 16383 2375101 42355950 210766920 420693273 408741333 216627840 67128490 12662650 1479478 106470 4550 105 1


plot(6:30, sapply(6:30, Stirling2, m=5), log = "y", type = "l") ## Abramowitz says something like S2(n,m) ~ m^n / m!

##                                 ------------------
nn <- 6:30; sapply(nn, Stirling2, m=5) / (5^nn / gamma(5+1)) ## 0.114 0.21 ...... 0.99033 0.992266 0.993812
-o<--o<-------------------------------------------------------------------

______________________________________________
R-help@stat.math.ethz.ch 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. Received on Thu Jul 20 00:23:12 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Thu 20 Jul 2006 - 02:16:33 EST.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.