蒙特卡洛（随机试验）法计算π

学校课程有道题要求计算π值，蒙特卡洛是利用随机试验求解问题的方法，课程要求随机10**8次

最开始的代码：

from random import random
from time import clock

times = 100000000
hits = 0
clock()
for i in range(1, times):
    x = random()
    y = random()
    dist = x ** 2 + y ** 2
    if dist <= 1.0:
        hits = hits + 1
pi =  4*(hits/times)
print('Pi is %s' % pi)
print('time is %-5.5ss' % clock())

结果：

速度很慢，以为是循环中多次调用random()太耗时，所以第二版：

import numpy as np
from time import clock

# this function use np.random.uniform(0,1,a),but seems it's more slower

times = 100000000
hits = 0
clock()
x = np.random.uniform(0, 1, times)
y = np.random.uniform(0, 1, times)
for i in range(1, times):
    dist = x[i]**2 + y[i]**2
    if dist <= 1.0 and dist>=-1.0:
        hits = hits + 1
pi =  4*(hits/times)
print('Pi is %s' % pi)
print('time is %-5.5ss' % clock())

使用np.random.uniform()直接生成一亿个随机数，结果：

更慢了（- -！），数组的读取会让速度更慢，考虑利用矩阵进行计算：

import numpy as np
from time import clock

# this function use np.where

times = 100000000
hits = 0
clock()
numbers=np.random.rand(2,times)
hits += np.sum(np.where((numbers[0]**2 + numbers[1]**2) < 1, 1, 0))
pi =  4*(hits/times)
print('Pi is %s' % pi)
print('time is %-5.5ss' % clock())

矩阵运算相当于并行运算，相比for循坏的串行，速度快了很多，结果：

np.where参考：https://www.cnblogs.com/massquantity/p/8908859.html