Выборочный метод Монте-Карло

Показывать лекцию целиком

Теоретическая часть

Сущность метода Монте-Карло состоит в следующем: требуется найти значение $$а$$ некоторой изучаемой величины. Для этого выбирают такую случайную величину $$Х$$, математическое ожидание которой равно $$а$$, т. е. $$М(Х) = а$$ [6].

Практически же поступают следующим образом: производят $$n$$ испытаний, в результате которых получают $$n$$ возможных значений величины $$Х$$ ; вычисляют их среднее арифметическое

$$\bar{x}=\frac{1}{n}\sum\limits_{i=1}^{n} x_i$$

и принимают $$\bar{x}$$ в качестве оценки (приближенного значения) $$а*$$ искомого числа $$а$$: $$a\approx a^{*} =\bar{x}$$.

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

1. Оценка погрешности метода Монте-Карло

Как отмечалось, для получения оценки математического ожидания случайной величины $$Х$$ необходимо произвести $$n$$ независимых испытаний и по ним найти выборочную среднюю, которая принимается в качестве искомой оценки. При каждой конечной серии испытаний будут получаться различные значения случайной величины и, следовательно, другая средняя, а значит, и другая оценка математического ожидания. Очевидно, что получить точную оценку математического ожидания невозможно. Поэтому возникает вопрос о допускаемой ошибке. Обычно ограничиваются отысканием лишь верхней границы $$\delta$$ допускаемой ошибки с заданной вероятностью $$\gamma$$, т. е.

$$P(|\overline{X}-a|\le\delta)=\gamma.$$

При этом возможны следующие случаи оценки числа испытаний:

Случайная величина $$Х$$ распределена нормально и ее среднее квадратическое отклонение (стандартное отклонение) $$\delta$$ известно. В этом случае с заданной вероятностью $$\gamma$$ верхняя граница ошибки определяется по формуле

$$\delta=\frac{x\sigma}{\sqrt{N}},$$

где:

$$N$$ — число испытаний (разыгранных значений случайной величины $$Х$$ );

$$х$$ — значение аргумента функции Лапласа $$Ф(х)$$ или интеграла вероятности, при котором она равна половине заданной вероятности;

$$\sigma$$ — известное среднее квадратическое отклонение.

Из формулы (4.1) может быть найдено число испытаний. Один из вариантов интеграла вероятностей (функции Лапласа) имеет вид [6]

$$Ф(х) = \frac{1}{\sqrt{2\pi}}\int\limits_{0}^{x}e^{-t^2 /2} dt.$$

Значения $$Ф(х)$$ табулированы и приведены в большинстве учебников по теории вероятностей и математической статистике.

В зарубежной литературе большое распространение получила так называемая функция ошибок ( $$error function$$ ) $$erf$$:

$$erf(х) = \frac{2}{\sqrt{\pi}}\int\limits_{0}^{x}e^{-t^2} dt.$$

Связь между функцией ошибок (4.3) и интегралом вероятностей (4.2) выражается в виде

$$Ф(х) = 0.5 erf(x/\sqrt{2}).$$

Случайная величина $$Х$$ распределена нормально, причем ее среднее квадратическое отклонение неизвестно. В этом случае с заданной вероятностью $$\gamma$$ верхняя граница ошибки вычисляется по формуле

$$\delta=\frac{t_{\gamma}s}{\sqrt{N}},$$

где:

$$N$$ — число испытаний;

$$s$$ — "исправленное" среднее квадратическое отклонение;

$$t_{\gamma}$$ находят по специальным таблицам, например, приведенной в [6].

Из формулы (4.5) может быть найдено число испытаний для определения верхней границы ошибки.

Случайная величина $$Х$$ распределена по закону, отличному от нормального. В этом случае при достаточно большом числе испытаний ( $$n > 30$$ ), с вероятностью, приближенно равной $$\gamma$$ (заданной вероятностью), верхняя граница ошибки может быть вычислена по формуле (4.1), если среднее квадратическое отклонение случайной величины известно; если же оно неизвестно, то можно подставить в формулу (4.1) его оценку — "исправленное" среднее квадратическое отклонение — либо воспользоваться формулой (4.5). При этом чем больше число испытаний, тем меньше различие между результатами, которые дают обе формулы.

2. Вычисление кратных интегралов методом Монте-Карло

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

Сначала рассмотрим вычисление простого определенного интеграла

$$I=\int\limits_{a}^{b}f(x)dx$$

где $$b \ne a$$.

Введем под знак интеграла постоянный множитель (равный единице):

