Как найти площадь фигуры методом монте карло

Печатать книгу

Расчет площади пересечения окружностей методом Монте-Карло

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

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

Обоснование

Задача расчета площади пересечения двух окружностей является тривиальной геометрической задачей (координаты центров окружностей и их радиусы нам известны). Площадь пересечения двух окружностей – это сумма площадей соответствующих сегментов этих окружностей. Есть решения для расчета площади пересечения двух, трех, четырех окружностей в различных частных случаях.

А вот решения общего случая для пересечения даже трех окружностей уже далеко не так тривиальны. В процессе поиска я нашел даже исследования по расчету площади пересечения N окружностей, однако они настолько же интересны, насколько и сложны.

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

Итак, алгоритм расчета площади любой фигуры методом Монте-Карло сводится к следующему:

1. Фигура вписывается в прямоугольник. Координаты сторон прямоугольника известны, значит, известна его площадь.
2. Псевдослучайным образом внутри прямоугольника генерируется большое количество точек. Для каждой точки определяется, попала ли точка внутрь исходной фигуры или нет.
3. В результате площадь исходной фигуры вычисляется исходя из обычной пропорции: отношение количества точек, попавших в фигуру, к общему количеству сгенерированных точек равно отношению площади фигуры к площади ограничивающего ее прямоугольника.

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

Реализация задачи на JavaScript

Рисование окружностей делалось средствами замечательной библиотеки D3.js. Алгоритм начального взаимного расположения окружностей выходит за рамки данной статьи, поэтому примем начальное расположение как данность.

var nodes = d3.selectAll("circle.node");

var squares = [];
var intersections = [];

nodes.each(function(node){
   // считаем радиус и площадь окружности
   var r = this.r.baseVal.value;
   var s = 3.14159*r*r;
   squares.push({node: node, square: s, r: r});

   // ищем пересечения пар окружностей
    nodes.each(function(node2){
        // расстояние между центрами и сумма радиусов
        var center_dist = Math.sqrt(Math.pow(node.x-node2.x, 2)+(Math.pow(node.y-node2.y, 2)));
        var radius_sum = r + this.r.baseVal.value;
        if(center_dist <= radius_sum && node.index != node2.index){
            // окружности пересекаются. проверить, что это пересечение найдено впервые
            node.r = r;
            node2.r = this.r.baseVal.value;
            if(isNewIntersection(intersections, node, node2)) 
                  intersections.push({node1: node, node2: node2, center_dist: center_dist});
        }
    });
});

var areaCalculator = {
  intersections: [], // массив пересечений, устанавливается снаружи
  frame: {}, // рамка вокруг фигуры
  circles: [], // массив окружностей
  figureArea: 0, // искомая площадь фигуры
  monteCarlo:
      function(p){
          // получаем массив окружностей из пересечения
          var circles = [];
          var x1_, y1_, x2_, y2_; // координаты сторон прямоугольника
          var inCirclesArr = function(node){
              for(var j=0; j<circles.length; j++){
                  if(circles[j].index==node.index){
                      return true;
                  }
              }
              return false;
          };
          for(var i=0; i<this.intersections.length; i++){
              if(!inCirclesArr(this.intersections[i].node1)){
                  circles.push(this.intersections[i].node1);
              }
              if(!inCirclesArr(this.intersections[i].node2)){
                  circles.push(this.intersections[i].node2);
              }
          }

          this.circles = circles;

          circles.sort(function(a,b){
              return a.x-a.r > b.x-b.r ? 1 : -1;
          });
          x1_ = circles[0].x-circles[0].r;

          circles.sort(function(a,b){
              return a.x+a.r < b.x+b.r ? 1 : -1;
          });
          x2_ = circles[0].x+circles[0].r;

          circles.sort(function(a,b){
              return a.y-a.r > b.y-b.r ? 1 : -1;
          });
          y1_ = circles[0].y-circles[0].r;

          circles.sort(function(a,b){
              return a.y+a.r < b.y+b.r ? 1 : -1;
          });
          y2_ = circles[0].y+circles[0].r;

          this.frame.x1 = x1_;
          this.frame.x2 = x2_;
          this.frame.y1 = y1_;
          this.frame.y2 = y2_;
          this.frame.area = (x2_-x1_)*(y2_-y1_);
          
          // рисуем прямоугольник
          paintRect(this.frame);          

          // p - количество генерируемых точек. В примере использовалось 100.000, чего хватило для приемлемой точности
          var p_positive = 0; // количество точек попавших в фигуру

          // генерируем p точек для определения площади фигуры
          for(var i=0; i<p; i++){
              var x_rand = Math.random()*(x2_-x1_)+x1_;
              var y_rand = Math.random()*(y2_-y1_)+y1_;

              var yes = false;
              for(var j=0; j<circles.length; j++) {
                  if(!yes && (
                        (circles[j].x-circles[j].r) <= x_rand &&
                        (circles[j].x+circles[j].r) >= x_rand &&
                        (circles[j].y-circles[j].r) <= y_rand &&
                        (circles[j].y+circles[j].r) >= y_rand )
                    ){
                     yes = true;
                     p_positive++;
                  }
              }
          }

          // площадь фигуры = площадь прямоугольника*кол-во точек внутри фигуры / общее кол-во точек
          this.figureArea = this.frame.area*p_positive/p;
      }
};

