Программирование и научные вычисления на языке 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 xrange(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 1*, ...,

*x*, тогда аппроксимация решения будет найдена как

n(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
```


Рассмотрим интегрирование методом Монте-Карло на примере простой линейной функции , пределы интегрирования — от 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* как .

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



```
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 в обычном виде с помощью цикла, но скорость его будет существенно меньшей.