スターリング数と多重分割を計算する

スターリング数の漸化式

前章では重複順列から段階を踏んで多重分割を生成し、包除原理を使いスターリング数を計算しました。ここでは、スターリング数の興味深い漸化式を使い、個数を計算し多重分割を生成します。

スターリング数を漸化式を使って使って計算する

SymPyライブラリのstirling関数による一覧表の作成

漸化式を考える上でのヒントにするため、SymPyライブラリの#5のstirling関数を使ってn=10までのスターリング数の一覧表を作成します。

表を作成する関数

  1. def print_table(array, w):
  2. print(f' ' * (w+2), end='')
  3. for k in range(len(array[0])):
  4. print(f'{k:^{w+2}}', end=' ')
  5. print('sigma')
  6. print(f'-' * (w+1), end='+')
  7. for k in range(len(array)+1):
  8. print('-' * (w+2), end='+')
  9. print()
  10. for n in range(len(array)):
  11. print(f'{n:^{w+1}}', end='|')
  12. sigma = 0
  13. for k in range(len(array[0])):
  14. print(f"{array[n][k]:>{w+2}}", end=' ')
  15. sigma += array[n][k]
  16. print(f"{sigma:>{w+2}}")
  18. n=10;k=10
  20. array = [[0]*(k+1) for _ in range(n+1)]
  21. for i in range(1,n+1):
  22. for j in range(1,k+1):
  23. array[i][j]=int(stirling(i,j))
  25. print_table(array, 3)
 重複順列の集合Kを区別なしにする

3.　あらかじめ10×11の配列を初期化しておきます。2次元の配列を初期化するときはリスト内包表記を使わないと、集計でエラーになってしまいます。

SymPyライブラリのstirling関数を使ってスターリング数を計算し、arrayを書き換えます。

スターリング数の一覧表から漸化式を導く

表を眺めていると、スターリング数には次3つの法則があることがわかります。

①　k=1の場合

学生がn人いてタクシー1台に乗る場合は、分乗する余地もなく1通りしかありません。タクシーの定員を考えなければ、nがどんなに大きくてもスターリング数は常に1になります。

②　n=kの場合

n人でn台のタクシーに分乗する場合は、1人1台になります。タクシーにラベルはつけない前提なので、組分けを考える余地はなくスターリング数は1になります。

③　①②以外の多くの場合

前例のようにS(6, 4)つまり、6人が4台のタクシーに分かれて市内観光をする予定にしていました。ところが急遽学生が１人増えることになり、これにともないタクシーを1台増やすと考え、スターリング数はS(7, 5)となります。

S(6, 4)とS(7, 5)の差は、新しく来た学生の扱いにより、次の2つのケースが考えられます。

１　学生6人がタクシー4台への分割はそのままにして、新しく来た学生は1人でタクシーに乗る。この場合、新しい分割の数はS(6, 4)と同じです。

２　タクシーが増える前提で予め6人を5台に分けておき、新しく来た学生は5台のうちいずれかに乗る。この場合、新しい分割の数はS(6, 5)×5になります。

1,2は背反で同時に起こることはありえないので、S(7,5)=S(6,4)+S(6,5) ×5=65+15×5=140となるはずです。

漸化式が正しいか、個別に検証する

上記の数式を、SymPyライブラリのstirling関数を使って検証します。

漸化式でスターリング数を計算

  1. print(f'S(6, 4)= {stirling(6, 4)} S(6, 5)={stirling(6, 5)}')
  2. print(f'S(7, 5)={stirling(7, 5)}')
  3. print(f'S(6, 4)+5×S(6, 5)={stirling(6, 4)+5*stirling(6, 5)}')
S(6, 4)= 65  S(6, 5)=15
S(7, 5)=140
S(6, 4)+5×S(6, 5)=140

3.　1.で計算したS(6, 4)とS(6, 5)と2.で計算したS(7, 5)を使い、漸化式が正しいことを確認します。

(1)の考え方が正しいこと確認できました。漸化式を一般化すると次の式になります。

スターリング数の漸化式

① k==1のとき1

② n==kのとき1

③ $S(n, k)=S(n-1, k-1)+k\; S(n-1, k)$

