<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Python数值计算库NumPy
—— 一切向量化计算的基础
1. NumPy初探
1.1 开始使用
Step1:  从三种“列表“定义开始
Step2:  IPython能够支持自动补全和帮助：
Step3:  论numpy的正确打开方式：少用原生语法、多用ufunc和broadcasting，少用数据转换
Step4:  为了比较性能，使用ipython的“魔法函数”timeit，或者datetime库的计时函数：
Step5:  1.2 Numpy数据类型理解
Step6:  int8
Step7:  Python不显示指定数据类型与精度，如何控制数据类型？Numpy自定义
Step8:  1.3 Numpy数据操作
数据切片（Slicing）变形（Reshaping）拼接（Stacking）与拆分（Splitting）
Step9:  2 NumPy数学函数
2.0 使用Random
使用其中的rand、randint等方法生成随机数来方便后续的测试。
随机数生成函数与统计分布将在统计基础课程中另行介绍。
Step10:  2.1 数学函数的测试
统计函数返回一个值：
Step11:  Ufunc类函数返回逐个计算值：
很多Ufunc即Universal Functional在Python中是以C语言级别实现的，效率高，建议使用。
Step12:  过滤值：类似SAS的Compress（其实是反过来）
Step13:  过滤值：类似R语法
Step14:  类似R的语法中，用于过滤条件的np.array是可以进行布尔操作的：
Step15:  Python的Clip，对值进行裁剪
Step16:  可以方便的进行时间序列操作——差分
Step17:  差分其实是一种最简单的卷积：
Step18:  差分的反向操作——累计加和（累积乘法类似）：
Step19:  排序时间到了：给出排序结果，本地排序，并保留原数组下标位置
Step20:  使用enumerate和sorted排序：
关于更多匿名函数（Lambda）使用，请参考Python函数式编程，也就是Functional Programming相关内容。
Step21:  3 二维数组与矩阵操作
3.1 基本定义
直接定义值：
Step22:  直接填充常用值：
Step23:  善于利用函数者，可以使用函数，传入下标作为参数，生成二维数组：
”首相卡梅隆“函数：
Step24:  习惯R者，使用broadcasting，生成二维数组，更具体的原理和应用在后续介绍：
9 * 1的二维数组和长度为9的一维数组broadcasting两次，互相补齐对方较短的维度：
Step25:  更进一步缩减代码，生成正交栅格用于广播的方法，支持灵活的运算：
状如切片，功能如range/linspace，实数第三参数代表步长，虚数第三参数代表总长度
Step26:  习惯数据库的笛卡尔积者，可以用outer关键字进行某些运算，支持传入iterable更能节省空间和时间：
Step27:  二维数组和矩阵有一些不同：
矩阵可以方便的求逆、直接做乘法，二维数组的＊只能作为np.multiply而不是np.dot：
Step28:  一般最小二乘法求解线性回归（OLS）
要进行OLS，核心是矩阵求逆：
$ y = X \beta + \beta_0 \rightarrow y ＝ X \beta + \beta_0 + \epsilon \rightarrow \hat{\beta} = (X^TX)^{-1}X^Ty $
随机制造数据来做个OLS的例子：
Step29:  给出模型参数的估计：$\hat{\beta} = (X^TX)^{-1}X^Ty $
Step30:  对照statsmodels的OLS结果：
Step31:  Page-Rank算法的核心
Step32:  4. 还有一件事……文件读写
Numpy的分析功能虽然大部分被Pandas替代，但文件读写仍然是一个科学计算库的基本功。
熟悉Matlab和R语言的人对下列语法不会陌生：
Step33:  保存多个变量时，等号前面是变量的key，从结果r里面可以取到。
如果不使用key来保存a，那么通过“arr_0”,"arr_1"来取到。
Step34:  最后，对txt，csv，tsv的读写：
csv意为comma separated values
tsv意为tab separated values
<ASSISTANT_TASK:>
Python Code:
from numpy import cos,sin #避免使用
import numpy as np #np.method()
Explanation: Python数值计算库NumPy
—— 一切向量化计算的基础
1. NumPy初探
1.1 开始使用
End of explanation
r1 = range(5)
r2 = np.arange(5)
r3 = xrange(5)
print r1,r2,r3
for i in r1:
    print i,
print '\n'
for i in r2:
    print i,
print '\n'
for i in r3:
    print i,
