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

全射ではなく制約なしとした場合の重複分割

制約なしの多重分割の生成

制約なしとする場合の具体例

これまでは、乗らないタクシーは無いという全射の前提で考えてきました。ここでは、空車も許される制約なし(unristricted)にした場合のグループ数ごとの分割を生成し個数を計算します。

#問題2

背番号0から5までの6人について、4つのグループに分けタクシーを使い市内観光を予定しています。グループ分けはタクシーを使うので、どの組み合わせがどこにはいるかは区別しないものとします。この場合に、考えられる組み合わせは何通りになりますか。なお、タクシーは4台すべてを使う必要はないものとします。

写像で表すと、左側が全射の多重分割ですが、制約なしの場合は右の(A, A, B, B, C, C)のようにDが選ばれないような分割も認められます。さらに極論ですが、C,Dが選ばれない分割や、極端な話、(A,A,A,A,A,A)のように、1つに集中するケースも対象になります。

制約なしの多重分割を作成するため、perm_unristricted_funcを使い制約なしの重複順列を作成し、inverse_funcにより終域である集合Kを基準に整理し、結果をリストinversed_unristricted_listに格納します。その後、増加写像の考え方を使い集合Kを区別なしに変換します。

kの要素から分類　部屋割りを見る

def perm_unristricted_func(k_list, n):
array = [()]
for cnt in range(n):
temp = array
array = []
for elm in temp:
for item in k_list:
array.append(elm + (item,))
return array
perm_unristricted_list = perm_unristricted_func(k_list, n)
inversed_unristricted_list = inverse_func(k_list,perm_unristricted_list)
print('件数=',len(inversed_unristricted_list))
pprint.pprint(inversed_unristricted_list)

件数= 4096
[((0, 1, 2, 3, 4, 5), (), (), ()),
 ((0, 1, 2, 3, 4), (5,), (), ()),
 ((0, 1, 2, 3, 4), (), (5,), ()),
 ((0, 1, 2, 3, 4), (), (), (5,)),
 ((0, 1, 2, 3, 5), (4,), (), ()),
 ((0, 1, 2, 3), (4, 5), (), ()),
 ((0, 1, 2, 3), (4,), (5,), ()),
 ((0, 1, 2, 3), (4,), (), (5,)),
 ((0, 1, 2, 3, 5), (), (4,), ()),
 ((0, 1, 2, 3), (5,), (4,), ())]
・
・

12. 重複順列を計算する関数を使い、perm_unristricted_list を作成します。$4^6=4096個になります。

13. 12.の結果に対しinverse_func集集合Kの個数ごとに集計し、リストinversed_unristricted_listを作成します。

重複組み合わせの全射の場合の組み合わせを生成し、 ()はそのタクシーにだれも乗らないことを示しています。

制約なしとした場合の値域の個数による分類

タクシーの台数k=4のうち乗車する台数をxとしたときの重複順列の個数を計算する方法を考えます。ここではx=3のケースを取り上げます。

図の2つの写像は、集合Kの要素4つのうち3つが集合Nに対応付けられており、左はA,B,C、右はA,B,Dに分乗していることを示します。これらは、集合Kをラベル無しの分割に変換するので最終的には1つの分割に集約されますが、現段階では別の重複順列としてカウントされています。k=4のうち3つが選ばれるので、その選び方は$ \ _4 C _3=3$通りあり、その各々についてn=6,k=3の包除原理により計算した全射の重複順列の個数との積になります。

n=6, k=4の重複分割を計算すると$ \ _4 C _3=3\times(\ _3 C _1\times3^6-\ _3 C _2\times2^6+\ _3 C _3\times1^6)=2160$個となります。

一般化すると、k個からx個を選ぶ場合の重複順列は、\ _k C\ _{x}通りに対し、集合Nを選ばれたx個に対する全射の重複順列になるので、包除原理の計算式のkをxに置き換えた数字を掛け合わせることにより個数を計算することができます。このため数式は次の通りです。

包除原理による全射の個数の計算

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

制約なしとした場合の個数

$\ _k C\ _{x}\sum\limits_{i=0}^{x}(-1)^i\cdot _{x} C _i\cdot(x-i)^n $

数式を使って、制約なしの場合の集合Kの個数ごとの個数を計算します。関数名は制約なしの順列を数えるという意味でperm_unstricted_count関数とします。

