Вычислите биномиальные коэффициенты точно

В этом примере показано, как получить точные значения для биномиальных коэффициентов и найти вероятности в бросающих монету экспериментах с помощью Symbolic Math Toolbox.

Задайте символьную функцию, P(n,k), это вычисляет вероятность для голов, чтобы подойти точно k времена из n броски.

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

(nk)2nnchoosek (n, k)/2^n

Предположим, вы бросаете монету 2000 раз. Вероятностью, что головы подходят в половине бросков, является P(n, n/2), где n = 2000. Результатом является рациональное выражение с большими количествами и в числителе и в знаменателе. Symbolic Math Toolbox возвращает точный результат.

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

320023691717107767864869141090753913186941388747093286534434787136210654094075586848270780341032844718978869737499691201939919030938206998752819929651447121650470765099112601116954925617580705043388419705530541025231105279169475241319583777521515740399490335296296986416422587436741374826165947451122678987232111208072306646349645444817081649045390224780439737117235599681944329390758068320042910611589088468568382563818033485652790542246479008106992991688671397929568442239221591755153140107136654000394609234559552383364790451240932289636703829891579702548735354741015983305611110777614681873617051793954211366022694113801876840128100034871409513586250746316776290259783425578615401030447369541046747571819748417910583511123376348523955353017744010395602173906080395504375010762174191250701116076984219741972574712741619474818186676828531882286780795390571221287481389759837587864244524002565968286448146002639202882164150037179450123657170327105882819203167448541028601906377066191895183769810676831353109303069033234715310287563158747705988305326397404720186258671215368588625611876280581509852855552819149745718992630449787803625851701801184123166018366180137512856918294030710215034138299203584sym('32002369171710776786486914109075391318694138874709328653443478713621065409407558684827078034103284471897886973749969120193991903093820699875281992965144712165047076509911260111695492561758070504338841970553054102523110527916947524131958377752151574039949033529629698641642258743674137482616594745112267898723211120807230664634964544481708164904539022478043973711723559968194432939075806832004291061158908846856838256381803348565279054224647900810699299168867139792956844223922159175515314010713665400039460923455955238336479045124093228963670382989157970254873535474101598330561111077761468187361705/1793954211366022694113801876840128100034871409513586250746316776290259783425578615401030447369541046747571819748417910583511123376348523955353017744010395602173906080395504375010762174191250701116076984219741972574712741619474818186676828531882286780795390571221287481389759837587864244524002565968286448146002639202882164150037179450123657170327105882819203167448541028601906377066191895183769810676831353109303069033234715310287563158747705988305326397404720186258671215368588625611876280581509852855552819149745718992630449787803625851701801184123166018366180137512856918294030710215034138299203584')

Аппроксимируйте это рациональное число 10-разрядной точностью с помощью digits и vpa.

previous_digits = digits(10);
vpa(central)
ans = 0.01783901115vpa ('0.01783901115')

Вычислите вероятность, что количество "голов" отличается от ожидаемого значения не больше, чем двумя стандартными отклонениями.

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

1368352466395056520289440683203474916723919590470093450966749985815252789703206921185078166194364368628241604751390276770010366566502905435584053306977090361373888807166934911669058359617267999440891473946928814790334474319243690058420836636313463906041227215461588805899687850446915983749011992904463517362389956795960668716644110185293844046419746984562010264700409766370663868950458143724613258018264237350867641134477342417401137040307238557519525647433371411711016304972400159139383774198454413097509455103606581317850616759661664811032314997464091016929165118720278926779247759469497371411534346514351633690928181552910415014721024800278971276108690005970534210322078267404628923208243578956328373980574557987343284668088987010788191642824141952083164817391248643164035000086097393530005608928615873757935780597701932955798545493414628255058294246363124569770299851118078700702913956192020527746291585168021113623057313200297435600989257362616847062553625339588328228815251016529535161470158485414650824874424552265877722482300505269981647906442611179237761490069369722948709004895010244652078822844422553197965751941043598302429006813614409472985328146929441100102855346352245681720273106393628672sym('13683524663950565202894406832034749167239195904700934509667499858152527897032069211850781661943643686282416047513902767700103665665029054355840533069770903613738888071669349116690583596172679994408914739469288147903344743192436900584208366363134639060412272154615888058996878504469159837490119929044635173623899567959606687166441101852938440464197469845620102647004097663706638689504581437246132580182642373508676411344773424174011370403072385575195256474333714117110163049724001591393837741984544130975094551036065813178506167596616648110323149974640910169291651187202789267792477594694973714115343465/14351633690928181552910415014721024800278971276108690005970534210322078267404628923208243578956328373980574557987343284668088987010788191642824141952083164817391248643164035000086097393530005608928615873757935780597701932955798545493414628255058294246363124569770299851118078700702913956192020527746291585168021113623057313200297435600989257362616847062553625339588328228815251016529535161470158485414650824874424552265877722482300505269981647906442611179237761490069369722948709004895010244652078822844422553197965751941043598302429006813614409472985328146929441100102855346352245681720273106393628672')

Аппроксимируйте результат числом с плавающей запятой.

vpa(withinTwoSigma)
ans = 0.9534471795vpa ('0.9534471795')

Сравните этот результат со следующими двумя оценками, выведенными из кумулятивной функции распределения (cdf) нормального распределения. Оказывается, что взятие средней точки между первым целым числом снаружи и первым целым числом в интервале 2D сигмы дает более точный результат, чем использование самого интервала 2D сигмы.

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

erf(10x-1000100)2+12erf ((sqrt (sym (10)) * (x - 1000))/100)/2 + sym (1/2)

estimate1 = vpa(cdf(n/2 + 2*sigma) - cdf(n/2 - 2*sigma))
estimate1 = 0.9544997361vpa ('0.9544997361')
estimate2 = vpa(cdf(floor(n/2 + 2*sigma) + 1/2) - ...
                cdf(ceil(n/2 - 2*sigma) - 1/2))
estimate2 = 0.9534201342vpa ('0.9534201342')
error1 = vpa(estimate1 - withinTwoSigma)
error1 = 0.001052556595vpa ('0.001052556595')
error2 = vpa(estimate2 - withinTwoSigma)
error2 = -0.00002704525904- vpa ('0.00002704525904')

Восстановите точность по умолчанию для расчетов с плавающей точкой.

digits(previous_digits);

Постройте это распределение для k в 2σИнтервал. Графиком является Кривая Гаусса.

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

Для просмотра документации необходимо авторизоваться на сайте