NOTES ON STATISTICS, PROBABILITY and MATHEMATICS


Bernoulli Numbers:


Bernoulli numbers are important in many fields. They can be generated by inverting Pascal’s triangle - or the matrix form of Pascal’s triangle, excluding the \(1\)’s on the right of triangle (and performing an alternate minus sign pre-“conditioning”, as explained here):


FINDING BERNOULLI NUMBERS IN THE INVERTED MATRIX OF PASCAL’S TRIANGLE:


The “second Pascal matrix” is the lower triangular form of the Pascal matrix without the diagonal of \(1\)’s:

library(MASS)

n = 11
second.Pascal = matrix(0,n,n)
for(i in 1:n){
  for(k in 0:(i-1)) 
    second.Pascal[i,k+1] = choose(i,k)
} 
(M = second.Pascal)
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
##  [1,]    1    0    0    0    0    0    0    0    0     0     0
##  [2,]    1    2    0    0    0    0    0    0    0     0     0
##  [3,]    1    3    3    0    0    0    0    0    0     0     0
##  [4,]    1    4    6    4    0    0    0    0    0     0     0
##  [5,]    1    5   10   10    5    0    0    0    0     0     0
##  [6,]    1    6   15   20   15    6    0    0    0     0     0
##  [7,]    1    7   21   35   35   21    7    0    0     0     0
##  [8,]    1    8   28   56   70   56   28    8    0     0     0
##  [9,]    1    9   36   84  126  126   84   36    9     0     0
## [10,]    1   10   45  120  210  252  210  120   45    10     0
## [11,]    1   11   55  165  330  462  462  330  165    55    11

The equation for the sums of powers of natural numbers can extracted from this second Pascal matrix with some sign changes in alternating diagonals:

m = M * (((row(M) - col(M))%%2==1) * -1) +
M * (((row(M) - col(M))%%2==0) * 1) 
m
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
##  [1,]    1    0    0    0    0    0    0    0    0     0     0
##  [2,]   -1    2    0    0    0    0    0    0    0     0     0
##  [3,]    1   -3    3    0    0    0    0    0    0     0     0
##  [4,]   -1    4   -6    4    0    0    0    0    0     0     0
##  [5,]    1   -5   10  -10    5    0    0    0    0     0     0
##  [6,]   -1    6  -15   20  -15    6    0    0    0     0     0
##  [7,]    1   -7   21  -35   35  -21    7    0    0     0     0
##  [8,]   -1    8  -28   56  -70   56  -28    8    0     0     0
##  [9,]    1   -9   36  -84  126 -126   84  -36    9     0     0
## [10,]   -1   10  -45  120 -210  252 -210  120  -45    10     0
## [11,]    1  -11   55 -165  330 -462  462 -330  165   -55    11

The Faulhaber matrix is the inverse:

(Faulhaber = fractions(solve(m)))
##       [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9]  [,10] [,11]
##  [1,]     1     0     0     0     0     0     0     0     0     0     0
##  [2,]   1/2   1/2     0     0     0     0     0     0     0     0     0
##  [3,]   1/6   1/2   1/3     0     0     0     0     0     0     0     0
##  [4,]     0   1/4   1/2   1/4     0     0     0     0     0     0     0
##  [5,] -1/30     0   1/3   1/2   1/5     0     0     0     0     0     0
##  [6,]     0 -1/12     0  5/12   1/2   1/6     0     0     0     0     0
##  [7,]  1/42     0  -1/6     0   1/2   1/2   1/7     0     0     0     0
##  [8,]     0  1/12     0 -7/24     0  7/12   1/2   1/8     0     0     0
##  [9,] -1/30     0   2/9     0 -7/15     0   2/3   1/2   1/9     0     0
## [10,]     0 -3/20     0   1/2     0 -7/10     0   3/4   1/2  1/10     0
## [11,]  5/66     0  -1/2     0     1     0    -1     0   5/6   1/2  1/11

Compare to the slide on Mathologer:

which is immediately apparent in the Faulhaber matrix, here with minimal rearrangements:

flipF = Faulhaber[, ncol(Faulhaber):1]

N = matrix(0,nrow(flipF),ncol(flipF))
for(j in 1:nrow(flipF)){
  for(i in 1:j)
    N[j,j -i+1] <- flipF[j, ncol(flipF) - i + 1] 
} 

rownames(N) <- c(paste('S',0:(nrow(N)-1)))
fractions(N)


SUM OF POWERS

For instance, in the matrix above, named the Faulhaber matrix, the sum of the first \(n\) natural numbers (to the power of \(1\)) would be given by:

