从蒙特卡洛方法计算pi值谈random模块

计算机模拟常常需要用到随机选择的数。本文从随机数的一个简单应用开始简要地介绍Python的random模块。

 

使用蒙特卡洛方法计算pi值

 

Links:该问题来自于pudure university（普渡大学）python课程中的problem set2

Monte Carlo methods are used to simulate complex physical and mathematical systems by repeated random sampling. In simple terms, given a probability, p, that an event will occur in certain conditions, a program generates those conditions repeatedly. The number of times the event occurs divided by the number of times the conditions are generated should be approximately equal to p.

蒙特卡洛方法通过重复随机取样用于复杂的物理和数学系统的仿真。简单地说，给定一个变量p,用于描述事件在一定条件下发生的概率，程序重复地生成这些条件。事件发生的次数除以总产生次数应该近似地等于概率p。

The Monte Carlo method can be used to generate an approximate value of pi. The figure below shows a unit square with a quarter of a circle inscribed. The area of the square is 1 and the area of the quarter circle is pi/4. Using a random number generator, imagine “throwing” random points at the square. The ratio between the number of points that fall inside the circle (red points) and the total number of points thrown (red and green points) gives an approximation to the value of pi/4. This process is a Monte Carlo simulation approximating pi.

蒙特卡洛方法可以用于产生接近pi的近似值。下图显示了一个带有1/4圆在内的正方形单元、落在圈内（红点）的点和总的投在正方形（红和绿点）上的点的比率给出了pi/4的近似值。这一过程称为使用蒙特卡洛方法来仿真逼近pi实际值。

 

Write a Python program that runs this simulation. Your program will run a simulation for a specified number of iterations. let n be the variable representing the number of iterations. For a given n, executes n iteration, with each iteration generating a random point (x,y) and determining whether the point lies inside the circle or not. Assume that (O,O) is the lower left corner of the square. After the n iterations are completed, the program outputs the approximation of pi using the ratio of the points inside the circle and outside the circle (need to multiply the ratio by 4). Also output the number of iterations executed and the value of math.pi.

编写Python程序来运行此仿真，运行你的程序来模拟对一定数量数目的迭代。设n是表示迭代数量的变量n。对于给定一个n，执行n次迭代，每次迭代产生一个随机的点(x,y)并且决定该点是否落在圈内。假设(0,0)是正方形左下角的坐标。n次迭代完成后，程序使用落在圈内和圈外的比值来输出pi的近似值（需要将ratio乘以4噢）

Your program should run for n = 100, 1000, 10000, 100000, and 1000000. Do not paste in the code five times, but use a for-loop around your simulation code that executes the simulation five times.

你的程序需要运行n = 100, 1000, 10000, 100000, and 1000000。不要复制代码5次，而是使用for循环来代替。 

Comments:

• A random point is generated by generating two random numbers, each between 0 and 1 (including 0, but not 1). See Lab 2 how to generate random numbers using library random.

• A point (x,y) lies inside the circle if sqrt(x^2 + y^2) <1.

• When testing your code, run it only once (not 5 times).

• Note that for the same value of n, you will see slightly different answers for pi as each execution uses a new set of random points.

　　说的通俗点，也就是在一个单位长度为1*1的正方形上随机地落点n次，比较落在内切1/4圆中的点和落在整个正方形上的点的比值。

　　取n= 100, 1000, 10000, 100000, 1000000，计算出得出的pi值与实际的pi 的误差。

　　那么我们就需要用到Python中的random模块的相关函数，用random()函数来随机地生成取值在（0,1）的点(x,y)。

　　给出相关代码：

import random

import math

 

def main():

# You didn't know how to do the below loop for this problem set, but you

# should know what it does now. It creates a list that contains all of

# the number of iterations you were supposed to test. Try it!

for n in [10, 100, 1000, 10000, 100000, 1000000]:

 total = 0

for i in xrange(n):

 x, y = random.random(), random.random()

if math.sqrt(x **2+ y **2) <1.0:

 total +=1

 mypi =4.0*total / n

print'Estimating pi with', n, 'iterations:', mypi

print'Value of math.pi is', math.pi

print'Error is', abs(math.pi - mypi) / math.pi

print

 

main()

　　

　　以上，我们介绍了random模块的一个简单应用，接下来是谈谈random模块的时候了：

random.random 

 

     random.random()用于生成一个0到1的随机符点数: 0 <= n < 1.0

 

>>> random.random()

0.72759757730421304

random.uniform

　　random.uniform的函数原型为：random.uniform(a, b)，用于生成一个指定范围内的随机符点数，两个参数其中一个是上限，一个是下限。如果a < b，则生成的随机数n: a <= n <= b。如果 a >b， 则 b <= n <= a。

 

>>> print random.uniform(10, 20)

17.5255098383

>>> print random.uniform(20, 10)

16.6823931069

 

random.randint

    random.randint()的函数原型为：random.randint(a, b)，用于生成一个指定范围内的整数。其中参数a是下限，参数b是上限，生成的随机数n: a <= n <= b

>>> print random.randint(10, 20)

18

>>> print random.randint(20, 20)

20

>>> print random.randint(30, 20)

 

Traceback (most recent call last):

  File "<pyshell#100>", line 1, in <module>

    print random.randint(30, 20)

  File "C:\Python26\lib\random.py", line 228, in randint

    return self.randrange(a, b+1)

  File "C:\Python26\lib\random.py", line 204, in randrange

    raise ValueError, "empty range for randrange() (%d,%d, %d)" % (istart, istop, width)

ValueError: empty range for randrange() (30,21, -9)

