Метод Монте-Карло

Вычисления по методу Монте-Карло с использованием языка программирования Python

Монте-Карло казино.

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

Определение орбиты по измерениям предполагает нахождение оценок векторов координат и скоростей объекта, исходя из некоторого их "нулевого" приближения и измерений (например, измерений дальности до объекта с Земли). В зависимости от того, насколько удачным было выбрано начальное приближение, зависело будет получена адекватная оценка параметров движения объекта или нет — иными словами — сойдется алгоритм оценивания или нет. Множество таких точек, вокруг неизвестного набора точных параметров, для которых алгоритм сходился достаточно близко к точному решению, представляет особый интерес — это область сходимости алгоритма.

Эта область располагалась в шестимерном пространстве — пространстве трехмерных координат и скоростей объекта, а ее двумерные сечения, найденные при помощи вычислительных экспериментов, свидетельствовали, что ее форма очень причудлива, и вряд ли может быть точно описана какой-либо аналитикой.

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

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

С теоретико-вероятностных позиций метод представлял оценку вероятности в рамках схемы Бернулли по результатам эксперимента. Исходом в этой схеме была вероятность сходимости для случайно выбранной точки в кубе допустимых значений параметров. $$P = \frac{V_{\mbox{сход}}}{V_{\mbox{куб}}}\approx\frac{M}{N}, $$ где \(M, N\) — число точек, в которых наблюдалась сходимость алгоритма и общее их число соответственно. Из этой формулы видно, что точность оценки неизвестной вероятности \(P\) определяет точность определения объема области сходимости.

Формулы для интервальных оценок вероятности по результатам испытаний в схеме Бернулли следующие: $$\begin{array}{l}
\Delta = \frac{\Phi_{0}^{-1}(\alpha/2+1/2)}{2\sqrt n}\\
P = \frac{k}{n}\pm\Delta\end{array}
$$

В этой формуле \(k\) — наблюдаемое число успехов в схеме Бернулли; \(n\) — общее число испытаний\наблюдаений; \(\Phi_0(x)\) — функция стандартного нормального распределения (т.е. с нулевым математическим ожданием и единичной дисперсией).

Таким образом, если \(V_1\) подлежащий оценке неизвестный объем, целиком содержащийся, например, в кубе известного объема \(V\), а на основе вычислительного эксперимента, заключающегося в проверке принадлежности случайно выбранных в объеме \(V\) точек объему \(V_1\) удается найти отношение \(\frac{k}{n}\), где \(k\) — число точек, выпавших в \(V_1\), \(n\) — общее число точек, то на его основе удается построить оценку объема \(V_1\): $$ V_1 = V\dfrac{k}{n}\pm V\Delta$$

Это и есть базовая схема метода Монте-Карло. Вообще, в современном статистическом анализе методами Монте-Карло называют разнообразные вычислительные подходы, опирающиеся на целеноправленно проводимые эксперименты по имитации случайности. Это достаточно расплывчатое определение, но оно отражает многообразие проявлений данного подхода. Такие подходы правильней представлять как составляющие единой концепции — Монте-Карло — имитации случайности. Название подхода (концепции) происходит от названия (ныне) района в княжестве Монако Monte-Carlo (означает гора Карло), где находится первое в Европе казино. Рулетка всегда считалась генератором случайности, отсюда и происхождение названия подхода.

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

Площадь пресечения треугольников на плоскости. Пример вычислений на Python.

Метод Монте-Карло для вычисления площади пересечения треугольников

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

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

 

# -*- encoding: utf-8 -*-

from __future__ import print_function
import numpy as np

def point_in_triangle(p1,p2,p3, p):
    '''Проверка принадлежности точки p треугольнику, 
     заданному точками (p1, p2, p3).
    '''
    e1 = p2 - p1
    e2 = p3 - p1
    # Точка внутри треугольника если сумма ее коррдинат в базисе двух сторон меньше 1
    res = np.dot(np.linalg.pinv(np.vstack([e1, e2]).T), p-p1)
    if (res >= -4.0e-16).all() and res.sum() <= 1:
        return True
    else:
        return False



def estimate_intersection_area(p1i, p2i, p3i, n=100000):
    cmax = np.max(np.vstack([p1i, p2i, p3i]), axis=0)
    cmin = np.min(np.vstack([p1i, p2i, p3i]), axis=0)
    a = np.random.uniform(cmin[0], cmax[0], size=n)
    b = np.random.uniform(cmin[1], cmax[1], size=n)
    res = 0
    # Пробегаем по всем пробным точкам
    for x, y in zip(a,b):
        # подсчитываем те, которые принадлежат всем треугольникам сразу
        if all([point_in_triangle(p1, p2, p3, np.array([x,y])) for p1, p2, p3 in zip(p1i, p2i, p3i)]):
           res+=1
    return (cmax[1]-cmin[1])*(cmax[0]-cmin[0])*res/float(n)



if __name__ == '__main__':
    p1i = np.array([[0.0, 0.5],
                    [0.0, 1.0],
                    [0.0, 2.0]])
    p2i = np.array([[1.0, 0.0],
                    [1.0, 0.0],
                    [1.0, 0.0]])
    p3i = np.array([[0.0, 0.0],
                    [0.0, 0.0],
                    [0.0, 0.0]])
    # В данном случае задан набор треугольников:
    # [0,0.5], [1,0], [0,0]
    # [0,1], [1,0], [0,0]
    # [0,2], [1,0], [0,0]
    print('Точное значение площади: 0.25', 'вычисленное:', estimate_intersection_area(p1i, p2i, p3i))
    

Представленная реализация метода Монте-Карло для вычисления площади пересечения треугольников позволяет рассчитать искомую площадь (0.25 - в примере) достаточно точно, но для этого требуется 100000 зондирующих точек. Время расчетов здесь весьма ощутимо, т.к. это пример — одна из самых неэффективных реализаций метода. Циклы на Python медленные, а здесь выполняется цикл по всем зондирующим точкам, и, как следствие, возрастает время вычислений. Код можно сделать быстрее, организовав его так, чтобы избежать длинных циклов на Python. Но это уже другая задача (ее решение возможно, либо на базе более аккуратной организации работы с массивами с использованием numpy, либо за счет реализации циклов на Си, например, применяя Cython).

blog comments powered by Disqus