Monte Carlo method - 倒水概率问题

http://www.cnblogs.com/EdwardLiu/p/6264968.html
http://www.stealthcopter.com/blog/2009/09/python-calculating-pi-using-random-numbers/

As we should know _pi_ is the ratio of circle’s radius to its circumference, which is conveniently the same as the ratio of a circle’s area to the square of its radius (wiki…)

So what we are going to be doing is picking lots of random coordinates in an x-y grid and calculating if they are within the circle or the square.

We will assign the radius to be 1, because that makes it easy to work with. By default a random number in python ( random() ) will return a floating point number between 0 and 1. To test if a point is within a circle we simply use Pythagoras.

So if the sqrt(a**2+b**2)<=1 then the point lies inside the circle’s radius. In the diagram above we see that point A lies within the circle, and point B lies outside the circle.

We can really don’t need to use the whole circle as it has symmetry, so we can just take a quartre, which makes the generating of random numbers easier as you only need to use a random number for x and y between 0 and 1, rather than -1 and 1. It will look like the diagram below.

Now for a confusing bit of maths. We are calculating the ratio of the area of a circle to the area of a square.

from random import *
from math import sqrt
inside=0
n=1000
for i in range(0,n):
 x=random()
 y=random()
 if sqrt(x*x+y*y)<=1:
  inside+=1
pi=4*inside/n
print pi

So we can see that the program quickly solves pi to about two decimal places, but it is a terribly inefficient method and will struggle to get much more accuracy than this.

倒水概率问题 | 书影博客
将 1L 水依次倒入 3 个杯子（倒完前 2 杯后，剩余水全部放入第 3 杯），求任意 1 杯水大于 0.5L 的概率。
假设杯子容量大于 1L ，倒出水的体积服从[0, V]上的均匀分布， V 为剩余水的体积。

任意 1 杯水大于 0.5L，可以分解为下面三种情况：

P1. 倒出第一杯水超过 0.5L 的概率为 0.5

P2. 倒出第二杯水超过 0.5L 的概率：

  设倒出第一杯水的体积为 x （ x < 0.5 ），则有：P2 = ∫（ 1 - x - 0.5) / (1 - x)   x ∈ [0, 0.5]

P3. 倒出第三杯水超过 0.5L 的概率：

  等价于求倒出前两杯水的和小于 0.5L 的概率

  设倒出第一杯水的体积为 x （ x < 0.5 ），则有：P3 = ∫（ 0.5 - x) / (1 - x)    x ∈ [0, 0.5]

P1 + P2 + P3 = 0.5 + ∫ (2 * x - 1) / (x - 1)    x ∈ [0, 0.5]

上式中 ∫ (2 * x - 1) / (x - 1) 的原函数为： 2x + ln |x - 1| + C

最终答案为P = 1.5 + ln(0.5) ≈ 0.807

class PourWater(object): def montecarlo(self, num): cnt = 0 for x in range(num): p1 = random.uniform(0, 1) p2 = random.uniform(0, 1 - p1) if max(p1, p2, 1 - p1 - p2) > 0.5: cnt += 1 return 1.0 * cnt / num print PourWater().montecarlo(100000)
Read full article from 倒水概率问题 | 书影博客