метод Монте-Карло для вычисления многомерного интеграла

[shwedka]

Здесь еще одно обстоятельство исходной задачи не учтено: Nobody хочет интегруровать не по квадрату какому-то, а по бесконечной области. Конечным числом точек бесконечную область не аппроксимировать. Потому нужна предварительная аналитическая работа. В зависимости от требуемой точности, следует обрезать область интегрирования так, чтобы интеграл по бесконечному обрезаемому куску был бы меньше, скажем, половины допускаемой погрешности.
А по ограниченной области считеть уже по какой квадратурной ('кубатурной') формуле.
Оценку интеграла по отбрасываемому куску проводить можно в зависимости от вида интегранда. Или оценивать интегранд чем-то простым, либо, если там осциляции, оценивать интеграл по стацфазе....

[Николай]

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

[Yuri Gendelman]

К сожалению физики и математики, применяющие Монте-Карлу, редко заглядывают во 2-й том "Искусства программирования" Д.Кнута. А в нем очень подробно рассматриваются такие темы как:
- методы оценки качества генераторов псевдо-случайных чисел;
- низкое качество некоторых популярных ГПСЧ, в т.ч. применяемых в физике и включенных в стандартные библиотеки;
- генерация случайных векторов (:!: некоторые ГПСЧ дают хорошие 1-мерные последовательности, но плохие, например, тройки);

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

Цитата:

...Равномерность распределения для вычисления интеграла важнее случайности...

Не приводит ли следование этому принципу, в пределе, к сеточным методам?

[maxal]

Руст писал(а):

Любая случайная величина, используемая в компьютере, не является случайной.

При отсутствии источников энтропии - да, но сейчас есть генераторы вовсю использующие внешнюю энропию. Вот неплохой обзор ГПСЧ.

[Руст]

Yuri Gendelman писал(а):

Не приводит ли следование этому принципу, в пределе, к сеточным методам?

Сеточные методы не удобны в многомерном случае из за плохого приближения в смысле равномерности на краях.
Поэтому, оптимальными для вычисления многомерных интегралов (по моему) являются методы основанные на теории чисел, изложенные в книге: К.М. Коробов "Теоретико числовые методы в приближённом анализе".

[Николай]

У меня есть книжка С.В.Поршнев "Компьютерное моделирование физических процессов в Matlab". Он заглядывал во 2-й том "Искусства программирования" Д.Кнута, о чем свидетельствует список литературы . Вот что Поршнев пишет про ГПСЧ:

Алгоритм генерации случайных чисел с равномерным законом распределения
Aлгоритм генерации случайных чисел с заданным законом распределения основан на использовании случайных чисел с равномерным законом распределения. Следовательно, точность вычисления интеграла зависит не только от количества точек, используемых для вычисления значения интеграла, но и от качества работы генератора случайных чисел.
Для применения в вычислительных экспериментах генератор должен обладать следующими характеристиками:
— хорошими вычислительными свойствами,
— эффективностью,
— большим периодом,
— воспроизводимостью.
Первыми и наиболее важными из рассматриваемых характеристик являются, конечно, статистические свойства. Известен пример, когда использование плохого генератора случайных чисел в методе молекулярной динамики привело к обнаружению особого поведения системы определенного размера. Результаты независимого исследования, проведенного другими авторами, не подтвердили этой особенности. Это позволило сделать вывод о том, что причиной обнаруженной особенности явилось плохое качество генератора случайных чисел.
Также очень важна вычислительная эффективность генератора. Программы для компьютерных экспериментов требуют огромного числа случайных чисел. Для получения порядка 10^10 чисел за ограниченное время одно случайное число должно вычисляться очень быстро. Кроме того, важно требование на размер необходимой для работы генератора памяти компьютера. Генератор должен быть достаточно быстрым и экономичным по памяти.
В действительности, генерируются не случайные числа, а последовательности псевдослучайных чисел, поэтому последовательность генерируемых чисел после определенного количества шагов повторяет саму себя. Ее период должен быть заведомо больше длины, необходимой для последовательности случайных чисел. Опасность представляет также исчерпание последовательности случайных чисел или даже приближение конца ее цикла. Это также приводит к неверным результатам.
Существуют различные способы генерации на компьютере случайных чисел. Далее мы рассмотрим, в первую очередь, наиболее широко распространенные в настоящее время методы, объединенные обидим названием «линейные конгруэнтные генераторы». Они производят определенным способом при задании начального целого числа последовательность случайных чисел с равномерным законом распределения. Несомненным достоинством генераторов данного типа является их воспроизводимость.
Линейный конгруэнтный генератор выдает последовательность произвольных целых чисел, которые можно привести к единичному интервалу [0,1]. Для этого достаточно разделить все члены последовательности на максимальное возможное целое число. Ненормированная последовательность целых чисел, однако, предпочтительнее по причине большей вычислительной эффективности алгоритма.

