Stirling numbers

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

Stirling numbers

Robin Hankin
Hi

anyone coded up Stirling numbers in R?

[I need unsigned Stirling numbers of the first kind]


cheers

Robin




--
Robin Hankin
Uncertainty Analyst
National Oceanography Centre, Southampton
European Way, Southampton SO14 3ZH, UK
  tel  023-8059-7743

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Stirling numbers

Martin Maechler
>>>>> "Robin" == Robin Hankin <[hidden email]>
>>>>>     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<-------------------------------------------------------------------

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Stirling numbers

davidr-2
In reply to this post by Robin Hankin
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: [hidden email]
[mailto:[hidden email]] On Behalf Of Martin Maechler
Sent: Wednesday, July 19, 2006 9:19 AM
To: Robin Hankin
Cc: [hidden email]
Subject: Re: [R] Stirling numbers

>>>>> "Robin" == Robin Hankin <[hidden email]>
>>>>>     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<-----------------------------------------------------------------
--

______________________________________________
[hidden email] 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.

______________________________________________
[hidden email] 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.