[R] Stirling numbers

[davidr at rhotrading.com]

Well, given Martin's 2nd kind function and the fact that the 1st and 2nd

kind matrices are inverses, you can get at least some of what you want
by:

... # fill a matrix S2 with second kind numbers and zeros
> S2
     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    1    0    0    0    0    0    0
[2,]    1    1    0    0    0    0    0
[3,]    1    3    1    0    0    0    0
[4,]    1    7    6    1    0    0    0
[5,]    1   15   25   10    1    0    0
[6,]    1   31   90   65   15    1    0
[7,]    1   63  301  350  140   21    1
> abs(round(solve(S2)))
     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    1    0    0    0    0    0    0
[2,]    1    1    0    0    0    0    0
[3,]    2    3    1    0    0    0    0
[4,]    6   11    6    1    0    0    0
[5,]   24   50   35   10    1    0    0
[6,]  120  274  225   85   15    1    0
[7,]  720 1764 1624  735  175   21    1 

Not sure how big arguments you need.

HTH,
David L. Reiner
Rho Trading Securities, LLC
Chicago  IL  60605
312-362-4963

-----Original Message-----
From: r-help-bounces at stat.math.ethz.ch
[mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Martin Maechler
Sent: Wednesday, July 19, 2006 9:19 AM
To: Robin Hankin
Cc: r-help at stat.math.ethz.ch
Subject: Re: [R] Stirling numbers

>>>>> "Robin" == Robin Hankin <r.hankin at 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 at 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.