スターリング数と集合分割

スターリング数について理解するため、次の問題を考えます。

問題1

出席番号0から5までの6人を4つのグループに分け、タクシーで市内を観光します。タクシーを区別しないとすると、考えられる組み合わせは何通りになりますか。なお、各タクシーには少なくとも1人は乗るものとします。

重複順列では学生の部屋割りを例に取り、学生のグループ分けが同じでも泊まる部屋が異なれば別の順列として扱いました。これに対し、スターリング数はタクシーを区別せず、だれとだれが同じ組になるかということだけに集中します。区別のある部屋はラベルあり（labeled）に対し、区別の無いタクシーはラベルなし(unlabeled)とし、その組み合わせの数をスターリング数、これにもとづく組み合わせをスターリング数に対応する集合分割（set partition）といいます。なお、「スターリング数に対応する集合分割」というのは長すぎるので、これ以降は、「集合分割」と呼ぶことにします。

全射の重複順列とスターリング数、集合分割の関係

写像から見たスターリング数と集合分割

重複順列、順列では、定義域の集合Nから値域の集合Kへの写像で、ともにラベルありで順列の生成と個数の計算を行いました。一方、重複組み合わせでは、集合Nの要素がラベルなしとし、集合Nの要素に対し増加関数の考え方を使いラベルを外しました。これに対し、スターリング数は、集合Nはラベルあり、集合Kはラベルなしにします。それでは、実施にどのようなものか、写像を使って考えます。

重複順列と重複組み合わせと集合分割との相違
重複順列と重複組み合わせと集合分割との相違

集合Kについてラベルなしにする前の重複分割の例を挙げると図表2のようになります。

集合Kを基準に集合の分割を整理
集合Kを基準に集合の分割を整理

2つの重複順列は別のものとして数えますが、集合Kのラベルを外すとA,B,C,Dの区別がなくなるので、(0, 1),(2,3),(4),(5)の組み合わせでそれぞれ1台のタクシーに乗るので同じ集合分割になります。とはいえ、重複順列の図を眺めていても、2つを同じ集合分割として捉えるのは容易ではありません。そこで、切り口を変え、A,B,C,Dの4つのタクシーにだれが乗るか、という視点でリストを変換するinverse_value_mapping関数を作成します。図表3で示すように、新たなリストnew_elmにk_listの要素の数だけ第２階層のリストを初期化し、A→0,B→1,C→2,D→3と順番をつけリストの添え字にします。そして、それぞれの重複順列の要素について、例えば出席番号2の学生がAの車に乗るとしたら、new_elm[0]に2を追加します。こうすることにより、どのタクシーにどの学生の組み合わせが乗るのかが明らかになります。

次に、増加関数の考え方を使い昇順に並んでいるものに絞り込みます。これについては以前作成したfilter_increasing_map関数を使用します。

集合Kについて区別なしの分割に変換
集合Kについて区別なしの分割に変換

右側の配列にすると、重複組み合わせのラベルを外すと全て同じ組み合わせになることわかります。k=4なので第2階層の要素数4のタプルについての順列、$k!=24$個が1つに集約されることがわかります。

ラベルありの全射の重複順列の生成

全射の重複順列はp90111で作成したgenerate_surjective_perm関数を使い、念のためcalc_surjective_perm関数を使い、全射の重複順列から抜き出したリストsurjective_perm_listを生成します。

全射の重複組み合わせを生成する関数（再掲）

  1. def generate_surjective_perm(k_list, n):
  2. perm_list = [()]
  3. for cnt in range(n):
  4. temp = perm_list
  5. perm_list = []
  6. for elm in temp:
  7. for item in k_list:
  8. if cnt < n - 1:
  9. perm_list.append(elm + (item,))
  10. else:
  11. if len(set(new_elm := elm + (item,))) == len(k_list):
  12. perm_list.append(new_elm)
  13. return perm_list
  15. k_list = ['A', 'B', 'C', 'D']
  16. n = 6
  17. surjective_perm_list = generate_surjective_perm(k_list, n)
  18. print(f'generate_surjective_perm({len(k_list)}, {n}) = {len(surjective_perm_list)}')
  19. surjective_perm_list
