-2

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

3 回答 3

2

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

于 2011-09-26T00:08:32.020 回答
1

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

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

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

于 2011-09-28T12:13:13.273 回答
0 
%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);
于 2011-09-29T19:28:36.227 回答