Монте-Карло интегрирования ехр (-х^2/2) от х = -бесконечности х = + бесконечность

Question

Я хочу, чтобы интегрировать

п (х) = ехр (-х^2/2)

от x = -инфекция до x = + бесконечность

с использованием метода Монте-Карло. Я использую функцию randn() для генерации всех x_i для функции f (x_i) = exp (-x_i^2/2). Я хочу интегрировать, чтобы затем вычислить среднее значение f ([x_1, .. x_n]). Моя проблема заключается в том, что результат зависит от того, какие значения я выбираю для своих границ x1 и x2 (см. Ниже). Мой результат удаляется от реального значения, увеличивая значение x1 и x2. На самом деле результат должен быть лучше и лучше, увеличивая x1 и x2.

Кто-нибудь видит мою ошибку?

Вот мой Matlab код

clear all; 
b=10;     % border 
x1 = -b;    % left border 
x2 = b;    % right border 
n = 10^6;    % number of random numbers 
x = randn(n,1); 
f = ones(n,1); 
g = exp(-(x.^2)/2); 
F = ((x2-x1)/n)*f'*g; 

Право значение должно быть ~ 2,5066.

Благодаря

Answer 1

Попробуйте это:

clear all; 
b=10;     % border 
x1 = -b;    % left border 
x2 = b;    % right border 
n = 10^6;    % number of random numbers 

x = sort(abs(x1 - x2) * rand(n,1) + x1); 
f = exp(-x.^2/2); 
F = trapz(x,f) 

F = 

    2.5066 

Answer 2

Ok, давайте начнем с написания общего случая интеграции MC:

I = S f(x) * p(x) dx, x in [a...b] 

S здесь является неотъемлемым признаком.

Обычно p(x) является функция нормализуется плотность вероятности, f(x) вы хотите интегрировать, и алгоритм очень прост:

• комплект аккумулятор s к нулю
• начало петли N событий
• образца x в случайном порядке от p(x)
• x, вычислить f(x) и добавить t О аккумулятор
• назад, чтобы начать цикл, если не сделано
• если сделано, разделите аккумулятор на N и вернуть его

В простейшем учебнике случае вы имеете

I = S f(x) dx, x in [a...b] 

, где это означает PDF равно к равномерно распределенному

p(x) = 1/(b-a) 

но что вы должны сумма фактически (b-a)*f(x), потому что ваш интеграл теперь выглядит

I = S (b-a)*f(x) 1/(b-a) dx, x in [a...b] 

В общем, если оба f(x) и p(x) могли бы служить в качестве PDF, то это вопрос выбора интеграции ли вы f(x) над p(x) или p(x) над f(x) ,Нет разницы! (Ну, за исключением, может быть, время вычисления)

Итак, вернемся к частным интегралом (который равен \sqrt{2\pi}, я верю)

I = S exp(-x^2/2) dx, x in [-infinity...infinity] 

Вы можете использовать более традиционный подход, как @Agriculturist и писать

I = S exp(-x^2/2)*(2a) 1/(2a) dx, x in [-a...a] 

и образец x из U (0,1) в [-a ... а] интервале, и для каждого x вычислительного ехра() и усреднить и получить результат

Из того, что я понимаю, вы хотите использовать exp(), как PDF, так что ваши неотъемлемые выглядит как

I = S D * exp(-x^2/2)/D dx, x in [-infinity...infinity] 

PDF быть нормированы так оно должно включать в себя коэффициент нормализации D, которая в точности равна \sqrt{2 \pi} от гауссова интеграла.

Сейчас f(x) всего лишь постоянная равная D. Это не зависит от x. Это означает, что вы за каждый выборку x должны добавить к аккумулятору значение ПОСТОЯННОЕ значение D. После запуска N образцов, в аккумуляторе у вас будет ровно N*D. Чтобы найти среднее значение, вы разделите на N, и в результате вы получите совершенный D, что составляет \sqrt{2 \pi}, что, в свою очередь, составляет 2.5066.

Слишком ржавый, чтобы писать MATLAB и С Новым годом в любом случае