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

全射としてのスターリング数

スターリング数を理解するために重複順列について振り返る

重複順列とスターリング数との関係

重複順列、順列では、定義域の集合Nから値域の集合Kへの写像で、両者とも区別ありのパターンで順列の生成と個数の計算を行いました。ここでは、集合Kについて区別をしない(unlabeled)ケースを考えます。

問題1

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

タクシーをチャーターするというのは贅沢なような気もしますが、部屋割りとは異なり、チーム分けさえ決まれば、どのチームがどのタクシーに乗るかは考えず、誰と誰が同じ組になるかに集中すればよいことになります。写像で表すと重複順列との違いが分かりやすくなります。

重複順列とスターリング数との違い
重複順列とスターリング数との違い

図表1の左の集合Kを区別なしにすると、右のようになり(0,1)(2,3)(4)(5)のように表し、(2,3)((4)(0,1)(5)のようにグループ分けが変わらないものは1つとして考えます。このような区別のないグループ分けを多重分割(partitions of multisets)といい、その分割の個数をスターリング数(stirling number)といいます。

全射の重複順列を計算する

多重分割やスターリング数を直接計算するのは難しいので、いったん全射の重複順列を作成したのち集合Kを区別なし（unlabeled）にします。このためいったんタクシーにA,B,C,Dという仮のラベルを付けたのち、増加写像の考え方を使いラベルを消す方法を考えます。

重複順列はp90110で作成したperm_surjective_func関数を使い、全射の重複順列から抜き出したリストperm_surjective_listを作成します。

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

  1. def perm_surjective_func(k_list, n):
  2. array = [()]
  3. for cnt in range(n):
  4. temp = array
  5. array = []
  6. for elm in temp:
  7. for item in k_list:
  8. if cnt < n-1:
  9. array.append(elm + (item,))
  10. else:
  11. if len(set(new_elm := elm + (item,))) == len(k_list):
  12. array.append(new_elm)
  13. return array
  15. k_list= ['A', 'B', 'C', 'D' ]
  16. n = 6
  17. perm_surjective_list = perm_surjective_func(k_list, n)
  18. print('件数=',len(perm_surjective_list))
  19. pprint.pprint(perm_surjective_list)
件数= 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', 'A', 'D', 'B', 'C'),
・
・

重複順列の数は1560通りになりますが、件数が多いので一部の表示にとどめます。

包除原理による全射の重複順列の計算

全射の重複順列の個数は、包除原理により次の式で計算することができます。

#包除原理による重複順列の

$\sum\limits_{i=0}^{k}(-1)^i\,\ _k C _i(k-i)^n$

さっそく、計算式にもとづき関数を使って計算します。以前作成したperm_surjective_count関数で、n=6,k=4の場合の組み合わせの個数を計算します。

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

  1. def perm_surjective_count(n, k):
  2. sigma = 0
  3. for i in range(k+1):
  4. sigma += ((-1)**i) * comb(k,i) * (k-i)**n
  5. return int(sigma)
  7. perm_surjective_count(n, len(k_list))
1560

1560個と結果は一致します。

集合Kのラベル付き分割をラベル無しの分割に変換

ラベルを消す方法

#1で求めた重複順列はだれがどのタクシーに乗るかというラベル付きの状態ですが、ここからラベル無しの分割に変換します。このことを表したのが図表2です。

集合Kを基準に多重分割を整理
集合Kを基準に多重分割を整理

例えば2つの写像は重複順列では別のものとして数えますが、ラベルなしにすると同じタクシーに乗る顔ぶれとしては同じものとなります。もっとも、写像の図を眺めていても、左右を同じものとして捉えるのは難しいので、視点を変えて図のように誰がどの車に誰が乗るかを表すように変換し、A,B,C,Dに対応する4つのタプルで表します。このことにより、増加写像の考え方を使い$4!=24$通りの重複順列を1つの配列で代表させるようにすることができます。このことを表したのが図表3です。

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

#1で作成したperm_surjective_listを集合K（乗るタクシーごとに）の配列に変換し、さらにはincreasing_func関数で増加写像の考え方を使い集約します。

集合Kから整理するプログラムの作成

#1で作成したperm_surjective_listに対して、集合Kを基準に分類するための関数inverse_funcを作成します。引数として整理前の配列と集合Kのリストk_listを渡します。

重複順列を変数elmとして順次読み込み、elmの要素1つ1つをタクシーごとの分類を示す新たに作成する配列new_elmに変換し、最終的にarrayに追加してまとめます。

このためk_listにあるタクシーにA→0,B→1,C→2,D→3と順番をつけ、要素が4つの配列new_elmの添え字にします。例えば背番号2がAの車に乗るとしたら、new_elm[0]に2を追加します。

集合Kの要素から重複順列を整理する

  1. def inverse_func(k_list, perm_list ):
  2. array = []
  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)].append(num)
  7. array.append(new_elm)
  8. return array
  10. inversed_surjective_list = inverse_func(k_list, perm_surjective_list )
  11. print('件数=',len(inversed_surjective_list))
  12. pprint.pprint(inversed_surjective_list)