\[S_1 = 1/2 \;n + 1/2 \; n^2= \frac{n\;(n+1)}2\]

whereas the sum of the first \(n\) powers of \(5\) would be given by

\[S_4= -1/30 \;n + 1/3 \; n^3 + 1/2 \; n^4 + 1/5 \; n^5\]

The general closed form for the sum of powers is derived below (*).

The Bernoulli numbers are in the first column, and the rest of the numbers in each diagonal depends on that first value:

\[1/2 \quad 1/2 \quad 1/2 \quad 1/2 \quad \dots\] or

\[2/12 \quad 3/12 \quad 4/12 \quad 5/12 \quad 6/12 \quad 7/12 \quad \dots\]

Therefore the Bernoulli numbers can be extracted as:

(Bernoulli = Faulhaber[,1])
##  [1]     1   1/2   1/6     0 -1/30     0  1/42     0 -1/30     0  5/66

The sum of each row in the Faulhaber matrix adds to \(1\):

rowSums(Faulhaber)
##  [1] 1 1 1 1 1 1 1 1 1 1 1

which is part of the method explained by John Conway in here:

The difference between each natural number \(x\) to the \(s\) power, e.g. \(x^s= 3^s\)

in

\[1^s + 2^s + \color{red}{3}^s + 4^s+ \cdots + n^{s}\]

and the preceding sum to an immediately lower power:

\[1^{s-1} + 2^{s-1} + \color{red}{3}^{s-1} + 4^{s-1}+ \cdots + n^{s-1}\]

is

\[S_s(x)-S_{s-1}(s)=x^s\]

Differentiating,

\[S_s'(x)-S'_{s-1}(s)=s\;x^{s-1}\]

Therefore, for each summand, moving from one formula of a sum of powers to the next one amounts to multiplying the previous line times \(s\) and integrating term by term. For instance, to go from the sum of powers of \(5\) to the powers of \(6\) first we must multiply the row on the formula for the powers of \(5\) by \(6:\)

\[5\times\left(-1/30\; x + 1/3 \; x^3 + 1/2 \; x^4 + 1/5\; x^5\right)\]

and then integrating:

\[-1/12\; x^2 + 5/12 \; x^4 + 1/2 \; x^5 + 1/6\; x^6\]

Now for the punch-line all the factors have to add to \(1\) because this equation must be valid for all sums up to a given natural number \(n.\) Therefore if \(n=1\) all the powers will be \(1.\) Therefore, the missing number, i.e. the Bernoulli number is

\[B_5 = 1 + 1/12 - 5/12 - 1/2 -1/6=0\]

which is indeed the case.


CLOSED FORM AND GENERATING FUNCTION:

\[B_n=\sum_{k=0}^n\frac1{k+1}\sum_{j=0}^k(-1)^j\binom{k}{j}j^n,\quad n\ge0\]

The exponential generating function is:

\[\frac{x}{e^x-1} = \sum_{k=0}^\infty B_k \frac{x^k}{k!}\] This is frightening, but using Wofram Alpha:

series expansion of f(x) = x/(exp(x) -1)

yieds the following monster:

\[1 - \frac 1 2 x + \frac 1{12} x^2 - \frac 1{720}x^4 + \frac 1 {30240}x^6 - \frac 1{1209600} x^8+ \frac 1 {47900160}x^{10}-\frac {691}{1307674368000}x^{12}+\cdots\]

Not the Bernoulli numbers! So how to get them? Multiply each term \(a_n\) by \(n!\):

\[\begin{align} 1 \times 0! &= 1\\ -1/2 \times 1! &= -1/2 \\ 1/12 \times 2! &= 1/6 \\ -1/720 \times 4! &= - 1/30 \\ 1/30240 \times 6! & = 1/42 \\ -1/1209600 \times 8! &= -1/30 \\ 1/47900160 \times 10! &= 5/66 \\ -691/1307674368000 \times 12! &= -691/2730 \end{align}\]

which are Bernoulli numbers:

This can also be seen as follows:

  1. Adding and subtracting in the numerator \(t+\frac{t^2}{2!}+\frac{t^3}{3!}+\frac{t^4}{4!}+...\) the \(t\) is eliminated:

