## Sunday, November 13, 2016

### 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.

```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]

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)