【数据挖掘导论】HW1

最新推荐文章于 2022-09-10 00:01:06 发布

我叭想起名啊

最新推荐文章于 2022-09-10 00:01:06 发布

阅读量515

点赞数

分类专栏：数据挖掘导论文章标签：蒙特卡洛方法

本文链接：https://blog.csdn.net/qq_36318771/article/details/88970983

版权

数据挖掘导论专栏收录该内容

3 篇文章 0 订阅

订阅专栏

本文介绍了如何使用蒙特卡洛方法求解积分问题，包括π的近似值、一元和二元积分。通过在不同采样点数下进行多次实验，分析了积分结果的均值和方差，展示了随着采样点增加，积分结果的精确度提高。

摘要由CSDN通过智能技术生成

蒙特卡洛方法求积分

Exercise 1.

蒙特卡洛方法可以用于产生接近 $\pi$ 的近似值。图1显示了一个带有1/4内切圆在内的边长为1的正方形。正方形的面积是1，该1/4圆的面积为 $\pi/4$ 。通过编程实现在这个正方形中产生均匀分布的点。落在圈内（红点）的点和总的投在正方形（红和绿点）上的点的比率给出了 $\pi/4$ 的近似值。这一过程称为使用蒙特卡洛方法来仿真逼近 $\pi$ 实际值。令N表示总的投在正方形的点。当投点个数分别是20, 50, 100, 200, 300, 500, 1000, 5000时， $\pi$ 值分别是多少？对于每个N，每次实验算出 $\pi$ 值，重复这个过程20次，并在表中记下均值和方差。
在这里插入图片描述

Exercise 2.

我们现在尝试通过蒙特卡洛的方法求解如下的积分：
$\int_0^1 x^3$
该积分的求解我们可以直接求解，即有 $\int_{x=0}^1 x^3=1/4$ 。如果你用蒙特卡洛的方法求解该积分，你认为x可以通过什么分布采样获得？如果采样次数是分别是N = 5, 10, 20, 30, 40, 50, 60, 70, 80, 100，积分结果有多好？对于每个采样次数N，重复蒙特卡洛过程100次，求出均值和方差，然后在表格中记录对应的均值和方差。

Exercise 3:

我们现在尝试通过蒙特卡洛的方法求解如下的更复杂的积分：
$\int_{x=2}^4\int_{y=-1}^1 f(x,y)=\frac{y^2*e^{-y^2}+x^4*e^{-x^2}}{x*e^{-x^2}}$

你能够通过公式直接求解上述的积分吗？如果你用蒙特卡洛的方法求解该积分，你认为(x, y)可以通过什么分布采样获得？如果点（x, y）的采样次数是分别是N = 10, 20, 30, 40, 50, 60, 70, 80, 100, 200, 500, 积分结果有多好？对于每个采样次数N，重复蒙特卡洛过程100次，求出均值和方差，然后在表格中记录对应的均值和方差。

解答

Monte Carlo Methods的两种表达：

Monte Carlo投点估计法。利用落在有效范围内的点和总的投点数的比率的方法来仿真逼近实际值。
Monte Carlo公式计算法。
- 一重积分 $F=\int_a^bf(x)dx$ ，利用公式 $<F^N>=(b-a)*\frac{1}{N}\sum_{i=0}^{N-1}f(X_i)$ 求解。
- 双重积分 $F=\int_{a}^b\int_{c}^df(x,y)dxdy$ ，利用公式 $<F^N>=(b-a)*(d-c)*\frac{1}{N}\sum_{i=0}^{N-1}f(X_i,Y_i)$ 求解。

第一题

Way1—投点估计法

function [mu,sigma_2,pis] = MonteCarloPI(N)
    round = 100;

    pis = zeros(1,round);
    for i=1:round
        x=rand(1,N);
        y=rand(1,N);

        dis_2=x.*x+y.*y;    % 采样点与远点的距离平方
        pis(i)=numel(find(dis_2<1))/N*4;    %pi=距离平方<r^2的点的个数/采样点总数*4

    end
    mu=mean(pis);       % 计算pi数组均值
    sigma_2=var(pis);   % 计算pi数组方差

每轮循环，通过rand函数产生服从均匀分布的N个x,y坐标样本
计算每点到原点的距离的平方
记录距离平方< $r^2$ 的点的个数与采样点总数的比值，即为 $\pi/4$
计算pis数组的均值和方差

测试表格：

