唯物是真 @Scaled_Wurm

プログラミング(主にPython2.7)とか機械学習とか

複数のサイコロの目の和がある値になる確率を求める

\(K\)面ダイス(サイコロ)を\(N\)個投げた時に、合計がある値になる確率を求めたい

合計がある値になる場合の数を求めて\(K^N\)で割って確率にする方がやりやすそうなので、場合の数を求める方法を考える

総当り

\(N\)個のダイスについてそれぞれ\(1\)から\(K\)までループするナイーブな総当りで和を求めていくと\(O(K^N)\)ぐらいの計算量が必要になってしまい、サイコロの数が増えると現実的な時間で終わらなくなってしまう

動的計画法

漸化式を立てて動的計画法で計算すると\(O(N^2K^2)\)ぐらいで計算できる
具体的には\(n\)個サイコロを使って合計が\(c\)になる場合の数を求めたいとき\(f_n(c) = \sum_{i=1}^K f_{n-1}(c-i)\)という漸化式が成り立つ

K = 6
N = 100

dp = [[0] * N * K for i in xrange(2)]

for i in xrange(K):
    dp[0][i] = 1

for i in xrange(1, N):
    for j in xrange(N * K):
        dp[i % 2][j] = 0
        for k in xrange(max(0, j - K), j):
            dp[i % 2][j] += dp[(i - 1) % 2][k]

for j in dp[N - 1 % 2]:
    print j,
print

動的計画法+累積和

動的計画法の一つのサイコロに対するループで、添え字に\(1\)から\(K\)足した範囲に同じ数を足していくので、累積和に関するアルゴリズムを使えばより高速化できる
\(0\)初期化した配列を用意して、足したい範囲の開始地点に足したい数を加えて範囲の終了地点の次の添字の場所に足したい数を引いた配列を作って、各地点までの要素の累積和の配列を計算すればよい
下記記事の図がわかりやすい

計算量は\(O(N^2K)\)になる

K = 6
N = 100

dp = [[0] * N * K for i in xrange(2)]

for i in xrange(K):
    dp[0][i] = 1

for i in xrange(1, N):
    dp[i % 2] = [0] * N * K
    for j in xrange(N * K - 1):
        dp[i % 2][j + 1] += dp[(i - 1) % 2][j]
        if j + K + 1 < N * K:
            dp[i % 2][j + K + 1] -= dp[(i - 1) % 2][j]
    cur = 0
    for j in xrange(N * K):
        cur += dp[i % 2][j]
        dp[i % 2][j] = cur

for j in dp[N - 1 % 2]:
    print j,
print

フーリエ変換

確率変数の和の分布は元の分布の畳み込みになるので、フーリエ変換して係数をかけあわせて逆フーリエ変換すると目的の値が得られる(はず)

www.slideshare.net
計算量は\(O(NK \log NK)\)
ただしサイコロの数を増やしたらオーバーフローしたり変な値になったりして正しい答えが得られなかった(´・ω・｀)
浮動小数点演算の誤差なのかなぁ……(？)

検証してくださった方がいたので追記(2015-09-09)


# -*- coding: utf-8 -*-
import math,string,itertools,fractions,heapq,collections,re,array,bisect

#http://www.slideshare.net/chokudai/fft-49066791
import cmath

def pow_2_at_least(x):
    n = 1
    while n < x:
        n *= 2
    return n

def convolution(g, N):
    Lg = len(g)
    n = pow_2_at_least(Lg * N)
    g = g + [0] * (n - Lg)
    zeta_arr = [cmath.exp(2j * cmath.pi / n) for n in xrange(1, n + 1)]
    gf = fft(g, n, zeta_arr)
    ff = [pow(gf[i], N) for i in xrange(n)]
    return ifft(ff, n, zeta_arr)


def _fft(f, n, zeta_arr, inv):
    m = n / 2
    if m == 1:
        f0 = [f[0]]
        f1 = [f[1]]
    else:
        f0 = [f[2 * i + 0] for i in xrange(m)]
        f1 = [f[2 * i + 1] for i in xrange(m)]
        f0 = _fft(f0, m, zeta_arr, inv)
        f1 = _fft(f1, m, zeta_arr, inv)
    
    sign = 1 - 2 * inv
    if inv:
        zeta = 1.0 / zeta_arr[n - 1]
    else:
        zeta = zeta_arr[n - 1]
    pow_zeta = 1
    for i in xrange(n):
        f[i] = f0[i % m] + pow_zeta * f1[i % m]
        pow_zeta *= zeta
    return f

def fft(f, n, zeta_arr):
    return _fft(f, n, zeta_arr, False)


def ifft(f, n, zeta_arr):
    return [i / n for i in _fft(f, n, zeta_arr, True)]


K = 6
N = 100

print '\n'.join([str(int(round(i.real))) for i in convolution([0] + [1.0] * K, N)][1:N * K + 1])

参考URL

いろいろな確率とかが載ってる

\(K\)面ダイスを\(N\)個投げた時に合計がある値\(C\)になる確率を求める以下の式の導出が載っている
$$P(N, K, C)=\frac{1}{K^N}\sum_{i=0}^{\lfloor \frac{C - N}{K} \rfloor}(-1)^i \binom Ni \binom{C-si-1}{N-1}$$

'); jQuery.noConflict(true);