統計数理（2009）
第 57 巻 第 1 号 159–178
c 2009 統計数理研究所
［原著論文］

変形バケットソートに現れる離散型確率分布とEulerian 数
1

土屋 高宏 ・中村 永友

2（受付 2008 年 6 月 23 日；改訂 2008 年 11 月 25 日）

要旨

バケットソートを変形したソーティング・アルゴリズムに現れる離散型確率分布を導出した．
このソーティングの過程で，ある規則性を持つ数が現れる．その数に関係する場合の数を漸化式で表現し，確率分布を与えた．また，確率分布の導出の際に現れる漸化式は Eulerian 数と呼ばれる数列を構成するが，導出過程において Eulerian 数との間に微妙な差異があることを示した．さらに数値実験によりこれを確認した．確率分布についてモーメントなどの基本統計量を求め，正規分布との比較，確率分布近似の精密化を行った．さらに，求めた確率分布と連続型
一様分布との関係について言及した．本研究が Eulerian 数を導出する新たな例となることを示した．
キーワード： Eulerian 分布，一様分布，キュムラント，漸近展開，正規分布．
1.

はじめに

バケットソートと呼ばれる計算時間が O(n) で高速に行えるソーティング・アルゴリズムがある．このアルゴリズムは並べ替えをしたいデータのとりうる値が k 通りあるとき，あらかじめ k 個のバケツを用意しておき，あるいは動的にバケツを増やしていきながら，各々の数字と対応するバケツにデータを入れていくアルゴリズムである．本論文は，n 個の連続した数字があり，あらかじめ用意するバケツ数が決められておらず，その数がデータの初期状態に依存する「変形バケットソート」を考え，そのバケツ数がどのような確率分布になるのかを考察する．
また，確率分布の導出の際に現れる漸化式は Eulerian 数（Euler, 1755; Graham et al., 1994）と呼ばれる数列を構成することを示す．また，導出過程において両者間に微妙な差異があるので，これについて言及する．ここで Eulerian 数とは，Eulerian numbers が原語であり，定ベキ化
係数，ベキ係数などの訳語もあるようであるが，定まった言い方がないため，この呼び方を用いる．
2 節では問題設定と変形バケットソートのアルゴリズムについて述べ， 3 節ではバケツ数の確率分布を導出するための考え方を示して確率分布の漸化式表現を与える．4 節ではバケツ数と Eulerian 数の関係について触れ，導出過程において両者間に微妙な差異があることに言及する．この確率分布について， 5 節ではモーメントを求め， 6 節では一様分布との関係について述べる．また，7 節では正規分布との比較と確率分布近似の精密化を行う．
1
2

城西大学 理学部：〒350–0295 埼玉県坂戸市けやき台 1–1
札幌学院大学 経済学部・総合教育センター：〒069–8555 北海道江別市文京台 11 番地

統計数理

160

2.

第 57 巻

第 1 号 2009

問題設定

簡単のため数字の書かれたカードを例に説明する．1 から n までの数字が書かれたよくシャッフルされている n 枚のカードが手元にある．カードを小さい順に並べ替えをするために，次の手順で行う．
（1）一番上のカードの数字が k のとき，数字 k ＋ 1 のカードがすでにテーブルに置かれていたらその上にカード k を載せる．
（2）もしテーブル上に k ＋ 1 が無ければ，どのカードの上にも載せず，テーブルにカード k
を置く．
（3）
（1）と
（2）を手元のカードが無くなるまで続ける．
（4）テーブルの上に m 個のカードの束ができる (1 ≤ m ≤ n)．
（5）各束の一番上のカードの数字を小さい順にまとめれば，並べ替えが完了する．
手元の n 枚のカードが 1 から n までの重複がないと通常のバケットソートとの相違点は，（1）
いう前提がある点，（2）バケツ数がいくつになるかわからないという点，（3）
数字 k ＋ 1 がすでに置かれていれば k を k ＋ 1 に載せる点，である．この方法による並べ替えでは，手順の（4）の段階で最大で n 個，最小で 1 個のカードの束ができることになる．このとき，カードをテーブ
ルに置き終わった時点におけるカードの束の数に関する確率分布がどのようなものになるのかを議論する．「束の数」とはバケットソートにおけるバケツ数に対応する．
数字（カード）の並びが与えられたとき，変形バケットソートにおけるバケツ数（束の数）
は以下のアルゴリズムで直接求めることができる．
アルゴリズム．n 個の数字（カード）の並びを ν1 , ν2 , . . . , νj , . . . , νn と表す．このとき，k < j(k = 1, . . . , j − 1; j = 2, . . . , n) でかつ νk = νj +1 を満足する総数を
m = #{k | νk = νj + 1, k < j}とすると，n − m がバケツ数である．ここで，記号 #{·} は，集合 {·} の要素の数を表す．
（カード）の並びを ν1 , ν2 , . . . , νj , . . . , νn として，バケツを n 個用意しておく．
解説．n 個の数字
j を固定して，νk = νj + 1 を満足する νk (k = 1, . . . , j − 1) があったときには，νk の入ったバケツに νj を入れることになるので，バケツ数が 1 個減ることになる．これを j = 2, . . . , n で動かしたとき，νk = νj + 1 を満足する総数が m 個あったとする．このとき m 個分だけバケツ数が減ることになるので，最終的なバケツ数は n − m となる．
3.

