当前位置：首页 > article >正文

蒙特卡洛方法(Monte Carlo，MC)

article 2024/11/26 13:39:46

1 序言

2 Monte Carlo法计算积分

3 最优化计算Monte Carlo法

1 序言

蒙特卡罗方法(Monte Carlo)是由冯诺依曼和乌拉姆等人发明的，“蒙特卡罗”这个名字是出自摩纳哥的蒙特卡罗赌场，这个方法是一类基于概率的方法的统称。是一种应用随机数来进行计算机模拟的方法，此方法随研究的系统进行随机观察抽样，通过对样本值的观察统计，求得所研究系统的某些参数。

2 Monte Carlo法计算积分

考虑二重积分

$I=\iint_{A}f(x,y)dxdy,\; \; f(x,y)\geq 0,\; \; \forall (x,y)\in A$

根据其几何意义，它是以f(x,y)为曲面顶，A为底的柱体C的体积。用下列简单思路求的近似值:假设C被包在几何体D的内部，D的体积已知，若在D内产生1个均分布的随机数，那么

P(随机数落在C内) $\approx$ C的体积/D的体积

现用Monte Carlo法计算： $I=\iint_{x^{2}+y^{2}\leq 1}\sqrt{1-x^{2}}dxdy$

% Monte Carlo Integration for f(x, y) = sqrt(1 - x^2) over x^2 + y^2 <= 1
clc; clear;

% Number of random points
N = 1e6; % You can increase this for better accuracy

% Initialize sum of function values
f_sum = 0;

% Loop to generate random points and calculate contributions
for i = 1:N
    % Generate random (x, y) within the bounding box [-1, 1] x [-1, 1]
    x = -1 + 2*rand(); % Random x in [-1, 1]
    y = -1 + 2*rand(); % Random y in [-1, 1]
    
    % Check if the point is inside the circle
    if x^2 + y^2 <= 1
        f_sum = f_sum + sqrt(1 - x^2); % Accumulate the function value
    end
end

% Calculate area of the bounding box
A_box = 4; % The bounding box [-1, 1] x [-1, 1]

% Calculate the integral estimate
integral_value = A_box * f_sum / N;

% Display result
fprintf('Estimated value of the integral: %.6f\n', integral_value);

2.1 代码解释：

1)随机点生成：

在 [−1,1]×[−1,1] 内均匀生成随机点。
使用条件 $x^{2}+y^{2}\leq 1$ 筛选落在单位圆内的点。

2)函数值累加：

对满足条件的点，计算 $\sqrt{1-x^{2}}$ 并累加到 f_sum。

3)积分估计公式：

估计积分值为：

这里的区域面积 Abox=4 是整个采样的矩形面积。

4)效率：

N越大，估计值越准确。
通过筛选 $x^{2}+y^{2}\leq 1$ ，只在实际目标区域内计算函数值。

2.2 运算结果

$N$	$I$
N = 1e2	2.655043
N = 1e4	2.685272
N = 1e6	2.666568
N = 1e8	2.666756

3 最优化计算Monte Carlo法

求下列函数的最大值：

$f(x)=(1-x^{3})sin(3x),\: \; \; \; \; -2pi<x<2pi$

为了方便理解，先绘制这个函数：

% Define the function f(x)
f = @(x) (1 - x.^3) .* sin(3 * x);

% Define the range for x
x = linspace(-2*pi, 2*pi, 1000); % Generate 1000 points in the range [-2*pi, 2*pi]

% Compute the function values
y = f(x);

% Plot the function
figure;
plot(x, y, 'b-', 'LineWidth', 1.5);
grid on;

% Add labels and title
xlabel('x');
ylabel('f(x)');
title('Plot of f(x) = (1 - x^3)sin(3x)');
legend('f(x) = (1 - x^3)sin(3x)', 'Location', 'Best');

matlab运行结果如下：

在给出计算代码：

% Optimization using Monte Carlo for f(x) = (1 - x^3) * sin(3x)
clc; clear;

% Number of random samples
N = 1e6; % Increase this for higher accuracy

% Define the function
f = @(x) (1 - x.^3) .* sin(3 * x);

% Generate random samples in the range [-2*pi, 2*pi]
x_samples = -2*pi + (2*pi - (-2*pi)) * rand(N, 1);

% Evaluate the function for each sample
f_values = f(x_samples);

% Find the maximum function value
f_max = max(f_values);

% Find the corresponding x value(s) for the maximum
x_max = x_samples(f_values == f_max);

% Display results
fprintf('Maximum value of f(x): %.6f\n', f_max);
fprintf('At x = %.6f (one of the possible values)\n', x_max(1));