二項係数を含む和を求めるプログラムの効率化を幾つかの方法で行う例を示します。今回のその１では、式を求めそれを変形して漸化式を使って$N=10^{10}$まで求めます。その２では全く違った観点でライプニッツの調和三角形から式を割り出してほぼ$N=\infty$まで求めます方法を紹介します。
【例題】2組のn枚のコインの表の数がkの時のパターン一致したら1/kポイントを貰えるゲームの期待値
ゲームは$n$ターン行います。$k$ターン目ではまず$n$個のコインを表が$k$個になるまで投げて、そのパターンを覚えておきます。もう一度表が$k$個になるまで投げて、そのパターンが最初と同じであれば$1/k$ポイントを得ます。各ターンのポイントの期待値を$E(n)$とするとき、トータルの期待値$S(N) = \sum_{n=1}^{N}E(n)$を１分以内でなるべく大きい$N$まで求めよ
期待値を求める式を求める
問題がややこしいのでまず具体的な例で考えます。
\begin{align}
&n = 5 \\
&k = 2 とすると\\
&一回目が\\
& （表裏表裏裏） \ \ \ だった時\\
& 二回目も同じパターンになる確率は\frac{1}{\binom{5}{2}}\\
& 取得ポイントの期待値は\frac{1}{\binom{5}{2}\cdot 2} \\
& したがって各ターンの期待値は \\
& e_n(k) = \frac{1}{\binom{n}{k}\cdot k}\\
& これをk=1..nで足せばよいので\\
& E(n) = \sum_{k=1}^{n}\frac{1}{\binom{n}{k}\cdot k} \cdots (1-1)\\
& S(n) = \sum_{n=1}^{N}\sum_{k=1}^{n}\frac{1}{\binom{n}{k}\cdot k} \cdots (1-2)\\
& となります\\
\end{align}
この式(1-2)の値をNを段々大きくしながら求めていきます。
1. sympy.binomialを使う
まずそのままsp.binomialを使ってn=10で求めます。1分以内で求められるのは$n=10^3$まで
import sympy as sp
import numpy as np
def E1(n):
B = np.array([sp.binomial(n,k) for k in range(1,n+1)])
W = np.array([k for k in range(1,n+1)])
return sum(np.ones(n)/B/W)
def S1(N):
return sum([E(n) for n in range(1,N+1)])
print(f"{S1(10):.09f}")
# 5.626190476
2. 二項係数をパスカルの三角形を使って求める
さすがに二項係数を毎回求めるのは無駄なのでパスカルの三角形を使って漸化的に求めます。だいぶ改善しましたが$n=10^4$まで。
B1 = [1.]
def E2(n):
global B1
B = [1.]+[B1[k]+B1[k+1] for k in range(n-1)]+[1.]
B1 = B
W = np.array([k for k in range(1,n+1)])
return sum(np.ones(n)/B[1:]/W)
def S2(N):
return sum([E2(n) for n in range(1,N+1)])
print(f"Answer = {S2(10**4):.09f}")
# Answer = 19.575012052
3. scipyのgammalnとlogsumexpを使う
二項係数は$n$が大きくなると指数的に大きくなるので、常套手段として対数を取ることを考えます。そこで階乗にはgammmalnを使い、足し算にはlogsumexpを使います。値のオーバーフローは抑えられますが、高速化の観点ではそれほど効果はありませんでした。$n=5 \times 10^4$まで。
from scipy.special import gammaln, logsumexp
def E3(n):
ks = np.arange(1, n + 1)
B = gammaln(n + 1) - gammaln(ks + 1) - gammaln(n - ks + 1)
W = np.log(ks)
return logsumexp(- B - W)
def S3(N):
return sum([np.exp(E3(n)) for n in range(1,N+1)])
print(f"Answer = {S3(5*10**4):.09f}")
# Answer = 22.793967898
4. 式の変形から漸化式を求める
そこで式(1)の変形を考えます。6. 式の変形の詳細で示した方法で以下のように変形できます。
\begin{align}
& E(n) = \sum_{k=1}^{n}\frac{1}{\binom{n}{k}\cdot k}=\frac{1}{2^{n}}\sum_{k=1}^{n}\frac{2^k}{k}&\cdots(3)\\
\end{align}
となり二項係数の逆数の式を$2^k$の式に変えることが出来きました。ここから漸化式を求めます。
\begin{align}
E(n) &=\frac{1}{2^{n}}\sum_{k=1}^{n}\frac{2^k}{k}\\
E(n-1) &=\frac{1}{2^{n-1}}\sum_{k=1}^{n-1}\frac{2^k}{k}\\
E(n) - \frac{1}{2}E(n-1) &= \frac{1}{2^{n}}(\sum_{k=1}^{n}\frac{2^k}{k}-\sum_{k=1}^{n-1}\frac{2^k}{k}) \\
&= \frac{1}{2^{n}} \frac{2^{n}}{n} = \frac{1}{n} \\
&従って、\\
E(n) &= \frac{1}{2}E(n-1)+ \frac{1}{n}&\cdots(3)\\
\end{align}
5. 漸化式を使ったコード
この漸化式を使ってで以下の非常にシンプルなコードになりました。$numba$を使えば$10^{10}$まで求めることができました。
from numba import njit
@njit
def S5(N):
e = res = 1.0
for n in range(2, N+1):
e = e/2 + 1./n
res += e
return res
n = 10**10
print(f"Answer = {S5(n):.09f}")
# Answer = 47.2061331898
6. 式の変形の詳細
\begin{align}
& 二項係数の公式より &\\
& \binom{r}{k} = \frac{r}{k}\binom{r-1}{k-1} &\cdots(a) \\
& これを(1)に代入すると& \\
& E(n) = \sum_{k=1}^{n}\frac{1}{n \cdot \binom{n-1}{k-1}} &\\
& = \frac{1}{n}\sum_{k=1}^{n}\frac{1}{\binom{n-1}{k-1}} \\
& = \frac{1}{n}\sum_{k=0}^{n-1}\frac{1}{\binom{n-1}{k}} &\cdots(2) \\
\end{align}
ここでCalculate sums of inverses of binomial coefficientsで紹介されている(b)を(2)に代入して$n$に$n-1$にすると
\begin{align}
&\sum_{k=0}^n\frac1{\binom{n}{k\vphantom{+1}}}
=\frac{n+1}{2^{n+1}}\sum_{k=1}^{n+1}\frac{2^k}{k}&\cdots(b)\\
& E(n) = \frac{1}{n}\frac{n}{2^{n}}\sum_{k=1}^{n}\frac{2^k}{k}\\
&=\frac{1}{2^{n}}\sum_{k=1}^{n}\frac{2^k}{k}&\cdots(3)\\
\end{align}