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 的小區間。先來測試下指數函數的積分。

\[\int_{1}^{4} e^{x}dx \]

CompIntegralbyladder(np.exp,1,4)

51.923094224367127

來看下正確答案。注意指數函數的不定積分還是指數函數本身。

from sympy.interactive import printing
printing.init_printing(use_latex=True)

\[\int_{1}^{4}e^{x}=e^{4}-e^{1} \]

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()

再看看計算個二次函數的積分。隨便寫個二次函數。

\[2x^{2}+3x+4 \]

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 的積分吧。

\[\int_{-5}^{5}2x^{2}+3x+4 dx \]

CompIntegralbyladder(Quadratic,-5,5)

206.69999999999825

下面來計算下正確的積分值。

因為：

\[\int2x^{2}+3x+4dx=\frac{2}{3}x^3+\frac{3}{2}x^2+4x \]

所以：

\[\int_{-5}^{5}2x^{2}+3x+4dx=[\frac{2}{3}x^3+\frac{3}{2}x^2+4x]_{-5}^{5} \]

def Integral(x):
    return (2*x**3)/3 + (3*x**2)/2 + 4*x

Integral(5)-Integral(-5)

207

可見，梯形法計算積分還是比較准確的。