束の数の確率分布

本節では関心のある確率分布の導出の基本的な考え方を示し，それに基づいて漸化式により確率分布を表現する．
3.1 束の数
カードの枚数 n が 2，3，4 の場合にカードの出方と束の数を列挙すると以下のようになる．
n=2:
n=3:



(1, 2) · · · 2(2, 1) · · · 1(1, 2, 3) · · · 3(1, 3, 2) · · · 2(2, 1, 3) · · · 2(2, 3, 1) · · · 2(3, 1, 2) · · · 2(3, 2, 1) · · · 1

変形バケットソートに現れる離散型確率分布と Eulerian 数

n=4:

⎧
⎪
(1, 2, 3, 4) · · · 4
⎪
⎪
⎪
⎪
⎪
⎪
(1, 2, 4, 3) · · · 3⎪⎪⎪⎪⎨ (1, 3, 2, 4) · · · 3⎪(1, 3, 4, 2) · · · 3
⎪
⎪
⎪
⎪
⎪
⎪
(1, 4, 2, 3) · · · 3⎪⎪⎪⎪⎩ (1, 4, 3, 2) · · · 2(2, 1, 3, 4) · · · 3(3, 1, 2, 4) · · · 3(4, 1, 2, 3) · · · 3(2, 1, 4, 3) · · · 2(3, 1, 4, 2) · · · 3(4, 1, 3, 2) · · · 2(2, 3, 1, 4) · · · 3(3, 2, 1, 4) · · · 2(4, 2, 1, 3) · · · 2(2, 3, 4, 1) · · · 3(3, 2, 4, 1) · · · 2(4, 2, 3, 1) · · · 2(2, 4, 1, 3) · · · 2(3, 4, 1, 2) · · · 3(4, 3, 1, 2) · · · 2(2, 4, 3, 1) · · · 2(3, 4, 2, 1) · · · 2(4, 3, 2, 1) · · · 1

161

カッコ内はカードの並びを，その右に束の数を示す．なお，カードはカッコ内の左から順にカードが現れるものとする．例えば，n = 3 のとき束の数が 2 になる場合は 4 通りあり，n = 4 のとき束の数が 2 または 3 になる場合はそれぞれ 11 通りあることがわかる．ちなみに，n 枚のカードがあるとき，カードの並べ方の総数は n! 通りある．
次に一般の場合について考える． n 枚のカードの束の数を求めるためには， n − 1 枚のカー
ドの並びにおける束の数を元に検証していく．例えばカードが 4 枚，つまり n = 4 の束の数は，3 枚のカードの並びに 4 のカードを 1 枚追加したときの束の数の変化を考えればよい． 3 枚のカードの並びの 1 つに (1, 3, 2) があるが，これは束の数が 2 になる．これにカードを 1 枚追加
したとき，4 のカードは (a 1 b 3 c 2 d) の 4 つの a, b, c, d のいずれかの位置に入る．このとき，4 が 3 の左側（a または b）に入るときは束の数は 2 のままで変化はない．一方，右側（c または d）に入るとき，束の数は 1 増加して 3 になる．
以上の考え方から一般に n 枚のカードにおける束の数に関して，次のことが言える． n − 1
枚のカードに数字 n を追加したとき，数字 n が数字 n − 1 の左側に入るとき束の数は変わらず，
右側に入るとき束の数は 1 増加する．したがって，n 枚のカードにおける束の数が i になる総数 Mn (i) は，n − 1 枚のカードにおいて，束の数が変わらない場合の数と，束の数が 1 増加して i になる場合の数の和として得られる．実際，次の漸化式が成り立つ．
(n ≥ 1),
Mn (1) = Mn (n) = 1(3.1)
Mn (i) = iMn−1 (i) + {n − (i − 1)}Mn−1 (i − 1)(n ≥ 3, i = 2, . . . , n − 1).
表 1 は，n = 12 までの束の数 Mn (i) の数表である（n = 13, 14, 15 は付録 A の表 4 ∼ 5 の最下行を参照）．以後，前後の文脈から明らかなときには，「束の数による分類の場合の数」を単に「束の数」と呼ぶことにする．
漸化式（3.1）は次節で述べる Eulerian 数を求める式と数学的に同値であり，これは 3.3 節で表 1. 束の数 Mn (i) (n = 1, 2, . . ., 12)．

統計数理

162

第 57 巻

第 1 号 2009示す．
3.2 Eulerian 数
Eulerian 数は Euler（1755）
によって提唱され，次のように定義されている（Knuth, 1998 など）
．


k

n+1
(3.2)
A(n, k) =(−1)j(n ≥ 0, k ≥ 0).
(k + 1 − j)njj=0

また，漸化式による表現は，
A(n, k) = (n − k)A(n − 1, k − 1) + (k + 1)A(n − 1, k),(3.3)
A(n, 0) = 1 (n ≥ 1), A(n, k) = 0 (k ≥ n ≥ 1)(n ≥ 1, n − 1 ≥ k ≥ 1),である．この漸化式の証明は右辺に（3.2）式を代入することにより示すことができる．
Eulerian 数は一般的に知られている Euler 数とは別のものである．これは主に代数学や離散数学の分野で研究されており，以下のように第 2 種スターリング数 S(n, k) やベルヌーイ数 Bnとの関係が知られている（例えば，Graham et al., 1994）．


