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

蒙特卡罗方法的代码实现——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次，并在表中记下均值和方差。

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

• 条形图

  • 均值
    

  • 方差
    

• 代码

  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:

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.
我们现在尝试通过蒙特卡洛的方法求解如下的积分：

该积分的求解我们可以直接求解，即有∫_(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

• 条形图

  • 均值
    

  • 方差
    

• 代码

  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:

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.
我们现在尝试通过蒙特卡洛的方法求解如下的更复杂的积分：

你能够通过公式直接求解上述的积分吗？如果你用蒙特卡洛的方法求解该积分，你认为(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

• 条形图

  • 均值
    

  • 方差
    

• 代码：

  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