Mathcad express

Пусть в аквариуме живут рыбки двух цветов.

Начнем с визуализации. Зададим число рыбок n=100 и договоримся что каждая из них имеет случайный цвет color 0 или 1, а также находится в случайной точке (x,y). Т.е. x, y, и color это три вектора длины n, а третью (z-) координату мы не рассматриваем.

%matplotlib inlineimport numpy as npimport matplotlib.pyplot as pltn = 100x, y = np.random.rand(n), np.random.rand(n)color = np.random.randint(0, 2, n)plt.scatter(x, y, c=color)

Модель аквариума

Те же три вектора в Mathcad Express можно сгенерировать так:

Генерация псевдослучайных векторов в Mathcad Express

Для векторов x,y используем генератор случайных чисел с равномерным распределением на отрезке (0, 1), а для вектора color биномиального распределения, моделирующего, в том числе, однократное подбрасывание монеты. В данном конкретном розыгрыше в Mathcad получилось 47 рыбок 1 и 100-47=53 рыбок нулевого цвета.

Теперь перейдем к главному вопросу: почему разные рыбки обладают цветом номер 0 или 1?

Все дело здесь в генетической информации, заключенной в биологическом материале, из которого состоит тело рыбки. Грубо говоря, живые существа, даже одного вида, различаются по химическому составу, причем в "длинных" молекулах биологических полимеров, из которых состоит все живое, зашифрована определенная информация, единицу которой обычно называют геном.

Совокупность генов каждого живого существа, называемая генотипом, определяет его внешний вид. В частности, генотип каждой рыбки определяет ее цвет. В простейшем случае, модель которого мы и будем дальше программировать, цвет рыбки (0 или 1) определяет один конкретный ген ("ген цвета"). Аналогичная ситуация (один ген определяет цвет) имела место в классических опытах Грегора Менделя, основоположника современной генетики, с растениями гороха, а в других случаях (например, для окраса собак-такс) на цвет отдельной особи влияет сочетание нескольких генов.

Как именно гены кодируют цвет? Прежде всего, ген цвета рыбки может существовать в двух вариантах, которые называются аллелями (я буду избегать этого и некоторых других важных для биолога терминов, чтобы сосредоточиться на информатике). В генетике их принято обозначать прописной и строчной буквами (например, А и а). Я буду чаще пользоваться цифрой А=0, а=1.

Каждое растение гороха в опытах Менделя и каждая рыбка из нашего аквариума несет в себе информацию о двух аллелях, т.е. описывается одной из трех возможных комбинаций:

• АА или 00

• Аа (то же самое, что аА) или 01

• аа или 11

Если рыбка несет ген цвета АА или Аа, то она серая, а если аа желтая. Т.е. тип гена А полностью блокирует тип а, поэтому вариант А (0) называют доминатным, а а (1) рецессивным.

Давайте реализуем расчет цвета рыбки по ее генетическому коду, который для i-й рыбки будем описывать вектором из двух компонентов, а генетический состав всей популяции соответственно, матрицей, которую назовем Аа. Пусть у нас есть три рыбки, каждая из которых с равной вероятностью несет генетический код 00, 10, 01 или 11. В первых трех случаях цвет рыбки будет серым, а в четвертом желтым.

n = 3Aa = np.zeros((n,2))color = np.zeros(n)for i in range(n):  Aa[i,0] = np.random.randint(0, 2, 1)  Aa[i,1] = np.random.randint(0, 2, 1)  color[i] = Aa[i,0] * Aa[i,1]print (Aa)print (color)

Результат работы программы вы видите на картинке:

Моделирование цвета рыбок через ее генотип

Очевидно, при такой модели, если рыбок будет много, то цветных будет около 25%. Аналогичные вычисления в Mathcad для популяции из n=1000 рыбок приведены на следующем рисунке.

Моделирование цвета рыбок через ее генотип (Mathcad)

Все эти расчеты я виду в режиме стримов на своем YouTube канале , т.е. фактически, данная статья это краткий дайджест первых двух роликов по вычислительной генетике. А в следующих сериях я планирую организовать пересчет генотипа рыбок от поколения к поколению, причем, видимо, мне придется перейти к ООП и описывать рыбок при помощи классов в Python, а Mathcad использовать, как средство вспомогательных расчетов. В качестве опорного материала я использую Википедию и вот этот курс на Степике.

В последнее время появляется все больше и больше аналитических обзоров результатов выборов, которые рассматривают их с точки зрения законов статистики и направлены, как правило, на изучение необычных явлений, сигнализирующих о возможных фальсификациях (см. "Гаусс против Чурова" и т.п. публикации). Думается, только что завершившиеся выборы президента США (подсчет голосов еще продолжается, причем беспрецедентно длительное время) дадут дополнительный толчок развитию этой "электоральной математике".

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

Так получилось, что, начиная с 5 ноября, я стал следить за промежуточными результатами подсчета голосов в штате Джорджия. Тогда было посчитано 95% голосов, соотношение было 49.6/49.1 в пользу Трампа, т.е. он был впереди на полпроцента. Исключительно из демонстрационных соображений (а я делаю разный образовательный контент и читаю лекции студентам), я решил оценить вероятность того, что штат "перевернется", т.е. окончательный результат будет в пользу Байдена. Собственно, 5-го ноября я начал стримить свои расчеты, но потом отвлекся и прекратил эфир, зато запись с исходными данными осталась. Когда я начал писать эту статью, в 9 утра ЕТ 6-го ноября, было посчитано 99% бюллетеней, соотношение голосов на этот момент вы видите на заставке, а текущие результаты вы можете сами посмотреть на сайте CNN, откуда я и беру данные. К слову, на момент публикации этой статьи, через двое суток, т.е. 9 утра ЕТ 8-го ноября, окончательных результатов все еще не было, а на табло по-прежнему отображались 99%, но с еще более убедительным результатом в пользу Байдена.