n−k

n−j
n−j−k
A(n, k) =S(n, j)(−1)
j!
,
k
j=1
n−1

k=0

(−1)k A(n, k) =2n+1 (2n+1 − 1)Bn+1.
n+1

が時系列に関する差の符号
Eulerian 数に対する確率分布の適用例としては， Kimber（1987）
検定や壺のモデルに適用している．
3.3 束の数と Eulerian 数の関係
束の数 Mn (i) と A(n, k) の 2 つの数の間には，k = i − 1 とおくと，Mn (i) = A(n, i − 1) という関係が成立している．例えば，M6 (1) = A(6, 0), M6 (2) = A(6, 1), M6 (3) = A(6, 2), . . . . である．
このように束の数と Eulerian 数が同値であることを通して，（3.1）
式が正しいことを以下に示す．

命題 3.1. 束の数と Eulerian 数の間に Mn (i) = A(n, i − 1)(i = 1, 2, . . . , n) が成り立つ．
証明. 命題 3.1 を示すためには，{1, 2, . . . , n} の置換全体の集合を，束の数によって分類する
ときの場合の数が，Eulerian 数に等しいことを証明すればよい．
自然数 n に対し，Ωn = {1, 2, . . . , n} とし，Ωn 上の置換全体の集合を Sn とする．


1
2
···
n
τ=
∈ Sn
τ (1) τ (2) · · · τ (n)に対して，a(τ ) = #{j | τ (j) < τ (j + 1), j = 1, 2, . . . , n − 1} を τ の上昇点の数とする．Graham etでは，Eulerian 数を
al.（1994）
A(n, k) = #{τ ∈ Sn | a(τ ) = k}と定義している．τ ∈ Sn に対して，τ の逆置換は，
τ (1) τ (2) · · ·
−1
τ =
1
2
···

τ (n)
n



変形バケットソートに現れる離散型確率分布と Eulerian 数

163となる．このとき，τ ∈ Sn に対して，τ (1), τ (2), . . . , τ (n) の順に Ωn の要素が現れるときの束の数を b(τ ) と表すと，b(τ ) = 1 + a(τ −1 )(3.4)

が成り立つ．これは以下のように証明することができる．j = 1, 2, . . . , n − 1 について，命題 Pj
を「τ (α) = j ，τ (β) = j + 1 となる α と β について，α < β 」と定義する．Ωk−1 上の置換に k(k = 2, 3, . . . , n) を挿入して Ωk 上の置換を構成する方法を繰り返して τ ∈ Sn を作ることを考える．このとき，最初に数字 1 のみが与えられると束の数は 1 である．j を増やしていき，各 Pjが成立するたびに束の数が 1 増加するので，b(τ ) = 1 +n−1
I(Pj )

j=1と表すことができる．ここで，I(Pj ) =−1

1(Pj が真)0(Pj が偽)−1である．命題 Pj は τ を用いて表すと，「τ (j) < τ −1 (j + 1)」となるので，b(τ ) = 1 + #{j | τ −1
−1(j) < τ (j + 1)} となり，
（3.4）式が成り立つ．
Sn から Sn への写像 ϕ : τ → τ −1 は全単射である．したがって，Sn の任意の部分集合 C を写像 ϕ で変換すると，#C = #ϕ(C) であるから，Mn (i) = #{τ ∈ Sn | b(τ ) = i} = #{τ ∈ Sn | b(τ −1 ) = i}となる．ここで，（3.4）
式より {τ ∈ Sn | b(τ −1 ) = i} = {τ ∈ Sn | a(τ ) = i − 1} となる．また，Eulerian数の定義から #{τ ∈ Sn | a(τ ) = i − 1} = A(n, i − 1) となるので，Mn (i) = A(n, i − 1) が成り立
つ．
命題 3.1 と Eulerian 数の漸化式（3.3）から，漸化式（3.1）はすべての自然数 n について成り立
つことがわかる．
3.4 確率分布
n 枚のカードが i 束になる確率 Pn (i) は，Pn (i) = Mn (i)/n! で定義され，次の漸化式で与えられる．
束の数の分布（Eulerian 分布）⎧
1⎪⎨ Pn (1) = Pn (n) =(n ≥ 1),n!
(3.5)⎪⎩ Pn (i) = i Pn−1 (i) + n − (i − 1) Pn−1 (i − 1)(n ≥ 3, i = 2, . . . , n − 1).
nn

図 1 に束の数 Mn (i) の確率分布 Pn (i) の例を示す．
4.

束の数と Eulerian 数のずれ

表 1 で示した束の数に関する漸化式により与えられる数列は Eulerian 数を構成する．本節ではこの数列と束の数との関係を考察する．
Eulerian 数とは 1 から n の順列に対して，上昇点の数（あるいは，それと同値な増加列・減少列の数）を表すものである．また，Stanley（1997）
は下降点の数による分類として Eulerian 数を定義している．上昇点の数とは n 個の数の並びにおいて，隣り合う 2 数の大小関係が < と

164

統計数理

図 1.

第 57 巻

第 1 号 2009

Pn (i) (n = 5, 10, 20, 30) の確率分布．

