Main Content

二項係数の厳密な計算

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

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

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

(nk)2n

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

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

320023691717107767864869141090753913186941388747093286534434787136210654094075586848270780341032844718978869737499691201939919030938206998752819929651447121650470765099112601116954925617580705043388419705530541025231105279169475241319583777521515740399490335296296986416422587436741374826165947451122678987232111208072306646349645444817081649045390224780439737117235599681944329390758068320042910611589088468568382563818033485652790542246479008106992991688671397929568442239221591755153140107136654000394609234559552383364790451240932289636703829891579702548735354741015983305611110777614681873617051793954211366022694113801876840128100034871409513586250746316776290259783425578615401030447369541046747571819748417910583511123376348523955353017744010395602173906080395504375010762174191250701116076984219741972574712741619474818186676828531882286780795390571221287481389759837587864244524002565968286448146002639202882164150037179450123657170327105882819203167448541028601906377066191895183769810676831353109303069033234715310287563158747705988305326397404720186258671215368588625611876280581509852855552819149745718992630449787803625851701801184123166018366180137512856918294030710215034138299203584

digitsvpa を使用して 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 = 

1368352466395056520289440683203474916723919590470093450966749985815252789703206921185078166194364368628241604751390276770010366566502905435584053306977090361373888807166934911669058359617267999440891473946928814790334474319243690058420836636313463906041227215461588805899687850446915983749011992904463517362389956795960668716644110185293844046419746984562010264700409766370663868950458143724613258018264237350867641134477342417401137040307238557519525647433371411711016304972400159139383774198454413097509455103606581317850616759661664811032314997464091016929165118720278926779247759469497371411534346514351633690928181552910415014721024800278971276108690005970534210322078267404628923208243578956328373980574557987343284668088987010788191642824141952083164817391248643164035000086097393530005608928615873757935780597701932955798545493414628255058294246363124569770299851118078700702913956192020527746291585168021113623057313200297435600989257362616847062553625339588328228815251016529535161470158485414650824874424552265877722482300505269981647906442611179237761490069369722948709004895010244652078822844422553197965751941043598302429006813614409472985328146929441100102855346352245681720273106393628672

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

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) = 

erf(10x-1000100)2+12

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 contains an object of type line.