r语言蒙特卡洛法求积分蒙特卡洛算法求积分

转载

mob64ca1414c613 2023-12-13 01:52:33

文章标签 r语言蒙特卡洛法求积分算法随机数定积分代码实现 文章分类 R语言后端开发

蒙特卡洛算法是一种基于随机数的算法，思想非常直观但是与其他算法相比精确度较低。其实我们在中学已经初步了解过蒙特卡洛算法的应用，例如通过在坐标轴中随机取点来估算圆周率pi的值。蒙特卡洛算法的其他应用本质上也和随机取点估值问题类似。

数值方法求定积分

对于定积分

$\int _a ^bf(x)dx$

，一般我们会找原函数然后通过牛顿莱布尼兹公式求解。但是遇到一些复杂的表达式可能就很难求出原函数。例如计算定积分：

$I=\int _{0.8}^{3.3}\frac{1}{1+\sin{x}+ln^2x}dx$

，很难通过一般的方法求解。这时候就要用到蒙特卡洛算法计算出一个近似的值。因为：

$\int _a^bf(x)dx=\displaystyle\lim_{n\rightarrow \infty}\frac{b-a}{n}\displaystyle\sum_{i=1}^nf(x_i)$

，我们可以在积分区间[a,b]上均匀地随机生成一列xi，之后求出积分的近似值。比如计算

$I=\int _{0.8}^{3.3}\frac{1}{1+\sin{x}+ln^2x}dx$

，就可以先在[0.8,3.3]上随机生成n个点，当n足够大时，便可以用近似值代替估算值。

Matlab代码如下：

clc
clear all;
close all;
syms x
f=1/(1+sin(x)+(log(x))^2)
n=1e3
x=unifrnd(0.8,3.3,[1,n]);
ans=0;
for i=1:n
    y(i)=1/(1+sin(x(i))+(log(x(i))^2))
    ans=ans+y(i)
end
ans=2.5*ans/n

最终运行结果是1.12~1.13左右，已经很接近真实值。但是显然需要n很大时才能保证更大的精度。

下雨问题

如果一周每天都有50%的可能下雨，求连续三天下雨的概率？

这里假设每一天是否下雨是相互独立的，记

$Y_i=1_{\{the\quad i^{th} day \quad rains\}}$

，则{Yi}独立。通过蒙特卡洛随机数的方法来模拟Yi的值：在[0,1]之间均匀分布任取一个数a，若a小于0.5则不下雨，否则下雨。之后再累加下雨天数，记录连续的下雨天数即可。对应的matlab程序如下：

n=1e4;ans=0;
for i=1:n
    c=0;y=0;p=0.5;
    for t=1:7
         r=unifrnd(0,1)
        if r<p
            c=c+1
        else
            c=0
        end
        if c>=3
            y=1
        end
    end
    ans=ans+y
end
ans=ans/n

n=1e4表示重复实验10000个星期，得出的结果为有3657个星期是有连续三天下雨的。

r语言蒙特卡洛法求积分蒙特卡洛算法求积分_算法_06