我需要编写 MATLAB 代码,使用 Monte Carlo 在 R^5 超立方体上进行积分。当我有通用函数时,我有一个基本算法。但我需要集成的功能是:
∫dA
A 是 R^5 的元素。
如果我有 ∫f(x)dA 那么我认为我的算法会起作用。
这是算法:
% Writen by Jerome W Lindsey III
clear;
n = 10000;
% Make a matrix of the same dimension
% as the problem. Each row is a dimension
A = rand(5,n);
% Vector to contain the solution
B = zeros(1,n);
for k = 1:n
% insert the integrand here
% I don't know how to enter a function {f(1,n), f(2,n), … f(5n)} that
% will give me the proper solution
% I threw in a function that will spit out 5!
% because that is the correct solution.
B(k) = 1 / (2 * 3 * 4 * 5);
end
mean(B)
无论如何,我想我明白这里的意图是什么,尽管它看起来确实有点做作。考虑一下尝试通过 MC 求圆面积的问题,如here所述。这里的样本是从一个单位正方形中抽取的,函数在圆内取值 1,在圆外取值 0。为了找到 R^5 中立方体的体积,我们可以从包含该立方体的其他物体中采样,并使用类似的过程来计算所需的体积。希望这足以让其余的实现变得简单。
我在这里猜测了一点,因为你给出的“正确”答案的数字与你陈述练习的方式不匹配(单位超立方体的体积为1)。
鉴于结果应为 1/120 - 您是否应该将 标准单纯形 集成到超立方体中?
你的职能会很清楚。 f(x) = 1 如果 sum(x) < 1; 0 otherwise
%Question 2, problem set 1
% Writen by Jerome W Lindsey III
clear;
n = 10000;
% Make a matrix of the same dimension
% as the problem. Each row is a dimension
A = rand(5,n);
% Vector to contain the solution
B = zeros(1,n);
for k = 1:n
% insert the integrand here
% this bit of code works as the integrand
if sum(A(:,k)) < 1
B(k) = 1;
end
end
clear k;
clear A;
% Begin error estimation calculations
std_mc = std(B);
clear n;
clear B;
% using the error I calculate a new random
% vector of corect length
N_new = round(std_mc ^ 2 * 3.291 ^ 2 * 1000000);
A_new = rand(5, N_new);
B_new = zeros(1,N_new);
clear std_mc;
for k = 1:N_new
if sum(A_new(:,k)) < 1
B_new(k) = 1;
end
end
clear k;
clear A_new;
% collect descriptive statisitics
M_new = mean(B_new);
std_new = std(B_new);
MC_new_error_999 = std_new * 3.921 / sqrt(N_new);
clear N_new;
clear B_new;
clear std_new;
% Display Results
disp('Integral in question #2 is');
disp(M_new);
disp(' ');
disp('Monte Carlo Error');
disp(MC_new_error_999);