perm_surjective_func(4, 6) =  1560
[('A', 'A', 'A', 'B', 'C', 'D'),
 ('A', 'A', 'A', 'B', 'D', 'C'),
 ('A', 'A', 'A', 'C', 'B', 'D'),
 ('A', 'A', 'A', 'C', 'D', 'B'),
・
('A', 'A', 'B', 'B', 'D', 'C'),
・
('B', 'B', 'A', 'A', 'D', 'C'),
・
('D', 'D', 'D', 'B', 'C', 'A'),
 ('D', 'D', 'D', 'C', 'A', 'B'),
 ('D', 'D', 'D', 'C', 'B', 'A')]

15. k_listには、4台のタクシーに仮に付けたラベルを定義します。nは集合Nの個数、つまり学生の人数を定義します。

重複順列の数は1560通りになりますが、件数が多いので一部の表示にとどめます。全射の重複順列の個数は、包除原理を使った次の式で計算することができます。

包除原理を使った全射の重複順列の個数

$PIE(n, k)=\sum\limits_{i=0}^{k}(-1)^i\,{}_kC_i(k-i)^n$

さっそく、計算式にもとづき関数を使って計算します。以前作成したcalc_surjective_perm関数を使い全射の重複順列の数を計算します。

重複順列の全射の公式を使って個数を計算（再掲）

  1. def calc_surjective_perm(n, k):
  2. sigma_surjective = 0
  3. for i in range(k + 1):
  4. sigma_surjective += ((-1)**i) * comb(k, i, exact=True) * (k - i)**n
  5. return sigma_surjective
  7. print(f'calc_surjective_perm({len(k_list)}, {n}) = {calc_surjective_perm(n, len(k_list))}')
count_surjective_perm(4, 6) =  1560

7.　k_listの要素数=4とn=6を引数として渡します。

1560個と#1の結果と一致することを確認することができました。

全射のラベルありの重複順列をラベルなしにする

全射の重複順列を生成する

全射の重複順列では集合Kを基準に増加関数の考え方を適用するのが困難なので、切り口変えてどのタクシーにどの学生が乗るかを示す集合分割を生成する関数を作成します。関数名は値の配列を逆にするという意味でinverse_value_mapping とし、生成するリストをinverse_value_mapping_listとします。

集合Kの要素を基準に重複順列を整理する

  1. def inverse_value_mapping(k_list, perm_list):
  2. converted_list = []
  3. for elm in perm_list:
  4. new_elm = [() for _ in range(len(k_list))]
  5. for num, k in enumerate(elm):
  6. new_elm[k_list.index(k)] += (num,)
  7. converted_list.append(tuple(new_elm))
  8. return converted_list
  10. inverse_value_mapping_list = convert_to_inverse_value_mapping(k_list, surjective_perm_list )
  11. print(f'inverse_value_mapping_list = {len(inverse_value_mapping_list)}')
  12. inverse_value_mapping_list
calc_surjective_perm(4, 6) = 1560
[((0, 1, 2), (3,), (4,), (5,)),
 ((0, 1, 2), (3,), (5,), (4,)),
 ((0, 1, 2), (4,), (3,), (5,)),
 ((0, 1, 2), (5,), (3,), (4,)),
・
((0, 1), (2, 3), (4,), (5,)),
・
((2, 3), (0, 1), (5,), (4,)),
・
((5,), (4,), (3,), (0, 1, 2))]

4.　new_elmで、AからDの4台のタクシーに誰が乗るかを表す組み合わせの1つを集計します。このため、k_listの数k=4の空のリストをあらかじめ定義しておきます。初期化をするだけなのでインデックスは必要ないので”for _ in”とします。

6.　k_list.index(k)]は乗る車のk_listの中での順番を表します。これをindexとして、new_elmにそのタクシーに乗る人の出席番号numを追加します。

配列を並べ替えただけなので、まだ、この段階では件数は1560件と変わりません。

集合Kに増加写像を適用してラベル付き分割をラベル無しの分割に変換

前節で部屋を区別の無いタクシーに変えるため、組み合わせと同じように増加写像の考え方を使い、各組み合わせにp90120作成したfilter_increasing_map関数を適用します。このことにより( (0, 1), (2, 3), (4), (5))、((0, 1), (2, 3), (5) , (4))、((2, 3), (0, 1), (4), (5))などは、小さい順番から並んでいる((0, 1), (2, 3), (4), (5))に代表させることにより、ラベル無しの分割に集約することができます。