Алгоритм линейного конгруэнтного генератора (Лемер, 1948)
Пусть m, а, с, х0 - причем m > х0, m > а, m > с и х0 > 0, а > 0, с > 0.
Тогда псевдослучайное число последовательности xi, получается из предшествующего ему числа x(i-1) как
xi=(a-xi+c)modm
(Здесь mod - функция деления по модулю, т. е. функция получения остатка от деления одного числа на другое. Например 19 mod 4 = 3.)
При соответствующем выборе чисел m, а, с, х0 алгоритм выдает последовательность случайных чисел. При этом генерируемая последовательность имеет повторяющийся цикл или период, не превышающий m чисел. Для получения достоверных результатов при моделировании физических систем необходимо использовать генераторы с максимально большим периодом. Максимальный период достигается при с<>0. Такой генератор называется смешанным. Однако, как показала практика, смешанные генераторы дают плохие результаты по сравнению с генераторами, у которых с=0 (мультипликативные генераторы). Отметим, что для обеспечения максимального периода мультипликативного конгруэнтного генератора модуль m и множитель а необходимо подбирать специально.
Алгоритм линейного конгруэнтного генератора в пакете MATLAB реализуется выполнением следующей последовательности команд:
» а=7;
» с=7;
» m=55;
» х(1)=20;
» N=100;
» for i=1:N
x(i+1)=mod(a*x(i)+c,m);
end;
» i=1:N+1; plot(i, x, ‘-kS’);
» axis([0 100 0 60]);

[Руст]

Как я уже писал, при вычислении интеграла случайность чисел не так важен, главное, чтобы искомоые точки (узлы) равномерно покрывали область интегрирования. В ногомерном случае количество точек на одну размерность $N^{1/n}$ оказывается небольшим числом даже когда N больше или порядка миллиона. Поэтому, сеточные методы работают не удовлетворительно. Методы теории чисел приводят выбору точек:
$a^{k_ji}(mod \ p),j=1,2,...,n, i=1,2,...,N.$
Здесь a натуральное число, имеющее больший период, чем N по модулю большого простого числа p. Соответственно, точки получаются разбросанными равномерно всюду (даже по краям) и результаты интегрирования лучше. Есть и другие методы, основанные на теории чисел. Можно проверить на тестах, когда результат известен (например вычисляются аналитический), и сравнить разные методы. На сколько мне известно, лучшего, чем методы теории чисел, для многомерного интегрирования пока не придумано.

[Николай]

Попробую встать на защиту Монте-Карло .
Да, сеточный метод действительно плохо подходит для случая функции $ n переменных. С методом Монте-Карло несколько другая ситуация и, как я помню, его преимущества проявляются именно в случае $ n-мерных интегралов. Более того, мне теперь очевидно, что програмку на предыдущей странице обобщить на случай 8-ми переменных совсем не сложно. Итак, если задана функция $\ f(x_1,...,x_n), то нужно сгенерировать $\ n+1 случайных массивов $ X_1, X_2,...,X_n,Y некоторой длины (ну, для определённости, 1000). Таким образом в пространстве $\mathbb{R}^{n+1}$ мы зададим $ 1000^{n+1} точек. И потом программа на каждом $\iota-ом "такте" своей работы будет вычислять значение функции $ f(x_1(\iota),..,x_n(\iota)) и сравнивать его с $ y(\iota) (Теперь понятно, что $ y(\iota) обязано быть случайной, равномерно-распределённой величиной). Таким образом, всё произойдёт в одном цикле. Объяснение получилось каким-то мудрёным , но суть преимущества, надеюсь, понятна.

P.S. Руст, я к сожалению и представления не имею о методе из теории чисел, который Вы хвалите. Поэтому скажите, позволяет ли он посчитать, например, 100-мерный интеграл (а такие, как известно, возникают)?

[photon]

Николай писал(а):

... Скажу, что этот интеграл я посчитал на Matlab и на Fortran. Так вот, когда я задал в Matlab 100000 "испытаний", он ушел в себя и н вернулся, пока я пил чай . А Fortran справился и с 10000000 испытаниями всего за 26 секунд (засекал секундомером ).

А Вы не могли бы привести код MatLAB - меня одолевают сомнения, что он построен "хорошо", если Вы умудрились отправить его "в себя".

[Николай]

Пожалуйста:
"Смерть" в виде пожизненной надписи "Busy" наступила уже на первом цикле при изменении N c 10000 на 100000.

Код писал(а):

>> %нижний "угол" прямоугольника
>> Xmin=0;
>> Ymin=0;
>> %верхний "угол" прямоугольника
>> Xmax=pi/2;
>> Ymax=1.5;
>> N=10000;
>> for i=1:N
x(i)=Xmin+(Xmax-Xmin)*rand(1);
end;
>> for i=1:N
y(i)=Ymin+(Ymax-Ymin)*rand(1);
end;
>> f=inline('sin(x)','x');
>> for i=1:N
if f(x(i))>=y(i)
s(i)=1;
else
s(i)=0;
end;
end;
>> S=sum(s);
>> A=(Xmax-Xmin)*(Ymax-Ymin);
>> A*S/N