件数= 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], [4], [5], [3]],
 [[0, 1, 2], [5], [4], [3]],
・
・
[[5], [3], [4], [0, 1, 2]],
 [[4], [5], [3], [0, 1, 2]],
 [[5], [4], [3], [0, 1, 2]]]

3.　#1で作成した重複順列をelmとして読み込みます。

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

5.　elmのindexつまり出席番号をnum、その学生が乗るタクシーをkとします。

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

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

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

次に、部屋を区別の無いタクシーに変えるため、組み合わせと同じように増加写像の考え方を使い、各組み合わせにp90120作成したincreasing_map_func関数を適用します。このことにより( (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 increasing_map_func(perm_list):
  2. array=[]
  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. array.append(elm)
  9. return array
  11. stirling_surjective_list = increasing_map_func(inversed_surjective_list)
  12. print('件数=',len(stirling_surjective_list))
  13. pprint.pprint(stirling_surjective_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], [2, 5], [3, 4]],
 [[0], [1], [2], [3, 4, 5]]]

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

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

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

スターリング数を計算する関数と一覧表の作成

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

SymPyライブラリのスターリング数を計算する関数

SymPyライブラリのstirling関数を使いスターリング数を計算することができます。stirling関数はNとKの個数を引数とします。

stirling関数によるスターリング数の計算

  1. from sympy.functions.combinatorial.numbers import stirling
  2. stirling(6, 4)
65

1.　stirling関数はSymPyライブラリのfunctions.combinatorial.numbers サブモジュールからimportします。

n=6,k=4で計算すると65と#4と計算結果は一致します。

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

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

スターリング数の計算

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

スターリング数を計算する関数を作成します。関数名は文字通りstirling_countとし、n,kを引数とします。式の通りの全射の重複組み合わせの個数をn!で割ることにより計算することができます。

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

  1. from math import factorial
  2. def stirling_count(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. stirling_count(n,len(k_list))
(65, 65)

7.　戻り値を返すときに、終域の件数(n!)で割ります。割り算の結果は小数点型になるので、int関数で整数に変換します。

SymPyライブラリで多重分割を生成する関数

SymPyライブラリのmultiset_partitions関数により多重分割を生成することができます。引数としてn,kを渡しますが、nは集合Nの要素をリストでも個数でも可能です。

スターリング数のリストを比較

  1. from sympy.utilities.iterables import multiset_partitions
  2. multiset_list = list(multiset_partitions(n, len(k_list)))
True

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

#4で生成した多重分割とmultiset_partitions関数で生成した多重分割を比べます。#4の計算結果が正しいことを確認することができました。

ictsr4 All Rights Reserved.