Numpy 练习题
1. 使用循环和向量化两种不同的方法来计算 100 以内的质数之和。
先定义个判断质数的函数。ps：纯手工打造，原生态，哈哈。
def checkprime(x):
if x<=1:
return False;
prime=True;
for i in range(2 , 1+x/2):
if x%i == 0:
prime = False;
break;
return prime;
使用循环方法来计算 100 以内的质数之和。
def sumprimebyiter(n=100):
primesum=0
for i in range(1, n+1):
if( True == checkprime(i)):
primesum += i
return primesum
%timeit sumprimebyiter(100)
10000 loops, best of 3: 138 µs per loop
使用向量化的方法来计算 100 以内的质数。ps：怎么将判断质数的函数应用到向量中的每一个元素，可是花了好几分钟来寻找，终于发现 map 可以实现这个功能。后来又发现 np.vectorize 也可以实现同样功能。
import numpy as np
def sumprimebyarr(n=100):
a = np.arange(1,n+1)
# return sum(a[np.array(map(CheckPrime, a))]) # 此处是之前用 Python 自带的 map 把函数应用到向量的每个元素
check_prime_vec = np.vectorize(CheckPrime) # 此处代码用到了 np.vectorize，可以把外置函数应用到向量的每个元素
return np.sum(a[check_prime_vec(a)])
%timeit sumprimebyarr(100)
10000 loops, best of 3: 204 µs per loop
上面两种方法都使用魔术函数 %timeit 计算了执行时间，意外的是，向量化的方法竟然没有循环快，一定是哪儿不对，待我好好检查下，再补充该题答案。
2. 模拟一个醉汉在二维空间上的随机漫步。
先试试一维空间上的随机漫步。既然本周学了 Numpy，这里就直接上 np.random 了。
%pylab inline
Populating the interactive namespace from numpy and matplotlib
nsteps = 1000
draws = np.random.randint(-1,2,size=nsteps)
walk = draws.cumsum()
plot(walk)
[<matplotlib.lines.Line2D at 0x7f52c3534250>]
再来看下二维的。
nsteps = 1000
draws = np.random.randint(-1,2,size=(2,nsteps))
walks = draws.cumsum(1)
plot(walks[0,:],walks[1,:])
[<matplotlib.lines.Line2D at 0x7f52c34cf3d0>]
先生成 1000 个随机漫步方向，方向是从 {-1, 0, 1} 中随机挑两个值（两个值也可相等）作为移动方向，所以每次移动有 3×3=9 种选择，初始位置也是 9 种选择，cumsum 函数是将每次的移动累加，最后通过 plot 画出来。
代码调通之前，我都没想到代码能这么少。Python 还是好用，这要是用 C++ 写，得多少代码啊。
3. 使用梯形法计算一个二次函数的数值积分。
梯形法计算数值积分，就是把自变量分成无数小段，每一小段的面积用一个梯形面积近似，当小段的个数无限多，小段的长度无限小时，所有的小梯形面积加起来，就近似等于该函数的数值积分。如下图所示：
原理很简单，代码也很简单：
import numpy as np
def CompIntegralbyladder(func,x0,x1):
wholearea = 0
step = 0.1
for i in np.arange(x0, x1, step):
wholearea += (func(i)+func(i+step))*step/2; # Compute the Trapezoidal area
return wholearea;
该函数可以计算任意函数的积分。函数写好了，都分隔成长度为 0.1 的小区间。先来测试下指数函数的积分。
CompIntegralbyladder(np.exp,1,4)
51.923094224367127
来看下正确答案。注意指数函数的不定积分还是指数函数本身。
from sympy.interactive import printing
printing.init_printing(use_latex=True)
np.exp(4)-np.exp(1)
51.879868204685188
附上指数函数的图形。
%pylab inline
Populating the interactive namespace from numpy and matplotlib
import numpy as np
x = np.linspace(-5, 5, num = 100)
y = np.exp(x)
import matplotlib.pyplot as plt
plt.plot(x,y)
plt.show()
再看看计算个二次函数的积分。随便写个二次函数。
def Quadratic(x):
return 2*x**2 + 3*x + 4
先来看看这个函数的图形。
import numpy as np
x = np.linspace(-5, 5, num = 100)
y = Quadratic(x)
import matplotlib.pyplot as plt
plt.plot(x,y)
plt.show()
我们就计算该二次函数从 -5 到 5 的积分吧。
CompIntegralbyladder(Quadratic,-5,5)
206.69999999999825
下面来计算下正确的积分值。
因为：
所以：
def Integral(x):
return (2*x**3)/3 + (3*x**2)/2 + 4*x
Integral(5)-Integral(-5)
207
可见，梯形法计算积分还是比较准确的。