\[\frac{t}{e^t-1}=\frac{t}{t+\frac{t^2}{2!}+\frac{t^3}{3!}+\frac{t^4}{4!}+...}=1+\frac{(1-1)t-\frac{t^2}{2!}-\frac{t^3}{3!}-\frac{t^4}{4!}-...}{t+\frac{t^2}{2!}+\frac{t^3}{3!}+\frac{t^4}{4!}+...}=1-\frac{+\frac{t^2}{2!}+\frac{t^3}{3!}+\frac{t^4}{4!}+...}{t+\frac{t^2}{2!}+\frac{t^3}{3!}+\frac{t^4}{4!}+...}\] 2. To now get the second Bernoulli coefficient, from the numerator we subtract from the numerator \(\frac {t} {2} \small \left(t+\frac{t^2}{2!}+\frac{t^3}{3!}+\frac{t^4}{4!}+... \right)\):

\[\frac{t}{e^t-1}=1-\frac{t}{2}+\frac{+(\frac{1}{2}-\frac{1}{2!})t^2+(\frac{1}{2.2!}-\frac{1}{3!})t^3+(\frac{1}{2.3!}-\frac{1}{4!})t^4+...}{t+\frac{t^2}{2!}+\frac{t^3}{3!}+\frac{t^4}{4!}+...}=1-\frac{t}{2}+\frac{(\frac{1}{2.2!}-\frac{1}{3!})t^3+(\frac{1}{2.3!}-\frac{t^4}{4!})t^4+...}{t+\frac{t^2}{2!}+\frac{t^3}{3!}+\frac{t^4}{4!}+...}\] 3. And so forth

\[\frac{t}{e^t-1}=1-\frac{1}{2}t+\frac{\frac{1}{12}t^3+\frac{1}{24}t^4+...}{t+\frac{t^2}{2!}+\frac{t^3}{3!}+\frac{t^4}{4!}+...}=1-\frac{1}{2}t+\frac{1}{12}t^2+\frac{(\frac{1}{24}-\frac{1}{2.12})t^4+...}{t+\frac{t^2}{2!}+\frac{t^3}{3!}+\frac{t^4}{4!}+...}\]


(*) Closed form derivation:

The first series are (see here):

\[\begin{align} S_0(n) &= n\\ S_1(n) &= 1/2 \, n^2 + 1/2 \, n\\ S_2(n) &= 1/3 \, n^3 + 1/2 \,n^2 + 1/6 \,n\\ S_3(n) &= 1/4 \, n^4 + 1/2 \,n^3 + 1/4 \, n^2\\ S_4(n) &= 1/5 \, n^5 + 1/2 \,n^4 + 1/3 \, n^3 - 1/30 \, n\\ S_5(n) &= 1/6 \, n^6 + 1/2 \,n^5 + 5/12\, n^4 - 1/12 \, n^2\\ S_6(n) &= 1/7 \, n^7 + 1/2 \,n^6 + 1/2 \, n^5 - 1/6 \, n^3 + 1/42 \,n) \end{align}\]

A common factor corresponding to the power can be extracted in front for each sum \(S_p(n) \mapsto 1/(p+1)\):

\[\begin{align} S_0(n) &= 1\,n\\ S_1(n) &= 1/2 \; (n^2 + n)\\ S_2(n) &= 1/3 \; (n^3 + 3/2 \,n^2 + 1/2 \,n)\\ S_3(n) &= 1/4 \; (n^4 + 2 \,n^3 + n^2)\\ S_4(n) &= 1/5 \; (n^5 + 5/2 \,n^4 + 5/3 \,n^3 - 1/6 \, n)\\ S_5(n) &= 1/6 \; (n^6 + 3 \,n^55 + 5/2\, n^4 - 1/2 \,^2)\\ S_6(n) &= 1/7 \; (n^7 + 7/2 \,n^6 + 7/2 \, n^5 - 7/6 \, n^3 + 1/6 \,n) \end{align}\]

Extracting common multiples by column:

\[\begin{align} S_0(n) &= 1\,n\\ S_1(n) &= 1/2 \; (n^2 +2 (1/2) \,n)\\ S_2(n) &= 1/3 \; (n^3 +3 (1/2) \,n^2 + 3 \,(1/6) \,n)\\ S_3(n) &= 1/4 \; (n^4 +4 (1/2) \,n^3 + 6 \,(1/6) \,n^2)\\ S_4(n) &= 1/5 \; (n^5 +5 (1/2) \,n^4 + 10 (1/6) \,n^3 - 5 (1/30) n)\\ S_5(n) &= 1/6 \; (n^6 +6 (1/2) \,n^5 + 15 (1/6) \,n^4 - 15 (1/30) n^2)\\ S_6(n) &= 1/7 \; (n^7 +7 (1/2) \,n^6 + 21 (1/6) \,n^5 - 35 (1/30) n^3 + 7(1/42)n) \end{align}\]

