6.00.2x学习笔记Week3 Part1
概率问题最早是有赌博引起的。
一个问题。🎲🎲 roll 24次，赌他们必有一次（6，6）是不是划算。
我们知道，roll一次，两个都是6的概率是(1/6)2, 既(1/36)。那么不是两个6的概率就是(35/36)。如果是roll24次，那么不是两个6的概率就是(35/36)24，既0.5085961238690966。用code来表示既：
import random
def rollDie():
"""returns a random int between 1 and 6"""
return random.choice([1,2,3,4,5,6])
def checkPascal(numTrials, roll):
yes = 0.0
for i in range(numTrials):
for j in range(24):
d1 = roll()
d2 = roll()
if d1 == 6 and d2 == 6:
yes += 1
break
print 'Probability of losing =', 1.0 - yes/numTrials
##checkPascal(10000, rollDie)
如果我们对🎲作弊
def rollLoadedDie():
# 因为如果不作弊，条件应该是random.random() < 1.0/6.0
# 因为有1/6机率🎲是6。分母变小，也就是结果变大，既是6的概率变大。
if random.random() < 1.0/5.5:
return 6
else:
return random.choice([1,2,3,4,5])
##checkPascal(10000, rollLoadedDie)
结果变成了0.4428，输的概率下降，赢的上升。
上述例子可以说明：simulation is a wonderful tool for exploring what if questions。
蒙特卡罗方法（英语：Monte Carlo method），也称统计模拟方法，是二十世纪四十年代中期由于科学技术的发展和电子计算机的发明，而被提出的一种以概率统计理论为指导的一类非常重要的数值计算方法。是指使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。
20世纪40年代，在冯·诺伊曼，斯塔尼斯拉夫·乌拉姆和尼古拉斯·梅特罗波利斯在洛斯阿拉莫斯国家实验室为核武器计划工作时，发明了蒙特卡罗方法。因为Ulam的叔叔经常在蒙特卡罗赌场输钱得名，而蒙特卡罗方法正是以概率为基础的方法。
与它对应的是确定性算法。
蒙特卡罗方法在金融工程学，宏观经济学，生物医学，计算物理学（如粒子输运计算、量子热力学计算、空气动力学计算）等领域应用广泛。
这个想法的得出是Ulam试图计算solitaire（纸牌）赢的概率时。发现计算起来实在太复杂。能不能用统计的方法来算出它的概率来。
A Monte Carlo simulation is a method of estimating the value of an unknown quantity using the principles of inferential statistics.
关键词是estimation，因为Monte Carlo simulation得出的不是准确结果。结果是inferential，是推断出的。
它的依据是：a random sample tends to exhibit the same properties as a population from which it is drawn.
例子：
- 如果你硬币丢100次，100次都是head，那你可以推断出，硬币两边都是head。
- 但如果head是52次，tail是48次。那么再丢100次的结果就不一定是头多尾少了。
def flip(numFlips):
heads = 0
for i in range(numFlips):
if random.random() < 0.5:
heads += 1
return heads/float(numFlips)
for i in range(5): #number of trials
print flip(10)
如果每轮扔10次，计算出的结果偏差很大。我们推断不出coin正反的概率。 每轮扔1,000,000次，计算结果就很接近。我们就此可以基本推断出硬币正反概率是各50％。那为什么呢？
Buffon-Laplace method of estimating pi by dropping needles
布冯投针求π。
在一个变长为2的正方形中画一个内切圆。 如果完全随机的情况下投针。
圆内针的数目：所有针的数目 ＝ 圆的面积：正方形的面积
因而：
圆的面积 ＝ （正方形面积＊圆内针数量）／ 所有针的数量。
π＊r2 = 4＊（圆内针数量／所有针数量）（内切圆r＝1）所以：
π ＝ 4＊（圆内针数量／所有针数量）
因为人不可能完全随机的丢针，但是计算机可以。 用代码表示：
def throwNeedles(numNeedles):
inCircle = 0
for Needles in xrange(1, numNeedles + 1, 1):
x = random.random()
y = random.random()
if (x*x + y*y)**0.5 <= 1.0:
inCircle += 1
return 4*(inCircle/float(numNeedles))
# x,y表示针的坐标，都 < 1。如果他们距原点的距离 <= 1。说明在圆内。
def getEst(numNeedles, numTrials):
estimates = []
for t in range(numTrials):
piGuess = throwNeedles(numNeedles)
estimates.append(piGuess)
sDev = stdDev(estimates)
curEst = sum(estimates)/len(estimates)
print 'Est. = ' + str(curEst) +\
', Std. dev. = ' + str(round(sDev, 6))\
+ ', Needles = ' + str(numNeedles)
return (curEst, sDev)
simulation一定要做多次，并且看他们标准方差，这样不仅仅是得到答案，还要知道答案有多准确。getEst() 这里给出了多次测试的算术平均curEst和标准方差sDev。
def estPi(precision, numTrials):
numNeedles = 1000
sDev = precision
while sDev >= precision/2.0:
curEst, sDev = getEst(numNeedles, numTrials)
numNeedles *= 2
return curEst
random.seed(0)
estPi(0.005, 100)
所以我们一直算到标准方差小于0.005的一半为止。如果标准方差大于它，就把numNeedle翻倍再算。
Est. = 3.14844, Std. dev. = 0.047886, Needles = 1000
Est. = 3.13918, Std. dev. = 0.035495, Needles = 2000
Est. = 3.14108, Std. dev. = 0.02713, Needles = 4000
Est. = 3.141435, Std. dev. = 0.016805, Needles = 8000
Est. = 3.141355, Std. dev. = 0.0137, Needles = 16000
Est. = 3.14131375, Std. dev. = 0.008476, Needles = 32000
Est. = 3.141171875, Std. dev. = 0.007028, Needles = 64000
Est. = 3.1415896875, Std. dev. = 0.004035, Needles = 128000
Est. = 3.14174140625, Std. dev. = 0.003536, Needles = 256000
Est. = 3.14155671875, Std. dev. = 0.002101, Needles = 512000
虽然计算结果64000的Est比32000的Est看上去偏差要大。但是如果把标准方差算进来，其实64000的范围比32000的范围要小。
It is not sufficient to produce a good answer. We have to have a valid reason to be confident that in fact, it really is a good answer.
但是good answer不能只看标准方差是不是随着测试基数增加而减小。因为这只是意味着，我们从一个集合中取的sample越多，我们的到的结果越接近。这并不代表数值正确。
要确保结果正确，必须做到以下三点：
- Our conceptual model is correct。从代数角度对π的算法正确
- Our implementation is correct。code里面可能有bug
- We have enough samples to actually believe the result。这个是标准方差做的事情。
最后，还要进行实际的测试，看看结果是否与simulation相符。
Read more: