并行与分布式计算-蒙特卡洛法

并行与分布式计算-蒙特卡洛法

1.蒙特卡罗模拟概念

蒙特卡罗方法又称 统计模拟法、 随机抽样技术，是一种随机模拟方法，以 概率和统计理论方法为基础的一种计算方法，是使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。将所求解的问题同一定的概率模型相联系，用电子计算机实现统计模拟或抽样，以获得问题的近似解。为象征性地表明这一方法的概率统计特征，故借用赌城蒙特卡罗命名。

1.1蒙特卡罗法优点

• 1.方法的误差与问题的维数无关。

• 2.对于具有统计性质问题可以直接进行解决。

• 3.对于连续性的问题不必进行离散化处理。

1.2蒙特卡罗法缺点

• 1.对于确定性问题需要转化成随机性问题。

• 2.误差是概率误差。

• 3.通常需要较多的计算步数N。

1.3蒙特卡罗法核心

蒙特卡罗方法是将确定的值通过抽样的方法转换成近似
计算。

• 人为地描述或者构造概率过程。

• 通过计算机产生已知概率分布的随机变量，常用的概率分布有均匀分布，正
  态分布，指数分布，泊松分布等。

[重点步骤:当不知道随机变量的概率模型服从那个分布时，可以使用均匀分布来构造；各种测量的误差、射击命中率、人的身高与体重等服从正态分布；指数分布可用在排队论与可靠性分析中；泊松分布可用在产品检验、排队系统、物理等领域中]

• 建立各种估计量。

• 

2.蒙特卡罗模拟案例

蒙特卡罗模拟最经典的案例就是近似计算π值。

已知：正方形长L=2，圆的半径r=1。所以正方形面积S正=22=4，圆的面积S圆=π®2= π，圆和正方向面积比：P= π /4;

通过蒙特卡罗模拟的方法模拟弓箭射出共m支箭，一共可以射中n支箭，那么射中的几率为P1=n/m；在m足够大的时候 P= P1。

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-heaDZxgA-1678246592744)(C:\Users\ASUS\AppData\Roaming\Typora\typora-user-images\image-20230308112115567.png)]

难点：怎么模拟射箭的场景？

在不考虑箭手熟练程度的问题下，以远点作为坐标中心，建立坐标系，那么x服从在[-1,1]上的均匀分布，y也是一样。

整个过程：

• 在服从在[-1,1]范围内均匀分布的总体中随机抽取2个值分别作为x和y的值；

• 通过验证x^2+y2≤1判断是与否；

• 将以上过程重复m（m是一个非常大的整数）次，其中一共有n次符合x^2+y2≤1 ；

• 将概率的理论值和验证值等式变化计算出π的值。

  import random
  import sys
  
  def one():
      x = random.uniform(-1.0,1.0)
      y = random.uniform(-1.0,1.0)
      if x**2 + y**2 <= 1:
          i = 1
      else:
          i = 0
      return i
  def m(num):
      result = 0
      for _ in range(num):
          i = one()
          result += i
      return result
  def pi_value(num):
      result = m(num)
      pi = 4*result/num
      return pi
  if __name__=='__main__':
      
      num = int(sys.argv[1])
      pi = pi_value(num)
      print('pi模拟值为{}'.format(pi))
  

3.蒙特卡罗模拟编程

思考：

在总体样本服从标准正太分布（N（0,1））的产品中随机抽取样本（样本之间独立），那么平均第一次抽到样本大于等于3或者小于等于-3的样本是第几次?

提示：

1. 从服从N（0,1）的总体中随机抽取一个样本；

2. 判断这个样本的大小，若这个样本的数值大小在(-3,3)之间证明这个样本是合格的。继续抽样，直到抽取的样本的数值在(- ∞ ,-3]U[3,∞)的范围内表示产品不合格停止抽样，记录总共抽取的样本总数量 𝑚𝑖(i 表示实验次数)，一次实验完成。

3. 重复n次步骤1~2，然后计算 𝑚𝑖 。

   import multiprocessing
   import random
   import time
   def task(_):
       x = random.uniform(-1,1)
       y = random.uniform(-1,1)
       if x**2 + y**2 < 1:
           output = 1
       else:
           output =0
       return output
   def main(num):
       task_list = [i for i in range(10**7)]
       pool = multiprocessing.Pool(processes = num)
       output_list = pool.map(task,task_list)
       true_num = sum(output_list)
       pi = 4*true_num/(10**7)
       return pi
   t1 = time.time()
   pi = main(12)
   t2 = time.time()
   print("takes :{}s".format(t2 -t1))
   print("pi value:{}".format(pi))
   
   

t1 = time.time()
pi = main(12)
t2 = time.time()
print(“takes :{}s”.format(t2 -t1))
print(“pi value:{}”.format(pi))

从总体为正态分布中抽取一个样本的代码。