ans =

1.0047

[Someone]

В связи с обсуждаемым вопросом: существует датчик псевдослучайных чисел (Mersenne Twister) с периодом $2^{19937}-1$, который обеспечивает равномерное распределение с точностью 32 бита в 623-мерном кубе. У него, однако, есть свои недостатки.

[Руст]

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

[photon]

Николай писал(а):

Пожалуйста:
"Смерть" в виде пожизненной надписи "Busy" наступила уже на первом цикле при изменении N c 10000 на 100000.

Ха . Я чуть-чуть модифицировал Ваш код:

clear all
tic
%нижний "угол" прямоугольника 
Xmin=0; 
Ymin=0; 
 %верхний "угол" прямоугольника 
Xmax=pi/2; 
Ymax=1.5; 

N=1e7; 

x=Xmin+(Xmax-Xmin)*rand(1,N); 
y=Ymin+(Ymax-Ymin)*rand(1,N); 

f=sin(x);

s=((f-y)./abs(f-y)+1)/2;

S=sum(s); 
A=(Xmax-Xmin)*(Ymax-Ymin); 
A*S/N
toc

В результате на машине Athlon 1700+ для $N=10^7$ получил ответ

ans =
    0.9996
elapsed_time =
    4.3760

хе-хе - меньше 4.5 секунд

Дописал: Существенные проблемы могут возникнуть, если не хватит памяти - подкачка осуществляется очень медленно. Так для $N=10^8$ я не смог дождаться - памяти не хватало на 4 (x,y,f,s) таких больших массива. Я не смог даже для $N=5\cdot10^7$. В два раза по памяти можно выгадать, если написать так, избавляясь от массивов s и f:

...
N=1e8; 

x=sin(Xmin+(Xmax-Xmin)*rand(1,N)); 
y=Ymin+(Ymax-Ymin)*rand(1,N); 

y=((x-y)./abs(x-y)+1)/2;

S=sum(y); 
...

Но для $N=5\cdot10^7$ моих 512 Мб все равно не хватило - требовалось почти Гб, а для $N=10^7$ в этом варианте по времени проиграл пару сотен миллисекунд. А для $N=1.5\cdot10^7$ - получилось (6.8с)

Но у Вас проблемы были не с памятью (кстати, сколько памяти в Фортране у Вас выделялось под один элемент массива? -я не знаю Фортран, в MatLAB - 8 байт).
Основные проблемы Вашего MatLAB-кода:
1) использование циклов - MatLAB очень удобен для работы с массивами - старайтесь избегать циклов - выигрыш во времени получите в 10 и более раз. Как видите, я не оставил ни одного цикла (Кстати, как Вам понравилось, как я обошел цикл с if-ом? Я прям собой горжусь . Хотя можно и по другому. Например: взять корень от разности, разделить на модуль корня разности и выделить реальную часть - только операция извлечения корня наверное будет работать медленнее);
2) переход к символьным вычислениям - это о-очень медленно работает.

[photon]

Подумал-подумал и решил, что вместо

y=((x-y)./abs(x-y)+1)/2;
S=sum(x);

будет лучше

y=x-y;
y=(y./abs(y)+1);
S=sum(x)/2;

проверил, но, правда, ощутимого выигрыша не обнаружил.

[Николай]

photon, Вы мне здорово помогли. Такой уровень мастерства програмирования мне не доступен (тем более, на Matlab, в котором я использую практически только pdeToolBox). Воспроизвёл у себя Ваш пример - действительно всё посчиталось, хоть и не так быстро (оперативы мало ). Кстати, чтобы сэкономить память, я вообще исключил массив S:

      DO I = 1,NMAX
        X(I)=XMIN+(XMAX-XMIN)*RANX(I) 
        Y(I)=YMIN+(YMAX-YMIN)*RANY(I) 
        IF (SIN(X(I)).GE.Y(I)) THEN
            SUM=SUM+1
        ENDIF
      END DO

Разумно?

На счет Фортрана - там тоже достаточно удобно работать с массивами. Сколько памяти выделяется под элемент массива? Не знаю - наверное, это зависит от ранга массива (В фортране максимум могут быть 7-мерные массивы), но если я правильно понял английскую справку - на элемент одномерного массива выделяется 4 байта. Удобно и то, что динамическую память можно освобождать, "отключая" массивы, если они уже не нужны (DEALLOCATE(...)).
Два слова про ГПСЧ, которым снабжён Фортран - для генерации чисел он использует (каким образом - тайна) текущие дату и время, что внушает некоторое доверие к результату.