Постановка задачи

Математическое моделирование

«Алгоритмы моделирования дискретных случайных величин»

Работу выполнил студент гр. 4084/21:

Кожаев Л.И.

Проверил доцент:

Хромов В.В.

Санкт-Петербург

Постановка задачи

Произвести исследование алгоритма моделирования дискретных случайных величин. За основу было взято биномиальное или распределение Ньютона.

Биномиальное распределение — дискретное распределение вероятностей случайной величины принимающей целочисленные значения с вероятностями:

Данное распределение характеризуется двумя параметрами: целым числом называемым числом испытаний, и вещественным числом называемом вероятностью успеха в одном испытании. Биномиальное распределение — одно из основных распределений вероятностей, связанных с последовательностью независимых испытаний. Если проводится серия из независимых испытаний, в каждом из которых может произойти "успех" с вероятностью то случайная величина, равная числу успехов во всей серии, имеет указанное распределение. Эта величина также может быть представлена в виде суммы независимых слагаемых, имеющих распределение Бернулли.

Характеристическая функция:

Моменты:

• Математическое ожидание:
• Дисперсия:
• Асимметрия: при распределение симметрично относительно центра

Тест проводился с использованием алгоритмов библиотеки Boost в языке С++.

Реализация задачи происходит на языке С++.

Листинг С++

Расчет вероятностей и квантилей для распределения Бернулли случайных величин, представляющих собой подбрасывания монеты.

Идеализированная монета состоит из круглого диска нулевой толщины, которая,
после подбрасывания в воздух и падения, будет возложена на одну из двух сторон лицевой стороной вверх ("Головы" H или "хвостов" T) с равной вероятностью.
Несмотря на небольшие различия между сторонами и ненулевой толщины фактического монеты, распределение их бросков имеет хорошее приближение к р == ½ распределения Бернулли.

Результат показывает биномиальное распределения для прогнозирования вероятности, когда бросали монеты, с головы или с хвоста.
Количество правильных ответов (скажем голов), X, распространяется в виде биномиальной случайной величины. Биномиальное распределение с параметрами числа испытаний (сальто) N = 10 и вероятности (success_fraction) получения головы р = 0,5 ("справедливой" монеты).

Наши монеты предполагаются справедливымы, но мы могли бы легко изменить параметр р success_fraction от 0,5 до другого значения для имитации несправедливой монеты.

#include <boost/math/distributions/binomial.hpp>

#include <iostream>

#include <iomanip>

using boost::math::binomial;

using std::cout; using std::endl; using std::left;

using std::setw;

int main()

