Программирование и научные вычисления на языке Python/§18

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

Принципы моделирования методом Монте-Карло

Представим, мы провели N экспериментов и результат каждого эксперимента случайное число. В этих N экспериментах некоторое событие случается M раз. Оценка вероятности тогда M/N и она становится более точной при стремлении числа экспериментов N к бесконечности (заметьте, что дробь при этом не обращается в ноль, поскольку и число M также стремится к бесконечности). Математический метод, заключающийся в использовании большого числа генерируемых случайных чисел, получил название метод Монте-Карло. Этот метод оказался чрезвычайно полезным в науке и промышленности, там где случайное поведение нельзя не учитывать, или, например, там где и не подразумевается наличие случайных чисел — в случаях сложного интегрирования.

Бросание костей

Пусть у нас спрашивают: какова вероятность при броске костей, чтобы по крайне мере на двух из них были шестерки? Итак, эксперимент в том, что кидаются четыре кости, нужное событие - не менее двух кубиков с шестью точками. Собственно, это нам и нужно записать еще раз перефразировав на Python:

import random as random_number

N = int(raw_input('Enter number of experiments: '))
M = 0
for i in range(N):
    six = 0
    r1 = random_number.randint(1, 6)
    if r1 == 6:
        six += 1
    r2 = random_number.randint(1, 6)
    if r2 == 6:
        six += 1
    r3 = random_number.randint(1, 6)
    if r3 == 6:
        six += 1
    r4 = random_number.randint(1, 6)
    if r4 == 6:
        six += 1
    if six >= 2:
        M += 1

p = float(M)/N
print 'probability:', p

Более общая задача — ведь мы можем задать и не только число экспериментов, а, например, число бросаемых костей (ndice) и минимальное число одинаковых кубиков (nsix). Таким образом, мы получаем более общий и даже более короткий код, а более общие решения гораздо проще модифицировать далее:

import random as random_number
import sys

N = int(raw_input('Number of experiments: '))
ndice = int(raw_input('Number of dice: '))
nsix = int(raw_input('Number of dice with six eyes: '))
M = 0
for i in range(N):
    six =0
    for j in range(ndice):
        r = random_number.randint(1, 6)
        if r == 6:
            six += 1
    if six >= nsix:
        M += 1

p = float(M)/N
print 'Probability:', p

Для малых вероятностей число M мало и аппроксимация M/N не очень хорошая аппроксимация, требуется большое число экспериментов, которые для стандартного модуля обходятся дорого, то есть долго. Конечно же, мы знаем, что при нехватке скорости надо перейти к более быстрым векторным операциям.

Векторизация состоит в том, что мы создаем двухмерный массив случайных чисел, в котором число строк определяется числом экспериментов (бросков), а число столбцов — числом испытаний (костей):

eyes = random.random_integers(1, 6, (N, ndice))

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

from numpy import random, sum

N = int(raw_input('Number of experiments: '))
ndice = int(raw_input('Number of dice: '))
nsix = int(raw_input('Number of dice with six eyes: '))

eyes = random.random_integers(1, 6, (N, ndice))
compare = eyes == 6
nthrows_with_6 = sum(compare, axis=1)  # суммирование по столбцам - элементам строки (axis = 1)
nsuccesses = nthrows_with_6 >= nsix
M = sum(nsuccesses)
p = float(M)/N
print 'probability:', p

И подсчет при этом для большого числа экспериментов, происходит существенно быстрее.

Шары из шляпы

В шляпе 12 шаров: четыре черных, четыре красных и четыре синих. И мы хотим написать программу, которая достает из шляпы три случайных шара. Шары у нас отличаются цветом, при этом эти цвета неизменны, значит их удобно представить в виде кортежа, а саму шляпу в качестве списка строк:

colors = 'black', 'red', 'blue'   # (кортеж строк)
hat = []
for color in colors:
    for i in range(4):
        hat.append(color)

Чтобы достать шар:

import random as random_number
color = random_number.choice(hat)
print color

Но нам нужно достать несколько шаров. При этом мы не возвращаем их обратно в шляпу, а это значит что элементы (шары) удаляются из списка hat. Три способа: 1) использовать hat.remove(color), то есть достали, например, синий шар - значит одним синим шаром в списке меньше, 2) выбирать шар по случайному индексу, а затем по нему же его удалять через del hat[index], 3) то же самое, но более коротким способом через hat.pop(index):

