Материал: Matlab. Практический подход. Самоучитель

Самоучитель Matlab

ствительное число в диапазоне от 0 до 1. Если в качестве аргументов функции rand() указать целые числа, то в качестве результата возвращается массив соответствующих размеров, заполненный случайными действительными числами в диапазоне значений от 0 до 1. Например, в результате вызова команды rand(5,10) будет создана матрица размерами 5 на 10 из случайных чисел. Если указать один аргумент, например rand(5), то будет создана квадратная матрица соответствующего размера (в данном случае 5 на 5) из случайных чисел.

На заметку

Функция rand() далеко не единственная функция Matlab для генерирования случайных чисел. Например, функция randn() позволяет генерировать случайные числа с функцией стандартного нормального распределения, а функция randi() позволяет генерировать равномерно распределенные целые числа.

Простой пример использования функции генерирования случайных чисел приведен в документе на рис. 8.21.

Рис. 8.21. Вычисление площади области методом Монте-Карло

В качестве иллюстрации мы попытаемся методами Монте-Карло вычис-

336

грала дают одно и то же значение 1 3 .

На заметку

Если быть более точным, то двойной интеграл сводится к повторному, а повторный – к однократному. Тем не менее, с точки зрения идеологии метода МонтеКарло, правильнее было бы интерпретировать наши вычисления как расчет двойного интеграла (который равен площади соответствующе области). Что касается непосредственно метода Монте-Карло, то в данном случае он реализован по следующей схеме. В пределах единичного квадрата, внутри которого находится область, площадь которой мы вычисляем, случайным образом выбираются точки. Их много. Мы подсчитываем отношение количества точек, которые попали внутрь области, к общему количеству точек. В граничном пределе (когда количество точек неограниченно возрастает) это отношение стремится

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

кобщему количеству точек принимаем как оценку для площади области.

Для выполнения вычислений создаем следующий программный код:

z=0:0.01:1;

y1=sqrt(z);

y2=z.^2; plot(z,y1,'r-','LineWidth',2); hold on; plot(z,y2,'r-','LineWidth',2); grid on;

title('Вычисление интеграла методом Монте-Карло'); n=1000;

s=0;

for i=1:n x=rand(); y=rand();

if and(y<sqrt(x),y>x^2) s=s+1/n;

end plot(x,y,'bo','LineWidth',1,'MarkerSize',2,'MarkerFaceColor','b'); end

hold off;

disp(['Количество точек: ',num2str(n)]); disp(['Результат вычислений: ',num2str(s)]);

337

Самоучитель Matlab

Вместе с вычислением площади области мы также выполняем некоторые графические построения. Вначале создается график с кривыми для зависимостей, которые ограничивают область. После этого командой n=1000 задаем количество точек, которые будут сгенерированы. В переменную s, которая инициализируется с начальным нулевым значением, будем записывать результаты вычислений площади. Вычисления выполняются в операторе цикла. Индексная переменная пробегает значения от 1 до n. Командами x=rand() и y=rand() генерируются случайные числа, которые служат координатами случайной точки. Затем в условном операторе проверяется условие and(y<sqrt(x),y>x^2), которое означает, что точка попадает внутрь области. Если это так, командой s=s+1/n увеличиваем значение переменной s. Но это еще не все. Сгенерированную точку мы отображаем на графике с помощью команды plot(x,y,'bo','LineWidth',1, 'MarkerSize',2,'MarkerFaceColor','b'). Точки отображаются синими кружками с синей заливкой маркера (опция MarkerFaceColor). В завершение кода выводится информация о количестве сгенерированных точек и оценке для площади области. В этом случае мы использовали команду disp(). Аргументом передается текстовая строка. Строка получается созданием массива с текстовыми элементами. Для перевода числа в текстовый формат использовалась функция num2str().

На заметку

Для объединения строк можно было бы использовать и функцию strcat(). Однако эта функция автоматически удаляет конечные пробелы в объединяемых строках, что в данном случае не очень приемлемо.

Приведенный выше код сохраняем в файле MCInt.m. После этого для выполнения комплекса вычислений в командном окне выполняем команду MCInt (рис. 8.22).

Рис. 8.22. Вычисление интеграла в командном окне

Здесь мы реализуем вычисления на основе 1000 точек. Это не очень много, поэтому результат от запуска к запуску может существенно меняться (на уровне сотых). График с точками, который создается в результате выполнения кода, представлен на рис. 8.23.

338

Глава 8. Обработка данных

Рис. 8.23. Графическое представление результата

На заметку

Даже из этого простого примера видно, что точность метода Монте-Карло оставляет желать лучшего. Чтобы добиться приемлемого результата, необходимо существенно увеличить количество тестов. А это требует значительных затрат времени (и ресурсов). Поэтому к методам вроде Монте-Карло обычно прибегают, когда другие методы неприменимы.

Выше мы рассматривали задачу, в которой нужно было генерировать случайные числа, равномерно распределенные в диапазоне от 0 до 1. На практике приходится иметь дело со случайными величинами, которые имеют самые разные законы распределения. В таких случаях можно либо воспользоваться встроенной функцией (если такая есть) для генерирования чисел с нужными характеристиками, либо создать собственную. Нас интересует последний случай. Пикантность ситуации состоит в том, что создать функцию, генерирующую случайные числа с практически любой приличной (не очень экзотической) функцией распределения, можно на основе функции, генерирующей равномерно распределенные числа. Другими словами, на основе функции rand() можно достаточно просто создать функцию, которая генерирует числа с заданным законом распределения.

339

Самоучитель Matlab

Полезным окажется то обстоятельство, что если случайная величина ξ распределена на интервале [0,1], то случайная величина η = F−1(ξ) имеет функцию распределения F(x). Здесь функция F(x) непрерывная и строго возрастающая, а через F−1(x) обозначена обратная функция к функции F(x). Таким образом, если нам нужно создать функцию для генерирования случайных чисел с функцией распределения F(x), то мы должны найти обратную функцию F−1(x), а затем вычислять случайные числа η в соответствии с соотношением η = F−1(ξ), где случайное число ξ генерируется, например, с помощью функции rand().

На заметку

Напомним, что по определению функцией распределения Fξ(x) случайной величины ξ называется вероятность P(ξ < x) того, что случайная величина принимает значение меньшее x (параметр −∞ < x < +∞). Плотностью

распределения случайной величины называется функция pξ(x) = dFξ(x) .

dx

Для равномерно распределенной на интервале от 0 до 1 случайной величины функция распределения имеет вид Fξ(x) = 0 при x < 0, Fξ(x) = x при

0 ≤ x ≤ 1 и Fξ(x) = 1 при x > 1. Плотность распределения для этого распределения pξ(x) = 1 при x [0,1] и pξ(x) = 0 при x [0,1].

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

На заметку

Случайная величина ξ распределена с показательным законом с параметром

α, если функция ее распределения равна Fξ(x) = 1 − exp(−αx) при x ≥ 0

иравна нулю в противном случае (то есть при x < 0).

Вданном случае функция распределения F(x) = 1 − exp(−αx). Обрат-

ная к ней функция может быть найдена в результате решения уравнения x = F(y) относительно y , то есть y = F−1(x) . Это нам позволит найти обратную функцию. Имеем y = 1 − exp(−αx) . После несложных преобразо-

этому если мы будем генерировать случайные числа ξ функцией rand(),

ηбудут иметь показательное распределение с параметром α . Обратимся

кдокументу на рис. 8.24.

340