KentのSmalltalkでの記事。何年ぶりなんだか。http://t.co/EN0kCYAE
@umejava 1月19日
経由で、
面白そうなのでちょっと趣向を変えて Squeak Smalltalk のジェネレーターを使って書いて動きを見てみようと思います。
サイコロより目の数がひとつ余計ですが、二番目のリンクの最後にあるのと同じ７つの目がそれぞれ 1/8, 1/5, 1/10, 1/4, 1/10, 1/10, 1/8 というふうに偏った確率で出る場合を考えます。
まず、ケント・ベックも冒頭で触れているエイリアス法を使わないナイーブな実装。各目の閾値（境界）を前述確率から算出しておき、0から1未満の乱数がそのどの範囲に収まるかで出る目を選んでいます。
| loadedDie bag | loadedDie := Generator on: [:g | | probs sum rand thresholds | probs := {1/8. 1/5. 1/10. 1/4. 1/10. 1/10. 1/8}. self assert: probs sum = 1. sum := 0. thresholds := probs collect: [:prob | sum := sum + prob]. rand := Random new. [ | index | index := rand next. g yield: (thresholds findFirst: [:threshold | index < threshold]) ] repeat]. bag := Bag new. 100000 timesRepeat: [bag add: loadedDie next]. bag sortedElements
=> {1->12466 . 2->20021 . 3->9944 . 4->25094 . 5->10003 . 6->10003 . 7->12469}
ちなみに、結果と比べやすくするため、それぞの目が出る確率として与えられた分数列を小数に直すとこんなふうになります。
{1/8. 1/5. 1/10. 1/4. 1/10. 1/10. 1/8} collect: #asFloat
=> #(0.125 0.2 0.1 0.25 0.1 0.1 0.125)
ついでにケント・ベックの #inject:into: を使った版も Squeak Smalltalk でジェネレーターを使って書いてみましょう。
| loadedDie bag | loadedDie := Generator on: [:g | | probs rand | probs := {1/8. 1/5. 1/10. 1/4. 1/10. 1/10. 1/8}. self assert: probs sum = 1. probs := probs collectWithIndex: [:prob :idx | idx -> prob]. rand := Random new. [ | index | index := rand next. [:exit | probs inject: 0 into: [:sum :each | | next | next := sum + each value. index < next ifTrue: [g yield: each key. exit value]. next]. ] valueWithExit. ] repeat]. bag := Bag new. 100000 timesRepeat: [bag add: loadedDie next]. bag sortedElements
=> {1->12465 . 2->20210 . 3->10025 . 4->24757 . 5->10012 . 6->9956 . 7->12575}
メソッドにしないこの書き方だと、リターン（^）が使えない（しかたがないので #valueWithExit 構造で代用）ので、うれしさ半減で冗長になってしまうようです。^^;
さて。本題のエイリアス法ですが、これは簡単には、偏りのある目の出方を「目の数と同じ枚数のコインから１枚のコインの抽出（等確率）」と「そのコインの偏った裏表の出方」という二段階の作業の結果として置き換える考え方のようです。使用する各コインの表には各目を割り振り、その裏には本来の 1/n より多く出る確率を持つ他の目を適切に分割・分配しておきます。この余剰分の「適切な分割・分配」は機械的にすることができて、二番目のリンクにある「Algorithm: Vose's Alias Method」の項がわかりやすいです。画付きの解説もあるので適宜参照して下さい。
これにほぼ従ったかたちで Squeak Smalltalk で書き下すとこんなかんじになります。
| loadedDie bag | loadedDie := Generator on: [:g | | probs n group small large coins rand | probs := {1/8. 1/5. 1/10. 1/4. 1/10. 1/10. 1/8}. self assert: probs sum = 1. n := probs size. probs := probs collectWithIndex: [:each :idx | idx -> (each * n)]. group := probs groupBy: [:each | each value < 1] having: #notEmpty. small := group at: true. large := group at: false. coins := OrderedCollection new. [small notEmpty] whileTrue: [ | heads tails | heads := small removeFirst. self assert: large notEmpty. tails := large removeFirst. tails value: tails value + heads value - 1. tails value >= 1 ifTrue: [large add: tails] ifFalse: [small add: tails]. heads value: heads value -> tails key. coins add: heads]. [large notEmpty] whileTrue: [ | heads | heads := large removeFirst. self assert: heads value = 1. heads value: heads value -> nil. coins add: heads]. rand := Random new. [ | coin | coin := coins atRandom. g yield: (rand next < coin value key ifTrue: [coin key] ifFalse: [coin value value]) ] repeat]. bag := Bag new. 100000 timesRepeat: [bag add: loadedDie next]. bag sortedElements
=> {1->12561 . 2->19976 . 3->10096 . 4->24892 . 5->9977 . 6->10031 . 7->12467}
このくらいの長さになったら面倒くさがらずにベックがやっているようにちゃんとクラスを作った方がよさそうですね。^^;
ちなみに上のそれぞれをブロックで括り timeToRun して速度を計測したところ、順に 1427ms, 1585ms, 622ms となりました（1.8GHz Core i7, Win 7, Squeak4.2 CogVM）。当たり前ですがそれなりに効果は出ているようです。
