Python：每日一题 56

ooxx7788

对于今天的题目，我原本只想写这样一句。
利用蒙特卡洛算法，计算pi.

但想想估计这样很多人都不知道从哪里开始，所以我把步骤稍微分解下。
先介绍下，如何根据蒙特卡洛算法来算出pi的。
基本思想： 利用圆与其外接正方形面积之比为pi/4的关系，通过产生大量均匀分布的二维点，计算落在单位圆和单位正方形的数量之比再乘以4便得到pi的近似值。样本点越多，计算出的数据将会越接近真识的pi（前提时样本是“真正的”随机分布）。

步骤：
1、产生一定数量的位于正方形内部的样本（样本越大越精确）。
2、计算样本中位于圆内的样本数量。
3、利用圆与其外接正方形面积之比为pi/4的关系，计算出pi值
4、和真实的pi值进行对比，获得误差。

网上已经有一些二维的例子，当然你也可以参考。

如果觉得练习不够，还可以试试下面这个三维的球体外接正方体计算pi，道理和步骤都是相似的。


冬雪给了个二维版的，我就给这个三维版的吧。

游客，如果您要查看本帖隐藏内容请回复

冬雪雪冬

先写一个，速度比较慢。

1. import math, random
   
2. a = 0
   
3. n = 10000000 #样本量
   
4. for i in range(n):
   
5.     x = random.random()
   
6.     y = random.random()
   
7.     if x * x + y * y <= 1:
   
8.         a += 1
   
9. pi = a / n * 4
   
10. print('计算出的pi值是:%f'%pi)
    
11. print('与内置pi值的相对误差是:%f%%'%((pi - math.pi) / math.pi * 100))
    
12. 
    

复制代码

SixPy

冬雪雪冬 发表于 2017-6-1 18:33
先写一个，速度比较慢。

蒙特卡洛算法

jerryxjr1220

貌似我的小练习中已经算过了
http://bbs.fishc.com/thread-80125-1-1.html

哨子1122

学习

solomonxian

一样的思想，用了4种写法，结果numpy快了好几个数量级

1. import random
   
2. import math
   
3. import time
   
4. import numpy
   
5. 
   
6. # 通过参数半径r及实验次数trial返回pi估计值，涉及平方和用整数可能会好算一点，集合可能会缩小概率，还是用列表
   
7. 
   
8. # 直接在推导式中判断，然后返回
   
9. def cal_pi_1(r,trial):
   
10.     return 4*len([1 for i in range(trial) if (random.randint(-r,r)**2+random.randint(-r,r)**2)<=r**2])/ trial
    
11. 
    
12. # 将生成和判断分开两半，瞅瞅有没有变化
    
13. def cal_pi_2(r,trial):
    
14.     q = [(random.randint(-r,r),random.randint(-r,r)) for i in range(trial)]
    
15.     w = [1 for i in q if i[0]**2+i[1]**2 <= r**2]
    
16.     return 4*len(w)/ trial
    
17. 
    
18. # 不用推导式，直接循环赋值
    
19. def cal_pi_3(r,trial):
    
20.     p=0
    
21.     for i in range(trial):
    
22.         a,b = random.randint(-r,r),random.randint(-r,r)
    
23.         if a**2+b**2 <= r**2:
    
24.             p +=1
    
25.     return (4*p/trial)
    
26. 
    
27. # 使用 numpy 库，数字比较大，没有定好dtype会出错
    
28. def cal_pi_4(r,trial):
    
29.     m = numpy.random.randint(-r,r+1,trial,dtype=numpy.int64)
    
30.     n = numpy.random.randint(-r,r+1,trial,dtype=numpy.int64)
    
31.     k = pow(m,2)+pow(n,2)
    
32.     return (4*len(numpy.nonzero(k<=r**2)[0])/trial)

复制代码

r = 100000,trial = 100000
结果

1. ##第1个函数耗时  1.0870623s -----:my pi is: 3.13892    ,error approaches to  0.0850732 %
   
2. ##第2个函数耗时  1.1000628s -----:my pi is: 3.14136    ,error approaches to  0.0074055 %
   
3. ##第3个函数耗时  1.0660610s -----:my pi is: 3.1406     ,error approaches to  0.0315971 %
   
4. ##第4个函数耗时  0.0110006s -----:my pi is: 3.14204    ,error approaches to  0.0142394 %

复制代码

shigure_takimi

1. def getPiBy(n=1000): #半径1000，样本数为2000*2000=4000000
   
2.     if n < 10:
   
3.         return '样本太小'
   
4.     else:
   
5.         count = 0
   
6.         total = 4*n**2
   
7.         for x in range(-n,n):
   
8.             for y in range(-n,n):
   
9.                 if x**2+y**2<=n**2:
   
10.                     count+=1
    
11.         pi = count*4/total
    
12.         return '%.7f'%pi
    
13. 
    
14. for i in (50, 1000, 2000, 3000):
    
15.     pi = getPiBy(i)
    
16.     print('采样数：{0}，pi近似值：{1},误差率：{2}%'.format(4*i**2, pi, '%.4f'%(100-float(pi)*(100/3.1415926))))
    
17. 
    
18. ##    采样数：10000，pi近似值：3.1372000,误差率：0.1398%
    
19. ##    采样数：4000000，pi近似值：3.1415470,误差率：0.0015%
    
20. ##    采样数：16000000，pi近似值：3.1415857,误差率：0.0002%
    
21. ##    采样数：36000000，pi近似值：3.1415772,误差率：0.0005%
    
22. 
    

复制代码

咕咕鸡鸽鸽

看看

yu123py

好像不难。

1. import random
   
2. counter_in = 0
   
3. counter_all = 0
   
4. for i in range(10000000):
   
5.     x = random.uniform(0, 2)
   
6.     y = random.uniform(0, 2)
   
7.     if (x-1)**2 + (y-1)**2 <= 1:
   
8.         counter_in += 1
   
9.     counter_all += 1
   
10.    
    
11. print(counter_in, counter_all)
    
12. pi = (counter_in * 4) / counter_all
    
13. print('pi = %f' % pi)

复制代码

holiday_python

学习

19971023

1

nononoyes

666

工口漫画家

学习

Yedada

nice

我就是来水一水

蒙特卡洛算法