$$\frac{b-a}{b-a},\\I=\int\limits_{a}^{b}f(x)\frac{b-a}{b-a}dx.$$

Вынесем из-под интеграла числитель дроби, получим

$$I=(b-a)\int\limits_{a}^{b}f(x)\frac{1}{b-a}dx.$$

Как известно, если случайная величина $$Х$$ распределена на заданном интервале (например, $$b – a$$ ) равномерно, то ее функция плотности $$p(x)$$ обратно пропорциональна длине интервала, т. е.

$$p(x)=\frac{1}{b-a}.$$

Кроме того, если известно распределение случайной величины $$Х$$, то функция от этой случайной величины $$f(x)$$ будет иметь тот же самый закон распределения. В этом случае математическое ожидание $$M[x]$$ непрерывной равномерно распределенной случайной величины рассчитывается по формуле

$$M[x]=\int\limits_{a}^{b}x\frac{1}{b-a}dx.$$

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

$$M[f(x)]=\int\limits_{a}^{b}f(x)\frac{1}{b-a}dx.$$

Сопоставляя (4.8) и (4.9), приходим к выводу, что определенный интеграл может быть рассчитан по формуле

$$I=(b-a)M[f(x)].$$

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

$$M[f(x)]\approx \frac{1}{n}\sum\limits_{i=1}^{n}f(x_i).$$

С учетом (4.11) получаем выражение для приближенного расчета определенного интеграла

$$I\approx \frac{b-a}{n}\sum\limits_{i=1}^{n}f(x_i).$$

Чем больше число испытаний $$n$$, тем точнее будет расчет математического ожидания (4.11) и, следовательно, определенного интеграла, вычисляемого по формуле (4.12).

Рассмотрим общий подход вычисления $$m$$ -кратного интеграла с помощью метода Монте-Карло [5].

Пусть задан $$m$$ -кратный интеграл вида

$$I=\mathop{\int\int ... \int}\limits_{D} f(x_1,x_2,...,x_m)dx_1dx_2...dx_m,$$

где подынтегральная функция $$f(x)$$ задана на замкнутой области $$D \subset R^m$$.

Погрузим область интегрирования $$D$$ в $$m$$ -мерный промежуток

$$D^* = [a_1, b_1]\times[a_2, b_2]\times ... [a_m, b_m],$$

имеющий меру

$$\mu(D^*)=\prod\limits_{j=1}^{m}(b_j-a_j).$$

Определим в промежутке (4.15) функцию

$$g(x)=\begin{cases} f(x),x\in D;\\ 0, x \in D^* \setminus D.\\ \end{cases}$$

Тогда в соответствии с (4.13) и (4.16) получим

$$I=\mathop{\int\int ... \int}\limits_{D^*} g(x_1,x_2,...,x_m)dx_1dx_2...dx_m,$$

Введем в рассмотрение $$m$$ -мерную случайную величину $$Х$$, имеющую в замкнутой области равномерное распределение вероятностей с дифференциальной функцией плотности

$$p(x)=\frac{1}{\mu (D^*)}$$

Функция плотности равномерного распределения есть величина постоянная, поэтому введем ее под знак интеграла (4.17) следующим образом:

$$I=\mathop{{\int\int ... \int} g(x_1,x_2,...,x_m)}\limits_{D^*}\frac{\mu (D^{*})}{\mu (D^{*})}dx_1dx_2...dx_m,$$

Вынесем числитель дроби за знак интеграла, т. е.

$$I=\mu (D^{*})\mathop{{\int\int ... \int} g(x_1,x_2,...,x_m)}\limits_{D^*}\frac{1}{\mu (D^{*})}dx_1dx_2...dx_m,$$

В (4.20) $$m$$ -кратный интеграл — это математическое ожидание от функции $$g(x)$$ случайной величины в предположении, что случайная величина $$Х$$ распределена равномерно с плотностью (4.18). Следовательно, можем записать

$$I=\mu (D^{*})M[g(x)],$$

где $$х = (х_{1}, х_{2}, ... , х_{m})$$.

В свою очередь математическое ожидание может быть оценено с помощью арифметического среднего. Тогда приближенное значение $$m$$ -кратного интеграла будет определяться приближенной формулой

$$I\approx \frac{\mu (D^{*})}{N}\sum\limits_{k=1}^{N}g(x^{(k)}),$$

где $$x^{(k)} \in D^{*}$$ — значение случайной величины $$Х$$ в $$k$$ -м испытании.