全射の集合Kの数から分類を計算する関数

def perm_unstricted_count(n, k, x):
sigma = 0
for i in range(x+1):
sigma += ((-1)**i) * comb(x, i) * (x-i)**n
return int(comb(k, x)*sigma)
perm_unstricted_count(6, 4, 3)

4.　$\sum\limits_{i=0}^{x}(-1)^i\cdot _{x} C _i\cdot(x-i)^n $の部分を足し込んでいきます。

5.　返り値を返すところで$\ _k C\ _{x}$を掛けます。このときcomb関数は小数点で結果を返すのでint関数を使い整数に変換します。

上記のプログラムと合致し、正しく計算できることがわかりました。

xの件数ごとに集計する

inversed_unristricted_listをxの個数ごとに集計し、perm_unstricted_count関数が正しく計算できているかを確認します。

全射のkの数から分類を計算する関数

cnt = [0 for _ in range(k+1)]
for elm in inversed_unristricted_list:
cnt[k-elm.count([])]+=1
sigma_func =[]
for i in range(k+1):
sigma_func.append(perm_unstricted_count(n, k, i))
print(f' i | cnt | func |')
print(f'-----+---------+---------|')
for i in range(1, 5):
print(f'{i:^5}|{cnt[i]:>9}|{sigma_func[i]:>9}|')
print(f'-----+---------+---------|')
print(f'sigma|{sum(cnt):>9}|{sum(sigma_func):>9}|')

  i  |   cnt   |  func   |
-----+---------+---------|
  1  |        4|        4|
  2  |      372|      372|
  3  |     2160|     2160|
  4  |     1560|     1560|
-----+---------+---------|
sigma|     4096|     4096|

1.　xごとに個数を集計するため配列cntを定義し0で初期化します。

2.　inversed_unristricted_listから重複順列をelmとして順次読み込みます。

3.　elmの中で空白のリスト[]の数だけ空きが出ているので、kからその数を差し引いた値をxとして、cntを足し上げていきます。

5.　集計する式を使い個数を計算し、その結果を配列sigma_funcに追加します。

重複順列4096通りの内訳の中で、xの個数毎に場合の数が一致していることがわかりました。

増加関数に変換する

前節のxの件数ごとにラベルありの件数を計算したperm_unstricted_count関数をラベル無しに変換する計算式を考え、inversed_unristricted_list に対してincreasing_map_func関数を使い集合Kのラベル無しの分割を生成します。

(3)で求めた制約なし集合Kのラベルなしの場合の数の計算式を求めます。まず、集合Kのラベルが無いということは、k個のうちどの組み合わせを選んでも1つして数えるので式の$\ _k C\ _{x}$の部分は不要となります。また、例えばk=3で (A,B,C)が選ばれた場合、ラベルをとると(A,B,C)の組み合わせが((0, 1), (2, 3), (4, 5))の例でみる通り、k!=3!で割る必要があります。つまりのx=3の重複順列2160通りを$2160/ \ _4 C _3/3!$で割って90通りになります。

このことから、制約なしとした場合のkの個数ごとに個数は次の式で計算することができます。

制約なしでKの区別が無い場合の数

$\sum\limits_{i=0}^{x}(-1)^i\cdot _{x} C _i\cdot(x-i)^n/x!$

この式をよく見るとスターリング数のkをxに置き換えた式になり、同じ計算をしています。そこで、inversed_unristricted_list に対してincreasing_map_func関数を使い集合Kのラベルを外した結果と、上記の数式とstirling数が合っているかを確認します。

増加関数にしてKの区別なしにする

stirling_unristricted_list = increasing_map_func(inversed_unristricted_list)
cnt = [0 for _ in range(k+1)]
for elm in stirling_unristricted_list:
cnt[k-elm.count([])]+=1
sigma_func =[]
sigma_stirling = []
for i in range(k+1):
sigma_func.append(int(perm_unstricted_count(n, k, i)/comb(k,i)/math.factorial(i)))
sigma_stirling.append(int(stirling(n, i)))
print(f' i | cnt | func | stirling|')
print(f'-----+---------+---------|---------|')
for i in range(1, 5):
print(f'{i:^5}|{cnt[i]:>9}|{sigma_func[i]:>9}|{sigma_stirling[i]:>9}|')
print(f'-----+---------+---------|---------|')
print(f'sigma|{sum(cnt):>9}|{sum(sigma_func):>9}|{sum(sigma_stirling):>9}|')

  i  |   cnt   |  func   | stirling|