漸化式を使てスターリング数を計算し、多重分割を生成する

漸化式を使ってスターリング数を計算する関数

この漸化式をもとに再帰関数を使いスターリング数を求めるstirling_r関数を作成し、S(7, 5)を計算します。

漸化式を使ったスターリング数の計算

  1. def stirling_r(n, k):
  2. if k==1 or n==k:
  3. return 1
  4. return stirling_r(n-1, k-1) + k*stirling_r(n-1, k)
  6. stirling_r(7, 5)
140

2.　漸化式の①、②に対応します。

4.　漸化式の③に対応します。stirling関数からstirling関数を引数をnやkを1つずつ減らしながら呼び出します。引数n,kの値はどんどん減っていくので最終的には2.の条件を満たし、1を返します。

このように簡単なプログラムで、複雑な計算ができるのは驚きです。s(7, 5)の結果は#9と同じく140となります。

さらに、漸化式を使いスターリング数の一覧を作成し、SymPyライブラリのstirling関数と比較し正しく計算できていることを確認します。漸化式を使った計算では、#10で作成したstirling_r関数を使うことも考えられますが、100通りの計算に対して1つ1つ関数を適用するのは相当な処理量になるので、漸化式の定義を使う方法でコンピュータへの負荷を減らすようにします。

漸化式の考え方を使ったスターリング数の一覧を関数と比較

  1. n=10;k=10
  2. recursive_array = [[0]*(k+1) for _ in range(n+1)]
  3. for i in range(1,n+1):
  4. recursive_array[i][1]=1
  5. recursive_array[i][i]=1
  6. for i in range(1,n+1):
  7. for j in range(2,i):
  8. recursive_array[i][j]=recursive_array[i-1][j-1]+j*recursive_array[i-1][j]
  9. recursive_array == stirling_sympy
  10. print_table(recursive_array, 4)
True

2.　漸化式による一覧を作成するため、あらかじめ10×10の配列を用意し、全ての要素に0で初期化します。

3.　2からi=1、n=kのときスターリング数は1になるので、該当する箇所を1に変更します。

6.　1からnまでをiでカウントアップし、各々のiについて2からi-1まで（iの場合は.3.で変更済）まで、漸化式の③を使いスターリング数を計算します。

漸化式を使い、多重分割を生成する

同じ漸化式を使い、多重分割をリストとし一覧を作成し、以前作成したリストと比較します。ここではS(6, 4)を生成するとして、

漸化式を使った多重分割の生成の考え方
漸化式を使った多重分割の生成の考え方

図のとおり、S(n-1, k-1)においては、新しく増えた5は単独の新しい多重分割になるので、リストとして追加します。一方S(n-1, k-1)に対しては、新しく増えた5はすでに4つある多重分割のいずれかに追加されるので、4つの重複分割に順次追加していきます。この考え方により、多重分割を生成する関数multiset_funcを作成します。

漸化式を使った重複分割の作成

  1. def multiset_func(n, k):
  2. if k == 1:
  3. return [[list(range(n))]]
  4. if n == k:
  5. return[[[i] for i in range(n) ]]
  6. array=[]
  7. for elm in multiset_func(n-1, k-1):
  8. array.append(elm+[[n-1]])
  9. for elm in multiset_func(n-1, k):
  10. for i in range(len(elm)):
  11. array.append(elm[:i] + [elm[i] + [n-1]] + elm[1+i:])
  12. return array
  14. sorted(stirling_surjective_list)==sorted(multiset_func(6, 4))
True

2～3　k=1の場合0からn-1までの要素が1つの分割になります。例えば、n=5の場合[[0, 1, 2, 3, 4]]となります。

4～5　k=nの場合、0からn-1までの要素が1つずつの分割になります。例えば、n=5の場合[[0], [1], [2], [3], [4]]となります。

7～8　漸化式③の左辺の処理をします。multiset_func(n-1, k-1)で返された多重分割について、新しいnは単独の分割になるので、これまでの分割に[[n-1]]を結合して返します。

9～11 漸化式③の左辺の処理をします。multiset_func(n-1, k)で返された多重分割は、分割のリストがk個あります。新しいn-1はこのk個の分割のいずれかに追加して返します。

ictsr4 All Rights Reserved.