N	20	50	100	200	300	500	1000	5000
$\mu$	3.1400	3.1296	3.1252	3.1466	3.1424	3.1497	3.1443	3.1416
$\sigma^2$	0.1232	0.0621	0.0257	0.0145	0.0085	0.0067	0.0023	0.0005

Way2—公式计算法
```
function [mu,sigma_2,pis] = MonteCarloMethodsPI(N)
    round = 100;

    pis = zeros(1,round);
    for i=1:round
        x=rand(1,N);        % 得到随机x样本
        y=sqrt(1-x.^2);     % 计算对应的y值,y=sqrt(1-x^2)
        F_N=mean(y);        % 计算y的均值

        pis(i)=F_N*4;       % 通过均值*x的宽度计算面积

    end
    mu=mean(pis);
    sigma_2=var(pis);
```
- 产生随机N个x样本
- 计算对应的函数值y， $y=\sqrt{1-x^2}$
- 计算随机样本对应的y的均值
- 通过近似均值*x的范围宽度模拟面积

第二题

Way1—投点估计法

function [mu,sigma_2,result] = MonteCarloIntegration(N)
    round = 100;
    
    result =zeros(1,round);
    for i=1:round
        x=rand(1,N);
        y=rand(1,N);
        dis_2=y-x.^3;   % 检测随机点是否在图像下方
        
        result(i)=numel(find(dis_2<0))/N; % 计算下方点的个数与样点总数的比值*样本区间面积1
        
    end
    mu=mean(result);
    sigma_2=var(result);

每轮循环，通过rand函数产生服从均匀分布的N个x,y坐标样本
检测随机点是否在图像下方
积分结果=下方点的个数/样点总数的比值*样本区间面积1
计算result数组的均值和方差

Way2—公式计算法

function [mu,sigma_2,result] = MonteCarloMethodsIntegration(N)
    round = 100;
    
    result =zeros(1,round);
    for i=1:round
        x=rand(1,N);    % 产生随机样本x
        y=x.^3;         % 计算样本x对应函数值
        
        result(i)=mean(y);  % 面积=y的均值*x区间范围（1-0）
        
    end
    mu=mean(result);
    sigma_2=var(result);

产生随机N个x样本
计算对应的函数值y， $y=x^3$
计算随机样本对应的y的均值
积分结果=近似均值*x的区间范围

测试表格：

N	10	20	30	40	50	60	70	80	100
$\mu$	0.2480	0.2552	0.2518	0.2541	0.2515	0.2460	0.2524	0.2517	0.2494
$\sigma ^2$	0.0083	0.0046	0.0027	0.0026	0.0016	0.0011	0.0013	0.0009	0.0008

第三题

定义积分函数：

function [f] = DoubleIntegrationFunc(x,y)

    numerator = y.^2.*exp(-y.^2)+x.^4.*exp(-x.^2);
    denominator = x.*exp(-x.^2);
    
    f = numerator./denominator;

Way1—投点估计法

function [mu,sigma_2,result] = MonteCarloDoubleIntegration(N)
    round = 100;
    zmax=1000000;   % 定义最大z值
    
    v=2*2*zmax;     % 测试体积=x区间宽度*y区间宽度*z区间宽度
    result =zeros(1,round);
    for i=1:round
        x=rand(1,N)*2+2;    % 产生随机x样本
        y=rand(1,N)*2-1;    % 产生随机y样本
        
        z=rand(1,N)*zmax;   % 产生随机z样本
        f=DoubleIntegrationFunc(x,y);   % 计算x,y对应函数值
        
        dis_2=z-f;  % 检测点的有效性
        
        result(i)=numel(find(dis_2<0))/N*v; % 积分体积=有效点个数/总个数*测试体积
        
    end
    mu=mean(result);
    sigma_2=var(result);

定义最大z值，计算测试体积=x区间宽度 * y区间宽度 * z区间宽度
产生x,y,z随机样本
计算x,y对应函数值，检测点的有效性
积分体积=有效点个数/总个数*测试体积
计算均值和方差

Way2—公式计算法

function [mu,sigma_2,result] = MonteCarloMethodsDoubleIntegration(N)
    round = 100;
    
    result =zeros(1,round);
    for i=1:round
        x=rand(1,N)*2+2;
        y=rand(1,N)*2-1;
        
        z=DoubleIntegrationFunc(x,y);   % 计算x,y对应函数值
                
        result(i)=2*2*mean(z);  % 积分体积=函数值的均值*x区间宽度*y区间宽度
        
    end
    mu=mean(result);
    sigma_2=var(result);