増加関数の考え方を使い集合Kをラベル無しの分割に変換

  1. def filter_increasing_map(perm_list):
  2. increasing_list = []
  3. for elm in perm_list:
  4. for i in range(len(elm) - 1):
  5. if elm[i] >= elm[i + 1]:
  6. break
  7. else:
  8. increasing_list.append(elm)
  9. return increasing_list
  11. unlabeled_partitions_list = filter_increasing_map(inverse_value_mapping_list)
  12. print(f'unlabeled_partitions_list = {len(unlabeled_partitions_list)}')
  13. unlabeled_partitions_list
unlabeled_partitions_list = 65
[((0, 1, 2), (3,), (4,), (5,)),
 ((0, 1, 3), (2,), (4,), (5,)),
 ((0, 1), (2, 3), (4,), (5,)),
 ((0, 1, 4), (2,), (3,), (5,)),
 ((0, 1), (2, 4), (3,), (5,)),
 ((0, 1), (2,), (3, 4), (5,)),
 ((0, 1, 5), (2,), (3,), (4,)),
 
・
・
((0,), (1, 5), (2,), (3, 4)),
 ((0,), (1,), (2, 5), (3, 4)),
 ((0,), (1,), (2,), (3, 4, 5))]

5.　変換された重複順列を1つずつelmとして読みこみ、要素の1番目と2番目、２番目と3番目・・というように次々と大きさを比べます。このとき、後の要素の方が小さい場合は増加写像の要件を満たさなくなるので、対象外となり6.でbreakします。

7.　5～6.に一度もあてはまらずループを抜けたものは増加写像となるので、surjective_listに追加します。

n=4,k=6のスターリング数は、包除原理により求めた全射の重複順列の1560通りをn!=4!=24で割った65通りに絞られました。

スターリング数の計算と集合分割を生成する関数

包除原理を使いスターリング数を計算する

スターリング数は、大文字のSを使いS(n, k)とあらわします。スターリング数は包除原理による全射の重複順列の個数をk!で割ることにより計算することができます。

スターリング数の計算

$S(n, k)=\sum\limits_{i=0}^{k}(-1)^i\,{}_kC_i(k-i)^n/k!$

スターリング数を計算する関数を作成します。関数名は文字通りcalc_stirling_countとし、n,kを引数とします。

スターリング数の計算式を関数化し、これまでの結果と一致しているか確認します。

スターリング数を計算する関数

  1. from math import factorial
  2. def calc_stirling_num(n, k):
  3. sigma = 0
  4. for i in range(k + 1):
  5. sigma += ((-1)**i) * comb(k, i) * (k - i)**n
  6. return int(sigma / factorial(k))
  8. print(f'calc_stirling_num({n}, {len(k_list)}) = {calc_stirling_num(n, len(k_list))}')
calc_stirling_num(6, 4) =  65

7.　#2 perm_surjective_counと同じ包除原理を使い、全射の重複順列の個数を計算し、戻り値を返すときに、集合Kの要素の個数kの階乗(n!)で割ります。割り算の結果は小数点型になるので、int関数で整数に変換します。

SymPyのstirling関数を使いスターリング数を計算する

Pythonの標準ライブラリを使いスターリング数と集合の分割を求め、前節での計算結果が正しいことを確認します。

スターリング数の計算にはSymPyライブラリのsympy.functions.combinatorial.numbersモジュールで提供される stirling関数を使用します。stirling関数は集合Nの個数nと集合Kの個数kを引数とします。

SymPyライブラリのstirling関数によるスターリング数の計算

  1. from sympy.functions.combinatorial.numbers import stirling
  2. print(f'sympy_stirling({n}, {len(k_list)}) = {stirling(n, len(k_list))}')
sympy_stirling(6, 4) =  65

2.　 sympy.functions.combinatorial.numbersからインポートします。

n=6,k=4で計算すると65と#4で導いた集合の分割の個数、及びSymPyライブラリのstirling関数の個数と一致していることがわかりました。

集合の分割を生成する関数