-----+---------+---------|---------|
  1  |        1|        1|        1|
  2  |       31|       31|       31|
  3  |       90|       90|       90|
  4  |       65|       65|       65|
-----+---------+---------|---------|
sigma|      187|      187|      187|

面白いことに3つの計算結果は一致しました。

制約なしでn=kとした場合のベル数

ベル数の例

これまで、6人を4つのタクシーに分乗することを考えました。これだけでも十分贅沢ですが、数学的に考え得る最大の6台に分乗させることも考えられます。つまり集合Kをチーム数の上限をn（つまり6人を1台から１人1台まで）個まで、ラベル無し、制約無しで分割するような場合の数をベル数といいます。ベル数はSymPyライブラリのbell関数で計算することができるので、perm_unstricted_countでの計算結果を比べます。

n=kの場合の全射の件数

from sympy.functions.combinatorial.numbers import bell
sigma = 0
for i in range(1,n+1):
sigma+=int(perm_unstricted_count(n, n, i)/comb(n,i)/math.factorial(i))
sigma ,bell(n)

(203, 203)

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

2.　nが1から6まで、perm_unstricted_count関数の結果をラベル無しにした値をsigmaに足し込みます。

4.　sigmaとbell関数と結果を比較します。

n=6のベル数は203となり答えは一致しました。これまでの結果を踏まえ、ベル数を計算する関数を作りたくなりますが、perm_unstricted_count 関数をｎ回繰り返しただけでは面白みもなく、nが大きくなると包除原理の増減を何度も繰り返す必要があり計算量が多くなってしまいます。そこで、次節ではベル数を導くための漸化式を検討します。

漸化式の考え方

ベル数については、ある特定の学生が何人の分割に入るかで分類すると興味深い漸化式をつくることができます。

例えば、単独で1つの分割を作る場合、残り5人で分割をつくるのでB(5)になります。次に誰か他の1人と2人で分割を構成する場合、同じ分割に入るのは5人のうち1人で5通り、その各々について他の4人でグループ分けをするのでB(4)通り、つまり5×B(4)になります。さらに他の2人と3人でグループを構成すると、5人のうち2人と同じグループになるので$ \ _5C_2$通り、その各々について残り3人でグループを作るので$ \ _5C_2×B(3)$になります。以下、同様になるのでまとめると次のとおりになります。

ここで、意外と厄介なのが6人で１つだけ分割をつくるときです。残りは0人なので1×B(0)となりますが、これも立派な多重分割で1通りと数える必要があるのでB(0)=1とします。このことも踏まえると次の漸化式になります。

ベル数の漸化式

$B(0)=1$

$B(n) = \sum_\limits{k=0}^{n-1}\ _{n-1} C_k B(k)$

漸化式を使ってベル数を計算する

漸化式を使い0から10までのベル数を計算します。ここではリストbell_arrayでは０から変数maxまでの順次ベル数を追加します。計算結果はB(0)からB(10)までbell_arrayというリストに累積します。

ベル数を漸化式で計算する

max = 10
bell_array = [1]
for n in range(1, max+1):
next_bell = 0
for k in range(n):
next_bell += int(comb(n-1,k) * bell_array[k])
bell_array.append(next_bell)
print(f' n | myfunc | sympy |')
print(f'---+---------+---------+')
for i, num in enumerate(bell_array):
print(f'{i:^3}|{num:>9}|{int(bell(i)):>9}|')

n |  myfunc |  sympy  |
---+---------+---------+
 0 |        1|        1|
 1 |        1|        1|
 2 |        2|        2|
 3 |        5|        5|
 4 |       15|       15|
 5 |       52|       52|
 6 |      203|      203|
 7 |      877|      877|
 8 |     4140|     4140|
 9 |    21147|    21147|
10 |   115975|   115975|

2.　0のベル数は1なので、bell_arrayに初期値1を設定します。

4.　nに格納した次のベル数をnext_bellで計算するため、漸化式に従いint(comb(n-1,k)*bell_array[k]) を次々に足し込みます。

0から6までのベル数をリストにすることができました。だいぶすっきりしてきましたが、nのベル数を計算するために、n-1以下のベル数をいちいち参照する必要があるので、もう一歩工夫が必要な気がします。