なる個数 m を表している．例えば，（1, 2, 3, 4）の並びは 1 < 2 < 3 < 4 であるから，上昇点の数
は 3 となり，（3, 2, 4, 1）の並びは 3 > 2 < 4 > 1 であるから，上昇点の数は 1 となる．
一方，束の数は m を用いて m + 1 と表され，束の数 = 上昇点の数 + 1の関係が成り立っているようにみえる．しかしながら，両者の間に食い違いが生じる場合があり，すべての並びについてこの等式が成り立つとは限らない．例えば，表 2 のように，番号 1と 2 の数字の並びでは両者が一致するが，番号 3 と 4 では ずれ が生じていることがわかる．
n = 4 のときにずれが生じるのは，この 2 組のみであり，束の数が 2 と 3 になる場合の総数
は「 上昇点の数 +1」と等しくなる．したがって，束の数は Eulerian 数を構成していることがわかる．
表 3 は n = 4, 5 のときの「束の数」と「 上昇点の数 +1」の一致と不一致を表すクロス表である．n = 4 の非対角要素，2 行 3 列と 3 行 2 列の 1 はそれぞれ「束の数」と「 上昇点の数+1」が 2 と 3 になるずれと 3 と 2 になるずれが 1 組ずつあることを示している．また，n = 5
表 2.

束の数と上昇点の数．

変形バケットソートに現れる離散型確率分布と Eulerian 数
表 3.

165

n = 4, 5 の「束の数」と「 上昇点の数 +1」の一致と不一致（ずれ）．

のときは 2 と 3 および 3 と 2 になるずれがそれぞれ 6 組あり，さらに 3 と 4 および 4 と 3 にな
るずれもそれぞれ 6 組あることがわかる．結果的に行と列の合計が等しくなり，ずれが相殺さ
れ，束の数が Eulerian 数になっている．こうしたずれが生じる理論的根拠は以下のように説明することができる．
τ ∈ Sn が与えられたとき，「 上昇点の数 +1」は 3.3 節の記号を用いると，a(τ ) + 1 であるが，束の数は a(τ ) + 1 ではなく，a(τ −1 ) + 1 であるから，両者の値が異なるような τ については，ずれが生じる．表 2 の番号 3 と 4 の置換は互いに逆置換になっていて，しかも a(τ ) = a(τ −1 )となっているから，ずれが生じることになる．逆置換を与える変換 ϕ は Sn 上の全単射であるから，対応する場合の数は同じになる．
数学的には束の数と Eulerian 数は同値であるので，ずれの議論は数学的に実りがないと思われるが，導出過程の相違を指摘したので，そのずれの様子を調べた．n = 6, . . . , 15 の一致・不一致の様子を付録 A に示す．表から，ずれの様子が主対角線対称および反主対角線対称であることがわかる．また，ずれの程度は主対角線要素の総和を n! で割った一致率で測ることができる（付録 A の表 5）．表を見る限り，n が大きくなるほど一致率は単調に減少している．導出過程に相違があるが，先の命題で束の数が Eulerian 数になることが示されたことは，新たなEulerian 数を導く一例を示したことになる．
以上の議論より，（3.5）式で与えられる束の数の確率分布 Pn (i) をこれ以降 Eulerian 分布と呼ぶことにする．Eulerian 分布は Eulerian 数を用いて，Pn (i) = A(n, i − 1)/n! と表される．
5.

確率分布のモーメント

本節では確率分布の期待値や分散をはじめとするモーメントを与える． Eulerian 分布の平均
まわりのモーメントの漸化式表現は Mann（1945）により導出されているが，ここでは Eulerian
分布に対して原点まわりのモーメントの漸化式表現を与える．この漸化式から期待値や分散の一般項を導出する．
定理 5.1. 確率分布
（3.5）にしたがう確率変数 X の r 次モーメント En (X r ) は，以下の漸化式で与えられる．

 


r

n+1 r
1 r+1
(5.1)
E1 (X r ) = 1, En (X r ) =−En−1 (X k ) (n ≥ 2).
n
n
k
k
k=0
証明. モーメントの定義より，n = 1 のとき E1 (X r ) = 1r P1 (1) = 1 である．n ≥ 3 のとき（3.5）
式の i を x に置き換え，そして xr をかけて x(x = 1, 2, . . . , n) について和をとると，次式が成り

統計数理

166

第 57 巻

第 1 号 2009

立つ．
En (X r ) =

(5.2)

1
1
En−1 (X r+1 ) + En−1 [(n − X)(X + 1)r ] .
nn（5.2）式の右辺第 2 項は1
En−1 [(n − X)(X + 1)r ]n1= En−1 [(X + 1)r ] − En−1 [X(X + 1)r ]
n
 
 
r
r−1

r
1 r
1
k
En−1 (X ) −En−1 (X k+1 ) − En−1 (X r+1 )
=
n
n
k
k
k=0
k=0
と変形できるので，（5.2）式は
 
 
r−1
r

r
1 r
(5.3) En (X r ) =En−1 (X k ) −En−1 (X k+1 )
n
k
k
k=0
k=0
 

 
 r 
r
r
 r

r
1  r+1kk
En−1 (X ) −En−1 (X ) −En−1 (X k )
=
n
k
k
k
k=0
k=0
k=0

 