Пара гвоздей в метод Бутстрапа

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

Заключение

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

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

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

Что такое метод Монте-Карло

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

Название этого метода выбрано в честь города Монте-Карло, известного своими казино. Рулетка в казино – это своего рода генератор случайных чисел. Этот метод впервые был описан в 1949 году в совместной работе математиков Николаса Метрополиса и Станислава Улама. Известно, что название придумал Н. Метрополис в честь одного из своих родственников, который был азартным игроком и частым посетителем казино.

Рис. 1. Рулетка.

Использование метода Монте Карло для вычисления площади фигуры

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

Рис. 2. Фигура сложной формы.

Чтобы вычислить площадь S такой фигуры, ее вписывают в квадрат стороной a. Затем, случайным образом выбирают точки внутри квадрата. Причем точки могут попадать как внутрь исследуемой фигуры сложной формы, так и вне ее.

Чем большую площадь занимает фигура внутри квадрата, тем чаще в нее будут попадать точки. Точки, попавшие внутрь фигуры, обозначают буквой m, а все выбранные точки внутри квадрата – буквой n.

Здесь проявляется закономерность: чем больше точек в эксперименте задействовано, тем с большей вероятностью можно утверждать, что процент точек, содержащихся в исследуемой фигуре стремится к отношению площади фигуры к площади квадрата. Математически это можно записать следующим образом: s=a*a*(m/n).

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

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

Алгоритм вычисления площади фигуры с помощью метода Монте Карло

Определим площадь фигуры, ограниченной уравнениями y = x^2, y = 0 и x = 1.

Составим модель: площадь искомой фигуры впишется в квадрат. Если значения координат x и y принадлежит диапазону [0;1], то точка принадлежит квадрату. Если выполняется система неравенств 0 < x < 1 и 0 < y – x^2, то точка принадлежит фигуре.

Так как длина стороны квадрата равна 1, то после подсчета количества точек, площадь фигуры определится по формуле: s=m/n

Рис. 3. Фигура, ограниченная уравнением y=x^2, и прямыми y=0 и x=1.

Алгоритм решения задачи будет выглядеть так:

• Задать n (определиться, сколько испытаний будет в задаче)
• m=0 (начальное значение точек, попавших в область фигуры. Переменная m будет выполнять функцию счетчика в алгоритме)
• Задать значение x путем генерации случайных чисел
• Задать значение y путем генерации случайных чисел
• Организовать повтор цикла (перебирать целочисленные значения от 0 до n с шагом 1). Тело цикла:
• проверить условие 0 < x < 1, если оно выполняется перейти к проверке второго условия, иначе увеличить счетчик цикла на 1
• проверить условие 0 < y – x^2, если условие верно, то организовать команду счетчика точек (m=m+1), иначе увеличить счетчик цикла на 1
• s=m/n
• Вывод результата s