def draw_ball(hat):
    color = random_number.choice(hat)
    hat.remove(color)
    return color, hat

# или:
def draw_ball(hat):
    index = random_number.randint(0, len(hat)-1)
    color = hat[index]
    del hat[index]
    return color, hat

# или:
def draw_ball(hat):
    index = random_number.randint(0, len(hat)-1)
    color = hat.pop(index)
    return color, hat

balls = []
for i in range(n):  # n - число доставаемых шаров
    color, hat = draw_ball(hat)
    balls.append(color)
print 'Got the balls', balls

Продлим наше решение на более конкретную задачу, чем просто доставание из шляпы шаров: какова вероятность достать два и более черных шара из этой шляпы. В таком случае мы можем проделать N экспериментов и подсчитать сколько M раз мы достали не менее двух черных шаров и найти вероятность такого события как M/N. Один эксперимент для нас состоит в создании нового списка hat (перемешивании шаров), доставании шаров и подсчете числа черных. Последнее может быть легко подсчитано с помощью метода списков count, то есть в нашем случае hat.count('black'). Тогда оставим одну понравившуюся нам функцию draw_ball(hat) и допишем нашу программу следующим образом:

import random as random_number

def draw_ball(hat):
    color = random_number.choice(hat)
    hat.remove(color)
    return color, hat

def new_hat():
    colors = 'black', 'red', 'blue'
    hat = []
    for color in colors:
        for i in range(4):
            hat.append(color)
    return hat

n = int(raw_input('How many balls are to be drawn? '))
N = int(raw_input('How many experiments? '))

M = 0
for e in range(N):
    hat = new_hat()
    balls = []
    for i in range(n):
        color, hat = draw_ball(hat)
        balls.append(color)
    if balls.count('black') >= 2:
        M += 1
print 'Probability:', float(M)/N

Ограничение роста населения

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

Представим, у нас есть n индивидуумов, назовем их parents, которым мы придаем некоторые равномерно распределенные случайные значения. Далее из статистических данных мы знаем долю мужчин (male_portion), и всех индивидуумов, у которых их случайное число меньше male_portion, считаем мужчинами и присваиваем значение 1, а у кого выше — женщинами и присваиваем значение 2. Чтобы код было удобнее читать, введем константы MALE=1 и FEMALE=2. Наша задача в том, чтобы посмотреть как меняется массив parents для каждого нового поколения при учете двух указанных политик. Итак, создаем исходный массив:

from numpy import random, zeros
r = random.random(n)
parents = zeros(n, int)  # массив нулей размером n
MALE = 1; FEMALE = 2
parents[r <  male_portion] = MALE
parents[r >= male_portion] = FEMALE

Число потенциально возможных супружеских пар это минимум от числа мужчин и числа женщин. Однако, даже если мы посчитаем, что все эти пары сложатся, то все равно только часть из них дадут детей. При учете политики «одна семья — один ребенок» получаем:

males = len(parents[parents==MALE])  # число мужчин
females = len(parents) - males       # число женщин
couples = min(males, females)        # число потенциальных супружеских пар
n = int(fertility*couples)   # пары, родившие ребенка, fertility - рождаемость

# следующее поколение
r = random.random(n)
children  =  zeros(n, int)
children[r <  male_portion] = MALE
children[r >= male_portion] = FEMALE

Код, написанный для рождении нового поколения будет нужен при каждой новой генерации. Поэтому эти инструкции более естественно заключить в соответствующую функцию:

def get_children(n, male_portion, fertility):
    n = int(fertility*n)
    r = random.random(n)
    children = zeros(n, int)
    children[r <  male_portion] = MALE
    children[r >= male_portion] = FEMALE
    return children

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

# сначала обычное
children = get_children(couples, male_portion, fertility)
# в случае, если дочь:
daughters = children[children == FEMALE]
while len(daughters) > 0:
    new_children = get_children(len(daughters), male_portion, fertility)
    children = concatenate((children, new_children))
    daughters = new_children[new_children == FEMALE]