r

n+1 r
1 r+1
=
−
En−1 (X k )
n
n
k
k
k=0
となる．n = 2 のとき（3.5）式から
1
E2 (X r ) = 1r P2 (1) + 2r P (2) = (1 + 2r )
2

(5.4)

を得る．ここで，（5.3）式で n = 2 とおいた式は（5.4）式と一致することがわかる．したがって（5.3）式は n ≥ 2 のときにも成り立つ．
この証明は Mann（1945）に基づく方法である．また，土屋・中村（2007）には離散型確率分布の期待値の定義式による証明がある．
次に，Eulerian 分布に対する期待値と分散は次の定理により得られる．
（3.5）にしたがう確率変数 X の期待値 En (X) と分散 Vn (X) は，それぞれ定理 5.2. 確率分布n+1
(n ≥ 1),
2
n+1
Vn (X) =(n ≥ 2)
12

En (X) =V1 (X) = 0,で与えられる．
証明.
(5.5)

r = 1 のとき
（5.1）式はE1 (X) = 1,

En (X) =n−1
En−1 (X) + 1
nとなり，さらに Bn (X) = nEn (X) とおくと
（5.5）式は(5.6)

B1 (X) = E1 (X) = 1,

Bn (X) = Bn−1 (X) + nとなる．
（5.6）式から n ≥ 2 に対して，(5.7)

Bn (X) = B1 (X) +

n−1

k=1

k=

n(n + 1)
2

変形バケットソートに現れる離散型確率分布と Eulerian 数

167

を得る．
（5.7）式は n = 1 のときも成り立つ．よって，期待値 En (X) の一般項は
En (X) =Bn (X) n + 1=
n2(n ≥ 1)で与えられる．
r = 2 のとき
（5.1）式はE1 (X 2 ) = 1,

En (X 2 ) =n−2
2n − 1
En−1 (X 2 ) +En−1 (X) + 1nnと表される．よって，分散 Vn (X) に関して
V1 (X) = E1 (X 2 ) − {E1 (X)}2 = 0,
Vn (X) = En (X 2 ) − {En (X)}2 =

(5.8)

n−2
1
Vn−1 (X) +
n
4

を得る．
（5.8）式がVn (X) − (an + b) =

(5.9)

n−2
[Vn−1 (X) − {a(n − 1) + b}]nを満たすように定数 a，b の値を求めると， a = b = 1/12 である．ここで，Cn (X) = Vn (X) −（5.9）式は Cn (X) = {(n − 2)/n}Cn−1 (X) となるので， n ≥ 2 に対して{(n + 1)/12} とおくと，Cn (X) = 0 を得る．したがって，Vn (X) =n+1
12(n ≥ 2)となる．
6.

Eulerian 分布と一様分布

Eulerian 分布は一様分布と以下のような関連がある． Hensley（1982）
は， Eulerian 分布が次の確率と等しいことをラプラス変換を用いて示した．
 n


(6.1)
Pn (k + 1) = Pr
Xj ∈ [k, k + 1] .
j=1

ここで Xj は (0, 1) 上の一様分布にしたがう確率変数である．ここでは，このことを区間 [k, k + 1]
で直接積分することにより示す．
一様分布にしたがう確率変数の和の確率密度関数は次のように表される（Mood et al., 1974,．
p.238）
 
 
n−1

n
n
1
(6.2)
f (s) =(s − 1)n−1 +(s − 2)n−1
sn−1 −
(n
−
1)!
1
2
j=0
 
n
− · · · + (−1)j(s − j)n−1 I(j,j+1] (s).
j
ただし，


IR (s) =である．

1(s ∈ R)

0

(s ∈
/ R)

統計数理

168

第 57 巻

第 1 号 2009

（6.2）式の積分は，区間 [k, k + 1] に対応する区間 (j, j + 1] のみで考えればよいから，和の記
号を除いて j を k で置きかえると以下のようになる．
k+1
kf (s)ds
 
 
 
k+1
n
n
n
(s − 1)n +(s − 2)n − · · · + (−1)k(s − k)n
1
2
k
k
 
 
 
n
n
n
1
=
kn +
(k − 1)n − · · · + (−1)k(k + 1)n −
n!
1
2
k
 
 
 
n
n
n
(k − 1)n +(k − 2)n − · · · + (−1)k− kn −
1
2
k






n+1
n+1
n+1
1
n
n
n
k
k +
(k − 1) − · · · + (−1)=(k + 1) −
n!
1
2
k


k
n+1
1
1 
(−1)j
=
(k + 1 − j)n = A(n, k) = Pn (k + 1).
n! j=0
n!
j
=

1
sn −
n!

したがって，（6.1）式が証明できた．
また，同時に次の別の結果を示すことができる．
（6.2）
式において n の代わりに n + 1 とおき，s = k + 1 とおくと，






n+1
n+1
n+1
1
n
n
n
k
f (s) =(s − 1) +(s − 2) − · · · + (−1)
s −(s − k)n
n!
1
2
k


k
n+1
1 
j
(−1)
=
(k + 1 − j)n = Pn (k + 1)
n! j=0jとなり，Eulerian 分布は n + 1 個の連続型一様確率変数の和の確率密度関数と一致する．つまりこれは，離散型確率分布の確率と連続型確率分布の確率密度という性質の違う 2 つの一致を意味し，大変興味深い結果である．
7.

