二項係数の厳密な計算

この例では、Symbolic Math Toolbox™ を使用して二項係数の正確な値を取得し、コイン投げ実験における確率を求める方法を説明します。

シンボリック関数 P(n,k) を定義します。これは n 回投げた中で厳密に k 回表が出る確率を計算します。

syms P(n,k)
P(n,k) = nchoosek(n,k)/2^n

P(n, k) = 
 $\frac{(\binom{n}{k})}{2^{n}}$

コインを 2000 回投げるとします。投げた回数の 2 分の 1 で表が出る確率は P(n, n/2) となります。ここでは、n = 2000 です。結果は分子と分母ともに大数の有理式となります。Symbolic Math Toolbox は厳密な結果を返します。

n = 2000;
central = P(n, n/2)

central = 
 $\frac{32002369171710776786486914109075391318694138874709328653443478713621065409407558684827078034103284471897886973749969120193991903093820699875281992965144712165047076509911260111695492561758070504338841970553054102523110527916947524131958377752151574039949033529629698641642258743674137482616594745112267898723211120807230664634964544481708164904539022478043973711723559968194432939075806832004291061158908846856838256381803348565279054224647900810699299168867139792956844223922159175515314010713665400039460923455955238336479045124093228963670382989157970254873535474101598330561111077761468187361705}{1793954211366022694113801876840128100034871409513586250746316776290259783425578615401030447369541046747571819748417910583511123376348523955353017744010395602173906080395504375010762174191250701116076984219741972574712741619474818186676828531882286780795390571221287481389759837587864244524002565968286448146002639202882164150037179450123657170327105882819203167448541028601906377066191895183769810676831353109303069033234715310287563158747705988305326397404720186258671215368588625611876280581509852855552819149745718992630449787803625851701801184123166018366180137512856918294030710215034138299203584}$

digits と vpa を使用して 10 桁の精度でこの有理数を近似します。

previous_digits = digits(10);
vpa(central)

ans =  $0.01783901115$

"表" の回数と期待値との差が 2 標準偏差以内となる確率を計算します。

sigma = sqrt(n/4);
withinTwoSigma = symsum(P(n, k), k, ceil(n/2 - 2*sigma), floor(n/2 + 2*sigma))

withinTwoSigma = 
 $\frac{13683524663950565202894406832034749167239195904700934509667499858152527897032069211850781661943643686282416047513902767700103665665029054355840533069770903613738888071669349116690583596172679994408914739469288147903344743192436900584208366363134639060412272154615888058996878504469159837490119929044635173623899567959606687166441101852938440464197469845620102647004097663706638689504581437246132580182642373508676411344773424174011370403072385575195256474333714117110163049724001591393837741984544130975094551036065813178506167596616648110323149974640910169291651187202789267792477594694973714115343465}{14351633690928181552910415014721024800278971276108690005970534210322078267404628923208243578956328373980574557987343284668088987010788191642824141952083164817391248643164035000086097393530005608928615873757935780597701932955798545493414628255058294246363124569770299851118078700702913956192020527746291585168021113623057313200297435600989257362616847062553625339588328228815251016529535161470158485414650824874424552265877722482300505269981647906442611179237761490069369722948709004895010244652078822844422553197965751941043598302429006813614409472985328146929441100102855346352245681720273106393628672}$

結果を浮動小数点数で近似します。

vpa(withinTwoSigma)

ans =  $0.9534471795$

この結果を正規分布の累積分布関数 (cdf) から導出した以下の 2 つの推定値と比較します。2 シグマ区間の外側にある最初の整数と区間の内側にある最初の整数との中点をとると、2 シグマ区間それ自体を使用するよりも精度の高い結果が得られることがわかります。

syms cdf(x)
cdf(x) = 1/2 * (1 + erf((x - n/2)/sqrt(sym(n/2))))

cdf(x) = 
 $\frac{\erf (\frac{\sqrt{10} (x - 1000)}{100})}{2} + \frac{1}{2}$

estimate1 = vpa(cdf(n/2 + 2*sigma) - cdf(n/2 - 2*sigma))

estimate1 =  $0.9544997361$

estimate2 = vpa(cdf(floor(n/2 + 2*sigma) + 1/2) - ...
                cdf(ceil(n/2 - 2*sigma) - 1/2))

estimate2 =  $0.9534201342$

error1 = vpa(estimate1 - withinTwoSigma)

error1 =  $0.001052556595$

error2 = vpa(estimate2 - withinTwoSigma)

error2 =  $- 0.00002704525904$

既定の精度に戻して浮動小数点計算を続けます。

digits(previous_digits);

$2 σ$ 区間内の k の分布をプロットします。プロットはガウス曲線になります。

k = n/2 + (-2*sigma:2*sigma);
plot(k, P(n,k), '-+');
title('2000 coin tosses');
xlabel('Number of heads');
ylabel('Probability');

Figure contains an axes object. The axes object with title 2000 coin tosses, xlabel Number of heads, ylabel Probability contains an object of type line.