您的位置: 首页  >  文章  >  蒙特卡罗方法的代码实现——python

蒙特卡罗方法的代码实现——python

分类: 文章 • 2024-09-08 18:39:40

蒙特卡罗方法的代码实现——python

简单抽样

Exercise 1.

The Monte Carlo method can be used to generate an approximate value of pi. The figure below shows a unit square with a quarter of a circle inscribed. The area of the square is 1 and the area of the quarter circle is pi/4. Write a script to generate random points that are distributed uniformly in the unit square. The ratio between the number of points that fall inside the circle (red points) and the total number of points thrown (red and green points) gives an approximation to the value of pi/4. This process is a Monte Carlo simulation approximating pi. Let N be the total number of points thrown. When N=50, 100, 200, 300, 500, 1000, 5000, what are the estimated pi values, respectively? For each N, repeat the throwing process 100 times, and report the mean and variance. Record the means and the corresponding variances in a table.
蒙特卡洛方法可以用于产生接近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次，并在表中记下均值和方差。
蒙特卡罗方法的代码实现——python
Figure 1 蒙特卡洛方法求解pi

answer

  • 结果

    (均值保留四位小数，方差保留四位小数)

    SIZE MEAN variance
    20 3.1000 0.1660
    50 3.1120 0.0524
    100 3.1940 0.0354
    200 3.1300 0.0135
    300 3.1320 0.0079
    500 3.1448 0.0053
    1000 3.1360 0.0029
    5000 3.1470 0.0005

  • 条形图

    • 均值
      蒙特卡罗方法的代码实现——python

    • 方差
      蒙特卡罗方法的代码实现——python

  • 代码

    import numpy as np


def is_in(location):
    distance = location[0] ** 2 + location[1] ** 2
    distance = distance ** 0.5
    if distance < 1:
        return True
    else:
        return False


sizes = [20, 50, 100, 200, 300, 500, 1000, 5000]
mean = [0.0] * 8
var = [0.0] * 8

for num in range(0, 8):
    size = sizes[num]
    result = [0.0] * 20

    for time in range(0, 20):
        points = np.random.rand(size, 2)

        for one in points:
            if is_in(one):
                result[time] += 1
        result[time] /= size
        result[time] *= 4

    results = np.array(result)
    mean[num] = results.mean()
    var[num] = results.var()

print mean
print var

求简单积分

Exercise 2

We are now trying to integrate the another function by Monte Carlo method:
蒙特卡罗方法的代码实现——python
A simple analytic solution exists here: ∫_(x=0)1▒x3 =1/4. If you compute this integration using Monte Carlo method, what distribution do you use to sample x? How good do you get when N = 5, 10, 20, 30, 40, 50, 60, 70, 80, 100, respectively? For each N, repeat the Monte Carlo process 20 times, and report the mean and variance of the integrate in a table.
我们现在尝试通过蒙特卡洛的方法求解如下的积分：
蒙特卡罗方法的代码实现——python
该积分的求解我们可以直接求解，即有∫_(x=0)1▒x3 =1/4。如果你用蒙特卡洛的方法求解该积分，你认为x可以通过什么分布采样获得？如果采样次数是分别是N = 5, 10, 20, 30, 40, 50, 60, 70, 80, 100，积分结果有多好？对于每个采样次数N，重复蒙特卡洛过程100次，求出均值和方差，然后在表格中记录对应的均值和方差。

answer

  • 结果

    （均值保留六位小数，方差保留七位小数）

    SIZE MEAN VARIANCE
    5 0.251013 0.0009399
    10 0.249390 0.0001036
    20 0.250202 0.0000165
    30 0.250305 0.0000052
    40 0.250036 0.0000026
    50 0.249897 0.0000013
    60 0.249802 0.0000008
    70 0.250031 0.0000004
    80 0.249972 0.0000002
    100 0.250015 0.0000001

  • 条形图

    • 均值
      蒙特卡罗方法的代码实现——python

    • 方差
      蒙特卡罗方法的代码实现——python

  • 代码

    import numpy as np

sizes = [5, 10, 20, 30, 40, 50, 60, 70, 80, 100]
mean = [0.0] * 10
var = [0.0] * 10

for num in range(0, 10):

    size = sizes[num]
    result = [0.0] * 100

    for time in range(0, 100):

        locations = np.random.rand(size, 1)

        for one in range(0, size):
            x = locations[one]/size + one*1.0/size
            result[time] += pow(x, 3)

        result[time] /= size

    results = np.array(result)
    mean[num] = results.mean()
    var[num] = results.var()


print mean
print var

求两重积分

Exercise 3

We are now trying to integrate a more difficult function by Monte Carlo method that may not be analytically computed:

蒙特卡罗方法的代码实现——python

Can you compute the above integration analytically? If you compute this integration using Monte Carlo method, what distribution do you use to sample (x,y)? How good do you get when the sample sizes are N = 5, 10, 20, 30, 40, 50, 60, 70, 80, 100, 200 respectively? For each N, repeat the Monte Carlo process 100 times, and report the mean and variance of the integrate.
我们现在尝试通过蒙特卡洛的方法求解如下的更复杂的积分：

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

answer

  • 结果：

    size mean variance
    10 113970 11186015413
    20 104968 5703306581
    30 111917 4574907015
    40 109778 3102348467
    50 118887 2743787666
    60 118691 2009142994
    70 110826 1862256602
    80 118461 1240587892
    100 113951 1251562500
    200 110962 666497079
    500 112985 282780621

  • 条形图

    • 均值
      蒙特卡罗方法的代码实现——python

    • 方差
      蒙特卡罗方法的代码实现——python

  • 代码：

    import numpy as np
import math


def calculate(_x, _y):
    cal_result = []
    for index in range(0, len(_x)):
        temp1 = _y[index] ** 2 * math.exp(-(_y[index] ** 2))
        temp2 = _x[index] ** 4 * math.exp(-(_x[index] ** 2))
        temp3 = _x[index] * math.exp(-(_x[index] ** 2))
        cal_result.append((temp1 + temp2) / temp3)

    return cal_result


sizes = [10, 20, 30, 40, 50, 60, 70, 80, 100, 200, 500]
mean = []
var = []
floor_space = 2.0 * 2.0

for num in range(0, len(sizes)):
    size = sizes[num]
    result = []

    for time in range(0, 100):

        # create the random data
        x = np.random.rand(size, 1)
        x *= 2
        x += 2

        y = np.random.rand(size, 1)
        y *= 2
        y -= 1
        temp = calculate(x, y)
        ones = np.array(temp)
        one = ones.mean() * floor_space
        result.append(one)

    results = np.array(result)
    # print results
    mean.append(results.mean())
    var.append(results.var())

print mean
print var

相关推荐