集合の分割には、SymPyライブラリのsympy.utilities.iterablesモジュールで提供されるmultiset_partitions関数を使用します。引数としてn,kを渡しますが、nは集合Nの要素をリストでも個数でも可能です。前節に作成した集合の分割と一致するか確認します。multiset_partitionsは集合を指定した数に分割するものですが、集合に同じ値の要素があることも許されます。このため、集合分割(setpartition)は多重分割の特殊なケースと考えらえます。

SymPyライブラリのmultiset_partitionsによる集合の分割の生成

  1. k_list= ['A', 'B', 'C', 'D' ]
  2. n = 6
  3. from sympy.utilities.iterables import multiset_partitions
  4. sympy_multiset_partitions_list = list(multiset_partitions(n, len(k_list)))
  5. sympy_multiset_partitions_list
[[[0, 1, 2], [3], [4], [5]],
 [[0, 1, 3], [2], [4], [5]],
 [[0, 1], [2, 3], [4], [5]],
 [[0, 1, 4], [2], [3], [5]],
 [[0, 1], [2, 4], [3], [5]],

[[0], [1], [2, 5], [3, 4]],
 [[0], [1], [2], [3, 4, 5]]]

1.　multiset_partitions関数はSymPyライブラリのutilities.iterablesサブモジュールからimportします。

生成した集合分割とSymPyライブラリでの計算結果を比較

これまでは、重複順列や組み合わせなど、itertoolsモジュールで組み合わせを生成してきました。ここでは一貫して、個々の集合分割を示す第2階層の配列はタプルを使うようにしているので、関数を自作する際もこれに合わせてきました。ところが、SymPyライブラリで生成する組み合わせは全てリストで生成するので、両者が一致しているかを確認するためには、どちらかに形式を合わせる必要があります。そこで、多重の配列をitertoolsに形式を合わせ、第1階層はリスト、第2階層以降はタプルに変換する関数を作成します。

関数名は、convert_to_tuple とします。ここでは、配列が第3階層以下にもわたっていることを想定して、再帰関数を使い、リストやタプルが現れない階層まで深堀りするようにします。

リストをタプルに変換する

  1. def convert_to_tuple(sequence, depth=1):
  2. if isinstance(sequence, (list,tuple)):
  3. temp_list = [convert_to_tuple(item, depth + 1) for item in sequence]
  4. if depth == 1:
  5. return temp_list
  6. else:
  7. return tuple(temp_list)
  8. else:
  9. return sequence
  11. print(convert_to_tuple([[1, 2], [[2], [3, 4]], [3]]))
  12. print(convert_to_tuple(((1, 2), ((2,), (3, 4)), (3,))))
  13. print(convert_to_tuple(([1, 2], [(2,), (3, 4)], (3,))))
[(1, 2), ((2,), (3, 4)), (3,)]
[(1, 2), ((2,), (3, 4)), (3,)]
[(1, 2), ((2,), (3, 4)), (3,)]

4.　引数がリストやタプルの場合には、それぞれに対してさらに　関数を適用して、その結果をリスト内包表記でまとめて仮の返り値temp_listにします。このとき、levelを1だけ増やして階層が1つ深くなることを表します。

5.　4.のtemp_listが第1階層の場合はリスト、第2階層以下の場合にはタプルを返り値にします。

8.　引数がリストやタプルでない場合、つまり文字列や数字の場合には、そのままそこに属するリスト、タプルの要素になればよいので、そのまま返り値とします。

さっそく、SymPyで作成したsympy_multiset_partitions_listと比較します。

リストをタプルに変換する

  1. if convert_to_tuple(sympy_multiset_partitions_list) == unlabeled_partitions_list:
  2. print('sympy_multiset_partitions_list == unlabeled_partitions_list')
sympy_multiset_partitions_list == unlabeled_partitions_list

2.　第1階層の各リストについて、タプルに変換します。

multiset_partitions関数で生成したsympy_multiset_partitions_listはconvert_to_tuple関数を作成し、unlabeled_partitions_listと同じ形式に変換します。

集合の分割を比べます。#8の計算結果が正しいことを確認することができました。次の節では、スターリング数を簡単に計算する漸化式を作成し、これをもとに集合の分割を生成する方法を考えます。

ictsr4 All Rights Reserved.