MATLAB 中 R^5 超立方体的蒙特卡罗积分

Question

我需要编写 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) 

Answer 1

无论如何，我想我明白这里的意图是什么，尽管它看起来确实有点做作。考虑一下尝试通过 MC 求圆面积的问题，如此处所述。这里的样本是从一个单位正方形中抽取的，函数在圆内取值 1，在圆外取值 0。为了求出 R^5 中立方体的体积，我们可以从包含该立方体的其他物体中进行采样，并使用类似的过程来计算所需的体积。希望这足以让其余的实现变得简单。

Answer 2

我在这里猜测一下，因为您作为“正确”答案给出的数字与您陈述练习的方式不匹配（单位超立方体的体积为 1）。

鉴于结果应该是 1/120 - 是否您应该集成 超立方体中的标准单纯形？

你的功能会很清楚。如果 sum(x) <，则 f(x) = 1 1;否则为 0

Answer 3

%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);