Чтобы смоделировать выборку $$m$$ -мерной случайной величины $$Х$$, равномерно распределенной в $$m$$ -мерном промежутке $$D^{*}$$, используются псевдослучайные числа. Для этого в каждом испытании с номером $$k$$ выбирают $$m$$ псевдослучайных чисел $$r_i^{(k)},\mbox{ }i=\overline{1,m}$$, и по ним определяют координаты случайной величины $$x_i^{k}=a_i+(b_i-a_i)r_i^{(k)},\mbox{ }i=\overline{1,m}$$, псевдослучайной точки $$x^{(k)}$$ [5].

Таким образом, техника применения метода Монте-Карло здесь будет заключаться в определении области $$D^{*}$$, генерировании в ней псевдослучайных чисел, подсчета числа попаданий этих чисел в область $$D$$ и применении формулы (4.22).

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

Практическая часть

Пример 1. Рассчитайте с помощью метода Монте-Карло площадь круга радиусом $$r = 5$$ и с координатами центра $$(1; 2)$$. Число испытаний примите 500. Выполните графические построения, поясняющие расчет.

Для определения площади плоской фигуры ее вписывают в соответствующую известную фигуру, площадь которой достаточно просто вычисляется. Вычисление площади круга может быть произведено через вычисление площади квадрата, в который вписывается круг. Выборочный метод Монте-Карло предполагает генерирование случайных (псевдослучайных) равномерно распределенных чисел. Этим числам сопоставляются координаты точек для рассматриваемой фигуры — квадрата, в которую вписывается круг. Если площадь квадрата $$S$$ подсчитана, то задается число испытаний $$N$$, из которых $$m$$ исходов могут оказаться внутри круга или на его границе, т. е. на окружности. Тогда площадь круга $$S0$$ будет определяться выражением

$$S0=\frac{m}{N}S,$$

где:

$$N$$ — число сгенерированных случайных чисел, соответствующих количеству точек, находящихся в данном квадрате:

$$m$$ — число случайных чисел (точек), которые попали в круг. Граница круга — это окружность, уравнение которой имеет вид

$$(x-x0)^2+(y-y0)^2=r^2,$$

где:

$$х0$$, $$y0$$ — координаты центра окружности;

$$х$$, $$y$$ — текущие значения переменных;

$$r$$ — радиус окружности.

Программный код решения примера:

clear all, clc,close all
%%%%%%% Расчет площади круга методом Монте-Карло
%% Параметры
r = 5;
x0 = 1;
y0 = 2;
%%% Построение окружности
t = 0 : r/500 : 2*pi;
x = r*cos(t) + x0;
y = r*sin(t) + y0;
line(x, y, 'linew', 2, 'color','r')
%%% Центр круга
line(x0,y0,'marker','o','markerfacecolor','r','color', 'r' )
 
%%% Построение квадрата, в который вписан круг
line([x0-r, x0+r],[y0+r, y0+r], 'lines','-.')
line([x0-r, x0+r],[y0-r, y0-r], 'lines','-.')
line([x0+r, x0+r],[y0-r, y0+r], 'lines','-.')
line([x0-r, x0-r],[y0-r, y0+r],  'lines','-.')
 
%%% Гененрирование случайных чисел для области D*
N = 500; %%% число испытаний
%%% Генерация чисел по горизонтальной стороне квадрата
rx = min(x) + (max(x) - min(x))*rand(N,1);
%%% Генерация чисел по вертикальной стороне квадрата
ry = min(y) + (max(y) - min(y))*rand(N,1);
%%%% Заполнение случайными числами квадрата с кругом
for J = 1 : N
line(rx(J),ry(J),'marker','o','markersize',2,'color', 'k',...
 'markerfacecolor', 'k' )
end
%%% Подсчет количества случайных чисел, попавших в круг
m = 0;
for J = 1 : N
    if (rx(J) - x0)^2 + (ry(J) - y0)^2 <= r^2
        m = m + 1;
    end
end
%%%%% Расчет площади квадрата
S = ((x0+r) - (x0-r))^2;
 
%%% Расчет площади круга методом Монте-Карло
S0 = m/N*S

%%% Проверка расчета 
Scontrol = pi*r^2

%%% Заголовок для диаграммы
str = sprintf('%s Радиус круга r = %g, координаты центра x_0 = %g,  y_0 = %g.%s%g%s%g. %s%g', '\bf',r,x0,y0, ... 
'\newline Теоретическая площадь круга S = ', Scontrol, ...
'\newline Площадь круга по методу Монте-Карло S0 = ', S0, ...
 'Число испытаний N = ', N);