Накопленный опыт мы можем представить в виде общей функции анализа популяции для обеих политик:

def advance_generation(parents, policy='one child',
                       male_portion=0.5, fertility=1.0,
                       law_breakers=0, wanted_children=4):

    males = len(parents[parents==MALE])
    females = len(parents) - males
    couples = min(males, females)

    if policy == 'one child':
        # у каждой пары только один ребенок:
        children = get_children(couples, male_portion, fertility)
        max_children = 1
    elif policy == 'one son':
        # каждая пара продолжает рожать детей до тех пор, пока у них не будет сына

        children = get_children(couples, male_portion, fertility)
        max_children = 1
        daughters = children[children == FEMALE]
        while len(daughters) > 0:
            new_children = get_children(len(daughters), male_portion, fertility)
            children = concatenate((children, new_children))
            daughters = new_children[new_children == FEMALE]
            max_children += 1

    # кроме того часть граждан нарушает закон и имеет столько детей, сколько хочет (wanted_children)
    illegals = get_children(int(len(children)*law_breakers)*wanted_children, male_portion, fertility=1.0)
    children = concatenate((children, illegals))
    return children, max_children

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

from numpy import random, concatenate, zeros
MALE = 1;  FEMALE = 2

def get_children(n, male_portion, fertility):
    n = int(fertility*n)
    r = random.random(n)
    children = zeros(n, int)
    children[r <  male_portion] = MALE
    children[r >= male_portion] = FEMALE
    return children

def advance_generation(parents, policy='one child',
                       male_portion=0.5, fertility=1.0,
                       law_breakers=0, wanted_children=4):

    males = len(parents[parents==MALE])
    females = len(parents) - males
    couples = min(males, females)

    if policy == 'one child':
        # у каждой пары только один ребенок:
        children = get_children(couples, male_portion, fertility)
        max_children = 1
    elif policy == 'one son':
        # каждая пара продолжает рожать детей до тех пор, пока у них не будет сына

        children = get_children(couples, male_portion, fertility)
        max_children = 1
        daughters = children[children == FEMALE]
        while len(daughters) > 0:
            new_children = get_children(len(daughters), male_portion, fertility)
            children = concatenate((children, new_children))
            daughters = new_children[new_children == FEMALE]
            max_children += 1

    # кроме того часть граждан нарушает закон и имеет столько детей, сколько хочет (wanted_children)
    illegals = get_children(int(len(children)*law_breakers)*wanted_children, male_portion, fertility=1.0)
    children = concatenate((children, illegals))
    return children, max_children

N = 1000000
male_portion = 0.51
fertility = 0.92
law_breakers = 0.06
wanted_children = 6

generations = 10
# начнем с "идеального поколения" родителей:
start_parents = get_children(N, male_portion=0.5, fertility=1.0)
parents = start_parents.copy()
print 'one child policy, start: %d' % len(parents)
for i in range(generations):
    parents, mc = advance_generation(parents, 'one child',
                                     male_portion, fertility,
                                     law_breakers, wanted_children)
    print '%3d: %d' % (i+1, len(parents))

parents = start_parents.copy()
print 'one son policy, start: %d' % len(parents)
for i in range(generations):
    parents, mc = advance_generation(parents, 'one son',
                                     male_portion, fertility,
                                     law_breakers, wanted_children)
    print '%3d: %d (max children in a family: %d)' % (i+1, len(parents), mc)

Интегрирование

Одним из самых ранних применений генерируемых случайных чисел было вычисление интегралов. Пусть мы генерируем равномерно распределенные между a и b случайные числа x₁, ..., x_n, тогда аппроксимация решения будет найдена как

$\int _{a}^{b}f(x)dx={\frac {b-a}{n}}\sum _{i=1}^{n}f(x_{i})$

(18.1)

Этот метод обычно и называется интегрирование по Монте-Карло. Мы можем представить выражение (18.1) в виде небольшой программы:

import random as random_number

def MCint(f, a, b, n):
    s = 0
    for i in range(n):
        x = random_number.uniform(a, b)
        s += f(x)
    I = (float(b-a)/n)*s
    return I

Обычно, достаточная точность метода задается большим числом n, поэтому векторизованная версия будет более удобной:

from numpy import *

