Вычисление интегралов методом Монте-Карло (стр. 2 из 3)

Таким образом, значение искомого интеграла можно выразить как произведение математического ожидания функции и объема n- мерного параллелепипеда

:

(8)

Следовательно, необходимо найти значение математического ожидания

. Его приближенное значение можно найти произведя n испытаний, получив, таким образом, выборку случайных векторов, имеющих равномерное распределение на . Обозначим и . Для оценки математического ожидания воспользуемся результатом

, (9)

где

,

,

- квантиль нормального распределения, соответствующей доверительной вероятности .

Умножив двойное неравенство из (9) на

получим интервал для I:

. (10)

Обозначим

точечную оценку . Получаем оценку (с надежностью ):

. (11)

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

:

. (12)

Если задана целевая абсолютная погрешность

, из (11) можно определить объем выборки, обеспечивающий заданную точность и надежность:

. (13)

Если задана целевая относительная погрешность, из (12) получаем аналогичное выражение для объема выборки:

. (14)

1.3 Сплайн – интерполяция.

В данном программном продукте реализована возможность задавать дополнительные ограничения области интегрирования двумя двумерными сплайн – поверхностями (для подынтегральной функции размерности 3). Для задания этих поверхностей используются двумерные сплайны типа гибкой пластинки \4\.

Под сплайном (от англ. spline - планка, рейка) обычно понимают агрегатную функцию, совпадающую с функциями более простой природы на каждом элементе разбиения своей области определения. Сплайн – функция имеет следующий вид:

. (15)

Исходные данные представляют собой

троек точек .

Коэффициенты

и определяются из системы:

, (16)

где

,

.

1.4 Алгоритм расчета интеграла

Реализованный алгоритм включает следующие шаги:

1) выбирается начальное значение

, разыгрываются случайные векторы из и определяются и ;

2) в зависимости от вида погрешности (абсолютная, относительная) определяется достигнутая погрешность; если она меньше целевой, вычисление прерывается;

3) по формулам (13) или (14) вычисляется новый объем выборки;

4) объем выборки увеличивается на 20%

5) переход к шагу 1;

6) конец.

2. ГЕНЕРАТОР ПСЕВДОСЛУЧАЙНЫХ ЧИСЕЛ

2.1 Генератор псевдослучайных чисел применительно к методу Монте – Карло.

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

2.2 Алгоритм генератора псевдослучайных чисел

В программе реализован конгруэнтный метод генерации псевдослучайных чисел \3\:

, (17)

где

=8192,

=67101323.

Авторский код, реализующий защиту от переполнения был, реализован на С++. Перед использование первые три числа последовательности удаляются. Для получении чисел из интервала (0,1) все числа делятся на

.

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

Проверка равномерности распределения псевдослучайных чисел проводилась с помощью стандартного критерия χ² \2\.

Были использованы 3 последовательности псевдослучайных чисел, определяемых стартовыми значениями 1, 1001, 1000000 длиной 300000.

Интервал (0,1) подразделялся на 50 равных интервалов и программно подсчитывались абсолютные частоты (рис. 1).

Рис. 1

Результаты проверки приведены в Таблице 1.

Таблица 1

	стартовое значение ГСЧ
	1	1001	1000000
хи-квадрат	44.0533333333333	45.007	48.618
df	50	50	50
p-значение	0.709735881642893	0.673522612551685	0.528941919633451

Следовательно, равномерность распределения не отвергается на уровне 5%.

ЗАКЛЮЧЕНИЕ

В заключение можно сказать, что поставленная задача была полностью выполнена. То есть на языке С++ были разработаны генератор псевдослучайных чисел, функция рассчитывающая интеграл методом Монте – Карло (Приложение 1); был проведен расчет тестовых многомерных интегралов (Приложение 2); в интегрированной среде разработки приложений Borland C++ Builder Enterprises 7.0 был создан программный продукт «CarloS», реализующий описанные выше алгоритмы (Приложение 3).

СПИСОК ИСПОЛЬЗОВАННЫХ ИСТОЧНИКОВ

1. Бережная Е. В., Бережной В. И. Математические методы моделирования экономических систем. – М.: Финансы и статистика, 2001. – 368 с.

2. Мюллер П., Нойман П., Шторм Р. Таблицы по математической статистике. – М.: Финансы и статистика, 1982. – 278 с.

3. Теннант-Смит Дж. Бейсик для статистиков. – М.: Мир, 1988. – 208 с.

4. Baranger J. Analyse numérique. Hermann, 1991.

5. Маделунг Э. Математический аппарат физики. Справочное руководство. М.: Наука, 1968., с.287.

6. В.Е. Гмурман Теория вероятностей и математическая статистика – М.: Высшая школа, 2003

ПРИЛОЖЕНИЕ 1

ЛИСТИНГИ ОСНОВНЫХ ФУНКЦИЙ

Листинг 1 Функция расчета интеграла

void integral ()

{

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

// размерность области интегрирования

unsigned d_int=fun_dim;

//----- 3 d график --------------------------------------------------------

// максимальное число троек

unsigned plot_dim_max=10000;

// матрица троек

pmatd xyz,xyz_tmp;

if (d_int==3) xyz=new matd(plot_dim_max,3);

//-------------------------------------------------------------------------

// индикатор относительной погрешности

mcres.relok=Read1double("error_type.txt");

// целевая погрешность

mcres.dlt_int=Read1double("error_value.txt");

// номер стандартного значения доверительной вероятности (начиная с 0)

int nome_int=Read1double("error_omega.txt");

// ГСЧ

unsigned long b=m_rng*m_rng-d_rng,c,r,i,PSChunk;

// "росток" ГСЧ

mcres.rng_seed=Read1double("rng_seed.txt");

pmatd fun_b, fun_A, con_b, con_A, con_U, con_v, \

a_int, b_int, ba_int, x_int, xyz_top, xyz_bottom;

unsigned j,ii,jj,con_ok;

struct date dat;

struct time tim;

pspl2d sp_top,sp_bottom;

// квантили нормального распределения

double omegas_int[6]={0.9,0.95,0.99,0.999,0.9999,0.99999};

double zs_int[6]={1.64485362695147,1.95996398454005,2.5758293035489, \

3.29052673149191, 3.89059188641317, 4.4171734134667};

mcres.omega_int=omegas_int[nome_int];

mcres.z_int=zs_int[nome_int];