Python：每日一题 56

ooxx7788

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

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

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

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

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


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

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

冬雪雪冬

先写一个，速度比较慢。

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)) 

SixPy

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

蒙特卡洛算法

jerryxjr1220

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

哨子1122

学习

solomonxian

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

import random
import math
import time
import numpy

# 通过参数半径r及实验次数trial返回pi估计值，涉及平方和用整数可能会好算一点，集合可能会缩小概率，还是用列表

# 直接在推导式中判断，然后返回
def cal_pi_1(r,trial):
    return 4*len([1 for i in range(trial) if (random.randint(-r,r)**2+random.randint(-r,r)**2)<=r**2])/ trial

# 将生成和判断分开两半，瞅瞅有没有变化
def cal_pi_2(r,trial):
    q = [(random.randint(-r,r),random.randint(-r,r)) for i in range(trial)]
    w = [1 for i in q if i[0]**2+i[1]**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)[0])/trial)

r = 100000,trial = 100000
结果

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

shigure_takimi

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%

咕咕鸡鸽鸽

看看

yu123py

好像不难。

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)

holiday_python

学习

19971023

1

nononoyes

666

工口漫画家

学习

Yedada

nice

我就是来水一水

蒙特卡洛算法