1个掷硬币问题,4个Python解法:读书笔记

1个掷硬币问题,4个Python解法:读书笔记

王勇王勇

关键词:统计,概率,机器学习,Pandas, Numpy, sympy scipy

预计阅读时间-10分钟

我在学习机器学习算法和玩Kaggle 比赛时候,不断地发现需要重新回顾概率、统计、矩阵、微积分等知识。如果按照机器学习的标准衡量自我水平,这些知识都需要重新梳理一遍。

网上或许有各种各样知识片断,却较难找到一本书将概率,统计、矩阵、微积分公式和Python结合起来。 要么是讲的比较浅显,要么跨度比较大。 最近看到一本书,恰好把上面的问题解决了。着重讲解Python for 概率,统计,机器学习. 相比较吴恩达教授的网上课程,各有千秋。

  • l 感觉吴教授的课程偏学术。 假设学生数学基础比较扎实(毕竟课堂受众是斯坦福大学学生,起点比较高),用矩阵推导公式,简洁明了。整个课程把绝大多数机器学习算法都教了一遍。
  • 这本书偏工程实践。 书中从多个层面来介绍经典算法。尤其是后期的泛化,正则化等章节。介绍的算法,但是每个算法都用2-5种python方法实现。例如: 
Python 循环或自带Itertools ((笛卡尔乘积,经典概率)


Python sympy(数学符号) (微积分公式推导和实现)


Python Pandas(分组计算) (程序员看得懂)


Python numpy (矩阵计算) (注:用矩阵计算,有速度飞起来的感觉)



Python scipy (科学计算库) (算法增强器)



个人感觉这本书比较适合我的学习目标。 也许也有人喜欢这样的书。书名和下载地址在文章最后面。

我先来翻译一段书中的一道期望计算题目,分享一下这种庖丁解牛和层次渐近的感觉。


题目:

三个硬币: 1角,2角,5角。 同时掷硬币,正面朝上的将面值加在一起求和。 只有两个硬币正面朝上的期望和是多少?  


  • Xi ∈ {0, 1} 注:硬币为Xi, (面值10,10,50,分别为X10,X20,X50,只有正面和反面,是服从二项分布(0,1)
  • ξ := 10X10 + 20X20+ 50X50 注: ξ 为正面朝上的硬币面值之和
  • η := X10X20(1 −X50) + (1 − X10)X20X50 + X10(1 − X20)X50 注: η 为只有两个硬币正面朝上的情况


这样此题就变成了计算 E(ξ |η) 在 η 条件为真时,ξ 的期望。我们首先需要找到一个函数 h(η)。 这个函数可以让残差最小化。

现在,计算两个硬币朝上的面值之和公式变成了如何定义h(η)函数。

注:η的结果是{0,1},所以h函数的只有两种输入值{0,1}。因此,正交内积条件为

因为我们只需要知道满足两个硬币朝上的情况,(即η =1 ),所以公式简化为:

两边积分计算求和


数学公式到此就结束了。本书中定义h(η) = αη,并求α。这样公式就变成了如下形式。


注:


公式推导完了,下面就看看Python的四种解法吧。


解法1 :Sympy数学符号方法

上述推导公式,直接可以用数学符号语言,在Sympy中计算。

计算结果精准 alpha = 160/3


E|η) = (160/3)*η

# Do it by sympy
import sympy as S
X10,X20,X50 = S.symbols("X10,X20,X50",real=True)

xi = 10*X10+20*X20+50*X50
print("ξ = ",xi)
eta = X10*X20*(1-X50)+X10*(1-X20)*(X50)+(1-X10)*X20*(X50)
print("η = ",eta)

num=S.summation(xi*eta,(X10,0,1),(X20,0,1),(X50,0,1))

den=S.summation(eta*eta,(X10,0,1),(X20,0,1),(X50,0,1))

alpha=num/den

print(alpha)


ξ = 10*X10 + 20*X20 + 50*X50
η = X10*X20*(-X50 + 1) + X10*X50*(-X20 + 1) + X20*X50*(-X10 + 1)
160/3

解法2 :用Pandas,分组和计算(程序员的方式)

注:做1000次试验(蒙特卡罗仿真),然后计算试验的均值。计算结果近似于推导结果。

# Do it by pandas
import pandas as pd
d = pd.DataFrame(columns=['X10','X20','X50'])
ntest = 10**6
d.X10 = np.random.randint(0,2,ntest)
d.X10 = np.random.randint(0,2,ntest)
d.X20 = np.random.randint(0,2,ntest)
d.X50 = np.random.randint(0,2,ntest)

grp=d.groupby(d.eval('X10+X20+X50'))
grp.get_group(2).eval('10*X10+20*X20+50*X50').mean()


53.340141339548644


解法3:用Numpy,矩阵计算(速度快,有飞起来的感觉)

# Do it by numpy
import numpy as np
from numpy import array
ntest = 10**6
x=np.random.randint(0,2,(3,ntest)) 

expectation =np.dot(x[:,x.sum(axis=0)==2].T,array([10,20,50])).mean()

print(expectation)


53.3542987559

解法4: 用笛卡尔笛卡尔乘积,过滤只有两个硬币朝上事件,计算期望

#do it by pure python
import itertools as it
Xi = list(it.product((0,1),(0,1),(0,1)))               #笛卡尔乘积函数 初始化 - 所有掷硬币的事件 8种情况,三个硬币,每个硬币只有正面和反面

Xi_conditioned = list(it.filterfalse(lambda i:sum(i)!=2,Xi)) #过滤 只有两个硬币正面朝上的事件,只有三个事件

results = list((map(lambda k:10*k[0]+20*k[1]+50*k[2],Xi_conditioned)))
dem = len(results)
num = sum(results)

print(Xi)
print(Xi_conditioned)
print(results)
print(num/dem)


[(0, 0, 0), (0, 0, 1), (0, 1, 0), (0, 1, 1), (1, 0, 0), (1, 0, 1), (1, 1, 0), (1, 1, 1)]
[(0, 1, 1), (1, 0, 1), (1, 1, 0)]
[70, 60, 30]
53.333333333333336


这道概率期望题,书中演示了4种方法: Sympy,Numpy, Pandas 和Itertools. 在科学计算和机器学习中,采用不同的实现方法可以有助于问题解决和交叉检查。


最后分享一下这本书的名字:

. Python for Probability, Statistics, and Machine Learning

网上有下载。


如果觉得本文对你有帮助,请随手点赞!


分享是最好的投资!


个人文章链接

Kaggle HousePrice : LB 0.11666(前15%), 用搭积木的方式(2.实践-特征工程部分)

Kaggle HousePrice : LB 0.11666(排名前15%), 用搭积木的方式(1.原理)

用Python分析指数: 11月16日热门指数Z值表

用Python分析指数: 10月18日指数高低Z值表。

沪深300进入低估区了吗? 没有。好像中证环保现在低估了。

一个韭菜用Python采集,清洗和分析中证指数

文章被以下专栏收录
2 条评论
推荐阅读