Python：每日一题 56
本帖最后由 ooxx7788 于 2017-6-2 08:55 编辑对于今天的题目，我原本只想写这样一句。
利用蒙特卡洛算法，计算pi.
但想想估计这样很多人都不知道从哪里开始，所以我把步骤稍微分解下。
先介绍下，如何根据蒙特卡洛算法来算出pi的。
基本思想： 利用圆与其外接正方形面积之比为pi/4的关系，通过产生大量均匀分布的二维点，计算落在单位圆和单位正方形的数量之比再乘以4便得到pi的近似值。样本点越多，计算出的数据将会越接近真识的pi（前提时样本是“真正的”随机分布）。
步骤：
1、产生一定数量的位于正方形内部的样本（样本越大越精确）。
2、计算样本中位于圆内的样本数量。
3、利用圆与其外接正方形面积之比为pi/4的关系，计算出pi值
4、和真实的pi值进行对比，获得误差。
网上已经有一些二维的例子，当然你也可以参考。
如果觉得练习不够，还可以试试下面这个三维的球体外接正方体计算pi，道理和步骤都是相似的。
冬雪给了个二维版的，我就给这个三维版的吧。
**** Hidden Message ***** 先写一个，速度比较慢。
import math, random
a = 0
n = 10000000 #样本量
for i in range(n):
x = random.random()
y = random.random()
if x * x + y * y <= 1:
a += 1
pi = a / n * 4
print('计算出的pi值是:%f'%pi)
print('与内置pi值的相对误差是:%f%%'%((pi - math.pi) / math.pi * 100))
冬雪雪冬 发表于 2017-6-1 18:33
先写一个，速度比较慢。
蒙特卡洛算法 本帖最后由 jerryxjr1220 于 2017-6-2 11:44 编辑
貌似我的小练习中已经算过了
http://bbs.fishc.com/thread-80125-1-1.html 学习 一样的思想，用了4种写法，结果numpy快了好几个数量级
import random
import math
import time
import numpy
# 通过参数半径r及实验次数trial返回pi估计值，涉及平方和用整数可能会好算一点，集合可能会缩小概率，还是用列表
# 直接在推导式中判断，然后返回
def cal_pi_1(r,trial):
return 4*len()/ trial
# 将生成和判断分开两半，瞅瞅有没有变化
def cal_pi_2(r,trial):
q = [(random.randint(-r,r),random.randint(-r,r)) for i in range(trial)]
w = **2+i**2 <= r**2]
return 4*len(w)/ trial
# 不用推导式，直接循环赋值
def cal_pi_3(r,trial):
p=0
for i in range(trial):
a,b = random.randint(-r,r),random.randint(-r,r)
if a**2+b**2 <= r**2:
p +=1
return (4*p/trial)
# 使用 numpy 库，数字比较大，没有定好dtype会出错
def cal_pi_4(r,trial):
m = numpy.random.randint(-r,r+1,trial,dtype=numpy.int64)
n = numpy.random.randint(-r,r+1,trial,dtype=numpy.int64)
k = pow(m,2)+pow(n,2)
return (4*len(numpy.nonzero(k<=r**2))/trial)
r = 100000,trial = 100000
结果
##第1个函数耗时1.0870623s -----:my pi is: 3.13892 ,error approaches to0.0850732 %
##第2个函数耗时1.1000628s -----:my pi is: 3.14136 ,error approaches to0.0074055 %
##第3个函数耗时1.0660610s -----:my pi is: 3.1406 ,error approaches to0.0315971 %
##第4个函数耗时0.0110006s -----:my pi is: 3.14204 ,error approaches to0.0142394 % def getPiBy(n=1000): #半径1000，样本数为2000*2000=4000000
if n < 10:
return '样本太小'
else:
count = 0
total = 4*n**2
for x in range(-n,n):
for y in range(-n,n):
if x**2+y**2<=n**2:
count+=1
pi = count*4/total
return '%.7f'%pi
for i in (50, 1000, 2000, 3000):
pi = getPiBy(i)
print('采样数：{0}，pi近似值：{1},误差率：{2}%'.format(4*i**2, pi, '%.4f'%(100-float(pi)*(100/3.1415926))))
## 采样数：10000，pi近似值：3.1372000,误差率：0.1398%
## 采样数：4000000，pi近似值：3.1415470,误差率：0.0015%
## 采样数：16000000，pi近似值：3.1415857,误差率：0.0002%
## 采样数：36000000，pi近似值：3.1415772,误差率：0.0005%
看看 好像不难。
import random
counter_in = 0
counter_all = 0
for i in range(10000000):
x = random.uniform(0, 2)
y = random.uniform(0, 2)
if (x-1)**2 + (y-1)**2 <= 1:
counter_in += 1
counter_all += 1
print(counter_in, counter_all)
pi = (counter_in * 4) / counter_all
print('pi = %f' % pi) 学习 1 666 nice
蒙特卡洛算法
页:
[1]