title(str)
grid on
xlabel('\bf - - - - - - - X  - - - - - - - ')
ylabel('\bf - - - - - - - Y  - - - - - - - ')
%%% Ограничения по осям
xlim([x0 - r - r/10, x0 + r + r/10])
ylim([y0 - r - r/10, y0 + r + r/10])
axis equal

Задание 1

Видоизмените программу так, чтобы выполнялось условие непревышения относительной погрешности заданной величины (например, 1%) и чтобы производился необходимый подсчет числа испытаний.

Предусмотрите, чтобы точки, попавшие в круг, были красного цвета, а не попавшие в него – синего цвета.

Напишите программу расчета числа $$\pi$$ по методу Монте-Карло. Расчет числа $$\pi$$ произведите с точностью в четыре знака после десятичной точки.

Пример 2. Рассчитайте с помощью метода Монте-Карло площадь фигуры, ограниченной эллипсом с большой осью $$a = 6$$ и малой осью $$b = 4$$. Координаты пересечения осей эллипса $$x0 = –2$$, $$y0 = 3$$, оси коллинеарны декартовым осям.

Теоретическая площадь $$S$$ эллипса вычисляется по формуле

$$S=\pi ab,$$

где:

$$a$$ — большая ось эллипса;

$$b$$ — малая ось эллипса.

Каноническое уравнение эллипса (относительно центра координат):

$$\frac{x^2}{a^2}+\frac{y^2}{b^2}=1.$$

Программный код решения примера:

clear all, clc, close all
%%%% Расчет площади эллипса методом Монте-Карло
%%% Параметры эллипса
a = 6;
b = 4;
x0 = -2;
y0 = 3;

%%% Построение эллипса
t = 0 : min([a, b])/500 : 2*pi;
x = a*cos(t) + x0;
y = b*sin(t) + y0;
line(x,y, 'linew', 2, 'color', [1,0,0])
%%% Центр пересечения осей эллипса
line(x0,y0,'marker','o','markerfacecolor','r','color','r' )
 
%%% Построение прямоугольника, в который вписан эллипс
line([x0-a, x0+a],[y0+b,y0+b], 'lines','-.')
line([x0-a, x0+a],[y0-b,y0-b], 'lines','-.')
line([x0-a, x0-a],[y0-b,y0+b], 'lines','-.')
line([x0+a, x0+a],[y0-b,y0+b], 'lines','-.')
 
%%% Гененрирование случайных чисел для области D* - области прямоугольника
N = 500; %%% число испытаний
%%% Генерация чисел по горизонтальной стороне прямоугольника
rx = (x0-a) + (x0+a - (x0-a))*rand(N,1);
%%% Генерация чисел по вертикальной стороне прямоугольника
ry = (y0-b) + (y0+b - (y0-b))*rand(N,1);

%%%% Заполнение случайными числами прямоугольника с эллипсом
for J = 1 : N
line(rx(J),ry(J),'marker','o','markersize',2,'color', 'k',...
 'markerfacecolor', 'k' )
end
 
%%% Подсчет количества случайных чисел, попавших в эллипс
m = 0;
for J = 1 : N
    if ((rx(J) - x0)/a)^2 + ((ry(J) - y0)/b)^2 <= 1
        m = m + 1;
    end
end
 
%%%%% Расчет площади прямоугольника
S = ((x0+a) - (x0-a))*((y0+b)-(y0-b));
 
%%% Расчет площади эллипса методом Монте-Карло
S0 = m/N*S
 
%%% Проверка расчета 
Scontrol = pi*a*b
 
%%% Заголовок для диаграммы
str = sprintf('%s Оси эллипса a = %g, b = %g, координаты центра x_0 = %g,  y_0 = %g.%s%g%s%g. %s%g', '\bf',a, b,x0,y0, ... 
'\newline Теоретическая площадь эллипса S = ', Scontrol, ...
'\newline Площадь эллипса по методу Монте-Карло S0 = ', S0, ...
 'Число испытаний N = ', N);
title(str)
grid off %% выключение сетки на диаграмме
xlabel('\bf - - - - - - - X  - - - - - - - ')
ylabel('\bf - - - - - - - Y  - - - - - - - ')
%%% Установление размеров и свойств графического окна
gfig = get(0, 'screensize');
set(gcf, 'color', 'w', 'position', [gfig(1) + 100, gfig(2) + 100, gfig(3)*0.8, gfig(4)*0.7])
%%% Ограничения по осям
xlim([x0-a-a/10, x0+a+a/10])
ylim([y0-b-b/10, y0+b+b/10])
axis equal

Задание 2