確率分布の近似

Eulerian 分布の正規分布への近似については， Mann（1945）や Hensley（1982）によって研究されている．Mann は Eulerian 数 A(n, k) のモーメントに関する漸化式から，A(n, k) を基準化した変数が n → ∞ のとき，標準正規分布のモーメントに収束することを示した．また，Hensley
は Eulerian 分布と一様確率変数の和の分布の関係を利用して，中心極限定理により正規分布への近似を行った．ここでは，高次のキュムラントを使って確率分布近似の精密化を行う．
7.1 近似の精密化
モーメントに関する漸化式から，より高次のキュムラントを導出することができる．6 次ま(n)でのキュムラント κr (r = 1, 2, . . . , 6) は次の定理で与えられる．

変形バケットソートに現れる離散型確率分布と Eulerian 数

169

定理 7.1. 確率分布
（3.5）の 6 次までのキュムラント κ(n)(r = 1, 2, . . . , 6) は，
rn+1(n ≥ 1),2n+1(n)(n ≥ 4),κ4 = −
120
(n)

κ1 =

(n)

κ2 =

n+1
(n ≥ 2),
12

(n)

(n)

κ3 = 0 (n ≥ 1),(n)

κ5 = 0 (n ≥ 1),

κ6 =

n+1
(n ≥ 6)
252

である．
証明.

付録 C 参照．

確率分布（3.5）は対称分布であるから，平均を除く奇数次のキュムラントは 0 となる．また，
1 次のキュムラントは (n + 1)(1 + B1 )，偶数次のキュムラントは (n + 1)Br /r (r ≤ 6) として表すことができる．8 次以上のキュムラントは，その式が非常に複雑であるためここでは示さない．ここで，Br は次の関係式を満たすベルヌーイ数である．
tet − 1

=

∞

Br r
t .
r!
r=0

r = 6 までのベルヌーイ数 Br は以下の通りである．
r

1

2

Br

1
−
2

1
6

3

4

5

6

0

1
−
30

0

1
42

これは，n + 1 個の一様分布にしたがう確率変数の和の分布のキュムラントと酷似している．
違いは，一様分布は任意の自然数 n について成り立つが，Eulerian 分布は n がキュムラント
の次数以上で成り立つという点である．n がキュムラントの次数より小さいときは一様分布のキュムラントと異なる．
したがって，確率分布（3.5）の r 次のキュムラント κ(n)は，（0,1）上の一様分布にした
r （r ≤ 6）
がう n + 1 個の互いに独立で同一な確率変数の和のキュムラントと一致する．
次に，確率分布の精密化は以下の定理により得られる．
定理 7.2. Tn を確率分布（3.5）にしたがう確率変数とするとき，√
N (Tn − N/2)
N Tn
1
√
=
(7.1)
Xn =
−
12 N
2
N/ 12
とおくと，確率関数 f (x) = Pr(Xn = x) (n ≥ 6) に対して，次のような 1/N 2 の項までの漸近展開式を得る．

1
1
f (x) ≈ φ(x) 1 − a1 (x) + 2 a2 (x) .
NN
ここで，x のとりうる値は Tn = i (i = 1, 2, . . . , n) を（7.1）式に代入した


1
N
i
−
(7.2)
x=
12 N
2
で与えられる．また，N = n + 1，φ(x) は標準正規密度関数，係数 a1 (x)，a2 (x) は，
1
1
1H4 (x),
H6 (x) +H8 (x)
a2 (x) =20
105
800である．ただし，Hj (x) (j = 1, 2, . . . ) は j 次のエルミート多項式である．
√
証明. N = n + 1, µ = N/2, σ = N/ 12 とすると，√
N (Tn − µ)Xn =σa1 (x) =

統計数理

170

第 57 巻

第 1 号 2009

は漸近的に標準正規分布 N (0, 1) に収束する．このとき，Xn の特性関数は次の式で与えられる．
 √
 
√

Nµ
N Tn
ψXn (t) = E(eitXn ) = exp −(7.3)
it E exp
it
σ
σ

√

2
 √
√
1
Nµ
N
N
(n)
(n)
it expit κ1 +
it κ2
= exp −
σ
σ
2
σ
3
4
√
√
1
1
N
N
(n)
(n)
it κ3 +
it κ4 + · · · .
+
3!
σ
4!
σ
(n)
2
ここで，κ(n)は r 次のキュムラントである．さらに（7.3）
式に κ(n)r
1 = µ = N/2，κ2 = σ = N/12，
(n)
(n)
(n)
(n)
κ3 = 0，κ4 = −N/120，κ5 = 0，κ6 = N/252 を代入すると，
 2 





1
1
t
1
1
1
(7.4)
ψXn (t) = exp −1−(it)4 + 2(it)6 +
(it)8 + · · ·
2
N 20
N
105
800

となる．いま，2
1(it)r e−t /2 e−itx dx = Hr (x)φ(x)2πを用いて（7.4）式を反転すると，Pr(Xn = x) に対して 1/N 2 の項までの式





1
1
1
1
1
H4 (x) + 2
H6 (x) +H8 (x)
f (x) ≈ φ(x) 1 −
N 20
N
105
800