Естественной моделью, которую я решил взять за основу, была модель независимых испытаний, т.е. классическая схема Бернулли. Проводя аналогию между выбором каждого из оставшихся 100 тыс. избирателей (за Трампа - за Байдена) и 100-тысячекратным бросанием монеты (орел - решка), я собираюсь оценить упомянутую вероятность победы Байдена в Джорджии. Недействительные бюллетени и голоса в поддержку других кандидатов учитывать не будем. Трампа в России любят больше, поэтому давайте условимся считать голос за него успехом (х=1), а голос против - неудачей (х=0). Сразу хочу сказать, что в этой статье мне не хочется обсуждать политику и делать какие-либо выводы о фальсификациях и т.п., а посвятить ее зарисовке моделирования схемы Бернулли в Mathcad Express.

Схема Бернулли: модель Монте-Карло

Итак, давайте считать, что голос каждого избирателя в Джорджии похож на бросание монеты, благо вероятность успеха (за Трампа) или неудачи (за Байдена), если судить по посчитанным 95% бюллетеней, практически равны. Давайте, следуя методам Монте-Карло, сгенерируем вектор х из N=1000 случайных чисел, каждое из которых (0 или1) добавляется в выборку с равной 50%-й вероятностью:

Модель Бернулли: однократное бросание "честной монеты"

В предпоследней строке выведены первые несколько компонент вектора х, а в последней - выборочное среднее (которое близко к математическому ожиданию 0.5).

Это пока "модель поведения одного избирателя". А если мы хотим добавить еще одного, т.е. промоделировать двукратное бросание монеты, то можно генерировать не одно, а два псевдослучайных числа: х1 и х2, а судить о выпадении того или иного количества успехов и неудач мы можем, просто вычисляя сумму этих чисел х=х1+х2. Например, вот так будут выглядеть векторы х1, х2 и х для N=10 двукратных бросаний:

Модель Бернулли: двукратное бросание монеты

Для того чтобы автоматизировать расчет частоты выпадения, к примеру, двух успехов (1+1), надо просто подсчитать количество двоек в векторе х. Сделать это можно при помощи условного оператора if, который доступен в бесплатном Mathcad Express.

Пример расчета количества исходов х=1+1

Многократное бросание монеты в Mathcad Express можно организовать проще, без привлечения переменных х1 и х2, сразу считая х при помощи той же встроенной функции rbinom(N,M,p), где M-число бросаний, а p=0.5 - вероятность успеха:

Пример расчета частоты исходов x<2 при двукратном бросании монеты

Моделируем выборы

Теперь все готово для решения основной задачи. Мы можем запросто смоделировать поведение М=250000 избирателей в штате Джорджия и посчитать искомую оценку вероятности, т.е. среднюю частоту событий x<A. Данные на начало дня 05.11 на сайте CNN выглядели так

Данные на начало 05.11 (95% голосов посчитано)

За Трампа было подано T=2 429 783 голоса, а за Байдена B=2 406 774, т.е. соответственно 50.24% и 49.76% соответственно (прочие голоса, ни за Трампа, ни за Байдена не считаем). Т.е. всего на тот момент было посчитано около 5 млн голосов, т.е. 95% бюллетеней. Поэтому можно считать, что осталось посчитать оставшиеся 5%, т.е. примерно 250 000 бюллетеней. Разрыв между Трампом и Байденом составлял 2 429 783 - 2 406 774 = 23 009. Иными словами, нам надо оценить вероятность того, что число успехов на выборке по М=250 000 бросаниям монеты будет на 23 тыс. меньше, нежели неудач. В наших обозначениях, надо установить с какой вероятностью х не превысит А = 250 000 / 2 - 11 500 (тогда разрыв между 125+11.5 тыс. и 125-11.5 тыс. будет как раз не меньше 23 тыс.).

Если взять объем выборки N=10 млн испытаний, то генератор псевдослучайных чисел Mathcad не дает ни одного случая x<A, т.е. оцениваемая вероятность победы Байдена (если верна модель Бернулли) получается ничтожной.

Если снизить (ради любопытства) А в 10 раз и попробовать посчитать среднюю частоту события x<A/10, то из миллиона испытаний оно выпадет всего 1-2 раза:

Модель Бернулли применительно к электоральной математике

Аналитические оценки

Модель Бернулли хороша тем, что для нее известны аналитические формулы вычисления подобных вероятностей. В частности, рассматриваемую задачу можно решить при помощи интегральной теоремы Муавра-Лапласа, которая дает тот же порядок вероятности (одна миллионная для условия x<A/10):

Расчет вероятности x<A/10 по теореме Муавра-Лапласа

Но об аналитических формулах для схемы Бернулли в следующей статье. А эту я завершаю в час ночи 09.11.2020 по московскому времени (соответствует 15 часам ET 08.11 в США). А читатель сам может оценить вероятность исхода, который в итоге реализуется, по скриншоту с сайта CNN. На нем цифры 99% показывают нам и всему миру, что "последнюю милю" никак не удается пройти, и несчастный 1 процент бюллетеней в Джорджии стоит насмерть и до последнего сопротивляется расчетам.