{

cout << "Using Binomial distribution to predict how many heads and tails." << endl;

try

{

const double success_fraction = 0.5; // = 50% = 1/2 for a 'fair' coin.

int flips = 10;

binomial flip(flips, success_fraction);

cout.precision(4);

cout << "From " << flips << " one can expect to get on average "

<< mean(flip) << " heads (or tails)." << endl;

cout << "Mode is " << mode(flip) << endl;

cout << "Standard deviation is " << standard_deviation(flip) << endl;

cout << "So about 2/3 will lie within 1 standard deviation and get between "

<< ceil(mean(flip) - standard_deviation(flip)) << " and "

<< floor(mean(flip) + standard_deviation(flip)) << " correct." << endl;

cout << "Skewness is " << skewness(flip) << endl;

// Skewness of binomial distributions is only zero (symmetrical)

// if success_fraction is exactly one half,

// for example, when flipping 'fair' coins.

cout << "Skewness if success_fraction is " << flip.success_fraction()

<< " is " << skewness(flip) << endl << endl; // Expect zero for a 'fair' coin.

/*Now we show a variety of predictions on the probability of heads:*/

cout << "For " << flip.trials() << " coin flips: " << endl;

cout << "Probability of getting no heads is " << pdf(flip, 0) << endl;

cout << "Probability of getting at least one head is " << 1. - pdf(flip, 0) << endl;

/*When we want to calculate the probability for a range or values we can sum the PDF's:*/

cout << "Probability of getting 0 or 1 heads is "

<< pdf(flip, 0) + pdf(flip, 1) << endl; // sum of exactly == probabilities

cout << "Probability of getting 0 or 1 (<= 1) heads is " << cdf(flip, 1) << endl;

cout << "Probability of getting 9 or 10 heads is " << pdf(flip, 9) + pdf(flip, 10) << endl;

cout << "Probability of getting 9 or 10 heads is " << 1. - cdf(flip, 8) << endl;

cout << "Probability of getting 9 or 10 heads is " << cdf(complement(flip, 8)) << endl;

cout << "Probability of between 4 and 6 heads (4 or 5 or 6) is "

// P(X == 4) + P(X == 5) + P(X == 6)

<< pdf(flip, 4) + pdf(flip, 5) + pdf(flip, 6) << endl;

cout << "Probability of between 4 and 6 heads (4 or 5 or 6) is "

// P(X <= 6) - P(X <= 3) == P(X < 4)

<< cdf(flip, 6) - cdf(flip, 3) << endl;

cout << "Probability of between 3 and 7 heads (3, 4, 5, 6 or 7) is "

// P(X <= 7) - P(X <= 2) == P(X < 3)

<< cdf(flip, 7) - cdf(flip, 2) << endl;

cout << endl;

/*Finally, print two tables of probability for the /exactly/ and /at least/ a number of heads.*/

// Print a table of probability for the exactly a number of heads.

cout << "Probability of getting exactly (==) heads" << endl;

for (int successes = 0; successes <= flips; successes++)

{

double probability = pdf(flip, successes);

cout << left << setw(2) << successes << " " << setw(10)

<< probability << " or 1 in " << 1. / probability

<< ", or " << probability * 100. << "%" << endl;

}

cout << endl;

// Tabulate the probability of getting between zero heads and 0 upto 10 heads.

cout << "Probability of getting upto (<=) heads" << endl;

for (int successes = 0; successes <= flips; successes++)

{

// Say success means getting a head

// (equally success could mean getting a tail).

double probability = cdf(flip, successes); // P(X <= heads)

cout << setw(2) << successes << " " << setw(10) << left

<< probability << " or 1 in " << 1. / probability << ", or "

<< probability * 100. << "%"<< endl;

}

/*The last (0 to 10 heads) must, of course, be 100% probability.*/

}

catch(const std::exception& e)

{

/*It is always essential to include try & catch blocks because

default policies are to throw exceptions on arguments that

are out of domain or cause errors like numeric-overflow.

Lacking try & catch blocks, the program will abort, whereas the

message below from the thrown exception will give some helpful

clues as to the cause of the problem.*/

std::cout << "\n""Message from thrown exception was:\n " << e.what() << std::endl;

}

return 0;

} // int main()

Результат для 10 бросаний:

Вероятность не получить значения (ребро) 0.0009766

Вероятность получения одной из сторон 0.999

Вероятность получения 0 или 1 heads is 0.01074

Вероятность получения 0 или 1 (<= 1) heads is 0.01074

Вероятность получения 9 или 10 heads is 0.01074

Вероятность между 4 и 6 heads (4 or 5 or 6) is 0.6562

Вероятность между 4 и 6 heads (4 or 5 or 6) is 0.6563

Вероятность между 3 и 7 heads (3, 4, 5, 6 or 7) is 0.8906

Вероятность получить именно (==) heads

0 0.0009766 or 1 in 1024, or 0.09766%

1 0.009766 or 1 in 102.4, or 0.9766%

2 0.04395 or 1 in 22.76, or 4.395%

3 0.1172 or 1 in 8.533, or 11.72%

4 0.2051 or 1 in 4.876, or 20.51%

5 0.2461 or 1 in 4.063, or 24.61%

6 0.2051 or 1 in 4.876, or 20.51%

7 0.1172 or 1 in 8.533, or 11.72%

8 0.04395 or 1 in 22.76, or 4.395%

9 0.009766 or 1 in 102.4, or 0.9766%

10 0.0009766 or 1 in 1024, or 0.09766%

Вероятность получения меньшего значения (<=) heads

0 0.0009766 or 1 in 1024, or 0.09766%

1 0.01074 or 1 in 93.09, or 1.074%

2 0.05469 or 1 in 18.29, or 5.469%

3 0.1719 or 1 in 5.818, or 17.19%

4 0.377 or 1 in 2.653, or 37.7%

5 0.623 or 1 in 1.605, or 62.3%

6 0.8281 or 1 in 1.208, or 82.81%

7 0.9453 or 1 in 1.058, or 94.53%

8 0.9893 or 1 in 1.011, or 98.93%

9 0.999 or 1 in 1.001, or 99.9%

10 1 or 1 in 1, or 100%