print type(r1),type(r2),type(r3)
print np.arange(0,5),np.arange(5,0),np.arange(5,0,-1)
Explanation: 从三种“列表“定义开始
End of explanation
np.arange #Tab
np.arange?
Explanation: IPython能够支持自动补全和帮助：
End of explanation
def vecadd(n):
    a = range(n)
    b = range(n)
    c = []
    for i in range(n):
        a[i] = i*2
        b[i] = i+3
        c.append(a[i]+b[i])
    return c
def np_vecadd(n):
    a = np.arange(n)*2
    b = np.arange(n)+3
    c = a+b
    return c
N = 10000000
Explanation: 论numpy的正确打开方式：少用原生语法、多用ufunc和broadcasting，少用数据转换
End of explanation
%timeit vecadd(N)
%timeit np_vecadd(N)
import datetime
t1 = datetime.datetime.now()
result = vecadd(N)
t2 = datetime.datetime.now()
t2-t1
np.arange(5)+3 #ufunc or broadcasting
np.arange(5)*2
a,b,c = 0,2,11
print np.linspace(a,b,c),'\n'*2,np.arange(c)/float(c-1)*(b-a)+a
Explanation: 为了比较性能，使用ipython的“魔法函数”timeit，或者datetime库的计时函数：
End of explanation
#使用np.ndarray类型的r2
print r2,type(r2),r2.dtype,r2.shape
#Immunity of Mutables
arr1 = np.array([np.arange(5),np.arange(5)])
arr2 = np.array([r2,r2])
arr3 = [r2,r2]
print arr1,'\n'*2,arr2,'\n'*2,arr3
r2[2] = 10
print arr1,'\n'*2,arr2,'\n'*2,arr3
r2[2] = 2
#不需要关心类型是否可变，在np.array创建过程中就已经建立独立拷贝
print arr1,arr1.shape
print arr1[0,0],arr1[1,2],arr1.dtype
print np.float64(arr1) #正确
arr1.dtype = np.float64 #错误
print arr1
arr1.dtype = np.int64
print arr1
# bool : 单个bit的布尔变量，仅存True或者False
# inti : 当前平台的整数类型，int32或int64
import platform
platform.architecture()
Explanation: 1.2 Numpy数据类型理解
End of explanation
i = np.float64(3)
print i
c = np.complex128(1+2j)
print c.real,c.imag
print i.dtype.itemsize,c.dtype.itemsize
Explanation: int8:     -128   ~ 127
int16:    -32768 ~ 32767
int32:    -$2^{31}$ ~ $2^{31}$-1
int64:    -$2^{63}$ ~ $2^{63}$-1
uint8:         0 ~ 255
uint16:        0 ~ 65535
uint32:        0 ~ $2^{32}$-1
uint64:        0 ~ $2^{64}$-1
float16:   符号＋5bit指数＋10bit小数
float32:   符号＋8bit指数＋23bit小数
float64:   符号＋11bit指数＋52bit小数
complex64:  单精度浮点复数
complex128: 双精度浮点复数
End of explanation
sku = np.dtype([('sku_id',np.int32),('desc',np.str_,50),('value',np.float64)])
print sku
sku2 = np.dtype({'names':['sku_id','desc','value'],'formats':['<i4','S50','<f8']})
print sku2
online_shop = np.array([(1,'apple',2.3),(2.1,5,5),(3,'banana',True)],dtype=sku)
print online_shop
Explanation: Python不显示指定数据类型与精度，如何控制数据类型？Numpy自定义
End of explanation
arr = np.arange(12)
print arr,arr.shape
print arr[3:5],arr[-5:-2],arr[3:28:5]
arr = arr.reshape(3,2,2)
print arr
arr.resize((3,4))
print arr,'\n',arr.ndim,'\n',arr.shape
print arr.T,'\n'*2,arr.transpose()
arrflat = arr.flat
print arrflat
print np.array([item for item in arrflat],dtype=np.int64)
print arr.flatten(order='C'),'\n'*2, arr.flatten(order='F'),'\n'*2,type(arr.flatten())
print arr.ravel(),'\n'*2,type(arr.ravel())
print np.vstack((arr,arr)),'\n'*2,np.row_stack((arr,arr)),'\n'*2, np.concatenate((arr,arr),axis=0)
print np.hstack((arr,arr)),'\n'*2,np.column_stack((arr,arr)),'\n'*2, np.concatenate((arr,arr),axis=1)
print np.dstack((arr,arr))
arr.resize((3,4,1))
print arr,'\n'*2,np.concatenate((arr,arr),axis=2)
print arr,'\n'*2,np.tile(arr,(1,1,3))
arr2 = np.arange(8).reshape(2,2,2)
print arr2
print np.vsplit(arr2,2), '\n'*2, np.split(arr2,2,axis=0)
print np.hsplit(arr2,2), '\n'*2, np.split(arr2,2,axis=1)
print np.dsplit(arr2,2), '\n'*2, np.split(arr2,2,axis=2)
Explanation: 1.3 Numpy数据操作
数据切片（Slicing）变形（Reshaping）拼接（Stacking）与拆分（Splitting）
End of explanation
arr_random = np.random.randint(1,6,15)
print arr_random
Explanation: 2 NumPy数学函数
2.0 使用Random
使用其中的rand、randint等方法生成随机数来方便后续的测试。
随机数生成函数与统计分布将在统计基础课程中另行介绍。
End of explanation
print arr_random.max(),arr_random.min(),arr_random.mean(),arr_random.var(),arr_random.std()
Explanation: 2.1 数学函数的测试
统计函数返回一个值：
End of explanation
print np.log(arr_random),'\n'*2,np.exp(arr_random),'\n'*2,arr_random>2,'\n'*2,np.mod(arr_random,2)
Explanation: Ufunc类函数返回逐个计算值：
很多Ufunc即Universal Functional在Python中是以C语言级别实现的，效率高，建议使用。
End of explanation
arr_random.compress(arr_random>2)
Explanation: 过滤值：类似SAS的Compress（其实是反过来）
End of explanation
arr_random[arr_random>2]
Explanation: 过滤值：类似R语法
End of explanation
print arr_random[(arr_random>1) & (arr_random<4)],
print arr_random[(arr_random>1) ^ (arr_random<4)]
Explanation: 类似R的语法中，用于过滤条件的np.array是可以进行布尔操作的：
End of explanation
arr_random.clip(2,4)
Explanation: Python的Clip，对值进行裁剪
End of explanation
arr_random = np.random.rand(10)
print arr_random,'\n'*2,np.diff(arr_random) 
Explanation: 可以方便的进行时间序列操作——差分
End of explanation
w = np.array([1,-1]) 
len = 2
print arr_random,'\n'*2,np.convolve(w,arr_random),'\n'*2,np.convolve(w,arr_random)[len-1:-len+1]
Explanation: 差分其实是一种最简单的卷积：
End of explanation
diff = np.convolve(w,arr_random)[0:-1]
print diff
print diff.cumsum(),'\n'*2,diff.cumprod()
Explanation: 差分的反向操作——累计加和（累积乘法类似）：
End of explanation
arr_random = np.random.rand(10)
print np.sort(arr_random)
print arr_random
arr_random.sort()
print arr_random
arr_random = np.random.rand(10)
Explanation: 排序时间到了：给出排序结果，本地排序，并保留原数组下标位置
End of explanation
s = sorted(enumerate(arr_random),key=lambda x:x[1],reverse=False)
print s
zip(*s)[0]
print type(enumerate(arr_random))
for idx,value in enumerate(arr_random):
    print idx,value