Для построения эллипса используйте его каноническое уравнение.

Предусмотрите, чтобы точки, попавшие в эллипс, были красного цвета, а не попавшие в него — синего цвета.

Координаты центра пересечения осей эллипса примите случайными при использовании функции $$randn$$.

Оси эллипса примите случайными, равномерно распределенными из интервала $$[Х; 5Х]$$, где $$Х$$ — номер компьютера, за которым выполняется лабораторная работа ( $$1, 2, 3, ..$$.).

В программу включите расчет общего числа испытаний, при котором достигается заданная относительная погрешность в 2%.

В программу включите построение графика относительной погрешности от числа испытаний начиная с $$N = 100$$ и более.

Пример 3. Найдите оценку значения тройного интеграла методом Монте-Карло

$$I=\mathop{{\int\int ... \int}}\limits_{V}zdxdydz,$$

где $$V$$ — область, ограниченная круговым конусом $$z^2=\frac{h^2}{R^2}(x^2+y^2)$$ и плоскостью $$z = h$$, $$h = 5$$, $$R = 2$$.

Программный код решения примера:

clear all, clc, close all
%%% Вычисление тройного интеграла методом Монте-Карло, 
R = 2;
h = 5;

%%% Для построения конуса
[x,y] = meshgrid(-1.2*R:0.02:1.2*R, -1.2*R:0.02:1.2*R);
z = h/R*sqrt(x.^2 + y.^2);
xmax = 1.2*R;
%% Радиус окружности в сечении конуса на высоте max(z)
rmax = sqrt(xmax^2 + xmax^2);
%%% Максимальное значение аппликаты
zmax = h/R*sqrt(xmax^2 + xmax^2);
%%% Тангенс угла фронтального сечения конуса
koef = (rmax/zmax);
%%% Радиус окружности в сечении конуса на высоте h
rh = h*(koef); 
 
%%% Круговой конус
mesh(x,y,z) %% функция построения пространственных фигур
hold on
%%% Плоскость на высоте h
line(x, y, h*ones(size(x)))
 
%%% Параллелепипед 
line([-rh, rh],[-rh,-rh],[0,0])
line([-rh, -rh],[-rh,-rh],[0,h])
line([rh, rh],[-rh,-rh],[0,h])
 
line([-rh, rh],[rh,rh],[0,0])
line([-rh, -rh],[rh,rh],[0,h])
line([rh, rh],[rh,rh],[0,h])
 
line([-rh, -rh],[-rh,rh],[0,0])
line([rh, rh],[-rh,rh],[0,0])
 
%%% Вычисление тройного интеграла по встроенным функциям int
syms x y z
In0 = int(int(int(z, z, h/R*sqrt(x^2+y^2), h), ...
    y, -sqrt(rh^2 - x^2), sqrt(rh^2 - x^2)), x, -rh, rh);
In = double(In0)
%%% Расчет тройного интеграла методом Монте-Карло
N = 1000; %%% Число испытаний
S = 0;
for J = 1 : N
    xr = -rh + (rh - (-rh))*rand;
    yr = -rh + (rh - (-rh))*rand ;
    zr = h*rand;
if  zr^2<=h^2/R^2*(xr^2 + yr^2)  h/R*sqrt(xr^2 + yr^2)>=0  ... 
         h/R*sqrt(xr^2 + yr^2) <= h 
 
        S = S + zr; 
      end
end
 
mD = [2*rh]*[2*rh]*[h]; %%% Объем параллелепипеда, мера
In_Car = (mD/N)*S
 
grid on
str = sprintf('%s%g;%s Точное значение тройного интеграла: %g',...
    '\bf Значение интеграла по методу Монте-Карло: ', In_Car, '\newline', In);
title(str)
 
view(-30,20)
axis equal

Задание 3

Формирование случайных переменных $$xr$$, $$yr$$, $$zr$$ вынесите за пределы цикла.

Объясните расстановку пределов интегрирования.

Вычислите по методу Монте-Карло объем тела для полученных пределов интегрирования. Сравните со значением, вычисленным аналитически.

Задание 4

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

Рассчитайте с помощью метода Монте-Карло площадь фигуры, ограниченной гипоциклоидой с параметрами $$R = 3$$ ; $$r = 2.7$$ ; $$m = r/R$$. Параметрическое уравнение гипоциклоиды:

$$x =(R-m*R)*cos(m*t)+m*R*cos(t-m*t),\\ y = (R-m*R)*sin(m*t)-m*R*sin(t-m*t).$$

Подберите массив значений параметра $$t$$.