Now the integers are all binomial coefficients corresponding to the row:

\[\small\begin{align} S_0(n) &= 1\; \; \left(\binom 1 0 \color{red} 1\,n\right)\\ S_1(n) &= \frac 1 2 \; \left(\binom 2 0 n^2 +\binom 2 1 \left(\color{red}{\frac 1 2} \right) \,n \right)\\ S_2(n) &=\frac 1 3 \; \left(\binom 3 0 n^3 +\binom 3 1 \left(\frac 1 2 \right) \,n^2 + \binom 3 2 \,\left(\color{red}{\frac 1 6} \right) \,n \right)\\ S_3(n) &= \frac 1 4 \; \left(\binom 4 0 n^4 + \binom 4 1 \left(\frac 1 2 \right) \,n^3 + \binom 4 2 \,\left(\frac 1 6 \right) \,n^2 + \binom 4 3 \,\color{red} 0 \,n \right)\\ S_4(n) &= \frac 1 5 \; \left(\binom 5 0 n^5 + \binom 5 1 \left(\frac 1 2 \right) \,n^4 + \binom 5 2 \left(\frac 1 6\right) \,n^3 - \binom 5 3 0 \,n^2+ \binom 5 4 \left(\color{red}{\frac 1 {30}} \right) n\right)\\ S_5(n) &= \frac 1 6 \; \left(\binom 6 0 n^6 +\binom 6 1 \left(\frac 1 2 \right) \,n^5 + \binom 6 2 \left(\frac 1 6 \right) \,n^4 + \binom 6 3 0\,n^3 - \binom 6 4 \left(\frac 1{30} \right) n^2 + \binom 6 5 \color{red} 0\, n\right)\\ S_6(n) &= \frac 1 7 \; \left(\binom 7 0 n^7 + \binom 7 1 \left(\frac 1 2 \right) \,n^6 + \binom 7 2 \left(\frac 1 6 \right) \,n^5 + \binom 7 3 0\,n^4- \binom 7 4 \left(\frac 1 {30} \right) n^3 + \binom 7 5 0\,n^2 + \binom 7 6 \left(\color{red}{\frac 1{42}} \right)n\right) \end{align}\]

leading to the final formula:

\[S_p(n) = \frac 1{p+1} \left( \sum_{k=1}^{p+1}\;\left(-1\right)^{\delta^k_p}\;\binom{p+1}k B^-_{p+1-k}\;n^k \right)\]

The actual proof is in this Michael Penn presentation.


Connections and Appearances of the Bernoulli Numbers:

Calculation of the zeta function for even multiples (Euler):

The Euler formula for even multiples of the the zeta function is:

\[\zeta(2 k)=\frac{(-1)^{k-1}(2 \pi)^{2 k} B_{2 k}}{2(2 k) !}\]

Riemann zeta functional equation with Bernoulli numbers:

Riemann formulated the following equation linking the zeta function and Bernoulli numbers:

\[\zeta(1-k)=-\frac{B_{k}}{k}\]

Generalized Bernoulli numbers:

For a Dirichlet character \(\chi,\) the functional Riemann equation right above can be extended to L-functions as:

\[L(1-k, \chi)=:-\frac{B_{k, \chi}}{k}\]

Dirichlet L-functions:

L-functions are infinite sums of the form of a Dirichlet series:

\[L(s, \chi)=\sum_{n=1}^\infty\frac{\chi_n}{n^s}\] whereas Bernoulli numbers appear in the sums of \(n\) powers, for \(k\) a positive whole number:

\[1^k + 2^k + 3^k + \cdots + n^k\]

Thinking of \(k\) above as a negative number, there has to be a connection to the Dirichlet series or the Dirichlet characters \(\chi_n,\) even though the L-functions sum to infinity, and the Bernoulli numbers appear in finite sums.

There is a relationship as in

\[\sum_{n \geq 0} B_{n, \chi} \frac{T^{n}}{n !}=\sum_{a=1}^{f} \frac{\chi(a) T e^{a T}}{e^{f T}-1}\]

in which \(f\) is the conductor of the Dirichlet character \(\chi,\) defined as:

The conductor of a Dirichlet character \(\chi\) modulo \(q\) is the least positive integer \(q_1\) dividing \(q\) for which \(\chi(n+kq_1)=\chi(n)\) for all \(n\) and \(n+kq_1\) coprime to \(q\). (From here).


Home Page

NOTE: These are tentative notes on different topics for personal use - expect mistakes and misunderstandings.