print zip(xrange(10),arr_random)
print np.lexsort((np.arange(10),arr_random)) #接受一个参数，传入元组
print np.argsort(arr_random) #相当于lexsort参数元组第一列给下标
Explanation: 使用enumerate和sorted排序：
关于更多匿名函数（Lambda）使用，请参考Python函数式编程，也就是Functional Programming相关内容。
End of explanation
arr2d_random = np.random.rand(10).reshape(5,2)
arr2d_zeros = np.zeros_like(arr2d_random)
arr2d_ones = np.ones_like(arr2d_random)
print arr2d_random,'\n'*2,arr2d_zeros,'\n'*2,arr2d_ones
Explanation: 3 二维数组与矩阵操作
3.1 基本定义
直接定义值：
End of explanation
print np.zeros((5,2)), '\n'*2 , np.ones((5,2)),'\n'*2,np.eye(3)
Explanation: 直接填充常用值：
End of explanation
def CameronThePrimeMinister(i,j):
    return (i+1)*(j+1)
print np.fromfunction(CameronThePrimeMinister,(9,9))
Explanation: 善于利用函数者，可以使用函数，传入下标作为参数，生成二维数组：
”首相卡梅隆“函数：
End of explanation
x = np.arange(1,10)
print x.reshape(-1,1),'\n'*2,x,'\n'*2,x.reshape(-1,1) * x
Explanation: 习惯R者，使用broadcasting，生成二维数组，更具体的原理和应用在后续介绍：
9 * 1的二维数组和长度为9的一维数组broadcasting两次，互相补齐对方较短的维度：
End of explanation
x,y = np.ogrid[0:1:6j,0:1:6j]
print x,y
print np.exp(-x**2-y**2)
Explanation: 更进一步缩减代码，生成正交栅格用于广播的方法，支持灵活的运算：
状如切片，功能如range/linspace，实数第三参数代表步长，虚数第三参数代表总长度
End of explanation
np.multiply.outer(xrange(1,10),xrange(1,10))
Explanation: 习惯数据库的笛卡尔积者，可以用outer关键字进行某些运算，支持传入iterable更能节省空间和时间：
End of explanation
cov =  np.cov(arr2d_random.T)
print cov
stdiag = np.diag(np.sqrt(np.diag(cov)))
print stdiag
invstdiag = np.array(np.mat(stdiag).I)
print invstdiag
#错误
invstdiag*cov*invstdiag
#正确
np.mat(stdiag).I*np.mat(cov)*np.mat(stdiag).I
#验证
np.corrcoef(arr2d_random.T)
Explanation: 二维数组和矩阵有一些不同：
矩阵可以方便的求逆、直接做乘法，二维数组的＊只能作为np.multiply而不是np.dot：
End of explanation
x = np.random.rand(30).reshape(15,2)
y = x[:,0]*0.7 - x[:,1]*0.2 + 0.1+0.1*np.random.rand(15)
print y,'\n'*2,x
Xcvr = np.mat(np.hstack((np.ones(15).reshape(15,1),x)))
print Xcvr
Explanation: 一般最小二乘法求解线性回归（OLS）
要进行OLS，核心是矩阵求逆：
$ y = X \beta + \beta_0 \rightarrow y ＝ X \beta + \beta_0 + \epsilon \rightarrow \hat{\beta} = (X^TX)^{-1}X^Ty $
随机制造数据来做个OLS的例子：
End of explanation
H = Xcvr*(Xcvr.T*Xcvr).I*Xcvr.T
betahats = np.dot((Xcvr.T*Xcvr).I*Xcvr.T,y)
preds = np.dot(H,y)
print betahats,'\n'*2,preds
Explanation: 给出模型参数的估计：$\hat{\beta} = (X^TX)^{-1}X^Ty $
End of explanation
import statsmodels.formula.api as sm
model = sm.OLS(y,Xcvr).fit()
model.summary()
Explanation: 对照statsmodels的OLS结果：
End of explanation
rand_relation = (np.random.rand(100)>0.7).reshape(10,10)
rand_relation.dtype=np.int8
rand_relation
rand_relation = rand_relation-np.diag(np.diag(rand_relation))
rand_relation
invrowsum = 1/(np.sum(rand_relation,axis=1)+1e-12)
print invrowsum
diag_invrowsum = np.diag(invrowsum)
tpm = np.dot(diag_invrowsum,rand_relation+1e-12/10)
tpm = (tpm>1e-12/10)*tpm
print tpm
import scipy.sparse as ssp
sp_tpm = ssp.csr_matrix(tpm)
sp_tpm
for i in xrange(10):
    sp_tpm = sp_tpm.dot(sp_tpm)
    print repr(sp_tpm)