Что мы узнали?

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

Предыдущая

ИнформатикаИстория развития ЭВМ кратко

Следующая

ИнформатикаКлассы логических элементов кратко

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

Метод Монте-Карло,
или метод статистических испытаний, —
это численный метод, основанный на
моделировании случайных величин и
построении статистических оценок для
искомых величин.

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

Такой метод
приближенного нахождения площадей
фигур и носит название метода Монте-Карло.

Пример. Вычисление
числа π
методом
Монте-Карло.

Постановка
задачи: для
вычисления числа π методом Монте-Карло
рассмотрим круг радиуса 1 с центром в
точке (1, 1). Круг вписан в квадрат, сторона
которого, а=2. Тогда площадь квадрата
S_{квадрата}=
a
²
= 2² =
4.

Решение.

Выбираем внутри
квадрата N
случайных точек. Выбрать точку означает
задать ее координаты – числа x
и y.

Обозначим N_круга
– число точек попавших при этом внутрь
круга.

Точка принадлежит
квадрату, если 0≤x≤2
и 0≤y≤2.

Если
(x-1)²
+(y-1)²
≤
1, то точка
попадает в круг, иначе она находится
вне круга. Геометрически очевидно, что

Отсюда

То есть для круга
единичного радиуса:

Но для круга
единичного радиуса
,
следовательно получаем:
.

Данная формула
дает оценку числа π. Чем больше N,
тем больше точность этой оценки. Следует
заметить, что данный метод вычисления
площади будет справедлив только тогда,
когда случайные точки будут не просто
случайными, а еще и равномерно разбросанными
по всему квадрату.

Для моделирования
равномерно распределенных случайных
чисел в интервале от 0 до 1 в языке
программирования Turbo
Pascal
используется датчик случайных чисел –
функция RANDOM,
которая выдает последовательность
случайных величин, равномерно
распределенных от 0 до 1.

Таким образом,
суть компьютерного эксперимента
заключается в обращении к функции RANDOM
для получения N
раз координат х
и у
точки. При этом определяется, попала ли
точка с координатами (х,у)
в круг единичного радиуса. В случае
попадания значение величины N_круга
увеличивается
на 1.

Программа:

Program
monte_karlo;

var

i, n, n1 : LongInt;
x, y, pi : real;
begin
Randomize;

WriteLn(‘Введите
количество точек n=’);

Readln(n);

for i:=1 to n do
begin
x:=2*Random;
y:=2*Random;

if sqr(x-1)+sqr(y-1)<=1 then n1:=n1+1;
end;

pi:=4*n1/n;
WriteLn(‘pi=’, pi:15:11);
end.

Соседние файлы в предмете [НЕСОРТИРОВАННОЕ]

• #
• #
• #
• #
• #
• #
• #
• #
• #
• #
• #

 В начале урока рассказываю идею метода.
Ребята, сегодня мы рассмотрим интересный метод
приближенного вычисления площадей фигур – метод
Монте – Карло. Пусть у нас есть какая – нибудь
фигура на плоскости, площадь которой ( S_fig )
нам необходимо найти. Ограничим ее другой
фигурой, площадь которой ( S_total ) мы можем
легко вычислить. Например, прямоугольником АСDB
со сторонами, параллельными координатным осям
(см. рис. 1). И пусть про любую точку прямоугольника
мы можем быстро узнать, попадает эта точка внутрь
фигуры, площадь которой мы ищем, или нет.