を得る．
7.2 数値近似（3.5）の近似誤差として，真の値の差 E と log10 |E| を示す（n = 6, 10, 20）．次
図 2 に確率分布
式の第 1 項（正規近似），1/N までの項，1/N 2 までの項の近似式がそれぞれ R，Z1，Z2 である．

(i−N/2)2
1
1
1
−
N/6
e
(7.5)
Pr (Tn = i) ≈ 
H4 (x)
1−
N 20
2πN/12
1
1
1
H6 (x) +H8 (x) .
+ 2
N
105
800

ここで，i = 1, 2, . . . , n であり，x は（7.2）式で与えられる．図中の R1 は（6.1）式（Hensley, 1982）に基づく正規分布近似である．
図 2 の左図を見ると，近似式 R でも誤差が ±0.003 の範囲内となっており，十分良い近似値を与えているが，近似式 Z1，Z2 はすべての値について近似精度が改善されていることがわかる．図 2 の右図では，提案方法が相対的に近似精度がよいことがわかる．
なお，離散型確率分布で行うような補正，例えば二項分布の補正のようなことはしていないが，提案した確率分布の近似の方法で十分な精度が得られたことを付け加えておく．
8.

おわりに

本論文は，変形バケットソートにおけるバケツ数の離散型確率分布を導出し，その漸化式表現を与えた．漸化式から導かれる数列は Eulerian 数を構成することを示した．また，束の数と Eulerian 数の導出過程にずれがあることを示し，数値的に検証した．以上より，束の数はEulerian 数を導出する新しい例を与えたことになる．また，Eulerian 分布の数学的特性や一様分布との関係を明らかにし，一様分布にしたがう確率変数の和の分布を利用して，従来の正規分布近似よりも精密な分布近似を行った．

変形バケットソートに現れる離散型確率分布と Eulerian 数

図 2.

171

正規分布近似と近似の精密化．左は n = 6, 10, 20 の確率分布の値を（7.5）式の第 1 項（正規近似），1/N までの項，1/N 2 までの項を用いて計算したときの誤差 E = 近似値− 真の値（それぞれ，R，Z1，Z2）と（6.1）式の正規分布近似による誤差 E ．右は誤差
E1 = log10 | 近似値 − 真の値 | である（縦軸は対数軸）．

今後の課題としては，（1）束の数に関する漸化式の直接的な証明，（2）n 枚の「ずれ」を表す各
クロス表の要素第 k 表の i 行 j 列が，k に関して別の意味のある数列を構成する可能性がある
ので，できる限り明らかにすること，などが挙げられる．
謝

辞

本研究を完成するにあたって，統計数理研究所 平野勝臣名誉教授と九州大学数理学研究院
小西貞則教授に貴重なコメントをいただきました．匿名査読者には命題 3.1 の証明を示していただくと同時に，有益な助言と不備の指摘をいただき，本論文の質を高めることができました．
ここに記して御礼を申し上げます．また，本研究は中村が 2008 年度札幌学院大学長期国内研究員の期間内に完成させました．
付 録
A. 「束の数」と「 上昇点の数 +1」の一致と不一致{1, 2, . . . , n} の n 個の数字
（カード）のすべての並べ替え（組み合わせ）{τ1 , τ2 , . . . , τn } から得られる「束の数」と「 上昇点の数 +1」の一致と不一致の総数，および一致率をまとめたものを
表 4 と表 5 にそれぞれ示す．任意の順列が与えられたときに，これらが一致するのか否かを一

172

統計数理

第 57 巻

第 1 号 2009

瞬で判断できる理論がないため，すべてについて確認した．行和および列和は Eulerian 数になる．また，主対角線に対して対称なので，行と列は「束の数」と「 上昇点の数 +1」のどちら
でも良い．
B.

束の数と Eulerian 数の計算時間についてここに示す表 6 は付録 A の表 4 を計算するときにかかった時間である．本論文で計算に用いた計算機の仕様は CPU: Intel(R) Core(TM) 2 Duo 3.00GHz である．確認した束の数（Eulerian数）のオーダー n，基礎となる計算量の n!，および実計算時間を示す． n = 15 までは計算確認済みであるが，明らかに組み合わせ爆発している n = 16 は予測時間を示した．

表 4.

束の数と上昇点の数の関係．

表 4．つづき

変形バケットソートに現れる離散型確率分布と Eulerian 数173

表 4．つづき

174

統計数理

第 57 巻

第 1 号 2009

変形バケットソートに現れる離散型確率分布と Eulerian 数

表 4．つづき

表 5.

表 6.

束の数と Eulerian 数の一致率．

付録 A の表を計算するのに要した時間．

175

統計数理

176

C.

第 57 巻

第 1 号 2009

定理 7.1 の証明

C.1 キュムラントとモーメントの関係
キュムラント κ(n)とモーメント En (X r ) の関係r∞∞ r
 tr (n)

t
exp
=
κr
En (X r )
r!
r!
r=1
r=0