random.randrange

 

　　random.randrange的函数原型为：random.randrange([start], stop[, step])，从指定范围内，按指定基数递增的集合中获取一个随机数。如：random.randrange(10, 100, 2)，结果相当于从[10, 12, 14, 16, ... 96, 98]序列中获取一个随机数。random.randrange(10, 100, 2)在结果上与 random.choice(range(10, 100, 2) 等效。

>>> random.randrange(10, 100, 2)

30

random.choice

　　random.choice从序列中获取一个随机元素。其函数原型为：random.choice(sequence)。参数sequence表示一个有序类型。这里要说明 一下：sequence在python不是一种特定的类型，而是泛指一系列的类型。list, tuple, 字符串都属于sequence。有关sequence可以查看python手册数据模型这一章。下面是使用choice的一些例子：

1. print random.choice("学习Python")   
2. print random.choice(["JGood", "is", "a", "handsome", "boy"])  
3. print random.choice(("Tuple", "List", "Dict"))   

 

第一个涉及到编码，在Python内部统一使用的是Unicode编码，我们得把这个字符串解码为Unicode：

>>> x='学习Python'

>>> y=x.decode('gbk')

>>> y

u'\u5b66\u4e60Python'

>>> print random.choice(y)

学

例子2和3的运行结果如下：

>>> print random.choice(["JGood", "is", "a", "handsome", "boy"])

handsome

>>> print random.choice(("Tuple", "List", "Dict"))

Tuple

random.shuffle

　　random.shuffle的函数原型为：random.shuffle(x[, random])，用于将一个列表中的元素打乱，注意该方法会改变原列表。如:

>>> p = ["Python", "is", "powerful", "simple", "and so on..."]

>>> random.shuffle(p)

>>> print p

['powerful', 'Python', 'simple', 'and so on...', 'is']

random.sample

　　random.sample的函数原型为：random.sample(sequence, k)，从指定序列中随机获取指定长度的片断。sample函数不会修改原有序列。

>>> list = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

>>> slice = random.sample(list, 5)

>>> print slice

[3, 6, 9, 5, 7]

>>> print list

[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

random.seed

伪随机数生成模块。如果不提供 seed，默认使用系统时间。使用相同的 seed，可以获得完全相同的随机数序列，常用于算法改进测试。

>>> a=random.Random()

>>> a.seed(1)

>>> [a.randint(1, 100) for i in range(20)]

[14, 85, 77, 26, 50, 45, 66, 79, 10, 3, 84, 44, 77, 1, 45, 73, 23, 95, 91, 4]

>>> b=random.Random()

>>> b.seed(1)

>>> [b.randint(1, 100) for i in range(20)]

[14, 85, 77, 26, 50, 45, 66, 79, 10, 3, 84, 44, 77, 1, 45, 73, 23, 95, 91, 4]

 

以下为用于计算三角、β分布、指数分布、伽马分布、高斯分布等专业的随机函数

参考来自：http://gchsuperman.blog.163.com/blog/static/16456746820101023112528365/

(1)random.triangular(low, high, mode) 三角

Return a random floating point number N such that low <= N <= high and with the specified mode between those bounds. The low and highbounds default to zero and one. The mode argument defaults to the midpoint between the bounds, giving a symmetric distribution.

 

(2)random.betavariate(alpha, beta)β分布

Beta distribution. Conditions on the parameters are alpha > 0 and beta > 0. Returned values range between 0 and 1.

 

(3)random.expovariate(lambd)指数分布

Exponential distribution. lambd is 1.0 divided by the desired mean. It should be nonzero. (The parameter would be called “lambda”, but that is a reserved word in Python.) Returned values range from 0 to positive infinity if lambd is positive, and from negative infinity to 0 if lambd is negative.

 

(4)random.gammavariate(alpha, beta)伽马分布

Gamma distribution. (Not the gamma function!) Conditions on the parameters are alpha > 0 and beta > 0.

 

(5)random.gauss(mu, sigma)高斯分布

Gaussian distribution. mu is the mean, and sigma is the standard deviation. This is slightly faster than the normalvariate() function defined below.

 

(6)random.lognormvariate(mu, sigma)对数正态分布

Log normal distribution. If you take the natural logarithm of this distribution, you’ll get a normal distribution with mean mu and standard deviation sigma. mu can have any value, and sigma must be greater than zero.

 

(7)random.normalvariate(mu, sigma)正态分布

Normal distribution. mu is the mean, and sigma is the standard deviation.

 

(8)random.vonmisesvariate(mu, kappa)  循环数据分布

mu is the mean angle, expressed in radians between 0 and 2*pi, and kappa is the concentration parameter, which must be greater than or equal to zero. 

If kappa is equal to zero, this distribution reduces to a uniform random angle over the range 0 to 2*pi.

 

(9) random.paretovariate(alpha)帕累托分布

Pareto distribution. alpha is the shape parameter.

 

(10) random.weibullvariate(alpha, beta) 威布尔分布

Weibull distribution. alpha is the scale parameter and beta is the shape parameter

 

OK，bla...bla说了一大堆，咱们再来看看

《Python核心编程》第5章的练习题

 

5-17.* 随机数。熟读随机数模块然后解下面的题： 

生成一个有 N 个元素的由随机数 n 组成的列表， 其中 N 和 n 的取值范围分别为： (1 <N<= 100), (0 <= n <= 2**31  -1)。然后再随机从这个列表中取 N (1 <= N <= 100)个随机数出来， 对它们排序，然后显示这个子集。

现在是不是可以轻松搞定了呢？