def MCint_vec(f, a, b, n):
    x = random.uniform(a, b, n)
    s = sum(f(x))
    I = (float(b-a)/n)*s
    return I

Рассмотрим интегрирование методом Монте-Карло на примере простой линейной функции $~f(x)=2x+3$, пределы интегрирования — от 1 до 2. Было бы интересно посмотреть как метод справляется с решением задачи для различных n. Оценку произведем следующим немного измененным MCint методом:

def MCint2(f, a, b, n):
    s = 0

    I = zeros(n)
    for k in range(1, n+1):
        x = random_number.uniform(a, b)
        s += f(x)
        I[k-1] = (float(b-a)/k)*s
    return I

Заметим, что k' изменяется от 1 до n, в то время как индексы I как и раньше идут от 0 до n-1. Поскольку n может быть очень большим, массив I может переполнить память. Поэтому следует записывать только каждое N-ое значение аппроксимации. Это возможно с помощью известной нам функции определения остатка:

for k in range(1, n+1):
    ...
    if k % N == 0:
        # store

Итак, каждый раз, когда k делится на N без остатка, мы записываем значение (в нашем случае каждое сотое). Соответствующая функция представлена ниже.

def MCint3(f, a, b, n, N=100):
    '''Cохраняет каждое N приближение интеграла
    в массив I и записываем соответствующее значение k'''
    s = 0

    I_values = []
    k_values = []
    for k in range(1, n+1):
        x = random_number.uniform(a, b)
        s += f(x)
        if k % N == 0:
            I = (float(b-a)/k)*s
            I_values.append(I)
            k_values.append(k)
    return k_values, I_values

Теперь у нас есть инструмент для того, чтобы посмотреть как изменяется ошибка в интегрировании методом Монте-Карло с ростом n. Законченная программа выглядит следующим образом, результат работы программы (может несколько отличаться ввиду случайности) представлен справа:

Результат работы программы

import random as random_number
import matplotlib.pyplot as plt
from numpy import array


def MCint3(f, a, b, n, N=100):
    '''Cохраняет каждое N приближение интеграла
    в массив I и записываем соответствующее значение k'''
    s = 0

    I_values = []
    k_values = []
    for k in range(1, n+1):
        x = random_number.uniform(a, b)
        s += f(x)
        if k % N == 0:
            I = (float(b-a)/k)*s
            I_values.append(I)
            k_values.append(k)
    return k_values, I_values

def f1(x):
    return 2 + 3*x

k, I = MCint3(f1, 1, 2, 1000000, N=10000)
error = 6.5 - array(I)

plt.title('Monte Carlo integration')
plt.xlabel('n')
plt.ylabel('error')
plt.plot(k, error)
plt.show()

Интегрирование через вычисление площади под кривой

Представим, у нас есть некоторая сложная фигура G , площадь которой мы хотим вычислить. Эту фигуру мы можем заключить в прямоугольник B определяемый координатами [X1, X2] x [Y1, Y2]. Далее в этот прямоугольник мы "бросаем" случайные числа из заданного массива [X, Y] и считаем сколько точек попало внутрь контура фигуры G. Далее, зная отношение числа точек, "упавших" в контур G к общему числу "брошенных" точек, и умножив это отношение на площадь простой фигуры B, мы получаем площадь сложной фигуры G.

Таким же образом мы могли бы находить интеграл от функции, поскольку он и определяет площадь под ней. Представим, мы знаем как выглядит функция, область интегрирования лежит в пределе от a до b, а максимум функции в указанном интервале равен m. Тогда аналогично, предыдущему рассмотрению, интеграл будет выражен через число "выброшенных" точек N и число упавших "под функцию" M как ${\frac {M}{N}}m(b-a)$.

Векторизованное решение может в простейшем случае выглядит так:

from numpy import random

def MCint_area_vec(f, a, b, m, N):
    x = random.uniform(a, b, N)
    y = random.uniform(0, m, N)
    M = y[y < f(x)].size   # в квадратных скобках условие, 
                           # size считает число ему удовлетворяющих элементов 
    area = M/float(N)*m*(b-a)
    return area

Решение можно написать и через стандартный random в обычном виде с помощью цикла, но скорость его будет существенно меньшей.