А теперь начнем опыт – будем бросать
на бумагу зерна случайным образом (вообще – то
это нелегко сделать, чтобы обеспечить
случайность). Когда нам покажется, что зерна
почти полностью покрыли бумагу, посчитаем,
сколько всего зерен на прямоугольнике (пусть их
число N_total ) и сколько из них на фигуре ( N_fig
). Ясно, что число зерен, попавших внутрь фигуры,
пропорционально ее площади: больше площадь –
больше зерен, меньше площадь – меньше зерен.
Поэтому, поделив количество зерен , попавших
внутрь фигуры , на количество всех зерен в
прямоугольнике, мы сможем найти, какую часть
площади прямоугольника занимает фигура:

N_fig / N_total » S_fig
/ S_total , отсюда S_fig » N_fig
/ N_total ? S_total

Промоделируем этот опыт на ЭВМ.
Предположим, нам надо найти площадь фигуры
ограниченной сверху кривой Y = F(X), а снизу – осью
абсцисс. Пусть Y = cos(X), а Х I [– p /2, p /2] (см. рис. 2).
Ограничим нашу фигуру прямоугольником АСDB, его
площадь равна p .

Из чего должен состоять алгоритм:

1. Бросание зерна – бросание случайной
точки, координаты X и Y которой случайны, причем Х
должна меняться от – p /2 до p /2, а Y – от 0 до 1. И в
этих интервалах X и Y должны появляться с
одинаковой вероятностью в любой точке этих
отрезков, т.е. X и Y должны быть равномерно
распределены по осям. Тут надо напомнить ребятам,
как получить равномерно распределенное
случайное вещественное число на интервале [А, В]:

Х = random* ( B – A) + A

2. Надо определить, куда попала точка –
под кривую или выше нее. И вести подсчет N_fig.
Условие попадания точки под кривую: Y ? sin (X).

3. Повторить пп.1 и 2 столько раз, чтобы
получить желаемую точность результатов.

Дети без труда напишут программу по
данному алгоритму (на практическом занятии). И
после этого наступает момент разочарования.
Сухие цифры…

Предложение учителя: давайте создадим
проект в Delphi, который наглядно демонстрировал бы
работу метода Монте – Карло. Разработаем
интерфейс программы. Тут ребята сами предлагают
несколько вариантов оформления проекта, один из
которых представлен на рис. 3.

Разместим на форме следующие
компоненты:

Edit – окно редактирования для ввода
общего количества испытаний (бросков зерен) – N_total;

Button – кнопка для запуска работы
метода Монте – Карло;

Panel – панель для вывода посчитанной
площади фигуры;

(все вышеперечисленные компоненты
расположены на вкладке Standart Палитры компонентов)

Image – для вывода точек, попавших в
искомую область (компонент расположен на вкладке
Additional Палитры компонентов).

Если необходимо, надо напомнить
ребятам некоторые свойства и методы выбранных
компонентов (это может быть и домашним заданием
предыдущего урока, чтобы не терять время на
“воспоминания”) .

Компонент Edit:

Свойство Text – содержит текст, который
пользователь набирает в окне Edit. Этот текст надо
преобразовать в число N_total. Для этого в Delphi
есть необходимая функция StrToInt:

N_total = StrToInt (Edit1.Text).

Компонент Panel:

Свойство Caption – содержит текст,
который выводится на панель. Чтобы вывести
полученное число S_fig на панель, мы должны
преобразовать его в строку S с помощью процедуры
Str:

Str(S_fig:10:4, S),

а потом вывести эту строку на панель
следующим образом:

Panel1.Caption := ‘Площадь фигуры = ‘ + S

Компонент Image:

Свойства Height и Width – соответственно
высота и ширина компонента;

Свойство Canvas – отводит канву (место)
для рисования на компоненте Image;

Методы Canvas:

FillRect(ClientRect) – закрашивает область
клиента компонента Image каким – либо цветом

(по – умолчанию – белым), т.е. стирает
предыдущую картинку;

MoveTo(X, Y) – перемещает перо в точку с
координатами X, Y без проведения линии (координаты
задаются в пикселях);

LineTo(X, Y) – проводит линию из текущей
точки в точку с координатами X, Y;