sp_tpm_a = sp_tpm.toarray()
print sp_tpm_a,'\n'*2,sp_tpm_a[0,:]
Explanation: Page-Rank算法的核心: 矩阵相乘
关系矩阵$R(i,j)$元素表示从第i个页面到第j个页面发生了跳转。
初始向量w反复乘以按行标准化后的概率转移矩阵（要保证行和全部为1）。
当内存空间足够时，直接用矩阵迭代乘法将矩阵与自己相乘，空间换时间。
End of explanation
a = np.random.rand(10).reshape(5,2)
print a
np.save('a.npy',a)
a_copy = np.load('a.npy')
print a_copy
Explanation: 4. 还有一件事……文件读写
Numpy的分析功能虽然大部分被Pandas替代，但文件读写仍然是一个科学计算库的基本功。
熟悉Matlab和R语言的人对下列语法不会陌生：
End of explanation
b = range(5)
np.savez('a_and_b.npz',a,bdata=b)
r = np.load('a_and_b.npz')
print r["arr_0"],'\n'*2,r["bdata"]
Explanation: 保存多个变量时，等号前面是变量的key，从结果r里面可以取到。
如果不使用key来保存a，那么通过“arr_0”,"arr_1"来取到。
End of explanation
np.savetxt("a.txt",a,delimiter = '\t')
a_copy2 = np.loadtxt("a.txt",delimiter = '\t')
print a_copy2
Explanation: 最后，对txt，csv，tsv的读写：
csv意为comma separated values
tsv意为tab separated values
End of explanation
<END_TASK>