から，キュムラントとモーメントの r = 6 までの関係は以下のように与えられる．ただし，(n)κr = κr ，En (X r ) = µr とする．
µ1 = κ1 ,µ2 = κ2 + κ21 ,µ3 = κ3 + 3κ2 κ1 + κ31 ,µ4 = κ4 + 4κ3 κ1 + 3κ22 + 6κ2 κ21 + κ41 ,µ5 = κ5 + 5κ4 κ1 + 10κ3 κ2 + 10κ3 κ21 + 15κ22 κ1 + 10κ2 κ31 + κ51 ,µ6 = κ6 + 6κ5 κ1 + 15κ4 κ2 + 15κ4 κ21 + 10κ23 + 60κ3 κ2 κ1 + 20κ3 κ31 + 15κ32+ 45κ22 κ21 + 15κ2 κ41 + κ61 .
C.2 1 次，2 次のキュムラント
定理 3.2 より明らか．
C.3 奇数次のキュムラント
X の確率分布は対称であるから，n ≥ 1 に対する r ≥ 3 次のキュムラントは 0 となる．
C.4 4 次のキュムラント
r = 4 のとき，
（5.1）式はE1 (X 4 ) = 1,n−44n − 66n − 4
En−1 (X 4 ) +En−1 (X 3 ) +En−1 (X 2 )En (X 4 ) =nnn
4n − 1+
En−1 (X) + 1nと表される．ここで，3 次以下のモーメントおよびキュムラントとモーメントの関係より，n ≥ 3
に対して
n − 4 (n−1)
1
(n)
κ4
(C.1)
κ4 =
−
n
24
(n)

を得る
（土屋・中村, 2007）．ここで，Cn (X) = κ4 + {(n + 1)/120} とおくと，（C.1）
式は Cn (X) ={(n − 4)/n}Cn−1 (X) となるから，n ≥ 4 に対して 4 次のキュムラントは(n)κ4 = −となる．
C.5 6 次のキュムラント
r = 6 のとき，
（5.1）式はE1 (X 6 ) = 1,

n+1
120

変形バケットソートに現れる離散型確率分布と Eulerian 数

En (X 6 ) =

177

n−6
6n − 1515n − 20
En−1 (X 6 ) +En−1 (X 5 ) +En−1 (X 4 )nnn
15n − 66n − 120n − 15
En−1 (X 3 ) +En−1 (X 2 ) +En−1 (X) + 1+nnnと表される．ここで，5 次以下のモーメントおよびキュムラントとモーメントの関係より，n ≥ 5
に対して
n − 6 (n−1)
1
(n)
κ4
(C.2)
κ6 =
+
n
36
を得る
（土屋・中村, 2007）．ここで，Cn (X) = κ(n)（C.2）
式は Cn (X) =4 − {(n + 1)/252} とおくと，{(n − 6)/n}Cn−1 (X) となるから，n ≥ 6 に対して 6 次のキュムラントは
(n)

κ6 =

n+1
252

となる．

参 考 文 献Euler, L.（1755）. Institutiones Calculi Diﬀerentialis.（John D. Blanton 訳（2000）．Foundations ofDiﬀerential Calculus, Springer, New York）
Graham, R. L., Knuth, D. E. and Patashnik, O.（1994）. Eulerian numbers, Concrete Mathematics:
A Foundation for Computer Science, 2nd ed., Section 6.2, 267–272, Addison-Wesley, Boston.
. Eulerian numbers and the unit cube, Fibonacci Quarterly, 20, 344–348.
Hensley, D.（1982）
Kimber, A. C.（1987）. Eulerian numbers and links with some statistical procedures, Utilitas Mathematica, 31, 57–65.
. Runs, The Art of Computer Programming, 2nd ed., Section 5.1.3, AddisonKnuth, D. E.（1998）
Wesley, Boston.
. On a test for randomness based on signs of diﬀerences, Annals of MathematicalMann, H. B.（1945）
Statistics, 16, 193–199.
Mood, A. M., Graybill, F. A. and Boes, D. C.（1974）. Introduction to the Theory of Statistics, 3rded., McGraw-Hill, Singapore.
. Enumerative Combinatorics, Vol. 1，Cambridge University Press, New York.
Stanley, R. P.（1997）
土屋高宏，中村永友（2007）. ある種の並べ替え算法における分布について，札幌学院大学商経論集，
24（2）, 31–47.

178

Proceedings of the Institute of Statistical Mathematics Vol. 57, No. 1, 159–178(2009)

Eulerian Numbers in Modified Bucket Sorting and Its Related
Distribution Theories
Takahiro Tsuchiya1 and Nagatomo Nakamura2
1 Department of Mathematics, Josai University
2 Department of Economics and General Education Center, Sapporo Gakuin University

This paper discusses the discrete probability distribution function induced from modiﬁed bucket sorting. The distribution is related to the unique sequence called “Eulerian
Numbers” proposed by Euler. The sequence number obtained by modiﬁed bucket sortingis expressed as a recurrence formula. It is mathematically equivalent to the recurrence
formula expression of Eulerian numbers. However, there are sensitive diﬀerences between
a derivation process of our proposal and Eulerian numbers expressed as a number of
permutation accents. The recurrence relation for the moments of the distribution is also
given. The mean, the variance and higher order cumulants are derived using the relation.
Interesting results of the relationship among the proposed, normal and uniform distributions are shown by asymptotic expansion. The eﬃciency of the proposed method is
illustrated through numerical experiments.

Key words: Asymptotic expansion, cumulant, Eulerian distribution, normal distribution, uniform
distribution.