Pixels[I, J] – содержит цвет точки с
координатами I, J.

У кнопки (компонент Button) мы будем
обрабатывать событие onClick (событие нажатие
кнопки). Т.е. вышеописанный алгоритм мы
программируем в процедуре Button1Click.

Тут учитель задает вопрос классу: “
Какая проблема возникает при выводе точки на
экран (на компонент Image)?” Ответ: расчетные
координаты очень малы (0 ? Y ? 1,меньше пикселя, -p /2 ? X
? p /2); а если взять другую кривую Y = F(X), они могут
оказаться слишком большими (больше размера
компонента Image). Поэтому, при выводе значений
функций (графиков) на экран монитора, необходимо
преобразовывать расчетные координаты в
графические с учетом дискретности растровой
сетки монитора, а также предусмотреть
возможность автоматического масштабирования
функции (графика) по осям координат. Для этого
желательно создать отдельную подпрограмму.

Для полного размещения функции
(графика) в расчетной области (это область
компонента Image) необходимо определить X_min, X_max, Y_min,
Y_max – минимальные и максимальные значения по X и
по Y соответственно. X_min =А, X_max = В. Как найти Y_min,
Y_max.

Коллективно обсуждается следующий
алгоритм:

1. Разобьем интервал [А, В] по Х на N
равных частей и определим массивы значений
аргумента и функции X[i] и Y[i] = F(X[i]), где I = 1..N;
2. Определяем наибольшее Y_max и наименьшее Y_min
значения функции в заданном интервале изменения
аргумента;
3. Находим коэффициенты масштабирования Kx, Ky при
построении графика в заданной области;
4. Т.к. коэффициенты масштабирования Kx, Ky могут
отличаться, то выводимый график может
искажаться. Устраняем искажения графика;
5. Преобразуем расчетные координаты точки X, Y в
графические Xg, Yg. С учетом того, необходимости
“переворота” оси Y, которая в координатах
монитора направлена сверху вниз.

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

1. Применить метод Монте – Карло для
приближенного вычисления площади фигуры,
ограниченной сверху кривой Y = sin (X), при Х I [ 0, p ];
2. Применить метод Монте – Карло для
приближенного вычисления числа p . Подсказка:
рассмотреть круг единичного радиуса с центром в
т. (1, 1). Его площадь и будет равна p .
3. Применить метод Монте – Карло для
приближенного вычисления площади фигуры,
ограниченной сверху кривой Y = sin (X), при Х I [ 0, 2p ];
4. Применить метод Монте – Карло для
приближенного вычисления площадей фигур,
представленных на рис. 5 – 7.



5. Доработать проект:

а) организавать проверку правильности ввода
информации в поле Edit (чтобы вводились только
целые числа);
б) разметить оси и подписать числовые значения.

На последующих уроках, на которых
предполагается изучение тем “Вычисление
площадей (интегралов) методом трапеций и методом
прямоугольников”, можно предложить ребятам
доработать проект, поместив на форму
дополнительные компоненты Image, Button, Edit, Panel (для
каждого численного метода – свои). В окно
компонента Edit пользователь будет вводить
количество разбиений интервала [А, В] по Х. Таким
образом, ребята смогут сравнить и наглядно
увидеть работу всех трех численных методов.

Листинг программы

{Нахождение максимального и
минимального значений функции и аргумента}

X_max:=B; X_min:=A;
Y_max:=Y[1]; Y_min:=Y[1];
for i:=2 to n do begin
if Y_max < Y[i] then Y_max:=Y[i];
if Y_min > Y[i] then Y_min:=Y[i];
end;

{Нахождение коэффициентов сжатия}

Kx:=right/(X_max-X_min);
Ky:=down/(Y_max-Y_min);
{Устранение искажения графика}
if Kx < Ky then begin
Ky:=Kx;
down:=round((Y_max-Y_min)*Ky);
end
else begin
Kx:=Ky;
right:=round((X_max-X_min)*Kx);
end;
end;