интеграция Монте-Карло на гиперкубе R ^ 5 в MATLAB - PullRequest
Новая
       44

интеграция Монте-Карло на гиперкубе R ^ 5 в MATLAB

0 голосов
/ 25 сентября 2011

Мне нужно написать код MATLAB, который будет интегрироваться через гиперкуб 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)

Ответы [ 3 ]

2 голосов
/ 26 сентября 2011

В любом случае, я думаю, я понимаю, что здесь задумано, хотя это кажется чем-то вроде надуманного упражнения.Рассмотрим проблему попытки найти площадь круга с помощью MC, как обсуждено здесь .Здесь образцы выбираются из единичного квадрата, и функция принимает значение 1 внутри круга и 0 снаружи.Чтобы найти объем куба в R ^ 5, мы могли бы взять образец чего-то еще, что содержит куб, и использовать аналогичную процедуру для вычисления желаемого объема.Надеюсь, этого достаточно, чтобы сделать остальную часть реализации простой.

1 голос
/ 28 сентября 2011

Я думаю, что здесь немного, поскольку числа, которые вы даете в качестве «правильного» ответа, не соответствуют тому, как вы формулируете упражнение (объем единичного гиперкуба равен 1).

Учитывая, что результат должен быть1/120 - может ли быть так, что вы должны интегрировать стандартный симплекс в гиперкубе?

Ваша функция будет понятна.f (x) = 1, если сумма (x) <1;0 в противном случае </p>

0 голосов
/ 29 сентября 2011 
%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);
...